A Combined Experimental and Mathematical Approach for Molecular-based Optimization of Irinotecan Circadian Delivery

Circadian timing largely modifies efficacy and toxicity of many anticancer drugs. Recent findings suggest that optimal circadian delivery patterns depend on the patient genetic background. We present here a combined experimental and mathematical approach for the design of chronomodulated administration schedules tailored to the patient molecular profile. As a proof of concept we optimized exposure of Caco-2 colon cancer cells to irinotecan (CPT11), a cytotoxic drug approved for the treatment of colorectal cancer. CPT11 was bioactivated into SN38 and its efflux was mediated by ATP-Binding-Cassette (ABC) transporters in Caco-2 cells. After cell synchronization with a serum shock defining Circadian Time (CT) 0, circadian rhythms with a period of 26 h 50 (SD 63 min) were observed in the mRNA expression of clock genes REV-ERBα, PER2, BMAL1, the drug target topoisomerase 1 (TOP1), the activation enzyme carboxylesterase 2 (CES2), the deactivation enzyme UDP-glucuronosyltransferase 1, polypeptide A1 (UGT1A1), and efflux transporters ABCB1, ABCC1, ABCC2 and ABCG2. DNA-bound TOP1 protein amount in presence of CPT11, a marker of the drug PD, also displayed circadian variations. A mathematical model of CPT11 molecular pharmacokinetics-pharmacodynamics (PK-PD) was designed and fitted to experimental data. It predicted that CPT11 bioactivation was the main determinant of CPT11 PD circadian rhythm. We then adopted the therapeutics strategy of maximizing efficacy in non-synchronized cells, considered as cancer cells, under a constraint of maximum toxicity in synchronized cells, representing healthy ones. We considered exposure schemes in the form of an initial concentration of CPT11 given at a particular CT, over a duration ranging from 1 to 27 h. For any dose of CPT11, optimal exposure durations varied from 3h40 to 7h10. Optimal schemes started between CT2h10 and CT2h30, a time interval corresponding to 1h30 to 1h50 before the nadir of CPT11 bioactivation rhythm in healthy cells.


The lactone-carboxylate model
A mathematical model was built to study the equilibrium between the lactone and the carboxylate forms of CPT11 with respect to the pH of the solution. The temperature was assumed to be constant and equal to 37 o C. We considered that reactions occurred in a closed system containing a buffer solution which kept the pH constant. Two inverse chemical reactions occurred: To model this problem we defined two variables: CP T l and CP T c respectively CPT11 lactone and carboxylate concentrations expressed in µM. As the system was closed, the total quantity of CPT-11 was conserved over time: CP T c + CP T l = CP T tot . Applying the law of mass action and using this conservation law gave the following equation for CP T l kinetics: dCP T l dt = −(k l 10 −pH + k c 10 pH−14 )CP T l + k l 10 −pH CP T tot As the pH remained constant because of the buffer solution we obtained a simple ODE that was solved analytically: with A = k l 10 −pH + k c 10 pH−14 and B = k l 10 −pH CP T tot At equilibrium the following holds: We used experimental data in Table I and Table II of [54] to estimate k c and k l by a least square approach using the CMAES algorithm for minimization. Estimated parameter values were: k l = 19300h −1 and k c = 8.54.10 −6 h −1 ( Figure S1, panel 1 and 2). We then validated our model by comparing it to the set of data contained in Figure 6 in [54] ( Figure S1, panel 3).
The mother solution of CPT11 is at pH 4.4. In our experimental conditions, cells were cultured in an extracellular medium at pH 7.8. We investigated the needed time to reach the lactone/carboxylate equilibrium after the mother solution was added to the culture medium. Thus we searched for t such that: CP T l (t) − CP T * l ≤ eps where eps was a tolerance parameter. At pH=4.4, all CPT-11 molecules were under their lactone form: [CP T l ](0) = [CP T total ] = 115µM . We ended with the analytical formula: Replacing numerical values in this formula and setting eps to 10 −3 gave : t ≥ 2.9 . Therefore CPT11 was added to the medium at least 3 hours before using it for cell exposure so that the lactone/carboxylate equilibrium was reached.  Figure S1: The lactone/carboxylate model. Panel 1: CPT11 carboxylate transformation into CPT11 lactone at pH 4 ( ), 5 (•) and 6 ( ). Panel 2: CPT11 lactone transformation into CPT11 carboxylate at pH 6 ( ), 7 ( ), 8 (•) and 9 ( ). Panel 3: Percentage of CPT11 lactone over CPT11 total quantity at equilibrium with respect to pH. See [54] for details about experimental results.

