Mathematical modeling of tumor growth: the MCF-7 breast cancer cell line

1.   Introduction

Breast cancer is the second leading cause of cancer death for women worldwide ^[1,2]. Approximately 75% of breast cancer cases are estrogen receptor positive ^[3]. The MCF-7 cell line is a common in vitro model used for studies on human breast cancer. The popularity of MCF-7 can be attributed to its estrogen responsiveness through expression of estrogen receptor, making it an idea model to study the effect of hormones on cell proliferation and protein synthesis ^[4,5,6].

White blood cells are the cells of the immune system that helps the body fight infectious disease and foreign invaders. All WBCs are produced in the bone marrow and found in the blood and lymphatic system ^[7]. The number of WBCs in the blood is often a marker of disease. A low WBC count may occur as a result of viral infections, cancer treatment, bone marrow condition, or certain diseases ^[8]. An individual with a low WBC count may increase frequency of infections and fevers and cause low immunity ^[9]. Shankar et al. ^[10] have found an association between WBC count and cancer mortality. Other clinical studies have reported that there is no significant association between WBC count and cancer risk or cancer mortality ^{[11,12,13,14]}.

Natural killer cells are specialized WBCs and form a very important component of the innate system. They are able to detect early-stage developing tumors and secrete immuno-regulatory cytokines and numerous chemokines that directly and effectively eliminate tumor cells. This process is sometimes called immune surveillance, and it is an important process to prevent neoplastic cells from developing into tumors. Cytotoxic T lymphocytes, which are produced in the bone marrow and migrate to the thymus for maturation, are key players of adaptive immune system. Upon activation of CTLs, the T-cell antigen receptors expressed on CTLs allow them to monitor all cells of the body, ready to destroy infected cells. Cytotoxicity is mediated through Fas or perforin pathways inducing apoptosis in target cells ^[15].

Despite the antitumor effect of NK cells and CTLs, there are still neoplastic cells that escape immune surveillance. However, even a fast growing tumor cannot grow beyond 1-2 mm in diameter without angiogenesis. At this stage, a balance exists between cell proliferation and apoptosis, and tumors may maintain dormant for years. Furthermore, clinical evidence has suggested that many of us have microscopic tumors that remain asymptomatic and never develop and progress to become large and lethal ^[16,17,18]. In fact, it has been estimated that more than one third of middle aged women who have never been diagnosed with breast cancer throughout their life were found at autopsy with in situ tumors in their breast ^[16]. When tumors develop their own vascularization, they can grow much larger and become lethal.

Mathematical modeling of tumor growth is a powerful tool to understand, predict, and improve the outcome of a treatment regimen. A mathematical model of tumor growth often includes tumor-immune interaction which is an important component for the model to produce clinical phenomena. Most mathematical models for cancer are based on differential equations ^[19]. While partial differential equations have been employed to model the spatio-temporal growth and spread of cancer cells, ordinary differential equations (ODEs) have been used to investigate the interactions among cancer cells, immune cells, and cytokines ^[20]. ODE systems with two or three cell populations, which usually contain tumor cells and effector cells, have been constructed to study the population dynamics ^{[21,22,23,24,25,26,27,28,29,30]}. These models are able to capture the essential patterns of tumor growth and decay: exponential/oscillatory growth or decay of tumor cells. Forys et al. ^[30] and de Vladar and Gonz$\acute{\mbox{a}}$lez ^[29] have reported the existence of stable small tumor equilibria in their models. These dynamics can describe dormant or persistent microscopic tumors. However these results are merely obtained by mathematical models ^[29,30] while the interaction terms and parameter values in the mathematical models have not been verified by experimental or clinical data.

More complex models involving more equations are devised to study particular aspects of tumor biology and have been able to provide significant insights into the mechanisms of tumor-immune interactions. These models inevitably involve more parameters and thus less well parameterized ^[20]. de Boer et al. ^[31] have considered an ODE system of eleven equations and reported that the cytotoxic T lymphocyte (CTL) count depends on the time of T helper cell activation. They have also pointed out that a tumor with strong antigenicity is destroyed by CD8+ T cells while a tumor with weak antigenicity is destroyed by macrophages. Numerical simulation using an ODE system of six equations by Cappuccio et al. ^[32] has suggested that IL-21 can be used to treat non-immunogenic tumors but not highly immunogenic tumors. Kronik st al. ^[33] have proposed a model, which includes gliomas cancer cells, CTLs, MHC class Ⅰ and Ⅱ molecules, and two cytokines, to study the effects of immunotherapy.

Jarrett et al. ^[34] have proposed a mathematical model of HER2+ breast cancer with immune response and studied the outcomes of trastuzumab treatment. Annan et al. ^[35] have constructed an ODE system describing the interactions between breast cancer cells and immune cells including CTLs, macrophages, NK cells, and helper T cells. Roe et al. ^[36] have presented a mathematical model to study the effect of chemotherapy on breast cancer. Most of the parameter values and the functional forms in these studies have not been justified by experimental and clinical data. The purpose of this paper is to build a mathematical model using ODEs which describe MCF-7 breast tumor growth, tumor-immune interaction, and tumor-estrogen interaction. Experimental data from several research studies are used to determine the functional responses and parameter values in the mathematical model. The dynamics produced by the mathematical model are to be compared with experimental and clinical phenomena. To provide details, the mathematical model is presented in Section 2. Numerical examples and discussion are given in Section 3. Finally, a brief conclusion is made in Section 4.

