Experimental Validation of a Mathematical Model to Describe the Drug Cytotoxicity of Leukemic Cells

: Chlorambucil (Chl), Melphalan (Mel), and Cytarabine (Cyt) are recognized drugs used in the chemotherapy of patients with advanced Chronic Lymphocytic Leukemia (CLL). The optimal treatment schedule and timing of Chl, Mel, and Cyt administration remains unknown and has traditionally been decided empirically and independently of preclinical in vitro efﬁcacy studies. As a ﬁrst step toward mathematical prediction of in vivo drug efﬁcacy from in vitro cytotoxicity studies, we used murine A20 leukemic cells as a test case of CLL. We ﬁrst found that logistic growth best described the proliferation of the cells in vitro. Then, we tested in vitro the cytotoxic efﬁcacy of Chl, Mel, and Cyt against A20 cells. On the basis of these experimental data, we found the parameters for cancer cell death rates that were dependent on the concentration of the respective drugs and developed a mathematical model involving nonlinear ordinary differential equations. For the proposed mathematical model, three equilibrium states were analyzed using the general method of Lyapunov, with only one equilibrium being stable. We obtained a very good symmetry between the experimental results and numerical simulations of the model. Our novel model can be used as a general tool to study the cytotoxic activity of various drugs with different doses and modes of action by appropriate adjustment of the values for the selected parameters.


Introduction
CLL is the most common leukemia in adults and is characterized by the uncontrolled growth of mature B lymphocytes (B cells) [1]. Advances in basic research and therapeutics over the last few years have significantly improved the treatment and clinical outcomes of patients with CLL [2][3][4][5]. Nonetheless, new mathematical models that could aid researchers and clinicians to decide how best to combine traditional and novel agents and to sequence these treatments might result in a reduction in the treatment side effects and the development of drug resistance, which would in turn lead to deep and long-term remissions.
Mathematical modeling of cancer growth has a long history and has been reviewed before [6][7][8][9]. Ordinary Differential Equations (ODE) are used in mathematical models to define the modification of continuous variables and to consider the change over time or space and may vary in complexity from a single equation to multivariable systems. Earlier efforts led to the creation of models with a simplified empirical structure, presented as ODEs, describing the growth of tumor size using exponential or sigmoidal functions [10]. Over the past decade, several mathematical models have described the dynamics of blood cancers under the influence of various drugs and/or immunotherapy [11][12][13][14]. With special regard to CLL, several studies [15][16][17] have modeled the interaction between peripheral blood lymphocytes and CLL cells. These studies produced dynamic systems and identified parameters that describe the growth rate of cancer cells. Kuznetsov et al. [18] was one of the first groups to provide a sensitivity analysis and to describe the observation that cancer cells remain in an apparently dormant state for a prolonged time before entering a rapid growth phase. More recently, Reference [19] used a compartmental model to separate circulating CLL cells from those in lymphoid tissues. This model tracked the production of new CLL cells and showed that the growth rate and death rate of CLL cells are similar. These examples demonstrate that mathematical modeling can have various applications in cancer therapy. However, a major limitation of these quantitative models is that they are not based on real-life experimental data. Furthermore, these models depict the relationship between CLL growth and the immune response to it, but they do not describe or predict the optimal use of drugs for this disease.
Previous studies have explored combining tumor growth kinetics and drug pharmacokinetics with computational machine learning to predict the treatment sequence [20], drug neurotoxicity, and efficacy in drug discovery [21,22]. However, machine learning requires large datasets, and the resulting algorithm is not readily available to the investigator. The model described in our work was based on that used by [23], who employed a system of ODEs to study chemotherapy and immunotherapy in CLL, but we introduced modifications and focused on the effect of chemotherapy on CLL cell growth.
We first grew A20 murine leukemic cells as a surrogate for CLL cells and described their growth rate mathematically. We then tested in vitro the cytotoxicity of three drugs, Chlorambucil (Chl), Melphalan (Mel), and Cytarabine (Cyt), against the cells, using the results to build a dynamic model that incorporates both cancer cell growth and death rates in relation to drug concentration. We then performed model validation and finally a stability analysis. Our results led to the development of a new model that allows the prediction of the optimal drug concentration that will produce the maximum death of cancer cells. Our novel model can be used as a general tool to study the cytotoxic activity of various drugs with different doses and modes of action by appropriate adjustment of the values for the selected parameters.