Normalization of HPLC measurements
CPT11 and SN38 extra-and intracellular concentrations were measured by HPLC. The direct results obtained from HPLC depended on the number of cells present in the Petri Dish and were therefore normalized as if all Petri dishes contained one million of cells. Control Petri dishes after 2, 12, 24 and 42 h of exposure were trypsinized and their cell number counted using the Malassez cell. Then we interpolated to obtain the number of cells for the other exposure durations.
HPLC measurements were normalized as follows. Let denote n the quantity of moles of drug and N c the measured number of cells in a Petri dish. () meas stood for the measured value from HPLC, () 1M for the normalized value to one million cells. For intracellular concentration, the following holds: (n in ) 1M = (n in )meas

Nc
. Therefore we get: For extracellular concentration normalization, we used the conservation law on the drug total quantity to compute (n out ) 1M : Finally we got:

Parameter estimation for circadian gene expressions in synchronized Caco-2 cells
The circadian expression of ten genes were measured in Caco-2 cells, namely three clock genes REV-ERBα, PER2, and BMAL1 ; three metabolism and target genes TOP1, UGT1A1, CES2 ; and four ABC transporters ABCB1, ABCC1, ABCC2 and ABCG2. We assumed that those genes oscillated with the same period T. Experimental data of mRNA expression were fitted to equation (1) (cf. Materials and Methods). In addition to the common period T, four parameters have to be estimated for each gene: R, λ, P and S. The 41 parameters were estimated simultaneously using a Bootstrap approach. Briefly, data were first normalized so that each mRNA expression had an average value equal to 0 and a standard deviation equal to 1. Then five hundreds datasets were generated by assuming that mRNA values at each time point followed a Gaussian law of average and standard deviation the ones of the corresponding data time point. For each generated dataset, parameters were estimated using a least square approach. The minimization of the cost function was performed using the Matlab function fmincon in which interval boundaries were specified for all parameters. Only few constraints of intervals were active: for PER2, ABCC1, ABCC2, ABCG2 the upper bound on λ was active (λ < 0.085); for TOP1 and ABCB1, the lower bound on R was active (R > −1.5), for UGT1A1 the upper bound on λ was active (λ < 0.07). Finally the mean value and standard deviation of the hundred parameter sets were computed and rescaled to obtain non-normalized values.

Circadian rhythms of mRNA amounts : variability between experiments
The mRNA amounts of CES2, TOP1, UGT and REV − ERBα were quantified in four independent experiments ( Figure S2). CES2 circadian rhythm was consistent in experiment 1, 2 and 4 with an amplitude of 15%, 29%, and 28% of their mean value respectively. On the contrary, no obvious circadian pattern of CES2 expression was noticed in experiment 3. The overshoot in REV − ERBα expression during the first cycle in experiments 1, 2 and 4 was not observed in experiment 3.

Parameter estimation for the CPT11 molecular PK-PD model
This section aims at giving more details about how the parameters of the CPT11 molecular PK-PD model were estimated. The first step consisted in determining reasonable search intervals for each parameter either based on literature or on unpublished data on our Caco-2 cell line. Then the second step used a bootstrap approach to fit the model to experimental data shown in Figures 2 and 4.