2.   Mathematical model

The model considered in this paper consists of five state variables: the MCF-7 tumor population $T$ (cells), the circulating level of estradiol $E$ (pmol/L), NK cell population $N$ (cells/L), CTL (or CD8+ T) cell population $L$ (cells/L), and WBC population $C$ (cells/L). The mathematical model is described as the following ODEs:

where $E(t)$ is a periodic function and $t$ is in days.

Tumor growth includes an initial exponential growth phase followed by a linear growth phase. As tumors grow, the availability of resources (nutrients, oxygen, and space) decreases resulting in deceleration of tumor growth, and the tumor size eventually reaches its maximum value ^[5,37]. The tumor population $T$ is assumed to follow a logistic or sigmoidal curve including three phases: an exponential phase, a linear phase, and a plateau phase. The parameter $a$ is the intrinsic growth rate, and $K$ is the carrying capacity of the tumor.

The proliferation of MCF-7 cells is known to be stimulated by E2. The term $I_{TE}(T, E)$ represents the interaction between MCF-7 cells and E2. The terms $I_{TN}(T, N)$ and $I_{TL}(T, L)$ represents tumor lysis by NK cells and CTLs, respectively. The term $I_{NT}(T, N)$ represents recruitment of NK cells in the presence of tumor cells. Because NK cells are a type of WBCs, the growth of NK cells, $eC$, depends on the concentration of WBCs. The parameter $f$ represents the death rate of NK cells, and the parameters $\alpha$ and $\beta$ represent the birth and death rate of WBCs, respectively. NK cells release perforin, which forms pores in the cell membrane of the target cell. They also release granzymes, which enter through these pores, inducing apoptosis. An activated NK cell can kill 4 target cells in 16 hours, and then it appeared to be exhausted due to significant reduction in perforin and granzyme B content ^[38,39]. The term $-p_2NT$, which is adopted from a model of tumor-immune interactions by de Pillis et al. ^[28], describes that NK cells become inactivated after some encounters with tumor cells, and the parameter $p_2 = 3.42\times 10^{-6}$ represents the mean inactivation rate.

A normal range for the WBC count is $4-12\times 10^9$ cells per liter ^[40]. Experimental data by ^[41] have shown that healthy young individuals have NK cell counts of $(3.87\pm 1.64)\times 10^8$ cells/L (Mean$\pm$SD) and healthy elderly individuals have NK cell counts of $(3.77\pm 2.21)\times 10^8$ cells/L. It is known that approximately 34% of WBCs are lymphocytes in human blood ^[42]. Pascal et al. have reported 2 to 31% of NK cells within blood lymphocytes in a group of healthy adults ^[43]. Thus an immune system with an equilibrium NK cell count of $4\times 10^8/L$ and WBC count of $5.7\times 10^9/L$ is considered to be a normal healthy immune system. The equilibrium concentration of WBC is $\alpha / \beta$ which is determined by Eq (2.4). The parameter value $\beta = 6.3\times 10^{-3}$ in a research study by de Pillis ^[44] is used in this paper. An equilibrium WBC concentration of $5.7\times 10^9$ cells per liter implies $\alpha = 5.7\times 10^9\times 6.3\times 10^{-3} = 3.6 \times 10^7$. Experimental data from a research study by Zhang et al. ^[41] have shown that the half-life of NK cells is about 10 days, and this implies $f = 0.0693$. An equilibrium NK cell count of $4\times 10^8/L$ in the absence of tumor cells implies $e = fN/C = 0.0693\times 4\times 10^8/5.7\times 10^9 = 0.00486$. Consider that the size of an MCF-7 cell is 15-25 $\mu$m in diameter ^[45,46] and a tumor of 2 cm in diameter is within the diagnostic capabilities ^[47]. Assuming that an MFC-7 cell is 20 $\mu$m in diameter, an MCF-7 tumor of 2 cm in diameter contains $10^9$ tumor cells. Several research studies ^{[27,28,48,49,50,51]} have also used $10^9$ cells for the carrying capacity $K$ in their mathematical models. In this model, the carrying capacity $K = 10^9$ is used.

The matured T cells are considered as naive T cells as they have not encounter appropriate antigen. The number of mature T cells increases progressively during young life and then remains relatively constant during adulthood ^[52]. An activation of naive CD8+ T cells requires two signals. The first signal occurs when a naive CD8+ T cell encounters and interacts with an antigen presenting cell, such as a dendritic cell and a macrophages, through the T cell receptor (TCR) binding to MHC molecules. The second signal, the co-stimulatory signal, is provided by interaction between CD28 and B7 on the membrane of the T cell and APC, respectively ^[15]. Once activated, CD8+ T cells undergoes clonal expansion with interleukin-2 (IL-2) stimulation ^[53]. At the end of immune response, CD8+ T cells undergo a contraction phase during which about 90% of CD8+ T cells die. The cells that survive apoptosis become memory cells for a rapid secondary response to antigen ^[54].