Cells and Reagents
A20 murine leukemic cells (obtained from the ATCC, VA, USA) were grown in RPMI 1640 (Thermo Fischer Scientific, Waltham, MA, USA). All media were supplemented with 10% Fetal Bovine Serum (FBS) (Thermo Fischer), 1% L-glutamine, and 0.33% Pen-Strep solution. Cells were maintained at 37°C and 5% CO 2 . Cell growth was measured with the XTT-based Cell Proliferation Kit (Biological Industries, Bet Haemek, Israel) according to the manufacturer's instructions.
To determine the growth rate, A20 cells were seeded at three concentrations: 5 × 10 3 ; 1 × 10 4 , and 5 × 10 4 mL. An aliquot was taken from the cultures every 12 h from the 48th until the 312th hour for total and viable cell counting. Experiments were repeated three times. The cells were stained using trypan blue, and live and dead cells were counted with a hemocytometer from duplicate samples of the same culture ( Figure 1).

Figure 1.
Hemocytometer cell counting. Using a microscope, the live (transparent circle) and dead (blue circle) cells were counted per milliliter separately. The average cell count from each of the sets of 16 corner squares was multiplied by 2 to correct for the 1:1 dilution from the trypan blue addition and then multiplied by 10,000.

Drug Cytotoxicity Assay
We chose to study the effect of three drugs known to induce cell death by interfering with DNA replication by different mechanisms. Chlorambucil is an alkylating agent that induces cross-links between DNA strands, resulting in interference with DNA replication [24]. Melphalan alkylates guanine, which leads to interstrand and intrastrand DNA adducts, resulting in inhibition of DNA and RNA synthesis and eventual cell death [25]. Cytarabine, also known as arabinosylcytosine (ARA-C), is a pyrimidine analog and competes with cytidine for incorporation into the DNA, leading to cessation of DNA replication and DNA damage response mechanisms [26].
Drug cytotoxicity was determined by culturing A20 cells in macroplate wells (Nunclon) at an initial concentration of 5 × 10 4 /well for 72 h in culture medium containing either Chl, Mel, or Cyt at concentrations ranging from 0-50 µM. Cell viability was assessed by using the XTT assay, which measures cellular ATP levels. Absorbance in the wells was measured at 450 nm and subtracted from the reference absorbance at 630 nm. Culture medium was used as the background control.

Validation of the Model
The data from the in vitro experiments and the parameters from Table 1 were used to validate our model. In addition, computer simulations were performed using fourthorder adaptive step Runge-Kutta integration, as implemented in the ODE45 subroutine of MATLAB.

Cell Growth and Death Dynamics
To determine cell growth and death rates, A20 cells were sampled every 12 h for total and viable cell counts ( Figure 2). The natural decrease in cell viability after the maximum peak was due to suboptimal growth conditions resulting from nutrient deprivation, decreased pH, and cell overcrowding, which induces cell death by necrosis and apoptosis [29].

Drug Cytotoxicity
Cells were cultured for 72 h in fresh medium containing either Chl, Mel, or Cut, and the cellular metabolic status was determined. The results are shown in Figure 3. As can be seen above (Figure 3), Chl was the least effective of the three drugs, producing only 56% cell growth inhibition at the maximum concentration used in the experiment, 50 µM, while Mel induced 92% inhibition at this dose. At the lowest dose of 1.56 µM, Chl induced only 7% growth inhibition, while Mel induced 49% inhibition. On the other hand, Cyt maintained a high plateau level of growth inhibition of between 96% and 92% inhibition down to 3.125 µM.
These experiments provided the basis for the development of our dynamic model.

