# Mathematical Modeling, Analysis, and Simulation of Tumor Dynamics with Drug Interventions

Authors: Pranav Unni; Padmanabhan Seshaiyer
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

Over the last few decades, there have been significant developments in theoretical, experimental, and clinical approaches to understand the dynamics of cancer cells and their interactions with the immune system. These have led to the development of important methods for cancer therapy including virotherapy, immunotherapy, chemotherapy, targeted drug therapy, and many others. Along with this, there have also been some developments on analytical and computational models to help provide insights into clinical observations. This work develops a new mathematical model that combines important interactions between tumor cells and cells in the immune systems including natural killer cells, dendritic cells, and cytotoxic CD8+ T cells combined with drug delivery to these cell sites. These interactions are described via a system of ordinary differential equations that are solved numerically. A stability analysis of this model is also performed to determine conditions for tumor-free equilibrium to be stable. We also study the influence of proliferation rates and drug interventions in the dynamics of all the cells involved. Another contribution is the development of a novel parameter estimation methodology to determine optimal parameters in the model that can reproduce a given dataset. Our results seem to suggest that the model employed is a robust candidate for studying the dynamics of tumor cells and it helps to provide the dynamic interactions between the tumor cells, immune system, and drug-response systems.


---

## Body

## 1. Introduction

Cancer is one of the leading causes of death in the world today. By 2030, over 13 million are estimated to harbor some form of the disease. While there have been many developments in cancer therapies including surgery, chemotherapy, immunotherapy, and radiotherapy, there is still a lot that is unknown about the dynamics of how cancer cells are created, propagated, and destroyed.Over the past few decades, there have been several experimental approaches and interventions developed that have helped us to understand the dynamics of tumor growth and its interactions with the immune system [1, 2]. This has also helped to inform how specific interventions such as immunotherapy can help strengthen our own ability to fight cancer by improving the effectiveness of the immune system [3–5]. While these developments have helped enhance our understanding about cancer dynamics, there are still several challenges in these experimental approaches to fully understand the interactions with the immune system.In the last two decades, there have also been several experimental advances in developing interventional therapies for cancer such as immunotherapy, virotherapy, targeted drug therapies, and chemotherapy. Along with these experimental developments, there have been some advances in scientific and engineering solutions to capture the dynamics of cancer. One of the promising approaches includes mathematical modeling [6, 7], which involves identifying the cells that play a role in cancer propagation, interactions between these bodies, and description of the dynamics of this interaction that has helped estimate parameters, perform stability analysis, and predict tumor dynamics [8–13]. These models have been able to demonstrate the importance of the presence of immune components for explaining clinically observed phenomena such as tumor dormancy [14], tumor size oscillations and regressions [15–18], nonspatial models of tumor and immune system interactions [19, 20], and tumor growth coupled with immunotherapy [8, 9, 12, 13, 21, 22].These mathematical models are often coupled system of governing differential equations that describe the dynamics of each of the interacting component cells. Specifically, the interactions between tumor growth and the immune system are often described using a system of coupled differential equations with prescribed initial conditions. These equations include nonlinear interactions and do not often admit an exact solution and therefore require computational methods to solve them. While these mathematical models have provided useful information regarding the importance of the immune system in controlling tumor growth, there is still a great need to continue to enhance existing models to incorporate new clinical developments and biological discoveries. For example, there have been studies suggesting the effectiveness of chemotherapy with immunotherapy and vaccine therapies [1, 13, 23]. The focus of this paper is to enhance existing models of tumor growth that incorporate tumor dynamics in conjunction with the immune system response and also study the effect of additional interventions including antitumor vaccination and immunotherapies along with chemotherapy.Of the many clinical approaches that are tested for cancer therapy, one of the popular approaches includes drug therapy to the tumor microenvironment. To understand the impact of the drugs delivered to the tumor cell site, it is important to include the effect of these drugs into the models as well. Towards this end, we develop a mathematical model that will combine essential interactions between growing tumor cells and cells of the innate and specific immune system coupled with models for drug delivery to these cell sites. Our goal is to use these models developed to study the effectiveness of anticancer drugs to reduce tumor growth.The growth of tumors has also been attributed to the dynamics of the cellular immune system within the human host. Two principal components of this immune system include the natural killer cells and cytotoxicCD8+ T cells which are known to kill tumor cells. Besides these, other important antigen-presenting cells include the dendritic cells that help stimulate, recruit, and activate the immune system. While research has been growing to discover potential mechanisms to describe immune system interactions with growing tumors, there is enough evidence that the dynamics of natural kills cells, cytotoxic CD8+ T-cells, and dendritic cells influence tumor dynamics. The model presented in this work will account for the influence of these as well.Over the years, there has been a lot of development in mathematical modeling of cancer. However, the mechanisms that are involved in the interactions of tumor cells with the immune system are still not clear. This paper attempts to make a new contribution in this direction by developing a coupled mathematical model that incorporates tumor dynamics and interactions between the dendritic cells, natural killer cells, andCD8+ T cells. Additionally, the model incorporates and studies the influence of various drug therapies including immunotherapy and chemotherapy. Finally, a new parameter estimation technique is proposed that helps to estimate parameters optimally for a given extrapolated dataset.