The term $T/(\alpha_5 +T)$ is used to model the phenomenon that CTLs are activated in the presence of tumor cells and undergo apoptosis at the end of immune response. The constant $\alpha_5$ is assumed to be a relatively small number compared with the carrying capacity, $K$, of tumor cells. In this paper, $\alpha_5 = 1000$ is assumed. From the clinical data by Gruber et al. ^[55] in a study of relationship between tumor cells and T lymphocytes in breast cancer patients, the patients had CTL counts of $365\pm 194$ cells/$\mu L$ (Mean$\pm$SD) and naive CTL counts of $237\pm 112$ cells/$\mu L$. The constant $L_N$ in Eq (2.3) represents the number of naive CTLs, and $L_N = 2.3\times 10^8$ cells/L is assumed in this paper. The constant $K_L$ represents the carrying capacity of CTLs and $K_L = 8\times 10^8$ cells/L, which is about mean+2SD, is assumed. Homann et al. ^[56] have reported that during primary immune response CTLs have a doubling time of 6-10 hours between days 5 and 7 and 24–36 hours between days 7 and 9. The number of CTL population reaches its peak on day 8. De Boer et al. ^[57] have reported that the half life of CTLs is approximately 41 hours (about 1.7 days) during the contraction phase. This gives a death rate of $d = ln2/1.7 = 0.41$ in Eq (2.3). The term $\frac{p_5I}{\alpha_4+I}L$ represents proliferation of CTLs, where $I$ represents the concentration of IL-2 (g/L). Assume that the average of doubling time during the expansion phase is 8 hours. This gives a growth rate of $\frac{p_5I}{\alpha_4+I} = 3ln2 = 2.07$. To estimate the parameter value for $p_4$, assume that $T/(\alpha_5 +T)\approx 1$ in the presence of tumor during the expansion phase, and the number of CTL population reaches a level at 90% of its carrying capacity $K_L$, which is $7.2\times 10^8$ on day 8. Equation (2.3) is reduced to

with conditions $L(0) = 0$ and $L(8) = 7.2\times 10^8$. Solving for $p_4$ in the above equation gives $p_4 = 9\times 10^{-5}$.

The interactions between immune cells and tumor cells are usually described by functional responses of the predator-prey type ^[20]. The typical forms of the predator-prey type are based on Lotka Volterra, Holling type Ⅱ, or Holling type Ⅲ ^[58]. These forms with predator interference can be sumarised as $I_{xy}(x, y) = \bar{p}x^{\lambda_1}y^{\lambda_2}/(1+\bar{\alpha} x^{\lambda_1}+\bar{\beta} y^{\lambda_2})$. For example, the interaction term $I_{xy}$ is Lotka Volterra functional response if $\lambda_1 = \lambda_2 = 1$ and $\bar{\alpha} = \bar{\beta} = 0$, and $I_{xy}$ is of Holling type Ⅱ with predator interference if $\lambda_1 = \lambda_2 = 1$. The functional responses with $\lambda_1 = 1$ or 2 and $\lambda_2 = 1$ or 2 will be considered in this paper.

2.1.   MCF-7 cells growth and the effect of $E2$

Let $T_d$ be the doubling time of the MCF-7 cell population. Then $a = ln2/T_d$. It has been found that the growth of MCF-7 cells in vitro varies over a large range in the absence of estrogens. Experimental data have shown that the population doubling time ranges between 30 hours and several days ^[59,60,61]. A doubling time of 30 hours is corresponding to $a = ln2/(30/24) = 0.55$ whereas a doubling time of 10 days is corresponding to $a = ln2/10 = 0.069$. A value falls in the range 0.069-0.55 is considered a reasonable parameter value for $a$. MCF-7 cells contain estrogen receptors and are estrogen responsive. Consider the following functional response for $I_{TE}(T, E)$:

The parameters $\lambda_1 = 1$ or 2, and $\lambda_2 = 1$ or 2. Parameters $\alpha_1$ and $\beta_1$ are positive constants. The MCF-7 growth in the presence of $E2$ becomes

Equation (2.7) is used to fit the experimental data by Nawata et al. ^[59], where the $E_2$ concentration is $10^{-8} M$. Figure 1 shows the experimental data and the fitting curve. The experimental data show that a population level of $4\times 10^6$ cells is at the log growth phase. Therefore, $K = 10^7$ and $E_2 = 10000$ (pmol/L) are used in this fitting. Let $x_i$, $i = 0, 1, \cdots, 4$ be the number of tumor cells from experimental data shown in Figure 1. The parameters $a$, $c$, $\alpha_1$, and $\beta_1$ are to be found for minimizing $\sqrt{ \sum_{i = 0}^4 ((T_i-x_i)/x_i)^2}$, where $T_i = T(3i)$ is the solution of Eq (2.7) with initial condition $T(0) = 2\times 10^4$. Different sets of $(\lambda_1, \lambda_2)$ values are fitted separatedly. The Nelder Mead simplex method is employed to solve the least square error of the fitting. The ODE solver is the Runge Kuta method. The Matlab build-in functions fminsearch and ode45 can also fulfill the tasks. Figure 1 shows that $\lambda_1 = 1$ and $\lambda_2 = 2$ can best fit Eq (2.7) to the data with the parameter values $a = 0.3$, $c = 1.93\times 10^{-6}$, $\alpha_1 = 0.507$, and $\beta_1 = 7.08\times 10^{-8}$.

Figure 1.  Curve fitting for $I_{TE}$ and parameter values $a$, $c$, $\alpha_1$, and $\beta_1$.

DownLoad: Full-Size Img PowerPoint