Formulation of the Model
Based on cell counts with a hemocytometer (Figures 1 and 2), A20 cells can be divided into two groups: live cells (transparent circle in Figure 1) and dead cells (blue circle in Figure 1), denoted as A and A d , respectively, in the mathematical model developed by us. Dead cells (A d ) are formed because of apoptosis and/or necrosis in culture after the death of living A cells. The death rate of living A cells depends on several conditions in the culture environment: the concentration of essential nutrients in the cells, the supply and uptake of oxygen, and the disposal of waste products such as ammonia and various acidic products [29][30][31]. The total number of living and dead cells cannot exceed a certain capacity (K) in the closed space of the culture. The dynamics of living and dead cells of type A20 in a closed environment is described using differential equations as follows: (1) (2) dA dt describes the dynamics of living A20 cells. It comprises two terms: one positive, corresponding to the logistic cancer growth characterized by the coefficient, r, which is limited by the maximal tumor cell number, K; one negative term, corresponding to living cells becoming dead at a rate of µ A . dA d dt describes the dynamic of dead A20 cells. It is comprised of two terms: the positive term is the death of A20 cells with a rate coefficient of µ A due to apoptosis or necrosis and depends on living A20 cells competing for survival (oxygen consumption and nutrition) in an enclosed space; the negative term corresponds to dissolution of dead cells at a rate of d. We extended the model by adding a new equation that represents the dynamics of chemotherapeutic drugs. Based on previous studies [23] and our experiments in vitro, we formulated an ODE model to mathematically explain the interaction between CLL cells and chemotherapeutic drugs: dC dt describes the first-order pharmacokinetics of a drug [32]. The drug was given only once at the beginning of the experiment, i.e., C(0) is a constant value and depends on the dose of the drug. The terms 1 and 2 of the Equations (3) and (4) represent the log-kill hypothesis [33], with a Michaelis-Menten drug saturation response [34], a + C; µ AC is the death rate resulting from the action of the drug on the cancer cells. The parameter µ AC changes depending on the drug dose and on the particular drug. We chose the parameter µ CA to be ten-times more than µ AC , assuming that there are 10-100 drug molecules attacking each cancer cell. Ultimately, the parameter µ CA does not play a significant role in the model, since there are many more drug molecules than cancer cells, and it only reflects the amount of available drug molecules. The parameter a represents the drug concentration that produces 50% of the maximum activity of the drug in each cell population [27].
We performed a mathematical analysis of our model by identifying fixed points and their stability. It was found that the system is characterized by three fixed points, one of which is stable asymptotically (Appendix A, Table A4).

Estimation of the Parameters of the Model
In this section, we evaluate the model parameters (Table 1) The number of drug molecules was calculated using the formula: The The cell doubling time was calculated using an exponential growth rate: where: Thus, according to Figure 2 the number of A20 cells double in less than 12 h. The parameters used in this study are summarized in Table 1.  To assess the closeness of the fit between the simulation and experimental curves, the root-mean-squared errors (RMSEs) were calculated using the formula:

Validation of the Cancer Cell Growth Dynamics Model
where Sim is the numerical simulation data (predicted value), Exp is the experimental data (observed value), and n is the number of observations (see Appendix A, Table A5). The RMSEs for the initial cell concentrations are: 5 × 10 4 cells/mL = 0.016; 1 × 10 4 cells/mL = 0.013; 5 × 10 3 cells/mL = 0.008. The data showed a strong correlation between the simulation and experimental results for both the time and cell concentration at both the maximum and final data points. The correlation held for all three starting concentrations (all RMSEs were < 0.1), but was most symmetric for the lowest starting concentration (5 × 10 3 cells/mL).

Validation of the Cancer Cell Drug Cytotoxicity Dynamics Model
By inserting the parameters from Table 1 into Equations (2) and (3), we could numerically simulate the effect of Chl, Mel, or Cyt on A20 cells, as depicted in Figure 5. After 72 h, the concentration of A20 growth without drug increased to A(72) = 1,646,260. The additional 50 µM of Chl (A) reduced the growth to A(72) = 717,369 cells/mL, which equals 56.4% growth inhibition; 50 µM of Mel induced 92% growth inhibition (B), and 6.25 µM of Cyt induced 94.8% growth inhibition (C). The complete set of calculated data for all drug concentrations is presented in Appendix A (Tables A1-A3).
We then compared the degree of A20 growth inhibition from the numerical simulation of the model (Equations (3)-(5)) ( Figure 5) with the output obtained by the in vitro experiments (Figure 3). The results are shown in Figure 6.  Figure 6 demonstrates that our simulation coincided with the experimental data for the high doses of each drug. For the intermediate doses, there appeared to be a slight deviation between the two sets of data. This deviation was artificial since we tried to maintain consistency in the decrease of the parameter µ AC in accordance with the decrease in the drug dose. This resulted in a discrepancy of 30% for each drug. However, it is important to note that maintaining this consistency is not necessary, and it is possible to juxtapose the experimental and simulation results with absolute accuracy by choosing the appropriate parameter µ AC .
The RMSE for each drug was calculated (see Appendix A, Table A6). The data showed a strong correlation between the simulation and experimental results for cell growth inhibition by various drug concentrations. The correlation held for all three drugs (all RMSEs were < 0.1), but was most apparent for Cyt (RMSE = 0.018).

