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

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 , tumor size oscillations and regressions [15-18], non-spatial 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 cytotoxic [CD8.sup.+] 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.sup.+] 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, and [CD8.sup.+] 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 (T(t)), natural killer cells (N (t)), dendritic cells (D(t)), and cytotoxic [CD8.sup.+] T cells denoted by (L(t)). 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:

[partial derivative][*]/[partial derivative]t + [nabla] x ([??][*])-[[delta].sub.D][[nabla].sup.2][*] = f(*)-g(*)- [K.sub.[*]]z(M)[*]. (1)

Here, the functions f (*) 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 through [K.sub.[*]]Z(M)[*]. The term z(M) = 1 - [e.sup.-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.sub.[*]] 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.sub.[*]] = 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 [??], we will consider only the temporal dynamics in this work and hence the associated ordinary differential equation:

d[*]/dt = f(*)-g(*)- [K.sub.[*]]Z(M)[*]. (2)

2.1. Modeling Tumor Cells. We begin with modeling tumor cells T which are assumed to have a proliferation rate that can be modeled by a logistic growth law aT(1 - 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.sup.+] T cells [26-28]. Denoting the corresponding competition rates as [j.bar], [c.sub.1],k, respectively, introduces the competition term as -([c.sub.1]N + j D + kL)T. The dynamics of tumor cells can then be described by the following ordinary differential equation:

dT/dt = aT(1-bT)-([c.sub.1]N+jD+kL)T - [K.sub.T]z(M)T. (3)

2.2. Modeling Natural Killer Cells. To model natural killer (NK) cells, we will assume that these cells have a constant source s1 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):

[g.sub.1] x [T.sup.2]/[h.sub.1]+[T.sup.2] x N, (4)

where [g.sub.1] denotes the maximum NK cell recruitment rate by tumor cells and [h.sub.1] denotes the steepness coefficient of the NK cell recruitment curve .

Next, the growth of NK cells will be impacted by two different interactions, namely, the interaction between NK cells and tumor cells  and the interaction between NK cells and dendritic cells [30-33]. We also introduce parameters [c.sub.2], [d.sub.1] 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

dN/dt = [s.sub.1] + [g.sub.1][NT.sup.2]/[h.sub.1]+[T.sup.2] - ([c.sub.2]T-[d.sub.1]D)N - [K.sub.N]z(M)N - eN. (5)

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 to [CD8.sup.+] 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.sup.+] T cells. There is, however, experimental evidence that dendritic cells play an important role in modeling tumor immunotherapy .

To study the dynamics of dendritic cells, we will assume [s.sub.2] to be a constant source of dendritic cells, [d.sub.2] to be the rate at which NK cells kill dendritic cells, [d.sub.3] to be a proliferation rate of dendritic cells due to tumor cells, [f.sub.1] to be the rate corresponding to the interaction of dendritic cells with [CD8.sup.+] T cells, and g to be the natural death rate of dendritic cells. We then have

dD/dt = [s.sub.2] - ([f.sub.1]L + [d.sub.2]N - [d.sub.3]T)D - [K.sub.D]z (M)D - gD. (6)

2.4. Modeling Cytotoxic [CD8.sup.+] T Cells. Among many factors that impact the growth of tumor cells, it is well known that [CD8.sup.+] T cells are an important component of the immune system that kills tumor cells. It has been seen that active tumor-specific [CD8.sup.+] 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 . It has been observed that mature [CD8.sup.+] T cells can remove dendritic cells [35, 36].

In our model, to describe the dynamics of the [CD8.sup.+] T-cells, we will consider [f.sub.2] to be the rate of interaction between dendritic cells and tumor cells to activate [CD8.sup.+] T cells; -hLT denotes the form of competition between [CD8.sup.+] T cells and tumor cells, and -iL denotes the natural death rate of [CD8.sup.+] T cells.

It has also been seen that [CD8.sup.+] T cells may be recruited by the debris from tumor cells lysed by NK cells . 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 [r.sub.1]NT. We will also need an additional term that helps describe the regulation and suppression of [CD8.sup.+] T-cell activity when there are high levels of activated [CD8.sup.+] T cells without responsiveness to cytokines in the system [38, 39]. This term is denoted by [uNL.sup.2]. We also include [CD8.sup.+] 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 :