2.2.   Tumor lysis by NK cells

The term $I_{TN}(T, N)$ representing cytolysis of NK cells is modeled as follows:

Parameter $p_1$, $\alpha_1$, and $\beta_1$ are positive constants to be determined by experimental data. The interaction between NK cells and tumor cells includes $I_{TN}(T, N)$ and $-p_2NT$, where $p_2 = 3.42\times 10^{-6}$. Consider the following equations:

The experimental data used in this fitting were provided by Caragine et al. ^[62], and 4-hour chromium-51 release assays were used in the cytolysis assays. Note that the experiment has been conducted to study MFC-7 tumor cell lysis by NK cells. Figure 2 shows the experimental data and the fitting curve. Let $x_i$, $i = 0, 1, 2, 3$ be the experimental data points in Figure 2 and $y_i = (T(0)-T(4/24))/T(0)$, where $T(4/24)$ is the solution to Eqs (2.9) and (2.10) with $T(0) = 10^6$ given in ^[62] and $N(0) = N_i$, where $N_0 = 0$, $N_1 = 25T(0)$, $N_2 = 50T(0)$, and $N_3 = 100T(0)$. The fitting process is formed using different sets of $(\lambda_1, \lambda_2)$. The parameter values for $p_1$, $\alpha_2$, and $\beta_2$ are to be determined for minimizing $\sqrt{ \sum_{i = 0}^3 (y_i-x_i)^2}$. It has been found that $\lambda_1 = 1$ and $\lambda_2 = 2$ can best fit Eqs. (2.9) and (2.10) to the data, with the parameter values $p_1 = 8.7\times 10^{-4}$, $\alpha_2 = 7\times 10^6$, $\beta_2 = 5.4\times 10^{-5}$.

Figure 2.  Curve fitting for $I_{TN}$ and parameter values $p_1$, $\alpha_2$, and $\beta_2$.

DownLoad: Full-Size Img PowerPoint

2.3.   Tumor lysis by CLTs

The term $I_{TL}(T, L)$ representing cytolysis by CTLs is modeled as follows:

The fitting process is similar to that in Section 2.2. Using the experimental data in ^[63], $(\lambda_1, \lambda_2) = (1, 1)$ or (2, 1) can better fit Eq (2.11) to the data than $(\lambda_1, \lambda_2) = (1, 2)$. Note that the experiment has been conducted to study MFC-7 tumor cell lysis by CLTs. The parameter values $p_6 = 2.36\times 10^{-4}$, $\alpha_6 = 0$, $\beta_6 = 6.04\times 10^{-5}$ are determined for $(\lambda_1, \lambda_2) = (1, 1)$ and $p_6 = 2.38\times 10^{-8}$, $\alpha_6 = 0$, $\beta_6 = 4.3\times 10^{-5}$ are determined for $(\lambda_1, \lambda_2) = (2, 1)$. Figure 3 shows the experimental data and the fitting curve.

Figure 3.  Curve fitting for $I_{TL}$ and parameter values $p_6$, $\alpha_6$, and $\beta_6$.

DownLoad: Full-Size Img PowerPoint

2.4.   Recruitment of NK cells

The term $I_{NT}(T, N)$ represents the recruitment of NK cells to the tumor site. To model the recruitment of NK cells, consider the following equations:

The form and parameter values of the functional response for recruitment of NK cells are to be determined using the experimental data from a research study by Chen et al. ^[64]. The experiment has been conducted to study and compare different treatment starting after 7 weeks of the presence of MCF-7 breast cancer in ovariectomized nude mice. Therefore, the term $I_{TE}$ is dropped out of Eq (2.12). The experimental data set used in this fitting is the one without treatment. It is assumed that the CTL and NK counts have reached high levels after 7 weeks of the presence of tumor cells, and thus the initial conditions of Eqs (2.12)–(2.14) are given by $T(0) = 5.7\times10^6$, $N(0) = 2\times10^8$, and $L(0) = 3\times10^8$. The proliferation of CTLs slows down to a doubling time of 24–36 hours after reaching the peak level ^[56]. Using an average doubling time of 30 hours, 1.25 days, gives $p_5I/(\alpha_4+I) = ln2/1.25 = 0.55$ in this fitting. Other experimental data have shown WBC counts of $(4.3\pm 0.8)\times 10^3$ cells/$\mu L$ in mice ^[65]. The concentration of WBC is assumed to remain an equilibrium concentration $C = 4.3\times 10^9$ cells/L in this fitting. The other parameter values have been determined in the previous subsections. The fitting process is similar to that in Section 2.1, where $p_3$, $\alpha_3$, and $\beta_3$ are unknown parameters. Figure 4 shows the experimental data and the fitting curve for all sets of $(\lambda_1, \lambda_2)$ values.

Figure 4.  Parameter values determined for $p_6$, $\alpha_6$, and $\beta_6$ in Section 2.3 poorly fit the data.

DownLoad: Full-Size Img PowerPoint

Figure 4 also shows that the model, Eqs (2.12)–(2.14), with the parameter values determined in the previous subsections gives strong cytotoxic function that cannot well fit the in vivo data. The fitting process is repeated using $(\lambda_1, \lambda_2) = (1, 1)$ and unknown parameters $p_3$, $\alpha_3$, $\beta_3$, $p_6$, $\alpha_6$, and $\beta_6$. It has been found that the result of the fitting is satisfactory, Figure 5 shows the experimental data and the fitting curve, and the parameter values $p_3 = 1.87\times 10^{-8}$, $\alpha_3 = 1.6\times 10^{-5}$, $\beta_3 = 3.27$, $p_6 = 2.04\times 10^{-3}$, $\alpha_6 = 0.268\times 10^{-6}$, and $\beta_6 = 4343$ are determined.