## 2. Models and Background

In this work, we will consider a model that consists of four main cell populations including tumor cells (Tt), natural killer cells (Nt), dendritic cells (Dt), and cytotoxic CD8+ T cells denoted by (Lt). The dynamics of these cells will include interactions between each other as well as dynamics generated by interaction with chemotherapy as well immunotherapy drug concentrations in the blood stream.For developing the model for each of the cell populations, a standard approach is to begin with applying conservation of mass with diffusion and activation. This would often yield the following equation for the dynamics of the various types of cells:(1)∂⋅∂t+∇⋅u→⋅−δD∇2⋅=f⋅−g⋅−K⋅zM⋅.Here, the functionsf⋅ and g⋅ will be based on proliferation rates, competition terms, and inhibition terms based on the respective roles of each type of cell.Also, note that, in all the models, we will consider the effect of a chemotherapy drug (dynamics described later) kill term throughK⋅zM⋅. The term zM=1−e−M is used to denote the fact that chemotherapeutic drugs (for example, doxorubicin) are only effective during certain phases of the cell cycle and pharmacokinetics. The values of the kill parameters K⋅ for the four cell population considered here are based on their ability to disrupt the process of division and growth [24, 25]. Note that if K⋅=0, the equation is not impacted by the drug kill term.While equation (1) includes both the diffusion term and the advection term due to blood velocity u→, we will consider only the temporal dynamics in this work and hence the associated ordinary differential equation:(2)d⋅dt=f⋅−g⋅−K⋅zM⋅.


### 2.1. Modeling Tumor Cells

We begin with modeling tumor cellsT which are assumed to have a proliferation rate that can be modeled by a logistic growth law aT1−bT, with parameters a and b denoting the per capita growth rates and reciprocal carrying capacities of the tumor cells [8, 9, 26]. Also, it is known that the growth of the tumor cells is impacted by three different competitive interactions including interactions between tumor cells and dendritic cells, interactions between tumor cells and natural killer cells, and interactions between tumor cells and CD8+ T cells [26–28]. Denoting the corresponding competition rates as j¯,c1,k, respectively, introduces the competition term as −c1N+jD+kLT. The dynamics of tumor cells can then be described by the following ordinary differential equation:(3)dTdt=aT1−bT−c1N+jD+kLT−KTzMT.


### 2.2. Modeling Natural Killer Cells

To model natural killer (NK) cells, we will assume that these cells have a constant sources1 as well as a NK cell recruitment term that can be represented through a modified Michaelis–Menten term (commonly used to govern cell-cell interactions):(4)g1⋅T2h1+T2⋅N,where g1 denotes the maximum NK cell recruitment rate by tumor cells and h1 denotes the steepness coefficient of the NK cell recruitment curve [8].Next, the growth of NK cells will be impacted by two different interactions, namely, the interaction between NK cells and tumor cells [29] and the interaction between NK cells and dendritic cells [30–33]. We also introduce parameters c2,d1 to be the rates of killing (due to tumor cells) and proliferation (due to dendritic cells) of NK cells, respectively. The governing differential equation for the dynamics of NK cells then can be described as(5)dNdt=s1+g1NT2h1+T2−c2T−d1DN−KNzMN−eN.Note that we have also included a natural death of NK cells through−eN.


### 2.3. Modeling Dendritic Cells

