Tumor growth and chemotherapy: Mathematical methods, computer simulations, and experimental foundations

Validation of a Mathematical Model Describing the Dynamics of Chemotherapy for Chronic Lymphocytic Leukemia In Vivo

In recent years, mathematical models have developed into an important tool for cancer research, combining quantitative analysis and natural processes. We have focused on Chronic Lymphocytic Leukemia (CLL), since it is one of the most common adult leukemias, which remains incurable. As the first step toward the mathematical prediction of in vivo drug efficacy, we first found that logistic growth best described the proliferation of fluorescently labeled murine A20 leukemic cells injected in immunocompetent Balb/c mice. Then, we tested the cytotoxic efficacy of Ibrutinib (Ibr) and Cytarabine (Cyt) in A20-bearing mice. The results afforded calculation of the killing rate of the A20 cells as a function of therapy. The experimental data were compared with the simulation model to validate the latter’s applicability. On the basis of these results, we developed a new ordinary differential equations (ODEs) model and provided its sensitivity and stability analysis. There was excellent accordance between numerical simulations of the model and results from in vivo experiments. We found that simulations of our model could predict that the combination of Cyt and Ibr would lead to approximately 95% killing of A20 cells. In its current format, the model can be used as a tool for mathematical prediction of in vivo drug efficacy, and could form the basis of software for prediction of personalized chemotherapy.

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

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.

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 efficacy studies. As a first step toward mathematical prediction of in vivo drug efficacy from in vitro cytotoxicity studies, we used murine A20 leukemic cells as a test case of CLL. We first found that logistic growth best described the proliferation of the cells in vitro. Then, we tested in vitro the cytotoxic efficacy 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.

Optimization of additive chemotherapy combinations for an in vitro cell cycle model with constant drug exposures

Proliferation of an in vitro population of cancer cells is described by a linear cell cycle model with n states, subject to provocation with m chemotherapeutic compounds. Minimization of a linear combination of constant drug exposures is considered, with stability of the system used as a constraint to ensure a stable or shrinking cell population. The main result concerns the identification of redundant compounds, and an explicit solution formula for the case where all exposures are nonzero. The orthogonal case, where each drug acts on a single and different stage of the cell cycle, leads to a version of the classic inequality between the arithmetic and geometric means. Moreover, it is shown how the general case can be solved by converting it to the orthogonal case using a linear invertible transformation. The results are illustrated with two examples corresponding to combination treatment with two and three compounds, respectively.

Measuring differences between phenomenological growth models applied to epidemiology

Phenomenological growth models (PGMs) provide a framework for characterizing epidemic trajectories, estimating key transmission parameters, gaining insight into the contribution of various transmission pathways, and providing long-term and short-term forecasts. Such models only require a small number of parameters to describe epidemic growth patterns. They can be expressed by an ordinary differential equation (ODE) of the type C^'(t)=f(t,C;Θ) for t>0, C(0)=C₀, where t is time, C(t) is the total size of the epidemic (the cumulative number of cases) at time t, C₀ is the initial number of cases, f is a model-specific incidence function, and Θ is a vector of parameters. The current COVID-19 pandemic is a scenario for which such models are of obvious importance. In Bürger et al. (2019) it is demonstrated that some PGMs are better at fitting data of specific epidemic outbreaks than others even when the models have the same number of parameters. This situation motivates the need to measure differences in the dynamics that two different models are capable of generating. The present work contributes to a systematic study of differences between PGMs and how these may explain the ability of certain models to provide a better fit to data than others. To this end a so-called empirical directed distance (EDD) is defined to describe the differences in the dynamics between different dynamic models. The EDD of one PGM from another one quantifies how well the former fits data generated by the latter. The concept of EDD is, however, not symmetric in the usual sense of metric spaces. The procedure of calculating EDDs is applied to synthetic data and real data from influenza, Ebola, and COVID-19 outbreaks.

Modeling iontophoretic drug delivery in a microfluidic device

