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.
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 cytotoxic 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 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 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 (), natural killer cells (), dendritic cells (), and cytotoxic T cells denoted by (). 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:
Here, the functions and 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 . The term 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 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 , 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:
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 , 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 T cells [26–28]. Denoting the corresponding competition rates as , respectively, introduces the competition term as . The dynamics of tumor cells can then be described by the following ordinary differential equation:
2.2. Modeling Natural Killer Cells
To model natural killer (NK) cells, we will assume that these cells have a constant source 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 denotes the maximum NK cell recruitment rate by tumor cells and 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 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
Note that we have also included a natural death of NK cells through .
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 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 assume to be a constant source of dendritic cells, to be the rate at which NK cells kill dendritic cells, to be a proliferation rate of dendritic cells due to tumor cells, to be the rate corresponding to the interaction of dendritic cells with CD8+ T cells, and to be the natural death rate of dendritic cells. We then have
2.4. Modeling Cytotoxic CD8+ T Cells
Among many factors that impact the growth of tumor cells, it is well known that T cells are an important component of the immune system that kills tumor cells. It has been seen that active tumor-specific 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 T cells can remove dendritic cells [35, 36].
In our model, to describe the dynamics of the T-cells, we will consider to be the rate of interaction between dendritic cells and tumor cells to activate T cells; denotes the form of competition between T cells and tumor cells, and denotes the natural death rate of T cells.
It has also been seen that 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 . We will also need an additional term that helps describe the regulation and suppression of T-cell activity when there are high levels of activated T cells without responsiveness to cytokines in the system [38, 39]. This term is denoted by . We also include 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. The model for the dynamics for the T cell growth population becomes
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 T-cells are promoted through antigen-specific cytolytic immune cells. We do this by adding the term in equation (8) and we have
To include the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream as follows:
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 and , 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. 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 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 and , both not equal to zero, there is no trivial equilibrium which implies
For tumor-free equilibrium, at an equilibrium point, we have which yieldswhich yields
Similarly, setting 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 and .
For these to have biological meaning, we need
These conditions suggest critical values for the death rate and the source term for NK cells to bein 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 bywhere
Evaluating these terms at the general tumor-free equilibrium point gives
To solve for the eigenvalues λ, we solve the equation orwhich yieldswhere 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 beand and 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
This implies that the eigenvalues and that are the roots of equation (32) have negative real parts.
In order for the tumor-free equilibrium to be stable, we require , which implies that if the tumor growth rate a is lesser than the critical value given by , 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 [9](ii)b: tumor-carrying capacity estimated as [9](iii): NK cell tumor cell kill rate estimated as [12](iv): NK cell inactivation rate by tumor cells estimated as [9](v): rate of dendritic cell priming NK cells estimated as [12](vi): NK cell dendritic cell kill rate estimated as [12](vii): rate of tumor cells priming dendritic cells estimated as Estimate(viii)e: death rate of NK cell estimated as [12](ix): CD8+ T cell dendritic cells kill rate estimated as [12](x): rate of dendritic cells priming CD8+ T cell estimated as [12](xi): death rate of dendritic cells estimated as [12](xii)h: T inactivation rate by tumor cells estimated as [12](xiii)i: death rate of CD8+ T cells estimated as [9](xiv)j: dendritic cell tumor cell kill rate estimated as [12](xv)k: NK cell tumor cell kill rate estimated as [12](xvi): source of NK cells estimated as [12](xvii): source of dendritic cell estimated as [12](xviii)u: regulatory function by Nk cells of CD8+ T cells estimated as [9]
For the first computation, we will assume there is no additional recruitment terms for T cells and NK cells (), removing some of the regulation, suppression. and activation of T cells (), not including the influence of drug kill terms (), no influence of drug and vaccine interventions () along with the corresponding death rates (). We will also assume 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 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.
(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 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 T cell dynamics also changes as the proliferation rate is doubled. For the rest of the simulations, we will include the effect of and assume the value to be .
(a)
(b)
(c)
(d)
Along with , we also wanted to study the influence of the source term on the dynamics of tumor, NK, and T cells. Figure 4 illustrates this behavior. When we increase the source term of dendritic cells, it increases NK cells and T cells. We also note that these have an effect on tumor growth as both NK and 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 T cells earlier in the tumor growth phase.
(a)
(b)
(c)
Next, to study the effect of TIL drug intervention term only for the 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 from 1 to . 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 . We set and study the influence of increasing 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 but noticed that even with inclusion of a T activation described via Michaelis–Menten interaction given bythere 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 and .
Next, we turn our attention to the effect of the nonlinear term introduced as an inactivation term, which describes the regulation and suppression of T-cell activity [38, 39]. Specifically, Figure 7 shows the effect of the nonlinear term in system (11) for parameters and , and the dynamics show that there is a drastic drop in the number of cytotoxic T cells while a small increase in tumor cell growth.
(a)
(b)
In summary, our computations seem to suggest that a combination of immunotherapy through TIL drug intervention along with chemotherapy through 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:
Figure 8 clearly shows that combined chemotherapy and immunotherapy drug intervention helps reduce the tumor growth for system (35). We have used with and 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, which is the competition rate that impacts the dynamics of tumor cells due to natural killer cells and which is the parameter-related proliferation of dendritic cells due to tumor cells. Recalling the influence of doubling the latter parameter 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 and 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 (see Figure 10).
Next, we make a guess for values of and 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 .
We then set up an error expression that is the sum of the squared differences between the computed values and the experimental 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 is within a user-prescribed tolerance , we stop and accept the values of and to be the optimal values. If not, we use the updated values of parameters and 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 bywhere , the least squared error, denotes the differences in the amount of computed tumor cells from the simulation from the observed data (Figure 9) over N observations.
Using same initial conditions with poor guesses for and , the optimization algorithm estimates the parameters to be and . 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 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 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.
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.