Dendritic cells play an important role in the immune system response and in controlling tumor growth. Also known as antigen-presenting cells, they update and present antigens toCD8+ T cells. Some of the earlier models [8, 13] in the literature have not incorporated the dynamics of dendritic cells in directly suppressing tumor growth, stimulating resting NK cells, and impacting the dynamics of CD8+ T cells. There is, however, experimental evidence that dendritic cells play an important role in modeling tumor immunotherapy [28].To study the dynamics of dendritic cells, we will assumes2 to be a constant source of dendritic cells, d2 to be the rate at which NK cells kill dendritic cells, d3 to be a proliferation rate of dendritic cells due to tumor cells, f1 to be the rate corresponding to the interaction of dendritic cells with CD8+ T cells, and g to be the natural death rate of dendritic cells. We then have(6)dDdt=s2−f1L+d2N−d3TD−KDzMD−gD.


### 2.4. Modeling Cytotoxic CD8+ T Cells

Among many factors that impact the growth of tumor cells, it is well known thatCD8+ T cells are an important component of the immune system that kills tumor cells. It has been seen that active tumor-specific CD8+ T cells are only present in large numbers when tumor cells are present [16, 34], and after some interactions with tumor cells, they become inactive [29]. It has been observed that mature CD8+ T cells can remove dendritic cells [35, 36].In our model, to describe the dynamics of theCD8+ T-cells, we will consider f2 to be the rate of interaction between dendritic cells and tumor cells to activate CD8+ T cells; −hLT denotes the form of competition between CD8+ T cells and tumor cells, and −iL denotes the natural death rate of CD8+ T cells.It has also been seen thatCD8+ T cells may be recruited by the debris from tumor cells lysed by NK cells [37]. To include this effect, we will incorporate in our model a recruitment term that is proportional to the number of cells killed, which is denoted as r1NT. We will also need an additional term that helps describe the regulation and suppression of CD8+ T-cell activity when there are high levels of activated CD8+ T cells without responsiveness to cytokines in the system [38, 39]. This term is denoted by uNL2. We also include CD8+ T activation by IL-2 immunotherapy which is in the form of a drug that influences the immune system’s efficacy and described via a Michaelis–Menten interaction term [16]:(7)pILIgI+I.Here,I refers to the immunotherapy drug concentration in the bloodstream. The model for the dynamics for the CD8+ T cell growth population becomes(8)dLdt=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL.


### 2.5. Modeling Drug and Vaccine Interventions

In this study, we incorporate a variety of external intervention treatment options including tumor-infiltrating lymphocyte (TIL) injections as well as chemotherapy and immunotherapy drugs. TIL drug intervention may be thought of as an immunotherapy approach in which theCD8+ T-cells are promoted through antigen-specific cytolytic immune cells. We do this by adding the term vL=vLt in equation (8) and we have(9)dLdt=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL+vL.To include the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream as follows:(10)dMdt=vMt−d4M,dIdt=vIt−d5I.The drug intervention terms in these equations reflect the amount of chemotherapy and immunotherapy drug given over time. Note that we assume that the chemotherapy and immunotherapy drugs will be eliminated from the body over time at a rate proportional to its concentration, and these are given byd4M and d5I, respectively.


### 2.6. Overall Model

From Sections2.1–2.5, we have the following overall model:(11)T˙=aT1−bT−c1N+jD+kLT−KTzMT,N˙=s1+g1NT2h1+T2−c2T−d1DN−KNzMN−eN,D˙=s2−f1L+d2N−d3TD−KDzMD−gD,L˙=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL+vLt,M˙=vMt−d4M,I˙=vIt−d5I.Figure1 illustrates the network of the dynamics for system (11). Sharp arrows represent reproduction or activation, and blocked arrows represent inhibition or killing. The blocked arrow in red represents the nonlinear interaction. Note that the suppressing effects of the drug are not shown explicitly in the figure but included in the model.Figure 1
Network of the dynamics for system (11).


## 2.1. Modeling Tumor Cells