Figure 5.  Curve fitting for $I_{NT}$ and parameter values $p_3$, $\alpha_3$, $\beta_3$, $p_6$, $\alpha_6$, and $\beta_6$.

DownLoad: Full-Size Img PowerPoint

2.5.   E2 concentration

Finally, the blood content of estradiol provided by Wu et al. ^[66] is used for the periodical function $E(t)$. Let $\tau = 29$ be the period of $E(t)$, and $\tilde{E}(t)$, $t\in[0, \tau)$ be the function whose graph is shown in Figure 6. Then $E(t) = \tilde{E}(t-n\tau)$ for $t\in[n\tau, (n+1)\tau)$.

Figure 6.  Estradiol levels across the menstrual cycle.

DownLoad: Full-Size Img PowerPoint

2.6.   ODE model and parameter values

Using the specific forms for Eqs (2.1)–(2.4) leads to the following equations:

A simulation study by Nikolopoulou et al. ^[67] have reported that the IL-2 density is not very dynamic and the equilibrium density of IL-2 is about $2.37\times 10^{-8} g/L$ which is also the value for half saturation constant $\alpha_4$. This assumption gives $p_5 = 4.14$ since $\frac{p_5I}{\alpha_4+I} = 2.07$. The parameter values are summarized and shown in Table 1.

Table 1.  Parameter values in Eqs (2.16)–(2.20).

Parameter	Value	Units	Reference
$a$	0.3	Day$^{-1}$	fitted with data [59]
$c$	$1.93\times 10^{-6}$	L Cell$^{-1}$ Day$^{-1}$ pmol$^{-1}$	fitted with data [59]
$\alpha_1$	0.507	L pmol$^{-1}$	fitted with data [59]
$\beta_1$	$7.08\times 10^{-8}$	Cell$^{-1}$	fitted with data [59]
$K$	$10^9$	Cell	[27,28,48,49,50,51]
$p_1$	$8.7\times 10^{-4}$	L$^2$ Cell$^{-2}$ Day$^{-1}$	fitted with data [62]
$\alpha_2$	$7\times 10^6$	Cell$^{-1}$	fitted with data [62]
$\beta_2$	$5.4\times 10^{-5}$	L$^2$ Cell$^{-2}$	fitted with data [62]
$\beta$	$6.3\times 10^{-3}$	Day$^{-1}$	[44]
$\alpha$	$3.6\times 10^7$	Cell L$^{-1}$ Day$^{-1}$	Estimated within bounds [40]
$e$	0.00486	Day$^{-1}$	Estimated within bounds [41]
$f$	0.0693	Day$^{-1}$	Estimated from data [41]
$p_2$	$3.42\times 10^{-6}$	Cell Day$^{-1}$	[28,100]
$p_3$	$1.87\times 10^{-8}$	Cell $^{-1}$ Day$^{-1}$	fitted with data [62]
$\alpha_3$	$1.6\times 10^{-5}$	Cell$^{-1}$	fitted with data [62]
$\beta_3$	3.27	L Cell$^{-1}$	fitted with data [62]
$p_6$	$2.04\times 10^{-3}$	L Cell$^{-2}$ Day$^{-1}$	fitted with data [62]
$\alpha_6$	0.268	Cell$^{-2}$	fitted with data [62]
$\beta_6$	4343	L Cell$^{-1}$	fitted with data [62]
$L_N$	$2.3\times 10^{8}$	Cell L$^{-1}$	Estimated from data [55]
$K_L$	$8\times 10^{8}$	Cell L$^{-1}$	Estimated from data [55]
$p_4$	$9\times 10^{-5}$	Day$^{-1}$	Estimated from data [56,55]
$I$	$2.3\times 10^{-11}$	g L$^{-1}$	[67,107]
$\alpha_4$	$2.3\times 10^{-11}$	g L$^{-1}$	[67,107]
$p_5$	$4.14\times 10^{-3}$	L Cell$^{-2}$ Day$^{-1}$	Estimated from data [56]
$d$	$0.41$	Day$^{-1}$	Estimated from data [57]
$\alpha_5$	1000	Cell	Assumed

---

 | Show Table

DownLoad: CSV

3.   Numerical simulations and discussion

In this section, the dynamics of the model, Eqs (2.16)–(2.20), with parameter values shown in Table 1 are to be studied. For the tumor free case, Eqs (2.16)–(2.20) with parameter values shown in Table 1 give that an immune system has equilibrium levels of NK cell population and WBC population are $4\times 10^8$ cells/L and $5.7\times 10^9$ cells/L, respectively.

3.1.   Normal immune system

