Differential Response to Cytotoxic Drugs Explains the Dynamics of Leukemic Cell Death: Insights from Experiments and Mathematical Modeling

Abstract

This study presents a framework whereby cancer chemotherapy could be improved through collaboration between mathematicians and experimentalists. Following on from our recently published model, we use A20 murine leukemic cells transfected with monomeric red fluorescent proteins cells (mCherry) to compare the simulated and experimental cytotoxicity of two Federal Drug Administration (FDA)-approved anticancer drugs, Cytarabine (Cyt) and Ibrutinib (Ibr) in an in vitro model system of Chronic Lymphocytic Leukemia (CLL). Maximum growth inhibition with Cyt (95%) was reached at an 8-fold lower drug concentration (6.25 μM) than for Ibr (97%, 50 μM). For the proposed ordinary differential equations (ODE) model, a multistep strategy was used to estimate the parameters relevant to the analysis of in vitro experiments testing the effects of different drug concentrations. The simulation results demonstrate that our model correctly predicts the effects of drugs on leukemic cells. To assess the closeness of the fit between the simulations and experimental data, RMSEs for both drugs were calculated (both RMSEs < 0.1). The numerical solutions of the model show a symmetrical dynamical evolution for two drugs with different modes of action. Simulations of the combinatorial effect of Cyt and Ibr showed that their synergism enhanced the cytotoxic effect by 40%. We suggest that this model could predict a more personalized drug dose based on the growth rate of an individual’s cancer cells.

1. Introduction

Chronic Lymphocytic Leukemia (CLL) is characterized by the accumulation of clonal B cells in peripheral blood [1]. Currently, patients with CLL are mostly treated with drugs that target the activity of mutated proteins critical to the survival of CLL cells, such as Bruton’s Tyrosine Kinase (BTK), B-cell lymphoma 2 (Bcl-2) and phosphoinositide 3-kinase (PI3K). Inhibitors of BTK, such as Ibrutinib (Ibr) or more recently acalabrutinib, are highly effective, even in poor-risk or chemo-refractory patients [2,3]. In addition, combinations of drugs, especially those with different mechanisms of action, may be effective in CLL cases [4,5]. However, while these drugs have led to significant improvements in the overall survival of most subgroups of patients with CLL, the disease remains incurable and novel approaches are required that will lead to more effective therapies [6]. Unfortunately, there has been little published collaboration between mathematicians and experimentalists to devise methods that might lead to improved treatment protocols for patients with CLL.

We suggest that numerical simulations of cancer chemotherapies can provide a means for researchers and clinicians to explore the possibility of predicting efficient, personalized treatment protocols based on mathematical parameters of cancer development in an individual. Ordinary differential equation (ODE)-based models have been used to define the modifications of continuous variables in mathematical models and have already been applied to several types of leukemia and lymphoma, including chronic myeloid leukemia [7], acute myeloid leukemia [8], T-cell lymphoma and periodic hematological diseases [9,10]. However, there are still no validated models for CLL. Experimental and clinical data are essential to the development and validation of these models [11].

In a previous study, we developed an ODE model based on a [12,13], to describe the drug cytotoxicity of three chemotherapeutic drugs, Chlorambucil, Melphalan, and Cytarabine (Cyt) for murine leukemic cells [14]. Cyt, which was the most toxic of the three drugs, is an analogue of cytidine, whose incorporation into DNA leads to interruption of DNA replication and induction of DNA damage repair pathways and is used in the treatment of several types of cancer [15]. Here, we extend the validation of the model by comparing the in vitro results of Cyt, with that of Ibr. Ibr covalently binds to cysteine-481 in the BTK enzyme, thus inhibiting several essential intracellular biochemical pathways that influence tumor cell proliferation and migration [16]. As these two drugs have different mechanisms of action, the availability of computational ODE models based on experimental data will assist in developing rational approaches to combination therapies that may overcome the outgrowth of drug resistant cancer cell clones, commonly the outcome of mono-drug therapies [17,18]. As an in vitro experimental model system, we used A20 murine leukemic cells transfected with monomeric red fluorescent proteins cells (mCherry) in order to compare the cytotoxic efficacy between Ibr and Cyt. The results validate that our model can be used as a generic tool to explore the cytotoxic activity of various drugs with different doses by proper adjustment of the values for certain parameters.

