Antibiotic cycling versus mixing: The difficulty of using mathematical models to definitively quantify their relative merits

1.   Introduction

The world-wide burden of cancer (malignant tumor) is one of the major health problems, with more than 8 million new cases and 5 million cancer and cancer-related deaths per year ^[1]. Tumor cells can lead to the appearance of specific antigens that are not found on normal cells and trigger reactions by both the innate and adaptive immune systems ^[2,3,4,5,6]. Recent studies show that the immune system plays an important role in recognizing and destroying tumor cells ^[7,8,9]. The immune response to tumor cells is usually cell-mediated by the natural killer (NK) cells and CD8+ cytotoxic T lymphocytes (CTLs).

As the first line of defense in host body, NK cells kill tumor cells in a variety of distinct ways ^{[10,11,12,13,14]}. In addition, NK cells do not express antigen receptors encoded by recombination-activating genes, but instead express a range of germ-line encoded receptors, including both inhibitory and activating proteins, to sense target cells. Upon the target cell recognition through the integration of signals provided by a combination of inhibitory and activating receptors, NK cells become activated. The ligands recognized by these distinct receptors include major histocompatibility complex (MHC) class I, stress-induced molecules, adhesion proteins and other molecules that are used by NK cells to identify cells to be killed. Now we know that NK cells have an innate mechanisms for target cells recognition, as well as an antibody-dependent mechanism ^[15,16].

As part of the adaptive immunity, CTLs play a crucial role in recognizing and destroying tumor cells ^{[17,18,19,20,21]}. When CTLs bind to epitopes from endogenous antigens bound to MHC on the surface of the tumor cells and lyse the cells with perforins, the activation of T cells occurs, and then the activated CTLs can directly kill tumor cells. In short, the CTLs induced apoptosis of tumor cells is mediated through two distinct pathways: (a) the release of lytic granules containing perforins and granzymes which enter the target cells; (b) the interaction of Fas ligand (FasL) with Fas on the CTLs and target cells. More and more studies have shown that the CTLs primed to combat the target cancerous cells contributes to the success of immune surveillance in controlling the tumor development and growth ^{[22,23,24,25]}.

Mathematical modeling is a powerful tool that helps researchers understand complex regulatory mechanisms. An early work in the modeling of the immune system came from Sercarz and Coons ^[26]. After that, a deterministic predator-prey model was proposed to show that the survival rate was increased if the immune system was stimulated and an increase of effector cells increased the chance of tumor survival in some cases ^[27]. Recently, the majority of existing mathematical models in tumor immunology are based on the sets of ordinary differential equations (ODEs). In 1994, Kuznetsov et al. defined an ODE model of the cytotoxic T lymphocytes response to the growth of an immunogenic tumor ^[28]. They predicted a threshold above which there was uncontrollable tumor growth, but below which the disease was attenuated with a recurrent profile with a 3-4 month cycle. Moreover, Tessi et al. developed a set of ODEs considering all compartments of the immune system related to the tumor response, which was based on experimental and clinical results ^[22]. The results of this paper suggested that for a given tumor growth rate, there was an optimal antigenicity maximizing the response of the immune system. In the work of Mahasa et al., the model was designed to illustrate how tumor evaded both arms of host immunity, which was the first mathematical study of the tumor cells escape and acquisition of immune resistance to immune system mediated by NK cells and CTLs ^[23].

Recent progress in tumor immunology and advances in immunotherapy suggests that the immune system plays a fundamental role in combating tumor cells ^{[2,3,4,14,16]}. Therefore, it is important to investigate the mechanisms underlying an effective immune response to tumor cells by using mathematical models. For example, Tessi et al. developed a model including the primary cell populations involved in effector T-cell mediated tumor killing: regulatory T cells, helper T cells, and dendritic cells ^[22]. On the contrary, De Pillis and colleagues considered NK cells, cluster of differentiation 8(CD8), circulating lymphocytes and interleukin 2(IL-2) to introduce a popular model ^[24]. Despite the fact that the major advances in the fight against tumor cells, previous studies have considered too many immune cells, so the interactions between tumor cells and the immune system are unclear.