We begin with modeling tumor cellsT which are assumed to have a proliferation rate that can be modeled by a logistic growth law aT1−bT, with parameters a and b denoting the per capita growth rates and reciprocal carrying capacities of the tumor cells [8, 9, 26]. Also, it is known that the growth of the tumor cells is impacted by three different competitive interactions including interactions between tumor cells and dendritic cells, interactions between tumor cells and natural killer cells, and interactions between tumor cells and CD8+ T cells [26–28]. Denoting the corresponding competition rates as j¯,c1,k, respectively, introduces the competition term as −c1N+jD+kLT. The dynamics of tumor cells can then be described by the following ordinary differential equation:(3)dTdt=aT1−bT−c1N+jD+kLT−KTzMT.


## 2.2. Modeling Natural Killer Cells

To model natural killer (NK) cells, we will assume that these cells have a constant sources1 as well as a NK cell recruitment term that can be represented through a modified Michaelis–Menten term (commonly used to govern cell-cell interactions):(4)g1⋅T2h1+T2⋅N,where g1 denotes the maximum NK cell recruitment rate by tumor cells and h1 denotes the steepness coefficient of the NK cell recruitment curve [8].Next, the growth of NK cells will be impacted by two different interactions, namely, the interaction between NK cells and tumor cells [29] and the interaction between NK cells and dendritic cells [30–33]. We also introduce parameters c2,d1 to be the rates of killing (due to tumor cells) and proliferation (due to dendritic cells) of NK cells, respectively. The governing differential equation for the dynamics of NK cells then can be described as(5)dNdt=s1+g1NT2h1+T2−c2T−d1DN−KNzMN−eN.Note that we have also included a natural death of NK cells through−eN.


## 2.3. Modeling Dendritic Cells

Dendritic cells play an important role in the immune system response and in controlling tumor growth. Also known as antigen-presenting cells, they update and present antigens toCD8+ T cells. Some of the earlier models [8, 13] in the literature have not incorporated the dynamics of dendritic cells in directly suppressing tumor growth, stimulating resting NK cells, and impacting the dynamics of CD8+ T cells. There is, however, experimental evidence that dendritic cells play an important role in modeling tumor immunotherapy [28].To study the dynamics of dendritic cells, we will assumes2 to be a constant source of dendritic cells, d2 to be the rate at which NK cells kill dendritic cells, d3 to be a proliferation rate of dendritic cells due to tumor cells, f1 to be the rate corresponding to the interaction of dendritic cells with CD8+ T cells, and g to be the natural death rate of dendritic cells. We then have(6)dDdt=s2−f1L+d2N−d3TD−KDzMD−gD.


## 2.4. Modeling Cytotoxic CD8+ T Cells

Among many factors that impact the growth of tumor cells, it is well known thatCD8+ T cells are an important component of the immune system that kills tumor cells. It has been seen that active tumor-specific CD8+ T cells are only present in large numbers when tumor cells are present [16, 34], and after some interactions with tumor cells, they become inactive [29]. It has been observed that mature CD8+ T cells can remove dendritic cells [35, 36].In our model, to describe the dynamics of theCD8+ T-cells, we will consider f2 to be the rate of interaction between dendritic cells and tumor cells to activate CD8+ T cells; −hLT denotes the form of competition between CD8+ T cells and tumor cells, and −iL denotes the natural death rate of CD8+ T cells.It has also been seen thatCD8+ T cells may be recruited by the debris from tumor cells lysed by NK cells [37]. To include this effect, we will incorporate in our model a recruitment term that is proportional to the number of cells killed, which is denoted as r1NT. We will also need an additional term that helps describe the regulation and suppression of CD8+ T-cell activity when there are high levels of activated CD8+ T cells without responsiveness to cytokines in the system [38, 39]. This term is denoted by uNL2. We also include CD8+ T activation by IL-2 immunotherapy which is in the form of a drug that influences the immune system’s efficacy and described via a Michaelis–Menten interaction term [16]:(7)pILIgI+I.Here,I refers to the immunotherapy drug concentration in the bloodstream. The model for the dynamics for the CD8+ T cell growth population becomes(8)dLdt=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL.


## 2.5. Modeling Drug and Vaccine Interventions