Iontophoresis employs low-intensity electrical voltage and continuous constant current to direct a charged drug into a tissue. Iontophoretic drug delivery has recently been used as a novel method for cancer treatment in vivo. There is an urgent need to precisely model the low-intensity electric fields in cell culture systems to optimize iontophoretic drug delivery to tumors. Here, we present an iontophoresis-on-chip (IOC) platform to precisely quantify carboplatin drug delivery and its corresponding anti-cancer efficacy under various voltages and currents. In this study, we use an in vitro heparin-based hydrogel microfluidic device to model the movement of a charged drug across an extracellular matrix (ECM) and in MDA-MB-231 triple-negative breast cancer (TNBC) cells. Transport of the drug through the hydrogel was modeled based on diffusion and electrophoresis of charged drug molecules in the direction of an oppositely charged electrode. The drug concentration in the tumor extracellular matrix was computed using finite element modeling of transient drug transport in the heparin-based hydrogel. The model predictions were then validated using the IOC platform by comparing the predicted concentration of a fluorescent cationic dye (Alexa Fluor 594®) to the actual concentration in the microfluidic device. Alexa Fluor 594® was used because it has a molecular weight close to paclitaxel, the gold standard drug for treating TNBC, and carboplatin. Our results demonstrated that a 50 mV DC electric field and a 3 mA electrical current significantly increased drug delivery and tumor cell death by 48.12% ± 14.33 and 39.13% ± 12.86, respectively (n = 3, p-value <0.05). The IOC platform and mathematical drug delivery model of iontophoresis are promising tools for precise delivery of chemotherapeutic drugs into solid tumors. Further improvements to the IOC platform can be made by adding a layer of epidermal cells to model the skin.

Mathematical model of brain tumour with glia-neuron interactions and chemotherapy treatment

In recent years, it became clear that a better understanding of the interactions among the main elements involved in the cancer network is necessary for the treatment of cancer and the suppression of cancer growth. In this work we propose a system of coupled differential equations that model brain tumour under treatment by chemotherapy, which considers interactions among the glial cells, the glioma, the neurons, and the chemotherapeutic agents. We study the conditions for the glioma growth to be eliminated, and identify values of the parameters for which the inhibition of the glioma growth is obtained with a minimal loss of healthy cells.

Towards more effective advanced drug delivery systems

This position paper discusses progress made and to be made with so-called advanced drug delivery systems, particularly but not exclusively those in the nanometre domain. The paper has resulted from discussions with a number of international experts in the field who shared their views on aspects of the subject, from the nomenclature used for such systems, the sometimes overwrought claims made in the era of nanotechnology, the complex nature of targeting delivery systems to specific destinations in vivo, the need for setting standards for the choice and characterisation of cell lines used in in vitro studies, to attention to the manufacturability, stability and analytical profiling of systems and more relevant studies on toxicology. The historical background to the development of many systems is emphasised. So too is the stochastic nature of many of the steps to successful access to and action in targets. A lacuna in the field is the lack of availability of data on a variety of carrier systems using the same models in vitro and in vivo using standard controls. The paper asserts that greater emphasis must also be paid to the effective levels of active attained in target organs, for without such crucial data it will be difficult for many experimental systems to enter the clinic. This means the use of diagnostic/imaging technologies to monitor targeted drug delivery and stratify patient groups, identifying patients with optimum chances for successful therapy. Last, but not least, the critical importance of the development of science bases for regulatory policies, scientific platforms overseeing the field and new paradigms of financing are discussed.

In silico modelling of treatment-induced tumour cell kill: developments and advances

Mathematical and stochastic computer (in silico) models of tumour growth and treatment response of the past and current eras are presented, outlining the aims of the models, model methodology, the key parameters used to describe the tumour system, and treatment modality applied, as well as reported outcomes from simulations. Fractionated radiotherapy, chemotherapy, and combined therapies are reviewed, providing a comprehensive overview of the modelling literature for current modellers and radiobiologists to ignite the interest of other computational scientists and health professionals of the ever evolving and clinically relevant field of tumour modelling.