Determining a search interval for each parameter
Search intervals for parameters were determined as follows: • CPT11 Uptake: CPT11 uptake was assumed to occur passively at the rate k upCP T .
It was studied in Caco-2 cell culture by measuring CPT11 intracellular concentration after a 10-minute exposure at concentrations from 10 to 100µM (data not shown). We assumed that in the first ten minutes CPT11 efflux could be neglected as the drug had only time to enter the cell ([20]). This allowed the computation of an approximation of the uptake parameter k upCP T = 5.9h −1 . A deviation of ±50% was allowed in parameter estimation.
• SN38 Uptake: SN38 uptake was assumed to occur passively at the rate k upSN . As no measurement of SN38 uptake was available in our Caco-2 cells, we used experimental data from [20] to determine SN38 uptake parameter and found: k upSN = 15.3h −1 . A deviation of ±100% was allowed in parameter estimation.
• CPT11 Efflux: CPT11 efflux followed Michaelis-Menten kinetics with parameters V ef f CP T and K ef f CP T . We used data from Figure 2 and unpublished data on CPT11 accumulation in Caco-2 cells at concentrations 10µM, 40µM to compute a guess for the efflux parameters. Indeed once the parameter for CPT11 uptake was approximately known, it was possible to deduce CPT11 efflux parameters from equation (2)  • SN38 efflux: SN38 efflux followed Michaelis-Menten kinetics with parameters V ef f SN and K ef f SN . Those parameters were directly inferred from experimental data of Figure 2 and 4.
• SN38 deactivation into SN38G by UGT1As: the glucuronidation of SN38 into SN38G catalyzed by enzymes UGT1As followed Michaelis-Menten kinetics with parameters • Formation and dissociation of DNA/TOP1 reversible complexes: we first determined an approximation of TOP1 total concentration by assuming that all SN38 molecules present in the intracellular medium in Figure 2 were trapped in complexes with TOP1. SN38 intracellular concentration was equal to 0.01µM in average and data from Figure 3 suggested that around 40% of TOP1 total quantity was trapped in complexes with the drug. Therefore TOP1 total concentration could be estimated to 0.025µM and was searched within the interval [0 0.1]. In the absence of drug, around one third of TOP1 was linked to the DNA in Caco-2 cells (data not shown) which led to the following initial values for the formation (k f 1 ) and dissociation (k d1 ) rates of DNA/TOP1 complexes: k f 1 = 2, k d1 = 300. A deviation of ±100% was allowed for those two parameters.
• Formation and dissociation of DNA/TOP1/SN38 reversible complexes: the association (k f 2 ) and dissociation (k r2 ) rates of SN38 with DNA/TOP1 complexes were directly inferred from experimental data of Figure 2 and 4.
• Formation of irreversible complexes: the irreversible complexes come from collision of reversible complexes with replication or transcription mechanisms which occur at the rate k Irr . We assumed that this rate was independent of time meaning that the amount and speed of mechanisms going along the DNA was constant. This assumption was realistic in quiescent cells but would be untrue in proliferating cells in which the number of replication forks strongly increases during the S-phase. No prior information was available concerning the parameter k Irr which was directly fitted to experimental data of Figures 2 and 4.
• DNA total quantity: entry sites where TOP1 can bind to the DNA were assumed to occur every k entry pairs of bases. Here we defined an entry site as k entry pairs of bases. When a molecule of TOP1 associated with an entry site, it was consumed and then became a Compl molecule. As a cell contains approximately 3.2.10 9 pairs of base, the number of moles of entry sites was 3.2.10 9 N A kentry where N A is the Avogadro number. Thus the concentration of free entry sites in a cell with a volume equal to V cell = 8.10 −12 L is approximately: DN A tot = 3.2.10 9 N A kentryV cell M . k entry was assumed to be between 1 and 100. Therefore DN A tot was searched within the interval [5 500].
• Protein amount circadian amplitudes: circadian amplitudes for CES, U GT , ABC CP T and ABC SN were searched in the interval [0, 1] so that the protein formation term in equation (12) remained positive.

A bootstrap approach for the second step of parameter estimation
The second step of our parameter estimation consisted in the simultaneous evaluation of the 21 parameters which were searched in the intervals previously determined. The mathematical model was fitted to experimental data on CPT11 and SN38 pharmacokinetics in the presence or absence of verapamil ( Figure 2) and CPT11 chrono-pharmacodynamics measured at CT14 and CT28 (Figure 4). Parameters were estimated using a Bootstrap approach. Briefly, fifty datasets were generated from the original one, by assuming that concentrations at each time point followed a Gaussian law of mean and standard deviation the ones measured at the corresponding data time point. For each generated dataset, parameters were estimated using a least square approach in which the minimization task was performed by the CMAES function. This algorithm was preferred to Matlab function fmincon because it handled better local minima. Thus we got 50 parameter sets from which we derived the mean value and standard deviation for each parameter (Table 2).
In order to start the 50 least-square estimations with initial values close to the optimal ones, we firstly determined a preliminary fit. We fitted averaged experimental values of Figures 2 and 4 by a least square approach which was iterated for different initial values of parameters. We thus computed the parameter set which fitted the best the average experimental data and used it as initial values for the bootstrap approach.