In this study, we incorporate a variety of external intervention treatment options including tumor-infiltrating lymphocyte (TIL) injections as well as chemotherapy and immunotherapy drugs. TIL drug intervention may be thought of as an immunotherapy approach in which theCD8+ T-cells are promoted through antigen-specific cytolytic immune cells. We do this by adding the term vL=vLt in equation (8) and we have(9)dLdt=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL+vL.To include the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream as follows:(10)dMdt=vMt−d4M,dIdt=vIt−d5I.The drug intervention terms in these equations reflect the amount of chemotherapy and immunotherapy drug given over time. Note that we assume that the chemotherapy and immunotherapy drugs will be eliminated from the body over time at a rate proportional to its concentration, and these are given byd4M and d5I, respectively.


## 2.6. Overall Model

From Sections2.1–2.5, we have the following overall model:(11)T˙=aT1−bT−c1N+jD+kLT−KTzMT,N˙=s1+g1NT2h1+T2−c2T−d1DN−KNzMN−eN,D˙=s2−f1L+d2N−d3TD−KDzMD−gD,L˙=f2DT−hLT−uNL2+r1NT+pILIgI+I−KLzML−iL+vLt,M˙=vMt−d4M,I˙=vIt−d5I.Figure1 illustrates the network of the dynamics for system (11). Sharp arrows represent reproduction or activation, and blocked arrows represent inhibition or killing. The blocked arrow in red represents the nonlinear interaction. Note that the suppressing effects of the drug are not shown explicitly in the figure but included in the model.Figure 1
Network of the dynamics for system (11).


## 3. Stability Analysis

In this section, we employ mathematical analysis to identify conditions that can help eliminate tumor cells. Also, we will determine conditions for when tumor-free equilibrium is unstable and the tumor grows without bound.We now consider the system of (11) in the absence of treatment. When we eliminate chemotherapy and immunotherapy, the system reduces to a four-population system of ODEs. Let ET∗,N∗,D∗,L∗ be an equilibrium point of the system described by the system without drug intervention.At an equilibrium point, we have(12)dTdt=dNdt=dDdt=dLdt=0.Since we assume there is constant recruitment through source termss1 and s2, both not equal to zero, there is no trivial equilibrium which implies(13)ET∗,N∗,D∗,L∗≠0,0,0,0.For tumor-free equilibrium, at an equilibrium point, we havedN/dt=0 which yields(14)s1+d1D∗N∗−eN∗=0,which yields(15)N∗=s1e−d1D∗.Similarly, settingdD/dt=0 at an equilibrium point yields(16)s2−d2N∗D∗−gD∗=0.Substituting (15) in (16), we get(17)gd1D∗2−s2d1+d2s1+egD∗+es2=0.Solving the quadratic equation yields(18)D1,2∗=d1s2+d2s1+eg±d1s2+d2s1+eg2−4ges22gd1.Hence, we have 2 tumor-free equilibrium given byE10,N∗,D1∗,0 and E20,N∗,D2∗,0.For these to have biological meaning, we need(19)e−d1D∗>0,(20)d1s2+d2s1+eg≥2ges2.These conditions suggest critical values for the death rate and the source term for NK cells to be(21)e=d1D∗,s1=2ges2−d1s2+egd2,in order for the tumor-free equilibrium point to be positive for biological significance.The Jacobian matrix for linearization of system (11) without drug intervention is given by(22)a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44,where(23)a11=a−2abT∗−c1N∗−jD∗−kL∗,a12=−c1T∗,a13=−jT∗,a14=−kT∗,a21=2g1N∗h1T∗h1+T2−c2N∗,a22=−c2T∗+d1D∗−e,a23=d1N∗,a24=0,a31=d3D∗,a32=−d2D∗,a33=−f1L∗+d2N∗−d3T∗+g,a34=−f1D∗,a41=f2D∗−hL∗+r1N∗,a42=−uL∗2+r1T∗,a43=f2T∗,a44=−hT∗+2uN∗L∗+i.Evaluating these terms at the general tumor-free equilibrium point gives(24)A=a−c1N∗−jD∗0002g1N∗h1T∗h1+T2−c2N∗d1D∗−ed1N∗0d3D∗−d2D∗−d2N∗+g−f1D∗f2D∗+r1N∗00−i.To solve for the eigenvaluesλ, we solve the equation detA−λI=0 or(25)deta−c1N∗−jD∗−λ0002g1N∗h1T∗h1+T2−c2N∗d1D∗−e−λd1N∗0d3D∗−d2D∗−d2N∗+g−λ−f1D∗f2D∗+r1N∗00−i−λ=0,which yields(26)a−c1N∗−jD∗−λ−i−λdetB−λI=0,where matrix B is given by(27)B=d1D∗−ed1N∗−d2D∗−d2N∗+g.It may be noted thattrace and determinant for matrix B can be computed to be(28)trB=d1D∗−e−d2N∗+g,(29)detB=d1d2N∗D∗−d1D∗−ed2N∗+g.Solving (26), we compute the eigenvalues to be(30)λ1=a−c1N∗−jD∗,λ2=−i,and λ3 and λ4 to be the roots of the equation:(31)λ2−λd1D∗−e−d2N∗+g+d1d2N∗D∗−d1D∗−ed2N∗+g=0.From (28) and (29), this can be rewritten as(32)λ2−trBλ+detB=0.Using (19) in equations (28) and (29), we can conclude that(33)trB<0,detB>0.This implies that the eigenvaluesλ3 and λ4 that are the roots of equation (32) have negative real parts.In order for the tumor-free equilibrium to be stable, we requireλ1<0, which implies that if the tumor growth rate a is lesser than the critical value given by c1N∗+JD∗, then the tumor population can be eliminated.


