 Research
 Open Access
 Published:
Estimation of feasible solution space using Cluster Newton Method: application to pharmacokinetic analysis of irinotecan with physiologicallybased pharmacokinetic models
BMC Systems Biology volume 7, Article number: S3 (2013)
Abstract
Background
To facilitate new drug development, physiologicallybased pharmacokinetic (PBPK) modeling methods receive growing attention as a tool to fully understand and predict complex pharmacokinetic phenomena. As the number of parameters to reproduce physiological functions tend to be large in PBPK models, efficient parameter estimation methods are essential. We have successfully applied a recently developed algorithm to estimate a feasible solution space, called Cluster Newton Method (CNM), to reveal the cause of irinotecan pharmacokinetic alterations in two cancer patient groups.
Results
After improvements in the original CNM algorithm to maintain parameter diversities, a feasible solution space was successfully estimated for 55 or 56 parameters in the irinotecan PBPK model, within ten iterations, 3000 virtual samples, and in 15 minutes (Intel Xeon E51620 3.60GHz × 1 or Intel Core i7870 2.93GHz × 1). Control parameters or parameter correlations were clarified after the parameter estimation processes. Possible causes in the irinotecan pharmacokinetic alterations were suggested, but they were not conclusive.
Conclusions
Application of CNM achieved a feasible solution space by solving inverse problems of a system containing ordinary differential equations (ODEs). This method may give us reliable insights into other complicated phenomena, which have a large number of parameters to estimate, under limited information. It is also helpful to design prospective studies for further investigation of phenomena of interest.
Background
Pharmacokinetics is a field of study that analyzes and predicts behaviors of drugs in organisms [1]. A major purpose of this study area is to predict pharmacokinetic properties of new drugs in humans, without performing clinical studies, in order to accelerate the efficiencies of new drug development processes. Another important purpose is to facilitate the proper use of not only newly developed drugs but also already existing drugs. There are numerous factors altering pharmacokinetics, such as drugdrug interactions (DDIs) [1], pharmacogenetics [2], or disease states [3, 4], which can cause large interindividual variability in drug responses. By studying these complicated phenomena, we may be able to explain and predict the alterations in clinical settings to administer drugs properly to each patient [5].
Physiologicallybased pharmacokinetic (PBPK) modeling and simulation are essential in understanding and predicting the abovementioned, complicated pharmacokinetic phenomena [6–8]. The basic idea of PBPK modeling is to reproduce physiological functions (absorption, distribution, metabolism, and elimination) in mathematical equations in order to understand the physiological phenomena extensively from the information obtained from in vivo studies and to predict unknown phenomena quantitatively. PBPK models are becoming more important because of the increase in complicated pharmacokinetic phenomena, such as DDIs involving multiple interaction sites [9] or the combined effect of multiple inhibitors [10]. Current draft guidance on DDI studies by the U.S. Food and Drug Administration [11] emphasizes the importance of PBPK simulation in deciding whether clinical DDI studies are required or not during new drug developments.
In analyzing complex pharmacokinetic phenomena with PBPK models, GaussNewton or its modified algorithm is often used for parameter estimation. However, these methods require feasible initial parameters beforehand, which are often difficult for the pathways with little or no prior information, such as pharmacokinetic parameters of metabolites, or of enterohepatic circulations (EHC). Optimized parameters may highly depend on the initial parameter settings as well. Additionally, the PBPK model contains a lot of parameters to estimate in nature, compared to the limited information available in clinical studies, making the accurate estimation of parameters more difficult. Since the accuracy of the parameter estimation process is quite important in PBPK analyses, these characteristics make the interpretations or extrapolations of the obtained results complicated.
Aoki et al. recently developed a new parameter estimation algorithm called CNM [12]. In this algorithm, we first prepare a group of virtual samples with random samplings from a certain initial range for each parameter to estimate. Then, linear approximations of a projection from parameter space into target values give the initial parameter values for the next iteration. In our experience, fewer than nine iterations of this process achieve the final, optimized parameters, which can reproduce clinically observed phenomena.
This algorithm has multiple advantages over conventional parameter estimating algorithms. The first advantage is the simplicity of the initial parameter settings. The new method only requires the designation of relatively broad parameter ranges as an initial setting, while the conventional algorithm requires the identification of feasible initial parameters. The second advantage is low computational costs owing to the deterministic nature of this algorithm, unlike other algorithms that use collective intelligence, such as genetic algorithms or particle swarm optimization. The last and the most important advantage is that the results are obtained as a group of optimized parameter sets. It allows us to interpret the phenomena with higher confidence and to extrapolate the obtained insights into new phenomena.
Aoki et al. previously applied CNM for analyzing a complicated pharmacokinetic phenomenon [5, 12] using a PBPK model [13], in which the pharmacokinetic properties of an anticancer drug, irinotecan (also known as CPT11), differs a lot between a bileduct cancer (BDC) patient and other cancer (OC) patients [14]. After intravenous administration, irinotecan is metabolized by CYP3A4 or carboxylesterase 2 (CES2) to form APC, NPC, SN38, and M4 [15, 16] (Figure 1). NPC and SN38 are further metabolized by CES2 and UGT1A to form SN38 and SN38G, respectively. Organic anion transporting polypeptide 1B1 (OATP1B1) is involved in the hepatic uptake of SN38, while OATP1B1 does not actively transport irinotecan or SN38G, according to an in vitro study [17]. SN38G is said to be deconjugated to form SN38 by βglucuronidase in intestinal microflora [18]. Possible involvement of enterohepatic circulation (EHC) in determining irinotecan pharmacokinetics [19] makes the estimation of parameters difficult because of the model structure and limited information about the feasible parameter values of EHC.
Since the previous report is mainly focused on the establishment of CNM, the PBPK model was not suitable from the pharmacokinetic viewpoint. Most importantly, the previous PBPK model did not contain EHC. Biliary drainage in the BDC group was not considered either. Furthermore, our preliminary investigation showed that the original CNM is not applicable to analyze the phenomena with EHC in a PBPK model.
In this paper, we improved both the PBPK model of irinotecan and the CNM algorithm itself. Firstly, we included EHC and biliary drainage in the PBPK model to properly interpret our obtained results. We also tried to improve the CNM algorithm itself to maintain the diversity of virtual samples during iterations, since the application of the original CNM algorithm failed using the new PBPK model in this study.
Results and discussion
Improvements in CNM for analyzing irinotecan accumulation profiles
In this study, we have applied and improved the CNM algorithm for the accumulation profiles of irinotecan and its metabolites by using the PBPK model shown in Figure 1 and in Additional File 1. Parameter estimation in the PBPK model was performed with 9 and 14 objective values for OC and BDC, respectively (Table 1c). Initial ranges of parameters were set as shown in Table 1b, which were thought to be large enough to contain the real value for each parameter.
Firstly, we performed the modified CNM with various dS values, ranging from 0 to 0.9. The modified CNM replicates a group of parameters towards the direction of the estimated feasible solution space with the ratio of dS after each iteration (see Methods section for details). We examined the effect of dS values on the convergence in this system, which was newly introduced in the algorithm to maintain parameter diversities as explained in the Methods section, and the calculation with dS of 0 is identical to the original CNM algorithm. By changing the dS values from 0 to 0.9, we performed CNM on the accumulation profiles in the BDC group with the same initial parameter ranges. As shown in Figure 2, the objective values converged well to the observed values when dS was no less than 0.5 (Figure 2pt). For the OC group, dS of at least 0.2 was needed to complete the process (data not shown).
When dS was smaller than 0.5, the ranges of the objective values remained similar to the initial range (Figure 2fj). Particularly, parameters got diverged after 2 iterations for dS of 0 to 0.2, or after 3 iterations for dS of 0.3 (Figure 2ad). The divergences of parameters after 3 iterations were smaller with higher dS values. It might be due to the capability to maintain parameter diversities with higher dS values. Since we have some parameters sensitive to the changes in dS values, such as a biliary clearance to a transit compartment of SN38 (CL_{bile}, parameter #12), we might be able to control the convergence by observing the behaviors of sensitive parameters. These results suggest the importance of maintaining parameter diversities in performing CNM using the dS algorithm, while further investigations should be needed to rationally determine feasible dS values in other cases.
Some pharmacokinetic parameters gradually converged to a feasible solution space along with the iterations of the parameter estimation process with high dS values. The residuals in objective values were decreased towards zero as well. Figure 3 illustrates the estimated parameters and the objective values after iterations with dS of 0.5, and the corresponding values are described in Additional File 2. While the ranges of parameter values seem stable only after the first iteration, the residuals in objective values continued to decrease until the third iterations. This observation may suggest that correlations among parameter values converged while maintaining the diversity of the absolute value in each parameter.
After 9 iterations, we observed high convergence of certain parameters, such as renal clearance (CL_{R}) of irinotecan (parameter #6) in both the OC and the BDC groups, and ratio of biliary clearance into Ttube to biliary clearance into a transit compartment (CL_{bile,T} / CL_{bile,transit}, parameter #56) in the BDC group (differences in convergence will be discussed later). In particular, the estimated values of irinotecan CL_{R} in OC (5.7 ± 1.2 ml/min/kg, geometric mean ± SD, Additional file 2) was close to the values calculated from the original report (2.8 ml/min/kg). Furthermore, the averages of estimated irinotecan CL_{R} values with two different initial distributions were within three fold of the values from the original report (data not shown). We do not have clear explanations of the convergences in clearance parameters without information on plasma concentrations, since these parameters are defined as the ratio between urinary, biliary, or fecal accumulations of drugs with area under the blood concentrationtime curve (AUC). We suspect that these parameters were determinants of the behavior of our PBPK model structures, and that the relationships with other parameters restricted the absolute values.
The calculation time required for the whole process was short; it took less than 15 minutes with ten iterations and 3000 virtual samples (Intel Xeon E51620 3.60GHz × 1 or Intel Core i78702.93GHz × 1; see Methods section for details).
After the application of CNM, we simulated the accumulation timeprofiles using estimated parameter sets to compare them with the reported values (Figure 4a, b). In both groups, some parameter sets reproduced the observed time profiles well, and the range of the estimated parameters became smaller when we selected samples which reproduced the observed time profile well (Figure 4c, d).
Differences in optimized parameters between OC and BDC
We compared the estimated parameter distributions between OC and BDC. When comparing the distribution of parameters with all samples, parameter distributions were similar in the two groups (Figure 5a). As mentioned in the previous section, we observed strong convergence in parameter #6, CL_{R} of irinotecan. Expansions of parameter ranges were also similar between the two groups. On the other hand, the distributions of the numbers of parameters were shifted when we considered accumulation timeprofiles (Figure 5b). Interestingly, these parameters are mostly involved in hepatobiliary or intestinal pathways (#1213, 19, 20, 35, 41). These results may suggest the importance of hepatobiliary and intestinal pharmacokinetic processes in determining the pharmacokinetics of irinotecan in the BDC patient, and are partly in a good agreement with the reports on reduced MRP2 expressions in liver with biliary cancer [20] or in the duodenums of patients with hepatic cholestasis [21], and altered pharmacokinetics of irinotecan and its metabolites with genetic polymorphisms of MRP2 [15, 16]. However, further investigations are needed to clarify the cause of the pharmacokinetic alterations, since most of the parameter ranges in BDC group are actually included in the parameter ranges in OC group.
Correlations among parameters
Finally, we have observed correlations among the estimated parameter values (Figure 6). We selected parameters with an r^{2} value of more than 0.64 against at least one other parameter in this figure. Tendencies in parameter correlations were similar between the two groups, while the number of parameter combinations having higher correlations was large in BDC group. The number of objective values might again be the cause of this observation.
Most of these relationships can be explained theoretically. For example, in both groups, parameters among #6, 16, and 18 have a strong linear relationship with each other to keep the ratio of parameters to maintain the ratio in the final outcome. Another example is the inverse proportional relationship of #56 to #12 and #14. When the amount of bile collected in the Ttube was increased (#56 increased), biliary clearances decreased (#12 and #14) to keep the amount of compounds excreted into the Ttubes. While these phenomena can be explained retrospectively, we cannot easily observe the actual parameter dependencies quantitatively without performing this kind of parameter estimation algorithms.
Conclusions
Efficient and reliable estimations of parameter values are important for the analysis of complex pharmacokinetic phenomena with PBPK models. In this study, we have successfully improved and applied CNM to estimate 55 or 56 parameters in the irinotecan PBPK model, by implementing dS function in the algorithm. The application of this method presented not only the control parameters in this PBPK model, but also correlations of parameters, which are important in determining the behaviors of the PBPK model.
Methods
Data source
Throughout thee study, we analyzed reported pharmacokinetic profiles of irinotecan and its metabolites [14]. This report described the accumulation of irinotecan and its metabolites into urine and feces for two groups of patients, with bileduct cancer (BDC, one patient) or with other cancers (OC, seven patient). Bile was also collected via biliary Ttube for a patient in the BDC group.
PBPK model and parameter settings
In order to analyze the causes of the different pharmacokinetic profiles from the two groups of patients, we developed a simplified PBPK model for irinotecan and its four metabolites (Figure 1), based on the reported metabolic pathways of irinotecan [15, 16]. The ordinary differential equations of this model are described in the Additional File 1.
Fixed parameter values and initial parameter ranges were shown in Table 1a and 1b, respectively. While these ranges were set arbitrarily, they may be large enough to contain real values, and may be applicable to many compounds having different kinetic constants. We used the same parameter for the ratio between CL_{bile,T} and CL_{bile} for irinotecan and the metabolites, since this parameter can be regarded as the ratio of bile flow between in Ttube and in bile duct after Ttube. In OC group, CL_{bile,T} was set to be equal to zero, since biliary Ttube was not placed for this group of patients. We have performed the whole parameter optimization processes in the log space of all the parameters for avoiding negative parameter values in the original space.
Objective values for the parameter estimation processes were shown in Table 1c. Since the deconjugation of SN38G to SN38 by βglucuronidase in intestinal microflora may affect the fecal elimination amount of SN38G and SN38, we combined the accumulation of these two compounds in objective values #7.
Materials
We solved the system of ODEs by using the MATLAB stiff ODE solver ODE15s [22]. We performed all calculations with MATLAB using a desktop computer (CPU: Core i7870 2.93GHz × 1, OS: Windows 7 SP1 32 bit, RAM: 4GB, MATLAB version: 8.0.0) or a workstation (CPU: XeonE51620 3.60GHz × 1, OS: CentOS 6.4 64 bit, RAM: 16GB, MATLAB version: 8.1.0).
CNM Method and its modification
CNM was constructed previously [12]. Briefly, a group of initial parameter sets (3000 virtual samples) was prepared with a random sampling from given parameter ranges. The linear approximations of the projections from one group of parameter sets into objective values generated the next group, and nine iterations of this process yielded a group of optimized parameter sets.
In each iteration process, we have newly implemented a calculation using a parameter called dS to maintain parameter diversities. In each iteration of the original CNM, we have parameter sets before (X_{b}) and after (X_{a}) the parameter estimation process. In our new strategy, we calculated internally dividing point X_{i} with the ratio of dS:(1dS), and applied the same inverse matrix to obtain new estimated parameters X_{a}'. Parameter sets for the next iteration were obtained by randomly selecting X_{a} or X_{a}' for each virtual sample.
After completing the estimations of parameters, the accumulation timeprofiles were compared with the observed profiles, using sum of squares of log residuals (SS_{log}):
where A_{e,simulated} and A_{e,observed} represent the simulated or observed amount of total irinotecan radioactivity at each time point.
Correlations of parameters
After the parameter estimation process, we have observed parameter correlations using PLOTMATRIX function in MATLAB. Parameters to be displayed in Figure 6 were those who have high correlation (r^{2}> 0.64) with at least one other parameter.
Abbreviations
 CL_{12}:

clearance from a rapid to a late equilibrium compartment
 CL_{R}:

renal clearance
 CL_{UGT}:

metabolic clearance of SN38 by UGT to form SN38G
 CNM:

Cluster Newton Method
 CPT11:

irinotecan
 k_{21}:

kinetic constant from a late to a rapid equilibrium compartment
 k_{a}:

absorption rate constant
 k_{feces}:

kinetic constant for the transit from large intestine to feces
 k_{LI}:

kinetic constants for the transit from small intestine to large intestine
 k_{transit}:

kinetic constant for the transit in bile compartments to small intestine
 Q_{H}:

blood flow rate in liver
 V_{H}:

volume of a liver
 V_{rapid}:

apparent volume of a rapid equilibrium compartment.
References
 1.
Buxton I, Benet L: Chapter 2. Pharmacokinetics: The Dynamics of Drug Absorption, Distribution, Metabolism, and Elimination. Goodman & Gilman's The Pharmacological Basis of Therapeutics. Edited by: Brunton L, Chabner B, Knollmann B. 2010, New York: McGrawHill, 1739. 12th
 2.
Relling M, Giacomini K: Chapter 7. Pharmacogenetics. Goodman & Gilman's The Pharmacological Basis of Therapeutics. Edited by: Brunton L, Chabner B, Knollmann B. 2010, New York: McGrawHill, 145168. 12th
 3.
Verbeeck R: Pharmacokinetics and dosage adjustment in patients with hepatic dysfunction. Eur J Clin Pharmacol. 2008, 64 (12): 11471161. 10.1007/s002280080553z.
 4.
Verbeeck R, Musuamba F: Pharmacokinetics and dosage adjustment in patients with renal dysfunction. Eur J Clin Pharmacol. 2009, 65 (8): 757773. 10.1007/s0022800906788.
 5.
Konagaya A: Towards an In Silico Approach to Personalized Pharmacokinetics. Molecular Interactions. 2012, Edited by Meghea A. Rijeka, Croatia: Intech, 263282.
 6.
Rowland M, Peck C, Tucker G: Physiologicallybased pharmacokinetics in drug development and regulatory science. Annu Rev Pharmacol Toxicol. 2011, 51: 4573. 10.1146/annurevpharmtox010510100540.
 7.
Zhao P, Zhang L, Grillo JA, Liu Q, Bullock JM, Moon YJ, Song P, Brar SS, Madabushi R, Wu TC: Applications of physiologically based pharmacokinetic (PBPK) modeling and simulation during regulatory review. Clin Pharmacol Ther. 2011, 89 (2): 259267. 10.1038/clpt.2010.298.
 8.
Fan J, Chen S, Chow EC, Pang KS: PBPK modeling of intestinal and liver enzymes and transporters in drug absorption and sequential metabolism. Curr Drug Metab. 2010, 11 (9): 743761. 10.2174/138920010794328931.
 9.
Yoshida K, Maeda K, Sugiyama Y: Hepatic and intestinal drug transporters: prediction of pharmacokinetic effects caused by drugdrug interactions and genetic polymorphisms. Annu Rev Pharmacol Toxicol. 2013, 53: 581612. 10.1146/annurevpharmtox011112140309.
 10.
Niemi M, Backman JT, Neuvonen M, Neuvonen PJ: Effects of gemfibrozil, itraconazole, and their combination on the pharmacokinetics and pharmacodynamics of repaglinide: potentially hazardous interaction between gemfibrozil and repaglinide. Diabetologia. 2003, 46 (3): 347351.
 11.
Draft Guidance/Guidance for Industry. Drug Interaction StudiesStudy Design, Data Analysis, Implications for Dosing, and Labeling Recommendations. [http://www.fda.gov/downloads/Drugs/GuidanceComplianceRegulatoryInformation/Guidances/UCM292362.pdf]
 12.
Aoki Y, Hayami K, Sterc HD, Konagaya A: Cluster Newton Method for Sampling Multiple. Solutions of an Underdetermined Inverse Problem: Parameter Identification for Pharmacokinetics. NII Technical Reports. 2011
 13.
Arikuma T, Yoshikawa S, Azuma R, Watanabe K, Matsumura K, Konagaya A: Drug interaction prediction using ontologydriven hypothetical assertion framework for pathway generation followed by numerical simulation. BMC bioinformatics. 2008, 9 (Suppl 6): S1110.1186/147121059S6S11.
 14.
Slatter JG, Schaaf LJ, Sams JP, Feenstra KL, Johnson MG, Bombardt PA, Cathcart KS, Verburg MT, Pearson LK, Compton LD: Pharmacokinetics, metabolism, and excretion of irinotecan (CPT11) following I.V. infusion of [(14)C]CPT11 in cancer patients. Drug Metab lism and Disposition: The Biological Fate of Chemicals. 2000, 28 (4): 423433.
 15.
Smith NF, Figg WD, Sparreboom A: Pharmacogenetics of irinotecan metabolism and transport: an update. Toxicol In Vitro. 2006, 20 (2): 163175. 10.1016/j.tiv.2005.06.045.
 16.
Innocenti F, Kroetz DL, Schuetz E, Dolan ME, Ramirez J, Relling M, Chen P, Das S, Rosner GL, Ratain MJ: Comprehensive pharmacogenetic analysis of irinotecan neutropenia and pharmacokinetics. J Clin Oncol. 2009, 27 (16): 26042614. 10.1200/JCO.2008.20.6300.
 17.
Nozawa T, Minami H, Sugiura S, Tsuji A, Tamai I: Role of organic anion transporter OATP1B1 (OATPC) in hepatic uptake of irinotecan and its active metabolite, 7ethyl10hydroxycamptothecin: in vitro evidence and effect of single nucleotide polymorphisms. Drug metabolism and disposition: the biological fate of chemicals. 2005, 33 (3): 434439.
 18.
Takasuna K, Hagiwara T, Hirohashi M, Kato M, Nomura M, Nagai E, Yokoi T, Kamataki T: Involvement of betaglucuronidase in intestinal microflora in the intestinal toxicity of the antitumor camptothecin derivative irinotecan hydrochloride (CPT11) in rats. Cancer Res. 1996, 56 (16): 37523757.
 19.
Younis IR, Malone S, Friedman HS, Schaaf LJ, Petros WP: Enterohepatic recirculation model of irinotecan (CPT11) and metabolite pharmacokinetics in patients with glioma. Cancer Chemother Pharmacol. 2009, 63 (3): 517524. 10.1007/s0028000807698.
 20.
Yamada T, Arai T, Nagino M, Oda K, Shoda J, Suzuki H, Sugiyama Y, Nimura Y: Impaired expression of hepatic multidrug resistance protein 2 is associated with posthepatectomy hyperbilirubinemia in patients with biliary cancer. Langenbecks Arch Surg. 2005, 390 (5): 421429. 10.1007/s0042300505645.
 21.
Dietrich CG, Geier A, Salein N, Lammert F, Roeb E, Oude Elferink RP, Matern S, Gartung C: Consequences of bile duct obstruction on intestinal expression and function of multidrug resistanceassociated protein 2. Gastroenterology. 2004, 126 (4): 10441053. 10.1053/j.gastro.2003.12.046.
 22.
Shampine L, Reichelt M: The matlab ode suite. SIAM journal on scientific computing. 1997, 18: 110.1137/S1064827594276424.
Acknowledgements
This work was supported partly by a GrantinAid for Scientific Research (S) (24229002) and GrantinAid for Scientific Research on Innovative Area Molecular Robotics (24104004) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.
Declarations
The publication costs for this article were funded by a GrantinAid for Scientific Research on Innovative Area Molecular Robotics from the MEXT, Japan.
This article has been published as part of BMC Systems Biology Volume 7 Supplement 3, 2013: Twelfth International Conference on Bioinformatics (InCoB2013): Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/7/S3.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
K.Y. wrote the manuscript, designed research, performed research, and analyzed data. K.M. wrote the manuscript and designed research. H.K. wrote the manuscript and designed research. A.K. wrote the manuscript, designed research, and analyzed data. All authors read and approved the final manuscript.
Rights and permissions
About this article
Published
DOI
Keywords
 Pharmacokinetics
 PBPK models
 Optimization