To investigate the interactions between immune cells and tumor cells more clearly, we propose a mathematical model that focus on those elements of the immune system that are known to be significant in controlling tumor growth. So, the model consists three important agents, namely NK cells, CTLs, and tumor cells. We first assume that both tumor cells and NK cells grow logistically. The logistic model gives the description of the saturated growth dynamics. If the growth dynamics is well under the saturation, it is a linear growth dynamics that has been used in the current modelling approaches. Thus, the logistical growth mechanism provides an alternative method to describe the growth dynamics of NK cells. In addition, the proposed model takes into account an important role of the immune response in the control of tumor cells. Our mathematical model of tumor-immune interactions mainly elucidate the different roles of NK cells and CTLs in suppressing tumor cells. This work represents an attempt to use a relatively simple model to capture the desired behaviors of the system dynamics.

2.   Methods

2.1.   Mathematical model

In this paper, we propose a model which consists a set of three ordinary differential equations to describe the immune response to tumor cells. Our proposed model are based on the following assumptions.

1) For the growth law terms of tumor cells, we consider several possible models including exponential growth, Gompertz growth, and logistic growth. It has been found in numerous studies both in vivo and in vitro that the growth of a tumor cell population is exponential for small quantities while growth is slowed at large population sizes. According to the data in ^[29], we assume that a tumor follows logistic growth in the absence of immune system.

2) As part of the innate immunity, NK cells are present and active in the host all the time, even in the absence of tumor cells.

3) As part of the adaptive immunity, CTLs only occur in large numbers when tumor cells are present in the host body.

4) Both NK cells and the CTLs are capable of killing tumor cells. But CTLs play a leading role in killing tumor, as part of the adaptive immunity.

5) When immune system encounter with tumor cells, some number of NK cells and CTLs become inactive, but have no damage on the cells.

Based on the above assumption, we develop a three-population model that tracks the numbers of tumor cells, NK cells and CTLs. Figure 1 demonstrates the detailed interactions between these cells. We use $N(t)$, $L(t)$ and $T(t)$ to represent the number of NK cells, CTLs, and tumor cells at time $t$, respectively. The model of ordinary differential equations describing the growth, death, and interactions of these populations is given by:

with initial condition $N(0) = N_0 \geq 0$, $L(0) = L_0 \geq 0$, and $T(0) = T_0 \geq 0$.

In our model, all the parameters $a$, $b$, $c$, $d$, $r$, $\mu$, $\alpha_1$, $\alpha_2$, $\beta_1$, $\beta_2$ are positive constants. The rate of the NK cell populations growth is denoted by $a$, which includes both cell multiplication and death, and the parameter $\frac{1}{b}$ represents the maximum carrying capacity of NK cells. We assume that in absence of immune systems, the tumor cell populations follow logistic growth $cT(t)(1-dT(t))$, with intrinsic growth rate $c$ and maximum carrying capacity $\frac{1}{d}$. Since we have assumed that CTLs are only present in the host immune system only when the tumor cells are present, the growth term for CTLs consists only of natural cell death $-\mu L(t)$. The recruitment term $rN(t)T(t)$ represents interactions between NK cells and tumor cells, through which we model the fact that the specific immune response of CTLs is activated only after the activation of the earlier response of innate immunity. Additionally, immune cells are inactivated through contact with tumor cells, the parameters $\alpha_1$, $\beta_1$ represent the respective proportions of NK cells and CTLs that are detached, when tumor cells bound with them. The parameters $\alpha_2$, $\beta_2$ represent the respective proportions of the fractional tumor cells which are killed by NK cells and CTLs.

2.2.   The reduced model

In order to quantitatively analyze the model, we attempt to determine the model baseline parameter values or at least a reasonable physiological range for them, and fix the values of state variables $c$, $d$, $\mu$, $\alpha_1$, $\beta_1$. A detailed listing of the parameters in model (2.1), along with their units, descriptions, estimated values and reference sources from which these values were taken, is provided in Table 1.

Then the model (2.1) is as follows:

For the sake of simplicity, we put in dimensionless form the model (2.2), i.e,

By rewriting $\bar{N}, \bar{L}, \bar{T}, \tau$ as $N, L, T$ and $t$, respectively, this lead to the following dimensionless model:

where $p = \frac{a}{\mu}$, $q = \frac{b\mu}{\alpha _2}$, $\delta = \frac{r\beta_2}{\alpha _1\alpha _2}$, $\frac{c}{\mu} = 2.57\times10$, $\frac{d\mu}{\alpha _1} = 2.04\times10^{-4}$, $\frac{\beta_1}{\alpha _1} = 3.42\times10^{-3}$.