Vacca et al. have reported that the immune system is able to eliminate tumors at initial stages ^[68], and it has been demonstrated that most tumor cells entering the circulation are destroyed during the first 24 hours ^[69,70,71]. In the first numerical experiment, an initial tumor burden of $10^6$ cells, $T(0) = 10^6$, is considered. In this case, it is assumed that the tumor has escaped immune surveillance and grown to reach a population level of $10^6$ cells. Note that a tumor of $10^6$ cells has a size about 2 mm in diameter and is clinically undetectable. Furthermore, assume that this is the first immune response against the cancer so no effector CTL is present initially. That is $L(0) = 0$. While the activation of CTLs requires some time, Figure 7(a) shows that NK cells provide the first line of defense against cancer and are able to quickly eliminate the tumor.

Figure 7.  A healthy normal immune system using parameter values listed in Table 1 with initial conditions (a) $(T(0), N(0), C(0), L(0)) = (10^6, 4\times 10^8, 4\times 10^9, 0)$ and (b) $(T(0), N(0), C(0), L(0)) = (2\times 10^6, 4\times 10^8, 4\times 10^9, 0)$, (c) $(T(0), N(0), C(0), L(0)) = (4\times 10^6, 4\times 10^8, 4\times 10^9, 0)$, and (d) $(T(0), N(0), C(0), L(0)) = (2\times 10^6, 4\times 10^8, 4\times 10^9, 4\times 10^7)$.

DownLoad: Full-Size Img PowerPoint

Suppose that the immune system fails to become activated against a tumor so the tumor grows to a population level of $2\times 10^6$ cells. A tumor of $2\times 10^6$ cells has a size of less than 3 mm in diameter. Figure 7(b) shows that a healthy immune system is able to control but not able to completely eliminate this tumor. The size of the tumor decreases to a population level of about $6 \times 10^5$ cells and has a size of less than 2 mm in diameter. Then the tumor remains small and stable. Evidence suggests that many microscopic tumors never progress to become invasive, a phenomenon termed "cancer without disease" ^{[16,17,18,72]}. At this stage, a balance between immune response and tumor growth could result in long-term persistence of tumor dormancy. Figure 7(c) shows that a healthy immune system is not able to control a tumor of $4\times 10^6$ cells, which has a size of 3.2 mm in diameter. The tumor grows to its carrying capacity eventually. The NK cell population remains at a low population level after the tumor reaches its carrying capacity.

Clinical and experimental observations have demonstrated the existence of immunosurveillance process that can destroy cancer cells or prevent primary tumor growth. The cancer immunoediting hypothesis has been developed based on these findings. Cancer immunoediting, which describes the relation between the tumor cells and the immune system, is composed of elimination, equilibrium, and escape phases ^[73,74]. Elimination is refered to as classical concept of immunosurveillance. Equilibrium occurs when the immune system fails to completely eliminate the tumor but is able to control the tumor leading to long term dormancy. Escape refers to the tumor outgrowth. Figs. 7(a)-(c) show that the system has three locally stable equilibria: the tumor free equilibrium, the equilibrium representing a microscopic tumor, and the equilibrium representing a large tumor.

During the contraction phase, 5-10% of CTLs differentiate into memory T cells ^[56]. Memory T cells can be immediately activated and trigger the second immune response when the individual is exposed to the same pathogen the next time. Figs. 7(b) and (c) shows that the CTL population reaches a high level of $6.4\times 10^8$ cells/L after the expansion phase. Assume $T(0) = 4\times 10^7$ which is about 6% of the CTL population at the peak level. Figure 7(d) shows that the second immune response is able to eliminate a tumor of $2\times 10^6$ cells which is larger than the tumor that an immune system can eliminate in the first response. A tumor of $4\times 10^6$ cells grows to its carrying capacity during the second response (the graph is not shown). Note that although the solutions shown in Figs. 7(a)-(c) appear to approach equilibria, these solutions oscillate with very small amplitudes since $E(t)$ is periodic. These periodic solutions are referred to as equilibria throughout the paper for convenience because their oscillation amplitudes are too small to be seen.

3.2.   The effect of CTL cytotoxicity

Recall that the parameter values fitted to the in vitro data, Figure 3, gives relatively strong cytotoxicity of CTLs against tumor cells compared to the in vivo data. Tumors may secrete molecules that weaken and inhibit CTL functions ^[75] or express molecules that directly inhibit the killing ability of CTLs such as immune checkpoint molecules CD80 or PD-L1 ^[76]. Adoptive cell transfer (ACT) therapy can produce highly active CTLs which give strong immune response ^[77]. Immune checkpoint therapy using checkpoint inhibitor antibodies, such as anti PD-1, anti PD-L1, or anti CTLA-4, can overcome tumor-induced immune dysfunction and enhance immune response against cancer ^[78].

The parameter $p_6$, which is related to the killing activity of CTLs, can be used as an indicator of the strength of cytotoxicity of CTLs. Let $p_6 = 4\times 10^{-3}$ which is twice the cytotoxic strength used in Section 3.1. Figure 8(a) shows that the immune system is able to eliminate a tumor of $4\times 10^6$ cells. When the immune system detects the presence of the tumor, the NK cells fight against the tumor cells during the activation phase, where CTLs are at a low population level. The tumor continues to grow during the activation and expansion phases. After the expansion phase, CTLs are able to control the growth of the tumor and then eliminate the tumor. It is followed by a contraction phase at the end of the immune response. Figure 8(b) shows this immune system is able to eliminate a tumor of $3\times 10^7$ cells, which has a size of 6 mm in diameter, when $p_6 = 2.04\times 10^{-2}$. The population level of NK cells decreases initially when the immune system is fighting against the tumor. After the tumor has been killed, the immune system starts to recover and strengthen.

