Estimation of feasible solution space using Cluster Newton Method: application to pharmacokinetic analysis of irinotecan with physiologically-based pharmacokinetic models
BMC Systems Biology volume 7, Article number: S3 (2013)
To facilitate new drug development, physiologically-based 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.
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 E5-1620 3.60GHz × 1 or Intel Core i7-870 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.
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.
Pharmacokinetics is a field of study that analyzes and predicts behaviors of drugs in organisms . 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 drug-drug interactions (DDIs) , pharmacogenetics , or disease states [3, 4], which can cause large inter-individual 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 .
Physiologically-based pharmacokinetic (PBPK) modeling and simulation are essential in understanding and predicting the above-mentioned, 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  or the combined effect of multiple inhibitors . Current draft guidance on DDI studies by the U.S. Food and Drug Administration  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, Gauss-Newton 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 . 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 , in which the pharmacokinetic properties of an anti-cancer drug, irinotecan (also known as CPT-11), differs a lot between a bile-duct cancer (BDC) patient and other cancer (OC) patients . After intravenous administration, irinotecan is metabolized by CYP3A4 or carboxylesterase 2 (CES2) to form APC, NPC, SN-38, and M4 [15, 16] (Figure 1). NPC and SN-38 are further metabolized by CES2 and UGT1A to form SN-38 and SN-38G, respectively. Organic anion transporting polypeptide 1B1 (OATP1B1) is involved in the hepatic uptake of SN-38, while OATP1B1 does not actively transport irinotecan or SN-38G, according to an in vitro study . SN-38G is said to be deconjugated to form SN-38 by β-glucuronidase in intestinal microflora . Possible involvement of enterohepatic circulation (EHC) in determining irinotecan pharmacokinetics  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 2p-t). 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 2f-j). Particularly, parameters got diverged after 2 iterations for dS of 0 to 0.2, or after 3 iterations for dS of 0.3 (Figure 2a-d). 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 SN-38 (CLbile, 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 (CLR) of irinotecan (parameter #6) in both the OC and the BDC groups, and ratio of biliary clearance into T-tube to biliary clearance into a transit compartment (CLbile,T / CLbile,transit, parameter #56) in the BDC group (differences in convergence will be discussed later). In particular, the estimated values of irinotecan CLR 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 CLR 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 concentration-time 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 E5-1620 3.60GHz × 1 or Intel Core i7-8702.93GHz × 1; see Methods section for details).
After the application of CNM, we simulated the accumulation time-profiles 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, CLR 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 time-profiles (Figure 5b). Interestingly, these parameters are mostly involved in hepatobiliary or intestinal pathways (#12-13, 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  or in the duodenums of patients with hepatic cholestasis , 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 r2 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 T-tube was increased (#56 increased), biliary clearances decreased (#12 and #14) to keep the amount of compounds excreted into the T-tubes. While these phenomena can be explained retrospectively, we cannot easily observe the actual parameter dependencies quantitatively without performing this kind of parameter estimation algorithms.
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.
Throughout thee study, we analyzed reported pharmacokinetic profiles of irinotecan and its metabolites . This report described the accumulation of irinotecan and its metabolites into urine and feces for two groups of patients, with bile-duct cancer (BDC, one patient) or with other cancers (OC, seven patient). Bile was also collected via biliary T-tube 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 CLbile,T and CLbile for irinotecan and the metabolites, since this parameter can be regarded as the ratio of bile flow between in T-tube and in bile duct after T-tube. In OC group, CLbile,T- was set to be equal to zero, since biliary T-tube 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 SN-38G to SN-38 by β-glucuronidase in intestinal microflora may affect the fecal elimination amount of SN-38G and SN-38, we combined the accumulation of these two compounds in objective values #7.
We solved the system of ODEs by using the MATLAB stiff ODE solver ODE15s . We performed all calculations with MATLAB using a desktop computer (CPU: Core i7-870 2.93GHz × 1, OS: Windows 7 SP1 32 bit, RAM: 4GB, MATLAB version: 8.0.0) or a workstation (CPU: XeonE5-1620 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 . 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 (Xb) and after (Xa) the parameter estimation process. In our new strategy, we calculated internally dividing point Xi with the ratio of dS:(1-dS), and applied the same inverse matrix to obtain new estimated parameters Xa'. Parameter sets for the next iteration were obtained by randomly selecting Xa or Xa' for each virtual sample.
After completing the estimations of parameters, the accumulation time-profiles were compared with the observed profiles, using sum of squares of log residuals (SSlog):
where Ae,simulated and Ae,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 (r2> 0.64) with at least one other parameter.
clearance from a rapid to a late equilibrium compartment
metabolic clearance of SN-38 by UGT to form SN-38G
Cluster Newton Method
kinetic constant from a late to a rapid equilibrium compartment
absorption rate constant
kinetic constant for the transit from large intestine to feces
kinetic constants for the transit from small intestine to large intestine
kinetic constant for the transit in bile compartments to small intestine
blood flow rate in liver
volume of a liver
apparent volume of a rapid equilibrium compartment.
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: McGraw-Hill, 17-39. 12th
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: McGraw-Hill, 145-168. 12th
Verbeeck R: Pharmacokinetics and dosage adjustment in patients with hepatic dysfunction. Eur J Clin Pharmacol. 2008, 64 (12): 1147-1161. 10.1007/s00228-008-0553-z.
Verbeeck R, Musuamba F: Pharmacokinetics and dosage adjustment in patients with renal dysfunction. Eur J Clin Pharmacol. 2009, 65 (8): 757-773. 10.1007/s00228-009-0678-8.
Konagaya A: Towards an In Silico Approach to Personalized Pharmacokinetics. Molecular Interactions. 2012, Edited by Meghea A. Rijeka, Croatia: Intech, 263-282.
Rowland M, Peck C, Tucker G: Physiologically-based pharmacokinetics in drug development and regulatory science. Annu Rev Pharmacol Toxicol. 2011, 51: 45-73. 10.1146/annurev-pharmtox-010510-100540.
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): 259-267. 10.1038/clpt.2010.298.
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): 743-761. 10.2174/138920010794328931.
Yoshida K, Maeda K, Sugiyama Y: Hepatic and intestinal drug transporters: prediction of pharmacokinetic effects caused by drug-drug interactions and genetic polymorphisms. Annu Rev Pharmacol Toxicol. 2013, 53: 581-612. 10.1146/annurev-pharmtox-011112-140309.
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): 347-351.
Draft Guidance/Guidance for Industry. Drug Interaction Studies--Study Design, Data Analysis, Implications for Dosing, and Labeling Recommendations. [http://www.fda.gov/downloads/Drugs/GuidanceComplianceRegulatoryInformation/Guidances/UCM292362.pdf]
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
Arikuma T, Yoshikawa S, Azuma R, Watanabe K, Matsumura K, Konagaya A: Drug interaction prediction using ontology-driven hypothetical assertion framework for pathway generation followed by numerical simulation. BMC bioinformatics. 2008, 9 (Suppl 6): S11-10.1186/1471-2105-9-S6-S11.
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 (CPT-11) following I.V. infusion of [(14)C]CPT-11 in cancer patients. Drug Metab lism and Disposition: The Biological Fate of Chemicals. 2000, 28 (4): 423-433.
Smith NF, Figg WD, Sparreboom A: Pharmacogenetics of irinotecan metabolism and transport: an update. Toxicol In Vitro. 2006, 20 (2): 163-175. 10.1016/j.tiv.2005.06.045.
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): 2604-2614. 10.1200/JCO.2008.20.6300.
Nozawa T, Minami H, Sugiura S, Tsuji A, Tamai I: Role of organic anion transporter OATP1B1 (OATP-C) in hepatic uptake of irinotecan and its active metabolite, 7-ethyl-10-hydroxycamptothecin: in vitro evidence and effect of single nucleotide polymorphisms. Drug metabolism and disposition: the biological fate of chemicals. 2005, 33 (3): 434-439.
Takasuna K, Hagiwara T, Hirohashi M, Kato M, Nomura M, Nagai E, Yokoi T, Kamataki T: Involvement of beta-glucuronidase in intestinal microflora in the intestinal toxicity of the antitumor camptothecin derivative irinotecan hydrochloride (CPT-11) in rats. Cancer Res. 1996, 56 (16): 3752-3757.
Younis IR, Malone S, Friedman HS, Schaaf LJ, Petros WP: Enterohepatic recirculation model of irinotecan (CPT-11) and metabolite pharmacokinetics in patients with glioma. Cancer Chemother Pharmacol. 2009, 63 (3): 517-524. 10.1007/s00280-008-0769-8.
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): 421-429. 10.1007/s00423-005-0564-5.
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 resistance-associated protein 2. Gastroenterology. 2004, 126 (4): 1044-1053. 10.1053/j.gastro.2003.12.046.
Shampine L, Reichelt M: The matlab ode suite. SIAM journal on scientific computing. 1997, 18: 1-10.1137/S1064827594276424.
This work was supported partly by a Grant-in-Aid for Scientific Research (S) (24229002) and Grant-in-Aid for Scientific Research on Innovative Area Molecular Robotics (24104004) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.
The publication costs for this article were funded by a Grant-in-Aid 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.
The authors declare that they have no competing interests.
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.
About this article
Cite this article
Yoshida, K., Maeda, K., Kusuhara, H. et al. Estimation of feasible solution space using Cluster Newton Method: application to pharmacokinetic analysis of irinotecan with physiologically-based pharmacokinetic models. BMC Syst Biol 7 (Suppl 3), S3 (2013). https://doi.org/10.1186/1752-0509-7-S3-S3