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

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.


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]. is 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][4][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][9][10][11][12][13]. ese 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][16][17][18], nonspatial models of tumor and immune system interactions [19,20], and tumor growth coupled with immunotherapy [8,9,12,13,21,22]. ese 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. ese 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]. e 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. e 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 + 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. e 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.
is 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 + 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.

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 + T cells denoted by (L(t)). e 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. is would often yield the following equation for the dynamics of the various types of cells: (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 [·] z(M) [·]. e term z(M) � 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. e 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.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 + T cells [26][27][28]. Denoting the corresponding competition rates as j, c 1 , k, respectively, introduces the competition term as − (c 1 N + j D + kL)T. e dynamics of tumor cells can then be described by the following ordinary differential equation:

Modeling Natural Killer Cells.
To model natural killer (NK) cells, we will assume that these cells have a constant source s 1 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): where g 1 denotes the maximum NK cell recruitment rate by tumor cells and h 1 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][31][32][33]. We also introduce parameters c 2 , d 1 to be the rates of killing (due to tumor cells) and proliferation (due to dendritic cells) of NK cells, respectively. e governing differential equation for the dynamics of NK cells then can be described as Note that we have also included a natural death of NK cells through − eN.

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 + 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.
ere 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 assume s 2 to be a constant source of dendritic cells, d 2 to be the rate at which NK cells kill dendritic cells, d 3 to be a proliferation rate of dendritic cells due to tumor cells, f 1 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

Modeling Cytotoxic CD8 + T Cells.
Among many factors that impact the growth of tumor cells, it is well known that CD8 + 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 the CD8 + T-cells, we will consider f 2 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 that CD8 + 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 r 1 NT. 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]. is term is denoted by uNL 2 . 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]: Here, I refers to the immunotherapy drug concentration in the bloodstream. e model for the dynamics for the CD8 + T cell growth population becomes

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

Computational and Mathematical Methods in Medicine
To include the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream as follows: e 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 4 M and d 5 I, respectively.
2.6. Overall Model. From Sections 2.1-2.5, we have the following overall model: Figure 1 illustrates the network of the dynamics for system (11). Sharp arrows represent reproduction or activation, and blocked arrows represent inhibition or killing. e 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.

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 * , N * , D * , L * ) be an equilibrium point of the system described by the system without drug intervention.
At an equilibrium point, we have Since we assume there is constant recruitment through source terms s 1 and s 2 , both not equal to zero, there is no trivial equilibrium which implies For tumor-free equilibrium, at an equilibrium point, we have dN/dt � 0 which yields which yields Similarly, setting dD/dt � 0 at an equilibrium point yields Substituting (15) in (16), we get Solving the quadratic equation yields Hence, we have 2 tumor-free equilibrium given by E 1 (0, N * , D * 1 , 0) and E 2 (0, N * , D * 2 , 0). For these to have biological meaning, we need ese conditions suggest critical values for the death rate and the source term for NK cells to be e � d 1 D * , in order for the tumor-free equilibrium point to be positive for biological significance. e Jacobian matrix for linearization of system (11) without drug intervention is given by a 12 a 13 a 14   a 21 a 22 a 23 a 24   a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44 where Evaluating these terms at the general tumor-free equilibrium point gives To solve for the eigenvalues λ, we solve the equation where matrix B is given by It may be noted that trace and determinant for matrix B can be computed to be Solving (26), we compute the eigenvalues to be and λ 3 and λ 4 to be the roots of the equation: From (28) and (29), this can be rewritten as Using (19) in equations (28) and (29), we can conclude that is 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 c 1 N * + JD * , then the tumor population can be eliminated.

Computational Experiments
In this section, we will consider the system of (11) which will be solved via the Runge-Kutta methods.
e parameter values, their units, and their estimated value are itemized next, which we use in our computations. For the first computation, we will assume there is no additional recruitment terms for CD8 + T cells and NK cells (r 1 � 0, g 1 � 0, h 1 � 0), removing some of the regulation, suppression. and activation of CD8 + T cells (p I � g I � u � 0), not including the influence of drug kill terms (K T � K N � K D � K L � 0), no influence of drug and vaccine interventions (v L � v M � v I � 0) along with the corresponding death rates (d 4 � d 5 � 0). We will also assume d 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 + T cell. Figure 2 shows the dynamics of each of these cells. e 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. is is the term d 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 + T cell dynamics also changes as the proliferation rate d 3 is doubled. For the rest of the simulations, we will include the effect of d 3 and assume the value to be 1 × 10 − 4 .
Along with d 3 , we also wanted to study the influence of the source term s 2 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.
is decrease is also seen in this figure. is 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.
Next, to study the effect of TIL drug intervention term only for the CD8 + 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 L from 1 to 10 6 . e 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 M . We set v M � 1 and study the influence of increasing K 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 I but noticed that even with inclusion of a CD8 + T activation described via Michaelis-Menten interaction given by 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 1 and r 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 + 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.
In summary, our computations seem to suggest that a combination of immunotherapy through TIL drug intervention v M along with chemotherapy through v 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: 6 Computational and Mathematical Methods in Medicine    Computational and Mathematical Methods in Medicine Figure 8 clearly shows that combined chemotherapy and immunotherapy drug intervention helps reduce the tumor growth for system (35). We have used K T � 9 × 10 − 2 , K D � K N � K L � 6 × 10 − 2 with v L � 10 6 and v M � 1 for our computations.

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 1 which is the competition rate that impacts the dynamics of tumor cells due to natural killer cells and d 3 which is the parameter-related proliferation of dendritic cells due to tumor cells. Recalling the influence of doubling the latter parameter d 3 is illustrated in Figure 3. e 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 1 � 3.5 × 10 − 6 and d 3 � 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 T data (see Figure 10).
Next, we make a guess for values of c 1 and d 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 1 , d 3 ).
We then set up an error expression E(c 1 , d 3 ) that is the sum of the squared differences between the computed values T(c 1 , d 3 ) and the experimental data T 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.
is direct search method attempts to minimize a function of real variables using only function evaluations without any derivatives. If the error E(c 1 , d 3 ) is within a user-prescribed tolerance TOL, we stop and accept the values of c 1 and d 3 to be the optimal values. If not, we use the updated values of parameters c 1 and d 3 as the new values and iterate them back to solve the ODE system (35). We continue this until convergence. is parameter estimation approach is summarized in Figure 10, and the minimized objective function is given by where E(c 1 , d 3 ), the least squared error, denotes the differences in the amount of computed tumor cells from the simulation T(c 1 , d 3 ) from the observed data T data (Figure 9) over N observations. Using same initial conditions with poor guesses for c 1 � 1 × 10 − 8 and d 3 � 1 × 10 − 2 , the optimization algorithm estimates the parameters to be c 1 � 3.5 × 10 − 6 and Tumor cells K T = 9 × 10 -3 , K N = K D = K L = 6 × 10 -3 K T = 9 × 10 -2 , K N = K D = K L = 6 × 10 -2 K T = 9 × 10 -1 , K N = K D = K L = 6 × 10 -1 K T = 9 × 10 -4 , K N = K D = K L = 6 × 10 -4

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 + T cells that influence the growth of tumors.
e 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 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. e 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. is will require the use of sophisticated numerical methods like the finite element methods to solve the associated system of PDEs. is 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. ese predictions can help develop control mechanisms such as drug therapy. is will also be a focus of our future work.

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