Figure 8.  A immune system parameter values listed in Table 1 and different $p_6$ values and initial conditions, (a) $p_6 = 4\times 10^{-3}$ and $(T(0), N(0), C(0), L(0)) = (4\times 10^6, 4\times 10^8, 4\times 10^9, 0)$ and (b) $p_6 = 2.04\times 10^{-2}$ and $(T(0), N(0), C(0), L(0)) = (3\times 10^7, 4\times 10^8, 4\times 10^9, 0)$.

DownLoad: Full-Size Img PowerPoint

3.3.   NK cell lymphopenia

The dynamical system, Eqs (2.16)–(2.20), with parameter values listed in Table 1 exhibits multistability. The asymptotic behavior of a solution depends on its initial conditions, which are the initial level of tumor cell population $T(0)$, the initial NK cell count $N(0)$, the initial NK cell count $N(0)$, and the initial WBC count, $C(0)$. Lymphopenia is the condition of having an abnormally low lymphocyte count in the bloodstream. Ebbo et al. ^[79] have studied the impact of low circulating NK cell counts and found a correlation between NK cell lymphopenia and invasive bacterial infections. Patients with severe NK cell lymphopenia have higher incidence of invasive bacterial infections. Although the experimental data in their study have shown negative association between cancer risk and the NK cell count, the correlation is not significant.

An individual with severe or mild NK cell lymphopenia has an NK cell count of $< 5\times 10^7 \mbox{cells}/L$ or $5-9\times 10^7 \mbox{cells}/L$, respectively ^[79]. In a cohort of 457 common variable immunodeficiency patients, 21.7% of the patients have severe NK cell lymphopenia, 25.8% of the patients have mild NK cell lymphopenia, and 52.5% of the patients have normal NK cell counts. In a cohort of 145 healthy adults, 5% of theses healthy controls have NK cell counts below $63\times 10^6$ NK cells/L (Figure S1 of ^[79]) indicating that a healthy individual may have severe NK cell lymphopenia. Suppose that a low NK cell count is temporary. The case of persistent lymphopenia will be discussed in the next subsection. Figure 9(a) shows that an individual with mild NK cell lymphopenia, where $N(0) = 8\times 10^7$, can rapidly eliminate a tumor of $2\times 10^5$ cells, which has a size of 1 mm in diameter. Figure 9(b) shows that the immune system is able to control a tumor of $10^6$ cells but is not able to completely eliminate the tumor which a heathy immune system can effectively eliminate.

Figure 9.  Immune systems with NK lymphopenia using parameter values listed in Table 1 with initial conditions (a) $(T(0), N(0), C(0), L(0)) = (2\times 10^5, 8\times 10^7, 4\times 10^9, 0)$, (b) $(T(0), N(0), C(0), L(0)) = (10^6, 8\times 10^7, 4\times 10^9, 0)$, and (c) $(T(0), N(0), C(0), L(0)) = (2\times 10^5, 2\times 10^7, 4\times 10^9, 0)$.

DownLoad: Full-Size Img PowerPoint

When an individual suffers from severe NK cell lymphopenia, where $N(0) = 2\times 10^7$, Figure 9(c) show that the immune system is able to eliminate a tumor of 1 mm in diameter after the CTL population is expanded. Figure 9(a)–(c) show that an individual with severe or mild NK cell lymphopenia can eliminate a tumor of $2 \times 10^5$ cells but is not able to completely eliminate a tumor of $10^6$ cells which a heathy immune system can effectively eliminate.

3.4.   Weakened immune system

Low WBC counts weaken the immune system and increase the risk of bacteria and virus infection. Some studies have suggested correlation between WBC count and coronary heart disease, stroke, and cancer mortality ^[80,81,82]. Other studies have reported that impact of the WBC count on cardiovascular disease is inconsistent ^[83,84]. Clinical studies have found that a weakened immune system may increase risk of new cases of cancer in immunocompromised populations ^{[11,85,86,87]}. Jee et al. ^[80] have shown that patients with WBC counts of $< 5\times 10^9$ have higher cancer mortality than patients with WBC counts of $5-8\times 10^9$ do. Nevertheless Kim et al. ^[11] have reported that the association between cancer risk and WBC count is inconsistent.

Conditions, such as infectious diseases, chemotherapy, drugs, anemia, splenic dysplasia, birth disorders, or damaging the bone marrow can lead to a low WBC count ^[88]. In this subsection, a system with an eqilibrium WBC count and an eqilibrium NK cell count lower than the normal ranges is considered. The parameter $\alpha$ is related to the production rate of while blood cells. A lower production rate of WBCs corresponds to a smaller $\alpha$ value. An individual with $\alpha = 1.575\times 10^7$ has an equilibrium WBC count of $\alpha/\beta = 2.5\times 10^9$, where clinical studies consider $2.5\times 10^9$ cells/L as a low WBC count ^{[89,90,91,92]}. Let $e = 8.32\times 10^{-4}$, and the equilibrium NK cell count in the absence of a tumor becomes $5\times 10^7$ cells/L which is lower than the normal values. Furthermore, let $p_4 = 3\times 10^{-5}$ indicating a lower activation in CTL immune response. All other parameter values are listed in Table 1. Figure 10 shows the ability of an immune system, with a low WBC count and a low NK cell count, to fight a cancer. The immune system is able to effectively eliminate a tumor of $2\times 10^5$ cells as shown in Figure 10(a). Figure 10(b) shows that the immune system is able to control a tumor with an initial population level of $10^6$ cells, and Figure 10(c) shows that a tumor with an initial population level of $2\times 10^6$ cells will eventually grow to its carrying capacity.