## 4. Computational Experiments

In this section, we will consider the system of (11) which will be solved via the Runge–Kutta methods. The parameter values, their units, and their estimated value are itemized next, which we use in our computations.(i)
a: tumor growth rate estimated as 4.31∗10−1day−1 [9](ii)
b: b−1 tumor-carrying capacity estimated as 2.17∗10−8cells−1 [9](iii)
c1: NK cell tumor cell kill rate estimated as 3.5∗10−6cells−1 [12](iv)
c2: NK cell inactivation rate by tumor cells estimated as 1.0∗10−7cells−1day−1 [9](v)
d1: rate of dendritic cell priming NK cells estimated as 1.0∗10−6cells−1 [12](vi)
d2: NK cell dendritic cell kill rate estimated as 4.0∗10−6cells−1 [12](vii)
d3: rate of tumor cells priming dendritic cells estimated as 1.0∗10−4 Estimate(viii)
e: death rate of NK cell estimated as 4.12∗10−2day−1 [12](ix)
f1: CD8+ T cell dendritic cells kill rate estimated as 1.0∗10−8cells−1 [12](x)
f2: rate of dendritic cells priming CD8+ T cell estimated as 0.01cells−1 [12](xi)
g: death rate of dendritic cells estimated as 2.4∗10−2cells−1 [12](xii)
h: CD8+ T inactivation rate by tumor cells estimated as 3.42∗10−10cells−1day−1 [12](xiii)
i: death rate of CD8+ T cells estimated as 2.0∗10−2day−1 [9](xiv)
j: dendritic cell tumor cell kill rate estimated as 1.0∗10−7cells−1 [12](xv)
k: NK cell tumor cell kill rate estimated as 1.0∗10−7cells−1 [12](xvi)
s1: source of NK cells estimated as 1.3∗104cells−1 [12](xvii)
s2: source of dendritic cell estimated as 4.8∗102cells−1 [12](xviii)
u: regulatory function by Nk cells of CD8+ T cells estimated as 1.80∗10−8cell−2day−1 [9]For the first computation, we will assume there is no additional recruitment terms forCD8+ T cells and NK cells (r1=0,g1=0,h1=0), removing some of the regulation, suppression. and activation of CD8+ T cells (pI=gI=u=0), not including the influence of drug kill terms (KT=KN=KD=KL=0), no influence of drug and vaccine interventions (vL=vM=vI=0) along with the corresponding death rates (d4=d5=0). We will also assume d3=0 which corresponds to the new term that has been added to the model to indicate the growth of dendritic cells being impacted by tumor cells. We will also assume for simplicity and illustration purposes of a weak immune system that the dynamics start with 100 tumor cells with one natural killer, one dendritic, and one CD8+ T cell. Figure 2 shows the dynamics of each of these cells. The tumor cells initially increase to a peak before a full immune clearance starts.Figure 2
Dynamics of (a) tumor, (b) NK, (c) dendritic, and (d)CD8+ T cells.
(a)
(b)
(c)
(d)Next, we consider the effect of one of the terms in system (11) corresponding to the dynamics of dendritic cells. This is the term d3TD which includes the influence of tumor growth on the dynamics of dendritic cells that all the previous studies have not considered. Figure 3 illustrates how not only the dendritic cells are impacted but the CD8+ T cell dynamics also changes as the proliferation rate d3 is doubled. For the rest of the simulations, we will include the effect of d3 and assume the value to be 1×10−4.Figure 3
Dynamics of (a) tumor, (b) NK, (c) dendritic, and (d)CD8+ T cells as the proliferation rate d3 is doubled.
(a)
(b)
(c)
(d)Along withd3, we also wanted to study the influence of the source term s2 on the dynamics of tumor, NK, and CD8+ T cells. Figure 4 illustrates this behavior. When we increase the source term of dendritic cells, it increases NK cells and CD8+ T cells. We also note that these have an effect on tumor growth as both NK and CD8+ T cells can lyse tumor cells. This decrease is also seen in this figure. This suggests that an external source term of dendritic cells has the potential to decrease tumor growth. We also note that dendritic cells play an important role in recruiting CD8+ T cells earlier in the tumor growth phase.Figure 4
Dynamics of (a) tumor, (b) NK, and (c)CD8+ T cells as function of source s2.
(a)
(b)
(c)Next, to study the effect of TIL drug intervention term only for theCD8+ T-cell population as an immunotherapy where the immune cell levels are boosted by the addition of antigen-specific cytolytic immune cells, we increase the value of vL from 1 to 106. The result is shown in Figure 5, and we notice that the effect of adding the drug in small doses does not have a big impact on tumor growth.Figure 5
Dynamics of tumor cells as the immunotherapy TIL drug intervention termvL is varied.Next, we consider the effect of the chemotherapy drug only that is introduced through the termvM. We set vM=1 and study the influence of increasing KT in the dynamics of tumor cells. Figure 6 illustrates how tumor cells can be reduced through this technique.Figure 6
Dynamics of tumor cells as immunotherapy drug intervention with constantvM and varying levels of chemotherapy drug kill terms.We want to point out that we also performed the study on the influence of immunotherapy drug interventionvI but noticed that even with inclusion of a CD8+ T activation described via Michaelis–Menten interaction given by(34)pILIgI+I,there was negligible effect on the dynamics of all four cells. We also noted that there was not much effect in the dynamics of NK cells through the recruitment terms involving g1 and r1.Next, we turn our attention to the effect of the nonlinear term introduced as an inactivation term, which describes the regulation and suppression ofCD8+ T-cell activity [38, 39]. Specifically, Figure 7 shows the effect of the nonlinear term in system (11) for parameters u=0 and u=3×10−10, and the dynamics show that there is a drastic drop in the number of cytotoxic CD8+ T cells while a small increase in tumor cell growth.Figure 7
Influence of nonlinearity on dynamics of (a) tumor and (b)CD8+ T cells.
(a)
(b)In summary, our computations seem to suggest that a combination of immunotherapy through TIL drug interventionvM along with chemotherapy through vL provides an optimal way to reduce tumor growth. Taking this into account and removing terms that had negligible effects, one can consider the following simplified system that captures most prominent features:(35)T˙=aT1−bT−c1N+jD+kLT−KTzMT,N˙=s1−c2T−d1DN−KNzMN−eN,D˙=s2−f1L+d2N−d3TD−KDzMD−gD,L˙=f2DT−hLT−uNL2−KLzML−iL+vLt,M˙=vMt−d4M.Figure8 clearly shows that combined chemotherapy and immunotherapy drug intervention helps reduce the tumor growth for system (35). We have used KT=9×10−2,KD=KN=KL=6×10−2 with vL=106 and vM=1 for our computations.Figure 8
Influence of no drug intervention, independent drug interventions, and combined drug interventions.