ON THE GOMPERTZ LIMIT OF THE MONOTONIC NEOCLASSICAL GROWTH MODEL

The burden of proof of any theory aiming to represent a physical or biological reality by demonstrating its unifying properties is applied in the present paper in relation to the Neoclassical growth model and its ability to reproduce Gompertz growth. The Neoclassical growth model derived from first biological and physical principles was shown to capture all qualitative features that were revealed experimentally, including the possibility of a Logarithmic Inflection Point (LIP), the possibility of a LAG, concave as well as convex curves on the phase diagram, the Logistic growth as a special case, growth followed by decay, as well as oscillations. In addition, quantitative validation demonstrated its ability to reproduce experimental data in a few tested cases. This paper demonstrates that the Neoclassical growth model can reproduce a Generalized version of Gompertz growth too.

Growth of a virtual tumour using probabilistic methods of cell generation

A study into treatment enhancement in combined chemo-radiotherapy for unresectable head and neck cancer has initiated the development of a computer model of tumour growth. The model is based on biological parameters, and characterises tumour growth prior to chemo-radiotherapy. Tumour growth starting from a single stem cell is modelled using the Monte Carlo method. The type of the cell function, their relative proportions on mitosis, their proliferative capacity, the duration of the four phases of the cell cycle, the mean cell cycle time, and the cell loss due to natural causes are the main parameters of the basic model. A Gaussian distribution function operates in establishing the cell cycle time, with a mean value of 33 hours, while the cell type is sampled from a uniform distribution. With the established model, the sensitivity of the developed tumour's cell population to the stem, proliferative and nonproliferative ratio at mitosis was assessed. The present model accurately reflects the exponential distribution of cells along the cell cycle (70% cells in GI phase, 15% in S, 10% in G2, 5% in M) of a developed tumour as described in the literature. The proportion of stem, finitely proliferating and resting cells during tumour growth is maintained within their biological limits (2% stem, 13% finitely proliferating, 85% nonproliferating cells). The ratio (R = 3) between the time necessary to develop a clinically detectable tumour (10(9) cells) and the further time to grow to its lethal size (10(12) cells) is in accordance with the biological data when tumour volume is compared for the two periods (30 doublings and 10 doublings respectively). In conclusion, computer simulation can illustrate the biological growth of a tumour and the cell distribution along the cell cycle. These distributions may then be used in the assessment of tumour response to radiotherapy and to specific chemotherapeutic agents.

Pharmacokinetic/pharmacodynamic modeling of antitumor agents encapsulated into liposomes

Pharmacokinetic/pharmacodynamic (PK/PD) modeling of antitumor agents has been developed for doxorubicin (DOX) in order to predict the optimum conditions for a drug carrier to maximize the antitumor effect. A PK model was constructed for free and liposomal doxorubicin using a hybrid model wherein the disposition in the whole body is described by compartment models, which were linked to the tumor compartment via the blood flow rate. The PD model for doxorubicin was described by a cell-kill kinetic model, which represents the number of tumor cells quantitatively, as a function of the free concentration of doxorubicin in the tumor compartment. The influence of each parameter on the antitumor effects was examined by sensitivity analysis based on the PK/PD model, which clearly showed the importance of optimizing the release rate of DOX from liposomes. The validity of the model has been tested using animal experiments. Preliminary simulations were also performed for humans after scaling up the PK/PD model from rodents to humans. The optimum conditions in the rate of drug release from liposomes were different for rodents vis-a-vis humans, which indicates the limitations involved in extrapolating optimum conditions for experimental animals to those for humans.

Applying computer modeling to examine complex dynamics and pattern formation of tissue growth

For research in areas such as developmental biology and cancer, understanding the formation and regrowth of tissue is of great importance. Since complex system behavior makes it difficult to interpret dynamics and pattern formation of tissue growth, it is helpful to have a way to simulate the cell systems and test hypotheses about the mechanisms by which the system is responding. Computer "modeling experiments" can serve this purpose, as we show with an example of the small intestine epithelial cells' response to cytosine arabinoside. This example demonstrates that nonhomogeneities in the cell population can play an important role and emphasizes the need to use a modeling approach, like our spatial modeling, that addresses the differences within the tissue structure. These types of "modeling experiments" can guide researchers with further experiments and provide clues as to how complex cell proliferation behavior is linked to underlying molecular phenomena.