Conclusions
We constructed a mathematical model (Equations (3)-(5)) of the cytotoxic efficacy of different CLL drugs Chl, Mel, and Cyt. We showed that Cyt is the most effective drug amongst those tested in killing A20 leukemic cells. We applied a multistep approach to the parameter estimation, which is appropriate for the analysis of experimental in vitro data with different concentrations.
Our simulation results indicated that the model correctly describes the interaction of drugs with leukemic cells. The data showed a strong correlation between the numerical simulations and experimental results for both the growth of the tumor cells and of their killing by the drugs, essentially due to the fact that we separated the dynamics of dead and living cells, since this plays a large role in in vitro experiments in closed spaces.
The advantage of our model is that it is flexible and can be adjusted to various drugs with different doses by changing the values of certain parameters. While the results in this work were generated using one cell line, clearly, this model is fundamental and can be applied to any other cancer cell line once the cell's growth rate has been calculated.
Already at this stage of development, the presented model is a tool for predicting the most effective drug in the treatment against CLL. In our further researches, we plan to expand the current model and to achieve the approach for in vivo experiments.
Although recent years have witnessed significant progress in the therapy of patients with CLL, the disease remains incurable, and new treatments are needed. This improved tool for CLL chemotherapy modeling will ultimately have basic research and therapeutic value. We believe that quantitative simulation of cancer chemotherapies, based on experimentally validated mathematical modeling, provides an opportunity for the researcher, and eventually the clinician, to address data for predicting possible scenarios of cancer developing in individual patients and establishing an efficient treatment protocol.  Tables A1-A3 present the data on the variations of parameter µ AC , the type and dose of drug, and number of cells at the 72nd hour of the numerical simulations (see Figure 5). By setting all equations for A * , A * d , and C * to zero, the fixed points of our model were derived. We chose only the non-negative equilibria, assuming all initial conditions to be positive.
From Equation (A1), it follows: From Equation (A2), it follows: Notice that when C * ≠ 0, from Equation (A3), it follows that A * = − µ C (a + C * ) µ CA , i.e., either Thereby, from Equations (A1)-(A3), there are three non-negative equilibria for consideration: when A * = 0 and C * = 0, and when A * d = 0 and C * = 0 ⇒ A * = K, To analyze the asymptotic stability of each equilibrium of the above nonlinear system, the eigenvalues of the Jacobian were calculated at a particular equilibrium wherē λ = [λ 1 , λ 2 , λ 3 ] were set at: If Λ < 0, then all real parts of the eigenvalues are negative, and we can determine that the equilibrium is asymptotically locally stable.
Steady states of the system in Equations (A1)-(A3): The system is characterized by the three following fixed points. Using the parameters from Table 1, we obtain the eigenvalues for every equilibrium. 1.
The eigenvalues of this Jacobian are: Thus, the fixed point E * 0 is unstable. Equilibrium E * 0 exists only if A(0) = A d (0) = C(0) = 0, which has no biological significance; The eigenvalues of this Jacobian are: Thus, there is an asymptotic stability ( Figure A1A); 3.
The eigenvalues of this Jacobian are: Thus, E * 2 is unstable. From a biological point of view, E * 2 means that the dose of chemotherapy was insufficient, which permitted the cancer cells to achieve the maximum growth capacity. However, considering that E * 2 is unstable, as can be seen in Figure A1B, after 60 h, the cancer cells experience a natural death.