## 5. Parameter Estimation

In this section, we focus on estimating some parameters used in system (35), based on the measurements of tumor cells. Our goal is to accurately describe the dynamics of tumor growth on an individual basis which is very important both for growth prediction and designing personalized, optimal therapy schemes (e.g., when using model predictive control). To demonstrate this, let us consider two of the parameters in the model, namely, c1 which is the competition rate that impacts the dynamics of tumor cells due to natural killer cells and d3 which is the parameter-related proliferation of dendritic cells due to tumor cells. Recalling the influence of doubling the latter parameter d3 is illustrated in Figure 3.The purpose of parameter estimation is to identify values of parameters for given experimental data. In this work, we will demonstrate how to check the reliability of a mathematical model to estimate parameters optimally. For this, we consider a discrete dataset for tumor dynamics corresponding to the values ofc1=3.5×10−6 and d3=1×10−4 as indicated from the literature (see Figure 9). We then introduce randomness in the data by adding some Gaussian noise to each value in the tumor dynamics. We refer to the latter as the experimental data Tdata (see Figure 10).Figure 9
Experimental data for tumor cells forc1=3.5×10−6 and d3=1×10−4 with some noise.Figure 10
Parameter estimation description.Next, we make a guess for values ofc1 and d3 which are very different from the values used to create the experimental data and try to solve the ODE system (35) to obtain the computed values of the tumor dynamics given by Tc1,d3.We then set up an error expressionEc1,d3 that is the sum of the squared differences between the computed values Tc1,d3 and the experimental data Tdata as shown in equation (36). Employing an unconstrained nonlinear optimization algorithm such as the Nelder–Mead algorithm, the minimization algorithm searches for a local minimum using a regression approach. This direct search method attempts to minimize a function of real variables using only function evaluations without any derivatives. If the error Ec1,d3 is within a user-prescribed tolerance TOL, we stop and accept the values of c1 and d3 to be the optimal values. If not, we use the updated values of parameters c1 and d3 as the new values and iterate them back to solve the ODE system (35). We continue this until convergence. This parameter estimation approach is summarized in Figure 10, and the minimized objective function is given by(36)Ec1,d3=∑i−1NTc1,d3−Tdata2,where Ec1,d3, the least squared error, denotes the differences in the amount of computed tumor cells from the simulation Tc1,d3 from the observed data Tdata (Figure 9) over N observations.Using same initial conditions with poor guesses forc1=1×10−8 and d3=1×10−2, the optimization algorithm estimates the parameters to be c1=3.5×10−6 and d3=0.45×10−3. Clearly, the algorithm estimates the values pretty close to the values used originally to create the experimental dataset, and the predicted dynamics of the tumor cells for these parameter values along with comparison to the experimental data is illustrated in Figure 11.Figure 11
Prediction of the dynamics of tumor cells through estimated values forc1=3.5×10−6 and d3=0.45×10−3 in solid line plotted with the data.