[p.sub.I]LI/[g.sub.I]+I. (7)

Here, I refers to the immunotherapy drug concentration in the bloodstream. The model for the dynamics for the [CD8.sup.+] T cell growth population becomes

dL/dt = [f.sub.2]DT - hLT - [uNL.sup.2] + [r.sub.1]NT + [p.sub.I]LI/[g.sub.I]+I - [K.sub.L]z(M)L - iL. (8)

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 the [CD8.sup.+] T-cells are promoted through antigen-specific cytolytic immune cells. We do this by adding the term [v.sub.L] = [v.sub.L](t) in equation (8) and we have

dL/dt = [f.sub.2]DT - hLT - [uNL.sup.2] + [r.sub.1]NT + -[p.sub.I]LI/[g.sub.I]+I - [K.sub.L]z(M)L - iL + [v.sub.L]. (9)

To include the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream as follows:

dM/dt = [v.sub.M](t) - [d.sub.4]M, dI/dt = [v.sub.I](t) - [d.sub.5]I. (10)

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 by [d.sub.4]M and [d.sub.5]I, respectively.

2.6. Overall Model. From Sections 2.1-2.5, we have the following overall model:

[mathematical expression not reproducible] (11)

Figure 1 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.

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 E([T.sup.*], [N.sup.*], [D.sup.*], [L.sup.*]) be an equilibrium point of the system described by the system without drug intervention.

At an equilibrium point, we have

dT/dt = dN/dt = dD/dt = dL/dt = 0. (12)

Since we assume there is constant recruitment through source terms [s.sub.1] and [s.sub.2], both not equal to zero, there is no trivial equilibrium which implies

E([T.sup.*],[N.sup.*],[D.sup.*],[L.sup.*]) [not equal to] (0,0,0,0). (13)

For tumor-free equilibrium, at an equilibrium point, we have dN/dt = 0 which yields

[s.sub.1] + [d.sub.1][D.sup.*][N.sup.*] - e[N.sup.*] = 0, (14)

which yields

[N.sup.*] = [s.sub.1]/e-[d.sub.1][D.sup.*]. (15)

Similarly, setting dD/dt = 0 at an equilibrium point yields

[s.sub.2] - [d.sub.2][N.sup.*][D.sup.*] - g[D.sup.*] = 0. (16)

Substituting (15) in (16), we get

[gd.sub.1][D.sup.*2] - ([s.sub.2][d.sub.1] + [d.sub.2][s.sub.1] + eg)[D.sup.*] + [es.sub.2] = 0. (17)

[D.sup.*.sub.1,2] = ([d.sub.1][S.sub.2] + [d.sub.2][s.sub.1] + eg) [+ or -] [square root of ([([d.sub.1][s.sub.2] + [d.sub.2][s.sub.1] + eg).sup.2] - 4[ges.sub.2])]/2[gd.sub.1] (18)

Hence, we have 2 tumor-free equilibrium given by [E.sub.1](0, [N.sup.*], [D.sup.*.sub.1], 0) and [E.sub.2](0,[ N.sup.*], [D.sup.*.sub.2], 0).

For these to have biological meaning, we need

e - [d.sub.1][D.sup.*] > 0, (19)

[d.sub.1][s.sub.2] + [d.sub.2][s.sub.1] + eg [greater than or equal to] 2[square root of [ges.sub.2]]. (20)

These conditions suggest critical values for the death rate and the source term for NK cells to be

e = [d.sup.1][D.sup.*], [s.sub.1] = 2[square root of [ges.sub.2]] - ([d.sub.1][s.sub.2] + eg)/[d.sub.2], (21)

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

[mathematical expression not reproducible] (22)

where

[mathematical expression not reproducible] (23)

Evaluating these terms at the general tumor-free equilibrium point gives

[mathematical expression not reproducible] (24)

To solve for the eigenvalues A, we solve the equation det (A - [lambda]I) = 0 or