2. Materials and Methods

In order to predict the efficacy of an anticancer drug, we first developed a mathematical model. To validate the applicability of this model, experiments were then carried out in vitro using two chemotherapeutic anticancer drugs. The first drug, Cyt, was chosen because our previous study demonstrated its superior cytotoxicity over Chlorambucil and Melphalan for A20 leukemic cells. As Cyt is not used clinically to treat CLL, we compared Cyt to the cytotoxicity of Ibr, an FDA-approved drug for front-line therapy of CLL. Finally, we compared the simulation model with the experimental data. The simplified scheme below (Figure 1) illustrates the whole process.

2.1. Formulation of the Model

Based on our previous study [14], we consider a model that uses three critical variables to describe the in vitro experimental system, live leukemic cancer cells—dead cancer cells, and drug molecules, denoted by A, $A_d$ and C, respectively.

$$\{\begin{array}{}\hfill \mathrm{(1)}& \frac{dA}{dt}=rA\left(1-\frac{A}{K}\right)-{\mu }_{A}A{A}_d-\frac{{\mu }_{AC}AC}{a+C},\hfill \mathrm{(2)}& \frac{d{A}_d}{dt}={\mu }_{A}A{A}_d-d{A}_d+\frac{{\mu }_{AC}AC}{a+C},\hfill \mathrm{(3)}& \frac{dC}{dt}=-{\mu }_{C}C-\frac{{\mu }_{CA}CA}{a+C},\end{array}$$

${\displaystyle \frac{dA}{dt}}$ is the dynamic of living cancer cells. It is made up of three components: one of these is positive and relates to the logistic cancer growth which is characterized by the coefficient r and is limited by the maximal tumor cell number, K; a negative term, corresponding to living cells to become dead with the rate ${\mu }_A$. The last negative term represents the log-kill hypothesis [12], with a Michaelis–Menten drug saturation response [19], $a+C$; ${\mu }_{AC}$ is the death rate resulting from the action of the drug on cancer cells. The parameter ${\mu }_{AC}$ changes depending on the particular drug and the dose.

${\displaystyle \frac{d{A}_d}{dt}}$ is the dynamic of dead A20 cells. The first and the last terms are positive, represent the death of cancer cells with a rate coefficient of ${\mu }_A$ due to apoptosis or necrosis and depend on living cancer cells competing for survival (oxygen consumption and nutrition) in an enclosed space, and those that changed from living to dead under the influence of medicine; the negative term corresponds to the dissolution of dead cells with the rate d.

${\displaystyle \frac{dC}{dt}}$ is the first-order pharmacokinetics of a drug [20]. $C\left(0\right)$ is a constant value and depends on the drug’s dose, since it was given only once at the beginning of the experiment. ${\mu }_C$ is the deactivation rate calculated by formula ${\mu }_C=\frac{ln\left(2\right)}{t_{1/2}}$, where $t_{1/2}$ is the in vitro elimination half-life, 1–3 h (biphasic) for Cyt and 4–6 h for Ibr (www.drugbank.ca, accessed on 17 June 2022). Assuming that there are 10–100 drug molecules attacking each cancer cell, we decided the parameter ${\mu }_{CA}$ should be ten times more than ${\mu }_{AC}$. The parameter a represents the drug concentration which produces 50% of the maximum activity of the drug in each cell population [12].

2.2. Estimation of the Parameters

The model (Equations (1)–(3)) contains nine parameters that should be determined. Some of the parameter’s values we took from the publicly available resources, but their values may vary in different types of cancers. Other values we obtained from the experiments and due data fitting to estimate the unknown parameters.

Importantly, all the model’s parameters’ values must be positive in order to describe a realistic biological dynamic. We specify the initial conditions at t = 0, as:

• $A\left(0\right)=5\times {10}^4$ (cells/mL)-the initial number of A20 mCherry cells;
• $A_d\left(0\right)=2500$ (cells/mL)-the initial number of dead A20 mCherry cells (cell cultures commonly consists of at least 5% of dead cells);
• $C\left(0\right)=$ dose ($\mathsf{\mu }$M) $\times 6\times {10}^{11}$ (number of drug molecules/mL)-the dose concentration of Ibr or Cyt (this number may vary depending on the drug, but not significantly since both drugs are related to the same type of small molecules).

Number of drug molecules was calculated using the expression: $\frac{m\times Na}{M},$ where

• m = the mass of drug in kg,
• $Na$ = avogadro number = $6.022\times {10}^{23}$ (constant),
• M = the molar mass of drug (Ibr 440.5 g/mol; Cyt 243.217 g/mol).

Thus, for 50 µM of Ibr, the Mass = Concentration 4.6 mM × Volume 10.85 µL/mL × Molecular Weight 440.5 g/mol = 22 µg = $2.2\times {10}^{-8}$ kg, the number of Ibr molecules will be:

$$\frac{2.2\times {10}^{-8}\times 6.022\times {10}^{23}}{440.5}=3\times {10}^{13}=50\times 6\times {10}^{11};$$

for 6.25 µM of Cyt, the Mass = Concentration 4 mM × Volume 1.56 µL/mL × Molecular Weight 243.22 g/mol = 1.52 µg = $1.52\times {10}^{-9}$ kg, the number of Cyt molecules will be:

$$\frac{1.52\times {10}^{-9}\times 6.022\times {10}^{23}}{243.22}=3.75\times {10}^{12}=6.25\times 6\times {10}^{11}.$$

The parameters used in this study are summarized in

Table 1. Table of parameters based on experimental results related to the model of Equations (1)–(3).

Parameter	Physical Interpretation (Units)	Estimated Value	Reference
$A\left(0\right)$	The initial number of A20 [cells/mL]	$5\times {10}^4$	Experimental data
$A_d\left(0\right)$	The initial number of dead A20 [cells/mL]	2500	Experimental data
$C\left(0\right)$	Number of drug [molecules/mL]	dose ($\mathsf{\mu }$M) $\times 6\times {10}^{11}$	Experimental data
t	Time of cell culture [h]	0–72	Experimental data
r	A20 growth rate [h${}^{-1}$]	0.07	Experimental data
K	Maximal tumor cell population [cells/mL]	$4\times {10}^6$	Experimental data
${\mu }_A$	Living cells become dead [h${}^{-1}$]	$3.7\times {10}^{-8}$	Simulation
a	Drug dose that produces 50% maximum effect [mL]	$2\times {10}^3$	[12]
${\mu }_{AC}$	Cytotoxicity rate in the presence of drug [h${}^{-1}$]	see Table A1 and Table A2	Simulation
${\mu }_{CA}$	Deactivation rate of drug due to killing of A20 cells [h${}^{-1}$]	${\mu }_{AC}\times 10$	Simulation
${\mu }_C$	Chemical deactivation rate of drug [h${}^{-1}$]	0.231-Cyt; 0.116-Ibr	[12]
d	Dissolution rate of dead A20 cells [h${}^{-1}$]	0.017	Simulation

.

2.3. Cells and Reagents

Murine A20 leukemic cells transfected with monomeric red fluorescent proteins cells (a kind gift from Prof. Michael Sherman, Ariel University) were cultured in RPMI 1640 (Thermo Fischer Scientific, Waltham, MA, USA) to which was added 10% Fetal Bovine Serum (FBS) (Thermo Fischer), 1% L-Glutamine and 0.33% Pen-Strep solution and maintained at 37° and 5% CO₂. Cell viability and growth were measured with the XTT-based Cell Proliferation Kit (Biological Industries, Bet Haemek, Israel) according to the manufacturer’s instructions. At the end of the culture period, the medium was replaced with fresh culture medium containing XTT reagent to measure cellular metabolism as a correlate of cell viability. After 2 h of cell incubation, absorbance in the wells was measured at 450 nm and subtracted from the reference absorbance at 630 nm. As a background control, the culture medium was used.

2.4. Drug Cytotoxicity Assay

Drug cytotoxicity was determined by culturing A20 mCherry cells at 37°, 5% CO₂ in 96-well microplates (Nunclon) at an initial concentration of $5\times {10}^4$/well for 72 h in culture medium (as described above) containing either Cyt or Ibr at concentrations ranging from 0 to 50 $\mathsf{\mu }$M. At the end of the culture period, cell viability was measured with the XTT reagent as described above.

3. Results

The outcome of our present study is the development of a tool, which is able:

• To simulate the impact of Cyt and Ibr drugs on killing A20 leukemic cells in silico;
• To predict with a high level of accuracy of the cytotoxic efficacy of Cyt and Ibr drugs for high doses in comparison with the results in vitro experiments.

Below we disclose documentation about the experiments, introduce simulation capabilities of Cyt and Ibr drugs, and predict the improvement of cytotoxic efficacy using combined therapy in silico.

3.1. Validation of the In Vitro A20 mCherry Cell Drug Cytotoxicity Dynamic Model

We first checked that A20 mCherry cells have the same growth rate as a surrogate for CLL cells, $r=0.07$ (Table 1 in Guzev et al. [14] and Table 1 in the current work). We then tested in vitro the cytotoxicity of Cyt and Ibr against A20 mCherry cells and compared the degree of A20 mCherry growth inhibition received from numerical simulation (Figure 2) with output obtained by in vitro experiments. The results are shown in Figure 3.

Using Equations (1)–(3), we could numerically simulate the effect of Cyt or Ibr on A20 mCherry cells by substituting necessary parameters with the values from the Table 1. Fourth-order adaptive step Runge–Kutta integration applied with the ODE45 subroutine of MATLAB [21] was used to carry out computer simulations. The finite-distance analysis error of the proposed model is presented in Appendix B.

After 72 h the concentration of A20 mCherry growth without drug increased to $A\left(72\right)$ = 2,636,950 cells/mL. The additional 50 µM of Ibr (A) reduced the growth to $A\left(72\right)$ = 50,953 cells/mL, which equals 98.1% growth inhibition and 6.25 µM of Cyt induced $94.8\%$ growth inhibition (B). For the calculated data points for all drug concentrations, see Appendix A (

Table A1. Simulated effect of Ibrutinib on A20 growth parameters.

Concentration of Ibr (µM)	Parameter ${\mathit{\mu }}_{\mathbf{AC}}$	Number of Cells at 72nd Hour (Cells/mL)	Cell Growth Inhibition (%)
0	-	2,636,950	0
50	0.066	50,953	98.1
25	0.046	229,415	91.3
$12.5$	0.032	424,548	83.9
$6.25$	0.022	743,620	71.8
$3.125$	0.015	1,152,347	56.3
$1.5625$	0.01	1,579,533	40.1
$0.78$	0.007	1,914,425	27.4
$0.39$	0.005	2,244,044	14.9
$0.195$	0.003	2,539,382	3.7

and

Table A2. Simulated effect of Cytarabine on A20 growth parameters.

Concentration of Cyt (µM)	Parameter ${\mathit{\mu }}_{\mathbf{AC}}$	Number of Cells at 72nd Hour (Cells/mL)	Cell Growth Inhibition (%)
0	-	2,636,950	0
$6.25$	0.051	137,121	94.8
$3.125$	0.036	268,969	89.8
$1.5625$	0.025	466,740	82.3
$0.78$	0.017	725,161	72.5
$0.39$	0.012	1,020,500	61.3
$0.195$	0.008	1,326,386	49.7
$0.098$	0.006	1,608,540	39
$0.049$	0.004	1,853,776	29.7
$0.024$	0.003	2,051,547	22.2
$0.012$	0.002	2,209,764	16.2
$0.006$	0.0014	2,328,427	11.7

).

As one can see from Figure 3, Ibr was more effective with a higher dose, 50 µM-97.5% of growth inhibition. On the other hand, Cyt maintained a high level, 95.2% of growth inhibition under a much lower dose concentration, 6.25 µM of the drug.

Numerical simulations corresponded with the trial information for the high portions of both medications: under 50 and 25 µM of Ibr, as seen in

Table 2. Comparison of cell growth inhibition (%) between experimental data and numerical simulations.

Concentration of Drug (µM)	Cytarabine		Ibrutinib
	Exp	Sim	Exp	Sim
50			97.5	98.1
25			91.5	91.3
12.5			74	83.9
6.25	95.2	94.8	42.5	71.8
3.125	93.5	89.8	29.1	56.3
1.5625	88.7	82.3	20.3	40.1
0.78	87	72.5	18.8	27.4
0.39	85.2	61.3	10.8	14.9
0.195	80.5	49.7	4.3	3.7
0.098	66.25	39
0.049	29.5	29.7
0.024	18.9	22.2
0.012	14.5	16.2
0.006	7.7	11.7
RMSE	0.018		0.032
MAPE	0.198		0.412

, the results from experiments and simulations were the same, around 98% and 91%, respectively, and under 6.25 µM of Cyt, around 95%. For the moderate dosages, there had all the earmarks of being a slight deviation between the two arrangements of information. This deviation was imitated since we attempted to keep up with consistency in the lessening of the parameter ${\mu }_{AC}$ as per the abatement in the medication portion. This brought in a discrepancy of 30% for each medication. Nonetheless, it is vital to take note that keeping up with this consistency isn’t required, and it is feasible to compare the experimental and simulation results with outright precision by picking the fitting parameter ${\mu }_{AC}$.

3.2. Goodness of Fitting Evolution

To assess the goodness of fit between the simulation and experimental data, the root mean square of errors (RMSE) and the mean absolute percentage error (MAPE) were calculated for Cyt and Ibr (Table 2). First of all, we have to normalize our data in order to measure how much error there is between the experimental data and numerical simulations:

$$\sigma {\left(z\right)}_i=\frac{z_i}{{\sum }_i^n{z}_i},$$

then, we calculate RMSE using the formula:

$$RMSE=\sqrt{\frac{{\sum }_i^n{(\sigma {\left(Sim\right)}_i-\sigma {\left(Exp\right)}_i)}^2}{n}},$$

where $Sim$ is a predicted value; $Exp$ is an observed value and n is the number of observations.

For calculation, MAPE used the formula:

$$MAPE=\frac{1}{n}\sum _i^n\frac{(Ex{p}_i-Si{m}_i)}{Ex{p}_i},$$

The result of this can be seen in Table 2.

The correlation held for both drugs (all RMSEs are < 0.1) but was more apparent for Cyt (for Ibr, 0.032; for Cyt, 0.018). The prediction accuracy for Cyt is about 80% (MAPE = 0.198) and for Ibr is about 60% (MAPE = 0.412).

3.3. Prediction of the Synergistic Effect of Drugs

In order to test the synergistic effects of Ibr combined with Cyt, we first ran a simulation using $L{D}_{50}$ concentrations of each drug (

Table 3. Simulated effect of the synergism.

Concentration of Drugs (µM)	Parameter ${\mathit{\mu }}_{\mathbf{AC}}$	Number of Cells at 72nd Hour (Cells/mL)	Cell Growth Inhibition (%)
0	-	2,636,950	0
Ibr 3.125	0.015	1,152,347	56.3
Cyt 0.195	0.008	1,326,386	49.7
Ibr 3.125 + Cyt 0.195	0.023	702,330	73.4
Ibr 6.25	0.022	743,620	71.8
Cyt 0.39	0.012	1,020,500	61.3
Ibr 6.25 + Cyt 0.39	0.034	375,766	85.7

), 3.125 µM of Ibr and 0.195 µM of Cyt and adjusted the calculated parameters into the model (Equations (1)–(3)): $C\left(0\right)=(3.125+0.195)\times 6\times {10}^{11}=3.32\times 6\times {10}^{11}$ drug molecules. Under these conditions, Ibr makes up about 94% of the molecules, and Cyt about 6%. Hence, ${\mu }_C$ will be: $(0.116\times 0.94)+(0.231\times 0.06)=0.124$; according to the Table A1 and Table A2, ${\mu }_{AC}$ will be $0.015+0.008=0.023$. This simulation predicted 73.4% of cell growth inhibition, which represents a 40% increase in killing efficiency. Doubling the dose of each drug ( 6.25 µM of Ibr (71.8% killing) and 0.39 µM of Cyt (61.3% killing)) resulted in a simulated increase of cell killing to 85.7% with relevant calculated parameters, ${\mu }_C=0.124$ and ${\mu }_{AC}=0.034$.

4. Discussion

In this study, we demonstrated experimentally in vitro the cytotoxic efficacy of Ibr compared to that of Cyt in killing A20 leukemic cells. This allowed us to build on our previous model [14] which describes the cytotoxic efficacy of different CLL drugs in killing A20 leukemic cells. In order to validate Equations (1)–(3) of the model, we carried out numerical simulations in accordance with the initial conditions (Table 1) taken from experiments (Figure 2). We found that Cyt (Figure 2B) is a more effective drug than Ibr (Figure 2A), since it requires a lower dose concentration (6.25 µM of Cyt against 50 µM of Ibr) to achieve almost the same maximum growth inhibition (about 95%) of A20 mCherry cells. As can be seen in Table 2, the median lethal dose, $L{D}_{50}$ for Ibr is between 12.5 and 6.25 µM, while the growth cell inhibition with 6.25 µM of Cyt is 93.5%, and the required $L{D}_{50}$ of Cyt is less than 0.098 µM, which Ibr does not react with.

Our simulation and experimental data for the high doses of both drugs were similar (50 and 25 µM of Ibr, around 98% and 91%, respectively, and for 6.25 µM of Cyt, around 95% (Figure 3)). For the moderate dosages, a small deviation between the two sets of data was found which was due to the maintenance of consistency in the decrease in parameter ${\mu }_{AC}$, correlating to the decrease in drug dose. However, it is possible to obtain the experimental and simulation results with absolute accuracy by selecting the appropriate parameter ${\mu }_{AC}$ for the drug used.

Using two different statistical analyses, the calculations of the differences among values displayed a strong symmetry between the numerical simulations and experimental results: for Cyt, RMSE = 0.018; MAPE = 0.198 and for Ibr, RMSE = 0.032; MAPE = 0.412 (Table 2). The prediction accuracy for Cyt is about 80% and for Ibr is about 60%. It should be noted that these values represent the accuracy across the entire calibration curve. However, regarding the effective drug doses in vitro the prediction was 100% (for 50 and 25 µM of Ibr and for 6.25 µM of Cyt).

The advantage of mathematical modeling is the prediction of the synergistic effect of drugs before conducting experiments. As these two drugs have different modes of action, their combination may lead to improved therapy for CLL patients. To test the synergistic effects of Ibr in combination with Cyt, we performed simulations with a combination of these drugs (Figure 4). The synergism of drugs, manifested in appropriate doses of $L{D}_{50}$, increases the effective cytotoxicity by another 40% as see in Table 3.

This in vitro model provides important information for drug development. First, it allows the assessment of a critical factor in drug therapy – the ability of the compound to directly influence cell viability. As we demonstrate, important variables of this interaction can be assessed. Another practical value of the model is to highlight additional applications of drugs. Our results predict that Cyt is even more cytotoxic than Ibr for leukemic cells. Even though Cyt is currently not used in the treatment of this cancer, our results suggest that this possibility should be investigated further. The combination of two and more drugs may lead to improved therapy for CLL patients and this can be verified by mathematical modeling.

CLL remains an incurable disease, and new approaches to improving therapy are needed. Mathematical models such as the one described in this study which allow quantitative simulations of cancer chemotherapies based on experimentally validated modeling should provide an opportunity to predict personalized drug therapy for individual patients and establish more efficient treatment protocols.

5. Conclusions

We constructed a mathematical model (Equations (1)–(3)) describing the cytotoxic efficiency of several drugs for leukemic cells. Surprisingly, we found that Cyt was a more effective drug than those currently used in the treatment of CLL.

The model is flexible and can be applied to various cancer cells, drugs and doses by using appropriate values for r, ${\mu }_{AC}$ and d.

The model is an auxiliary tool for predicting the most effective drug in the treatment of CLL.