Figure 10.  A weakened immune system using parameter values $\alpha = 1.575 \times 10^7$, $e = 8.32 \times 10^{-4}$, $p_4 = 3\times 10^{-5}$, and all other parameter values listed in Table 1 with initial conditions (a) $(T(0), N(0), C(0), T(0), L(0)) = (2\times 10^5, 5\times 10^7, 4\times 10^9, 0)$, (b) $(T(0), N(0), C(0), L(0)) = (10^6, 5\times 10^7, 2.5\times 10^9, 0)$, and (c) $(T(0), N(0), C(0), L(0)) = (2\times 10^6, 5\times 10^7, 2.5\times 10^9, 0)$.

DownLoad: Full-Size Img PowerPoint

3.5.   Effect of $E2$

Experimental and clinical data support that circulating estradiol is positively associated with breast cancer risk in postmenopausal women ^{[93,94,95,96]}. However, in premenopausal women, most prospective studies found no significant association between breast cancer risk and estrogens ^{[96,97,98,99]}. The numerical simulation using parameter values in Table 1 with a 10-fold increase in the circulating estradiol level shows that the immune system is able to effectively eliminate a tumor of $10^6$ cells and control a tumor of $2 \times 10^6$ cells. A tumor of $4 \times 10^6$ cells grows to its carrying capacity. In the case of a 100-fold increase in the circulating estradiol level, the numerical simulation shows that the immune response against cancer is similar to the case of a 10-fold increase. Figures of numerical simulation in this subsection are not shown because the figures resemble those that appear in the previous subsections. The numerical results suggest that an increase in the circulating estradiol level does not show a significant increase in cancer risk.

4.   Conclusion

A review of previously published mathematical models of tumor-immune interactions ^[20] has pointed out that most of these studies present the results derived from mathematical models and very few experimental or clinical results have been compared with these results obtained with mathematical models. Moreover, many of these models have not considered a specific cancer type or justified the functional forms using clinical or experimental data. In this paper, a mathematical model for the growth of breast tumors of MCF-7 cell line is constructed. All the parameter values (except for $p_2$, $\beta$, and $\alpha_5$) as well as the forms of functional responses for the effect of estradiol, tumor lysis by NK cells, and recruitment for NK cells have been estimated using experimental data ^{[40,41,45,46,55,56,57,59,60,61,62,63,64,66]}. The parameter values for the mean inactivation rate for NK cells $p_2$ and the death rate of WBCs $\beta$ can be found in the study by de Pillis et al. ^[28,44,100]. The results of the fitting procedures are shown in Figures 1–6.

Numerical simulation shows that a healthy normal immune system exhibits multistability where the tumor free, microscopic tumor, and large tumor equilibria are stable. The phenomena resemble the 3 E's of cancer immunoediting ^[73,74]. The immune system is able to rapidly eliminate a tumor with a restricted initial size and the critical size of a tumor that the immune system can eliminate is about 2 mm in diameter (Figure 7). Several experimental studies have demonstrated that NK cells can rapidly eliminate a tumor at initial stages ^{[68,101,102,103]}. A healthy normal immune system may also produce long-term dormancy (Figure 7(b)). The tumor does not grow to become large or lethal and remains at an equilibrium size of less than 2 mm in diameter. This phenomenon has been widely observed in clinical studies ^{[16,17,18,72]}. Clinical and experimental studies have suggested that an immune system may increase its antitumor effect by enhancing the cytotoxicity of CTLs using ACT therapy, tumor vaccines, or checkpoint antibodies ^{[104,105,106]}. Numerical simulation shows that an immune system with a stronger killing activity, a larger $p_6$ value, can eliminate a larger tumor (Figure 8).

An individual suffering from temporary NK cell lymphopenia has a low NK cell count. Numerical simulation demonstrates that there is a positive association between the size of a tumor which an immune system is able to effectively eliminate and the initial NK cell count. Although an immune system with NK cell lymphopenia is not able to effectively eliminate a tumor that a healthy normal immune system is able to (Figure 7(b)), it can rapidly eliminate a smaller tumor (Figure 9(a)). This also suggests that a lower NK cell count may increase cancer risk but the relationship is not significant. Clinical data ^[79] have shown a similar result that the group with lower NK cell counts has higher tumor incidence rate but the association is not significant.

Immune system strength is important in determining whether or not an immune system is able to eliminate a small tumor. Leukopenia, a low WBC count, is associated with several clinical disorders ^[88]. Numerical simulation shows that a weakened immune system with persistently low NK cell and WBC counts and a slower activation of CTL immune response may fail to control the growth of a tumor which a healthy immune system is able to control (Figures 7(b) and 10(c)). Finally, numerical simulation also shows that the effect of the estradiol concentration $E2$ on formation of a tumor is not significant which agrees with observations by previous research studies ^{[96,97,98,99]}.

Acknowledgments

This research work was supported by the Ministry of Science and Technology under the grant MOST107-2115-M-035-006.

Conflict of interest

The author declare there is no conflict of interest.