[mathematical expression not reproducible] (25)

which yields

(a - [c.sub.1][N.sup.*] - j[D.sup.*] - [lambda]) (-i - [lambda])det (B - [lambda]I) = 0, (26)

where matrix B is given by

[mathematical expression not reproducible] (27)

It may be noted that trace and determinant for matrix B can be computed to be

tr(B) = ([d.sup.1][D.sup.*] - e) - ([d.sub.2][N.sup.*] + g), (28)

det (B) = [d.sub.1][d.sub.2][N.sup.*][D.sup.*] - ([d.sub.1][D.sup.*] - e) ([d.sub.2][N.sup.*] + g). (29)

Solving (26), we compute the eigenvalues to be

[[lambda].sub.1] = a - [c.sub.1][N.sup.*] - [jD.sup.*], [[lambda].sub.2] = -i, (30)

and [[lambda].sub.3] and [[lambda].sub.4] to be the roots of the equation:

[[lambda].sup.2] - A[([d.sub.1][D.sup.*] - e) - ([d.sup.2][N.sup.*] + g)] + [d.sub.1][d.sub.2][N.sup.*][D.sup.*] -([d.sub.1][D.sup.*] - e)([d.sub.2][N.sup.*] + g) = 0.

From (28) and (29), this can be rewritten as

[[lambda].sup.2] - tr(B)[lambda] + det(B) = 0. (32)

Using (19) in equations (28) and (29), we can conclude that

tr(B) < 0, det(B) > 0. (33)

This implies that the eigenvalues [[lambda].sub.3] and [[lambda].sub.4] that are the roots of equation (32) have negative real parts.

In order for the tumor-free equilibrium to be stable, we require [[lambda].sub.1] < 0, which implies that if the tumor growth rate a is lesser than the critical value given by [c.sub.1][N.sup.*] + [JD.sup.*], 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.sup.-1] [day.sup.-1] 

(ii) b: [b.sup.-1] tumor-carrying capacity estimated as 2.17 * [10.sup.-8] [cells.sup.-1] 

(iii) [c.sub.1]: NK cell tumor cell kill rate estimated as 3.5 * [10.sup.-6] [cells.sup.-1] 

(iv) [c.sub.2]: NK cell inactivation rate by tumor cells estimated as 1.0 * [10.sup.-7] [cells.sup.-1] [day.sup.-1] 

(v) [d.sub.1]: rate of dendritic cell priming NK cells estimated as 1.0 * [10.sup.-6] [cells.sup.-1] 

(vi) [d.sub.2]: NK cell dendritic cell kill rate estimated as 4.0 * [10.sup.-6] [cells.sup.-1] 

(vii) [d.sub.3]: rate of tumor cells priming dendritic cells estimated as 1.0 * [10.sup.-4] Estimate

(viii) e: death rate of NK cell estimated as 4.12 * [10.sup.-2] [day.sup.-1] 

(ix) [f.sub.1]: [CD8.sup.+] T cell dendritic cells kill rate estimated as 1.0 * [10.sup.-8] [cells.sup.-1] 

(x) [f.sub.2]: rate of dendritic cells priming [CD8.sup.+] T cell estimated as 0.01 [cells.sup.-1] 

(xi) g: death rate of dendritic cells estimated as 2.4 * [10.sup.-2] [cells.sup.-1] 

(xii) h: [CD8.sup.+] T inactivation rate by tumor cells estimated as 3.42 * [10.sup.-10] [cells.sup.-1] [day.sup.-1] 

(xiii) i: death rate of [CD8.sup.+] T cells estimated as 2.0 * [10.sup.-2] [day.sup.-1] 

(xiv) j: dendritic cell tumor cell kill rate estimated as 1.0 * [10.sup.-7] [cells.sup.-1] 

(xv) k: NK cell tumor cell kill rate estimated as 1.0 * [10.sup.-7] [cells.sup.-1] 

(xvi) [s.sub.1]: source of NK cells estimated as 1.3 * [10.sup.4] [cells.sup.-1] 

(xvii) [s.sub.2]: source of dendritic cell estimated as 4.8 * [10.sup.2] [cells.sup.-1] 