A mathematical model of periodically pulsed chemotherapy: tumor recurrence and metastasis in a competitive environment

A competition model describing tumor-normal cell interaction with the added effects of periodically pulsed chemotherapy is discussed. The model describes parameter conditions needed to prevent relapse following attempts to remove the tumor or tumor metastasis. The effects of resistant tumor subpopulations are also investigated and recurrence prevention strategies are explored.

A mathematical model for the growth and classification of a solid tumor: a new approach via nonlinear elasticity theory using strain-energy functions

Medically, tumors are classified into two important classes--benign and malignant. Generally speaking, the two classes display different behaviour with regard to their rate and manner of growth and subsequent possible spread. In this paper, we formulate a new approach to tumor growth using results and techniques from nonlinear elasticity theory. A mathematical model is given for the growth of a solid tumor using membrane and thick-shell theory. A central feature of the model is the characterization of the material composition of the tumor through the use of a strain energy function, thus permitting a mathematical description of the degree of differentiation of the tumor explicitly in the model. Conditions are given in terms of the strain energy function for the processes of invasion and metastasis occurring in a tumor, being interpreted as the bifurcation modes of the spherical shell, which the tumor is essentially modeled as. Our results are compared with actual medical experimental results and with the general behavior shown by benign and malignant tumors. Finally, we use these results in conjunction with aspects of surface morphogenesis of tumors (in particular, the Gaussian and mean curvatures of the surface of a solid tumor) in an attempt to produce a mathematical formulation and description of the important medical processes of staging and grading cancers. We hope that this approach may form the basis of a practical application.

Role of optimal control theory in cancer chemotherapy

This paper presents a review of the ways in which optimal control theory interacts with cancer chemotherapy. There are three broad areas of investigation. One involves miscellaneous growth kinetic models, the second involves cell cycle models, and the third is a classification of "other models." Both normal and tumor cell populations are included in a number of the models. The concepts of deterministic optimal control theory are applied to each model in such a way as to present a cohesive picture. There are applications to both experimental and clinical tumors. Suggestions for designing better chemotherapy strategies are presented.

A mathematical model of cancer chemotherapy with an optimal selection of parameters

An optimal parameter selection model of cancer chemotherapy is presented which describes the treatment of a tumor over a fixed period of time by the repeated administration of a single drug. The drug is delivered at evenly spaced intervals over the treatment period at doses to be selected by the model. The model constructs a regimen that both minimizes the tumor population at the end of the treatment and satisfies constraints on the drug toxicity and intermediate tumor size. Numerical solutions show that an optimal regimen withholds the bulk of the doses until the end of the treatment period. When a drug used is of either moderate or low effectiveness, an optimal regimen is superior to a schedule that delivers all of the drug at the beginning of the treatment. This study questions whether the current method for the administration of chemotherapy is optimal and suggests that alternative regimens should be considered.

A model of proliferating cell populations with inherited cycle length

A mathematical model of cell population growth introduced by J.L. Lebowitz and S.I. Rubinow is analyzed. Individual cells are distinguished by age and cell cycle length. The cell cycle length is viewed as an inherited property determined at birth. The density of the population satisfies a first order linear partial differential equation with initial and boundary conditions. The boundary condition models the process of cell division of mother cells and the inheritance of cycle length by daughter cells. The mathematical analysis of the model employs the theory of operator semigroups and the spectral theory of linear operators. It is proved that the solutions exhibit the property of asynchronous exponential growth.

Mathematical model for human myeloma relating growth kinetics and drug resistance