## 6. Discussion and Conclusion

In this work, we developed a mathematical model that incorporated the dynamics of four coupled cell populations including tumor cells, natural killer cells, dendritic cells, and cytotoxicCD8+ T cells that influence the growth of tumors. The novelty in the model was how it combined important interactions between growing tumor cells and cells of the innate and specific immune system coupled with models for drug delivery to these cell sites. A detailed stability analysis of the associate ordinary differential equation system was performed. A variety of computational experiments were simulated to study the influence of the dynamics of the four cell populations on various parameters. For example, we noted a significant effect of the influence of proliferation rate d3 on the dynamics of the cell populations which was neglected in past studies. We noted that the dynamics are a function of the source terms included in the model. The effect of TIL drug intervention as an immunotherapy approach had significant impact on tumor growth. Similar results were observed for a chemotherapy approach as well. Clearly, a combined chemotherapy and immunotherapy drug intervention approach was seen to reduce tumor growth greatly. Another feature of the work is the application of a parameter estimation algorithm to accurately predict proliferation parameters for a given set of data of tumor cell growth. Similar algorithms have been used in the past to characterize properties of soft tissues [40]. We were not only able to capture the data accurately but also were able to accurately quantify the parameters related to the data.While this paper investigates a system of ordinary differential equations, one must computationally study the corresponding partial differential equations (PDES) presented as we developed the model along with fluid equations that help the drugs to move towards the cancer cells. This will require the use of sophisticated numerical methods like the finite element methods to solve the associated system of PDEs. This will be the focus of a forthcoming paper.Also, we hope to extend our work to apply the models and validate them against actual experimental or laboratory data along with applying machine learning type algorithms to predict behaviors of the growth of the tumor cells. These predictions can help develop control mechanisms such as drug therapy. This will also be a focus of our future work.



---

*Source: 4079298-2019-10-08.xml*