3.   Results

3.1.   Equilibria of the reduced model

In this section, we first identify the positively invariant set of the model (2.3), and then study the existence and stability of equilibrium points. Regarding the positively invariant set, we have the following results.

Lemma 3.1. The following set is a positively invariant set of model (2.3) $R_{+}^{3} = \{(N(t), L(t), T(t))|N(t)\geq 0, L(t)\geq 0, T(t)\geq 0\}$.

Proof. Let $X(t) = (N(t), L(t), T(t))$ be any positive solution of model (2.3), with initial condition $X(0)\in R_{+}^{3}$. By the equations of the model (2.3), we have that

That is, if $N_0\geq 0, L_0\geq 0, T_0\geq 0$, then $N(t)\geq 0, L(t)\geq 0, T(t)\geq 0$. Thus, $R_{+}^{3}$ is a positively invariant set of model (2.3).

The equilibria of the model are determined by setting $N^{'}(t) = L^{'}(t) = T^{'}(t) = 0$, and solving the resulting algebraic equations:

Then the model (2.3) has two non-negative boundary equilibria: $O(0, 0, 0)$, $E_0(\frac{1}{q}, 0, 0)$ and endemic equilibria $E_1(N_{1}^{*}, L_{1}^{*}, T_{1}^{*})$, $E_2(N_{2}^{*}, L_{2}^{*}, T_{2}^{*})$. Here

and $T_{i}^{*}$ $(i = 1, 2)$ are solutions of the following equation, given by

It is clear that $N_{i}^{*} > 0$, $L_{i}^{*} > 0$, when $0 < T_{i}^{*} < p$. Thus, endemic equilibria exist if and only if $0 < T_{i}^{*} < p$. By using the Vieta's theorem, we have results about the existence of endemic equilibria of model (2.3) that are presented in Supplementary information Table S1.

From the above discussion, model (2.3) has three kinds of equilibrium points: $O(0, 0, 0)$, $E_0(\frac{1}{q}, 0, 0)$, $E_i(N_{i}^{*}, L_{i}^{*}, T_{i}^{*})$. $O(0, 0, 0)$ means that there is no tumor and no immune cells in the host, which can be ignored. The equilibrium point $E_0(\frac{1}{q}, 0, 0)$ shows that there are only NK cells, which means that the tumor cells should be completely eliminated by treatment and other means. However, positive equilibrium points $E_i(N_{i}^{*}, L_{i}^{*}, T_{i}^{*})$ occur where a small tumor mass might coexist with a large number of immune cells.

3.2.   Stability of the reduced model

We now analyze the local geometric properties of the non-negative equilibria of model (2.3). The first equilibrium $O(0, 0, 0)$ is the trivial state, where all the populations are zero. The eigenvalues of Jacobian matrix for $O(0, 0, 0)$ are $p$, $-1$, $25.7$. Therefore, $O(0, 0, 0)$ is always a locally unstable saddle. For practical and biological reasons, we do not take into account zero equilibrium $O(0, 0, 0)$.

The Jacobian matrix of the equilibrium at $E_0(\frac{1}{q}, 0, 0)$ is given by

The characteristic equation of $J(E_0)$ is

So the eigenvalues of characteristic equation are $-p$, $-1$, $25.7-\frac{1}{q}$. Hence, $E_0(\frac{1}{q}, 0, 0)$ is a saddle for $q > q_{3}^{*} = \frac{1}{25.7}$, or $E_0(\frac{1}{q}, 0, 0)$ is locally asymptotically stable for $0 < q < q_{3}^{*}$.

The Jacobian matrix of the equilibrium at $E_i(N_{i}^{*}, L_{i}^{*}, T_{i}^{*})$, $(i = 1, 2)$ is given by

The characteristic equation of $J(E_i)$ is

where

According to Routh-Hurwitz theroem ^[34], the local geometric properties of $E_i(N_{i}^{*}, L_{i}^{*}, T_{i}^{*})$, $(i = 1, 2)$ are provided in Supplementary information Table S2.

Therefore, we can obtain the stability and type of all equilibrium points of the model (2.3), and the meaningful results are listed in Table 2, detailed in Supplementary information Table S2.

The equilibrium $E_0$, $E_i$ mean no tumor burden and large tumor burden, respectively. In cases 1–3, $E_0$ is unstable, and $E_1$ is locally asymptotically stable. These show that the model allows for the possibility of a large tumor mass, and tumor cells coexist with immune cells. However, cases 4 and 5 present different outcomes, namely $E_0$ is locally asymptotically stable, and $E_2$ is unstable. The most ideal case is that $E_0$ is stable, which we hope to achieve through treatment. The stability of $E_i$ indicates that tumor and immune cells in host body are in balanced. Tables 2 and 4 indicate that the stability of the equilibrium points of the model (2.3) is mainly related to the parameters $p$, $q$ and $\delta$. Therefore, we can change the stability of the equilibrium points of the model by changing the values of the parameters $p$, $q$ and $\delta$.

3.3.   Numerical simulations of the reduced model

In this subsection, we show that how the number and stability properties of the equilibrium points are affected by the changes of parameters $p$, $q$, and $\delta$. In all the simulations in Figure 2, we use the parameter values in Table 1 but specify the different values of $p$, $q$, and $\delta$. For the selection of the parameters $p$, $q$, and $\delta$, the detailed algorithm are shown in section 2 of the Supplementary information. The numerical solutions of our model (2.3) along with the initial condition, are carried out using MATLAB ode45.

Firstly, for Case 1 in Table 2, we choose $p = 7.3\times10^{3}$, $q = 5.1\times10^{-1}$, $\delta = 7\times10^{-5}$. Simulated results for this case, with initial condition $1\times10^{4}$ NK cells, $1\times10^{2}$ CTLs and $1\times10^{6}$ tumor cells, are pictured in Figure 2(a). As shown in Figure 2(a), the tumor cell populations and immune cell populations are non-zero for long days, and they tend to a constant value, that is, positive equilibrium $E_1(N_{1}^{*}, L_{1}^{*}, T_{1}^{*})$ is stable in this case.

Then, we simulate Case 3 in Table 2. In this case, $p = 5\times10^{3}$, $q = 6.2$, $\delta = 1.08$ are used in simulation. The initial condition is $N(0) = 1\times10^{4}$, $L(0) = 1\times10^{2}$, $T(0) = 1\times10^{6}$. Figure 2(b) shows the number of immune and tumor cells. Compared with the simulation of Case 1, the number of tumor cells is much smaller than that in Case 1.

Finally, we simulate Cases 4 and 5 in Table 2. In both cases we use $q = 3\times10^{-2}$, but use $p = 4.9\times10^{3}$ and $\delta = 3\times10^{-2}$ in Case 4, while $p = 5\times10^{3}$, $\delta = 4\times10^{-3}$ in Case 5. The initial values for these two cases are the same as that in Case 1. Results in Figure 2(c), (d) show that, after a period of time, the number of NK cells approaches a constant value, while the numbers of tumor cells and CTLs tend to 0. The number of CTLs tends to 0 because we assume that CTLs only occur when tumor cells are present in the system. These results indicate that the equilibrium $E_0(\frac{1}{q}, 0, 0)$ is stable. In addition. although the tumor cells are eliminated and the number of NK cells tend to be a constant in Cases 4 and 5, the time of tumor elimination is different in Figure 2(c), (d).

The results in Figure 2 suggest that the parameters $p$, $q$, $\delta$ are important parameters to determine the stability of the equilibrium points. The parameter space are within a region where we have two equilibria: one with zero tumor burden ($E_0$), and one with a large tumor burden ($E_i$). In medicine, the goal of treatment is to eliminate tumor cells, which is reflected in mathematical model to drive the system toward the zero tumor burden equilibrium $E_0$. It is known from the above numerical simulations that the NK cells environmental maximum capacity $\frac{1}{q}$ and the CTLs-induced tumor death rate $\delta$ should be enhanced through treatment and other methods to achieve $E_0$ stability. However, based on the practical significance, NK cells environmental maximum capacity cannot be changed. Therefore, we can increase the source of CTLs through treatment to improve the lethality of CTLs.

3.4.   Tumor sensitivity of the parameter q

Despite some parameters the model (2.3), we note that it is significantly more sensitive to variations to a few parameters. A sensitivity analysis of the parameters can highlight those model parameters that have the greatest effect on model outcome. A standard approach to performing this analysis is to fix all the parameters of model but one value, and then to increase or decrease one parameter by a certain percentage, and examine the effect on the model outcome. After parameter sensitivity analysis, we find that the parameter with the highest influence on tumor cells is the maximum environmental capacity of NK cells $\frac{1}{q}$. In model (2.3), we set $p = 6\times10^{3}$, $\delta = 6.8\times10^{-5}$, and $q$ takes a value of 0.01 to 1, with intervals of 0.01. Then we plot the number of tumor cells at $t = 1000$, with initial condition $1\times10^{4}$ NK cells, $1\times10^{2}$ CTLs and $1\times10^{6}$ tumor cells. As can be seen from the Figure 3, a small change in parameter $q$ determines whether the number of tumor cells tends to be large or small or zero. This finding hints at the direction we should follow to treat tumor cells. Since $q = \frac{b\mu}{\alpha _2}$, we can search for effective treatments to eliminate tumor cells by changing the parameters $b$, $\mu$, $\alpha _2$. Given the human body's ability to tolerate, it is impossible to increase the maximum environmental capacity of NK cells $\frac{1}{b}$. Therefore, we can search for effective treatments to eliminate tumor cells by decreasing death rate of CTLs and increasing rate of NK-induced tumor death.

4.   Discussion

In this work, we develop a novel mathematical model to describe how tumor cells evolve and survive the brief encounter with the immune system, which is in conformity with experimental data ^[29]. Compared with other models developed previously, main novelty of this model is to integrate the dynamics of important parts of immune system. The main contributions of this model are the use of two immune cell populations to reduce the tumor burden, and understand the interactions between tumor cells and immune system clearly. Through an analysis of model (2.3), we get the equilibrium points of the model, as well as the conditions for stability or instability. The above results demonstrate the fact that an effective immune surveillance, mediated by NK cells and CTLs, is capable of controlling small tumor cells growth. More importantly, the model simulations (see Figure 2) indicate that the very small changes in parameters $p$, $q$, $\delta$ can have drastic consequences on the outcome of the tumor cells. This means that NK cells and CTLs play a crucial role in tumor immune surveillance. What's more, the immune system alone is not fully effective against tumor progression.

The simulations of model, depicted in Figure 2(c), (d), support the importance of CTLs in immune surveillance. What's more, an increment in the number of CTLs or rate of CTLs-induced tumor death results in depleted number of tumor cells. Thus, the CTLs based therapy can be used to enhance immune surveillance against developing tumors. More importantly, the parameter $q$ plays an important role in the state of tumor cells. As shown in the Figure 3, we can see that small change in parameter $q$ can change the state of tumor cells greatly. In fact, weak immunogenicity is suggested to be a major route of immune response failure under which adaptive immunity ultimately disappears without sufficient antigenic stimulations.

Through the discussion in this work, we know that the immune system is not able to completely control the development of tumors or destroy tumors. So, it is vital to consider treatment in the model. Past studies and clinical results have shown that surgery, chemotherapy, radiation therapy, and combination therapy are among the effective treatments for tumor patients ^{[24,35,36,30]}. Although these treatments have played key roles in tumor treatment, there are some side effects in many cases. Immunotherapy is being developed and tested as a promising new approach to treatment for tumor patients ^[2,3,4,9]. Immunotherapy for tumor operate via several different ways. The goal of immunotherapy is to boost the system's own immune defense against the tumor. Some therapies aim to promote the response from cytotoxic T-lymphocytes, such as CTLs, by attaching immune-stimulating adjuvants such as a virus or bacteria to a patient's own irradiated tumor cells^[31]. Based on the results of this work, we can consider adding a drug to the model that directly stimulates CTLs activity without side effects.

In collaboration with mathematicians, biologists and physicians, we hope to use mathematical models to better and more scientifically plan the schedules of the therapies. In the future work, we will propose a mathematical model to investigate the interactions between immunotherapy and tumor cells in more depth, understand dose-response profiles and then optimize therapeutic strategies to reduce tumor cells. This could be a nontrivial task for giving optimal therapy to take into account such important physiological processes as drug diffusion, nutrient consumption and the development in tumor cells. Despite these complexities and limitations, we still believe that the mathematical models are useful in aiding our understanding of tumor cells biology and treatment. And the models will be invaluable for testing hypotheses, as well as identifying biological mechanisms.