We present a computer-based mathematical model that can simulate characteristic features of the clinical time course of human myeloma. It asserts that therapy resistance in myeloma cells is an inherited trait associated with the longer inter-mitotic times of some cells and that the strength of this trait affects tumour growth characteristics. These kinetic differences within the malignant cell clone may also influence therapeutic efficacy. In the model, the same total therapy, administered in different time-dose fractions, could be 'curative' or 'minimally effective' depending on kinetic properties. For example, as others have shown, in myeloma pulsed intermittent therapy is often more effective than low dose continuous therapy. According to our model this finding is compatible with a high coefficient of inheritability of resistance from one cell generation to the next. The model also suggests that if there are subclones of varying resistance, a therapy must have some effect on each of them if it is to be employed in a curative fashion. While many aspects of the model are not yet clinically testable, exploration of its concepts might increase knowledge about fundamental neoplastic mechanisms.

On the equivalence of mathematical models for cell proliferation kinetics

Various types of mathematical models, such as partial differential equations, ordinary differential equations and difference equations, are available in the literature to describe the kinetics of cell proliferation, and different studies of cell kinetic phenomena have been conducted using these models. This paper discusses the equivalence between the different models identifying the conditions and approximations under which one type of models may be derived from another. Such an equivalence study is highly useful for an integration of the diverse results that have been obtained using different models in order to gain a more complete understanding of cell kinetic phenomena.

Modeling cellular systems and aging processes: I. Mathematics of cell system models-a review

A review of the literature on mathematical models of populations of cellular systems is presented. Continuous, discrete, and stochastic models are presented in a semihistorical manner as a prelude to answering the question of how to model an asynchronously dividing cellular system. This analysis is then broadened, in an attempt to broach the more general question of modeling the distribution of a set or collection of cell properties through an asynchronously dividing cellular system. Such properties might be cell motility, cell cycle length, time to mitosis, or number of epigenetic particles. It is shown that one fruitful approach to this modeling question is a coupled continuous-probabilistic model. The ramifications of this type of formalism are discussed.

Determination of steady state generation time distributions from labelled mitosis experimental data

The proper application of detailed deterministic cell kinetic models depends on the way in which cells are assigned their generation times. A method is presented for the determination of population generation time distributions from labelled mitoses experiments. The model assumes that the generation time of each new cell is a function of both the steady-state generation time distribution function of the population, and also the generation time frequency-function of the previous generation of cells. This approach is applied to two different cell types to successfully simulate extended labelled mitoses curves using a population balance model with constant maturation rates.

Tumor growth patterns in multiple myeloma

Serial changes in tumor mass were evaluated in 61 patients with multiple myeloma who had received intermittent courses of melphalan-prednisone until death. The variations in the kinetics of tumor reduction and relapse could be explained by a mathematical model based on two cell populations, one sensitive to and one resistant to chemotherapy. For all responding patients, the median tumor halving-time was 1.3 months and the median doubling time was 2.9 months. The duration of a constant tumor mass during remission was brief in most patients. A larger fraction of resistant cells prior to therapy was associated with a slower tumor doubling-time during relapse. With a constant fractional reduction of sensitive cells and a tumor halving-time of one month or less, all cells sensitive to alkylating agents would be eliminated with 3 years of uninterrupted intermittent therapy.

Kinetic analysis of cell size and DNA content distributions during tumor cell proliferation: Ehrlich ascites tumor study

In order to study the growth dynamics of proliferating and non-proliferating cells utilizing discrete-time state equations, the cell cycle was divided into a finite number of age compartments. In analysing tumor growth, the kinetic parameters associated with a retardation in the growth rate of tumors were characterized by computer simulation in which the simulated results of the growth curve, the growth fraction, and the mean generation time were adjusted to fit the experimental data. The cell age distibution during the period of growth was obtained and by a linear transformation of the state transition matrices, was employed to specify the cell size and DNA content distributions. In an application of the model, the time-course behavior of cell cycle parameters of Ehrlich ascites tumor is illustrated, and the parameters important for the transition of cells in the proliferating compartment to the non-proliferating compartment are discussed, particularly in relation to the G1-G0 and G2-G0 transitions of non-cycling cells as revealed by the variation of cell size distribution.