(xviii) u: regulatory function by Nk cells of [CD8.sup.+] T cells estimated as 1.80 * [10.sup.-8] [cell.sup.-2] [day.sup.-1] 

For the first computation, we will assume there is no additional recruitment terms for [CD8.sup.+] T cells and NK cells ([r.sub.1] = 0,[g.sub.1] = 0,[h.sub.1] = 0), removing some of the regulation, suppression. and activation of [CD8.sup.+] T cells ([p.sub.I] = [g.sub.I] = u = 0), not including the influence of drug kill terms ([K.sub.T] = [K.sub.N] = [K.sub.D] = [K.sub.L] = 0), no influence of drug and vaccine interventions ([v.sub.L] = [v.sub.M] = [v.sub.I] = 0) along with the corresponding death rates ([d.sub.4] = [d.sub.5] = 0). We will also assume [d.sub.3] = 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.sup.+] 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.

Next, we consider the effect of one of the terms in system (11) corresponding to the dynamics of dendritic cells. This is the term [d.sub.3]TD 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.sup.+] T cell dynamics also changes as the proliferation rate [d.sub.3] is doubled. For the rest of the simulations, we will include the effect of [d.sub.3] and assume the value to be 1 x [10.sup.-4].

Along with [d.sub.3], we also wanted to study the influence of the source term [s.sub.2] on the dynamics of tumor, NK, and [CD8.sup.+] T cells. Figure 4 illustrates this behavior. When we increase the source term of dendritic cells, it increases NK cells and [CD8.sup.+] T cells. We also note that these have an effect on tumor growth as both NK and [CD8.sup.+] 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.sup.+] T cells earlier in the tumor growth phase.

Next, to study the effect of TIL drug intervention term only for the [CD8.sup.+] 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 [v.sub.L] from 1 to [10.sup.6]. 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.

Next, we consider the effect of the chemotherapy drug only that is introduced through the term [v.sub.M]. We set [v.sub.M] = 1 and study the influence of increasing [K.sub.T] in the dynamics of tumor cells. Figure 6 illustrates how tumor cells can be reduced through this technique.

We want to point out that we also performed the study on the influence of immunotherapy drug intervention [v.sub.I] but noticed that even with inclusion of a [CD8.sup.+] T activation described via Michaelis-Menten interaction given by

[p.sub.I]LI/[g.sub.I]+I, (34)

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 [g.sub.1] and [r.sub.1].

Next, we turn our attention to the effect of the nonlinear term introduced as an inactivation term, which describes the regulation and suppression of [CD8.sup.+] 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 x [10.sup.-10], and the dynamics show that there is a drastic drop in the number of cytotoxic [CD8.sup.+] T cells while a small increase in tumor cell growth.

In summary, our computations seem to suggest that a combination of immunotherapy through TIL drug intervention vM along with chemotherapy through [v.sub.L] 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:

[mathematical expression not reproducible] (35)

Figure 8 clearly shows that combined chemotherapy and immunotherapy drug intervention helps reduce the tumor growth for system (35). We have used [K.sub.T] = 9 x [10.sup.-2], [K.sub.D] = [K.sub.N] = [K.sub.L] = 6 x [10.sup.-2] with [v.sub.L] = [10.sup.6] and [v.sub.M] = 1 for our computations.

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, [c.sub.1] which is the competition rate that impacts the dynamics of tumor cells due to natural killer cells and [d.sub.3] which is the parameter-related proliferation of dendritic cells due to tumor cells. Recalling the influence of doubling the latter parameter [d.sub.3] 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 of [c.sub.1] = 3.5 x [10.sup.-6] and [d.sub.3] = 1 x [10.sup.-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 [T.sub.data] (see Figure 10).

Next, we make a guess for values of [c.sub.1] and [d.sub.3] 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 T ([c.sub.1], [d.sub.3]).

We then set up an error expression E([c.sub.1],[d.sub.3]) that is the sum of the squared differences between the computed values T([c.sub.1],[d.sub.3]) and the experimental data [T.sub.data] 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 E([c.sub.1],[d.sub.3]) is within a user-prescribed tolerance TOL, we stop and accept the values of [c.sub.1] and [d.sub.3] to be the optimal values. If not, we use the updated values of parameters [c.sub.1] and [d.sub.3] 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

E([c.sub.1][d.sub.3]) = [N.summation over (i-1)][(T([c.sub.1],[d.sub.3]) - [T.sub.data]).sup.2], (36)

where E([c.sub.1],[d.sub.3]), the least squared error, denotes the differences in the amount of computed tumor cells from the simulation T ([c.sub.1],[d.sub.3]) from the observed data [T.sub.data] (Figure 9) over N observations.

Using same initial conditions with poor guesses for [c.sub.1] = 1 x [10.sup.-8] and [d.sub.3] = 1 x [10.sup.-2], the optimization algorithm estimates the parameters to be [c.sub.1] = 3.5 x [10.sup.-6] and [d.sub.3] = 0.45 x [10.sup.-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.

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 cytotoxic [CD8.sup.+] 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 [d.sub.3] 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 . 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.

https://doi.org/10.1155/2019/4079298

Data Availability

All data supporting the results reported have sources provided in the manuscript through references or generated during the study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

 J. N. Blattman and P. D. Greenberg, "Cancer immunotherapy: a treatment for the masses," Science, vol. 305, no. 5681, pp. 200-205, 2004.

 C. R. Parish, "Cancer immunotherapy: the past, the present and the future," Immunology and Cell Biology, vol. 81, no. 2, pp. 106-113, 2003.

 K. J. O'Byrne, A. G. Dalgleish, M. J. Browning, W. P. Steward, and A. L. Harris, "The relationship between angiogenesis and the immune response in carcinogenesis and the progression of malignant disease," European Journal of Cancer, vol. 36, no. 2, pp. 151-169, 2000.

 A. Ribas and J. D. Wolchok, "Cancer immunotherapy using checkpoint blockade," Science, vol. 359, no. 6382, pp. 1350-1355, 2018.

 T. H. Stewart, "Immune mechanisms and tumor dormancy," Medicina (Buenos Aires), vol. 56, no. 1, pp. 74-82, 1996.

 B. Ribba, C. Boetsch, T. Nayak et al., "Prediction of the optimal dosing regimen using a mathematical model of tumor uptake for immunocytokine-based cancer immunotherapy," Clinical Cancer Research, vol. 24, no. 14, pp. 3325-3333, 2018.

 P. Veeresha, D. G. Prakasha, and H. M. Baskonus, "New numerical surfaces to the mathematical model of cancer chemotherapy effect in Caputo fractional derivatives," Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 29, no. 1, Article ID 013119, 2019.

 L. G. de Pillis and A. Radunskaya, "A mathematical model of immune response to tumor invasion," in Computational Fluid and Solid Mechanics, pp. 1661-1668, Elsevier Science, Amsterdam, Netherlands, 2003.

 L. G. de Pillis, W. Gu, and A. E. Radunskaya, "Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations," Journal of Theoretical Biology, vol. 238, no. 4, pp. 841-862, 2006.

 P. Bonnotte, A. Palladini, F. Pappalardo, and S. Motta, "Predictive models in tumor immunology," Selected Topicsin Cancer Modeling, E. Angelis, M. A. J. Chaplain, and N. Bellomo, Eds., Springer Science+Business Media, LLC, Berlin, Germany, 2008.

 A. Matzavinos, M. A. J. Chaplain, and V. A. Kuznetsov, "Mathematical modelling of the spatio-temporal response of cytotoxic T-lymphocytes to a solid tumour," Mathematical Medicine and Biology, vol. 21, no. 1, pp. 1-34, 2004.

 T. Trisilowati, S. McCue, and D. Mallet, "Numerical solution of an optimal control model of dendritic cell treatment of a growing tumour," ANZIAM Journal, vol. 54, pp. 664-680, 2013.

 S. Wilson and D. Levy, "A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy," Bulletin of Mathematical Biology, vol. 74, no. 7, pp. 1485-1500, 2012.

 R. J. de Boer and P. Hogeweg, "Interactions between macrophages and T-lymphocytes: tumor sneaking through intrinsic to helper T cell dynamics," Journal of Theoretical Biology, vol. 120, no. 3, pp. 331-351, 1986.

 R. R. Gomis and S. Gawrzak, "Tumor cell dormancy," Molecular Oncology, vol. 11, no. 1, pp. 62-78, 2017.

 D. Kirschner and J. C. Panetta, "Modeling immunotherapy of the tumor--immune interaction," Journal of Mathematical Biology, vol. 37, no. 3, pp. 235-252, 1998.

 V. Kuznetsov and I. Makalkin, "Bifurcation-analysis of mathematical-model of interactions between cytotoxic lymphocytes and tumor-cells--effect of immunological amplification of tumorgrowth and its connection with other phenomena of oncoimmunology," Biofizika, vol. 37, no. 6, pp. 1063-1070, 1992.

 V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, and A. S. Perelson, "Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis," Bulletin of Mathematical Biology, vol. 56, no. 2, pp. 295-321,1994.

 F. Ansarizadeh, M. Singh, and D. Richards, "Modelling of tumor cells regression in response to chemotherapeutic treatment," Applied Mathematical Modelling, vol. 48, pp. 96-112, 2017.

 R. Eftimie, J. L. Bramson, and D. J. D. Earn, "Interactions between the immune system and cancer: a brief review of nonspatial mathematical models," Bulletin of Mathematical Biology, vol. 73, no. 1, pp. 2-32, 2011.

 S. Kato, A. Goodman, V. Walavalkar, D. A. Barkauskas, A. Sharabi, and R. Kurzrock, "Hyperprogressors after immunotherapy: analysis of genomic alterations associated with accelerated growth rate," Clinical Cancer Research, vol. 23, no. 15, pp. 4242-4250, 2017.

 F. Nani and H. I. Freedman, "A mathematical model of cancer treatment by immunotherapy," Mathematical Biosciences, vol. 163, no. 2, pp. 159-199, 2000.

 C. J. Wheeler, D. Asha, L. Gentao, J. S. Yu, and K. L. Black, "Clinical responsiveness of glioblastoma multiforme to chemotherapy after vaccination," Clinical Cancer Research, vol. 10, no. 16, pp. 5316-5326, 2004.

 National Cancer Institute, "Understanding chemotherapy," May 2005, http://www.nci.nih.gov/cancertopics/ chemotherapy-and-you/,page2.

 M. C. Perry, The Chemotherapy Source Book, Lippincott Williams & Wilkins, Philadelphia, PA, USA, 3rd edition, 2001.

 A. Diefenbach, E. R. Jensen, A. M. Jamieson, and D. H. Raulet, "Rae1 and H60 ligands of the NKG2D receptor stimulate tumour immunity," Nature, vol. 413, no. 6852, pp. 165-171, 2001.

 Y. Kawarda, R. Ganss, N. Garbi, T. Sacher, B. Arnold, and G. T. Hammerling, "NK- and [CD8.sup.+] T cells-mediate eradication of established tumors by peritumoral injection of CpG-containing oligodeoxynucleotides," Journal of Immunology, vol. 167, no. 1, pp. 5247-5253, 2001.

 N. Larmonier, J. Fraszczak, D. Lakomy, B. Bonnotte, and E. Katsanis, "Killer dendritic cells and their potential for cancer immunotherapy," Cancer Immunology, Immunotherapy, vol. 59, no. 1, pp. 1-11, 2010.

 J. A. Adam and N. Bellomo, "Basic models of tumor immune system interactions-identification, analysis and predictions," in A Survey of Models for Tumor-Immune System Dynamics, Birkhauser, Basel, Switzerland, 1997.

 M. Cooper, T. A. Fehniger, A. Fuchs, M. Colonna, and M. A. Caligiuri, "NK cell and DC interactions," Trends in Immunology, vol. 25, no. 1, pp. 47-52, 2004.

 G. Ferlazzo, M. L. Tsang, L. Moretta, G. Melioli, R. M. Steinman, and C. Muanz, "Human dendritic cells activate resting natural killer (NK) cells and are recognized via the NKp30 receptor by activated NK cells," The Journal of Experimental Medicine, vol. 195, no. 3, pp. 343-351, 2002.

 L. Fong and E. G. Engleman, "Dendritic cells in cancer immunotherapy," Annual Review of Immunology, vol. 18, no. 1, pp. 245-273, 2000.

 A. Moretta, "Natural killer cells and dendritic cells: rendezvous in abused tissues," Nature Reviews Immunology, vol. 2, no. 12, pp. 957-965, 2002.

 I. M. Tessier, J. Brostoff, and D. K. Male, Immunology, Mosby, St. Louis, MO, USA, 1993.

 F. Castiglione and B. Piccoli, "Optimal control in a model of dendritic cell transfection cancer immunotherapy," in Proceedings of the 2004 43rd IEEE Conference on Decision and Control (CDC), Maui, HI, USA, December 2004.

 Y. U. Wu, L. Xia, M. Zhang, and X. Zhao, "Immunodominance analysis through interactions of [CD8.sup.+] T cells and DCs in lymph nodes," Mathematical Biosciences, vol. 225, no. 1, pp. 53-58, 2010.

 A. Huang, P. Golumbek, M. Ahmadzadeh, E. Jaffee, D. Pardoll, and H. Levitsky, "Role of bone marrow-derived cells in presenting MHC class I-restricted tumor antigens," Science, vol. 264, no. 5161, pp. 961-965, 1994.

 S. M. Gilbertson, P. D. Shah, and D. A. Rowley, "NK cells suppress the generation of Lyt-2+ cytolytic T cells by suppressing or eliminating dendritic cells," Journal of Immunology, vol. 136, no. 10, pp. 3567-3571, 1986.

 A. V. Gett, F. Sallusto, A. Lanzavecchia, and J. Geginat, "T cell fitness determined by signal strength," Nature Immunology, vol. 4, no. 4, pp. 355-360, 2003.

 P. Seshaiyer and J. D. Humphrey, "A sub-domain inverse finite element characterization of hyperelastic membranes including soft tissues," Journal of Biomechanical Engineering, vol. 125, no. 3, pp. 363-371, 2003.

Pranav Unni [ID], (1) and Padmanabhan Seshaiyer [ID], (2)

(1) American International School Chennai, Chennai, Tamilnadu, India

(2) George Mason University, Fairfax, Virginia, USA

Received 25 April 2019; Revised 7 August 2019; Accepted 5 September 2019; Published 8 October 2019

Caption: Figure 1: Network of the dynamics for system (11).

Caption: Figure 2: Dynamics of (a) tumor, (b) NK, (c) dendritic, and (d) [CD8.sup.+] T cells.

Caption: Figure 3: Dynamics of (a) tumor, (b) NK, (c) dendritic, and (d) [CD8.sup.+] T cells as the proliferation rate [d.sub.3] is doubled.

Caption: Figure 4: Dynamics of (a) tumor, (b) NK, and (c) [CD8.sup.+] T cells as function of source [s.sub.2].

Caption: Figure 5: Dynamics of tumor cells as the immunotherapy TIL drug intervention term [v.sub.L] is varied.

Caption: Figure 6: Dynamics of tumor cells as immunotherapy drug intervention with constant [v.sub.M] and varying levels of chemotherapy drug kill terms.

Caption: Figure 7: Influence of nonlinearity on dynamics of (a) tumor and (b) [CD8.sup.+] T cells.

Caption: Figure 8: Influence of no drug intervention, independent drug interventions, and combined drug interventions.

Caption: Figure 9: Experimental data for tumor cells for [c.sub.1] = 3.5 x [10.sup.-6] and [d.sub.3] = 1 x [10.sup.-4] with some noise.

Caption: Figure 10: Parameter estimation description.

Caption: Figure 11: Prediction of the dynamics of tumor cells through estimated values for [c.sub.1] = 3.5 x [10.sup.-6] and [d.sub.3] = 0.45 x [10.sup.-3] in solid line plotted with the data.