Bio-algorithms for the modeling and simulation of cancer cells and the immune response

Introduction

Cancer can be distinguished as a set of genetic diseases. It is caused by different mutations in specific genes and leads to uncontrolled cell proliferation and genomic instability. Breast cancer is the second leading cause of death worldwide in cancer-associated women [1], [2]. In spite of the latest diagnostic and therapeutic methods for breast cancer, it is developing in different regions and countries at a different speed. According to Peng Ji [3], there were 1,960,682 cases of breast cancer and 611,625 deaths globally in 2017. This number is increasing every year and breast cancer has become a major public health burden worldwide.

There are different types and stages of breast cancer in which Lobular Carcinoma in Situ (LCIS), Ductal Carcinoma in Situ (DCIS), Inflammatory Breast Cancer (IBC) and Metastatic cancer are common [4]. LCIS is a pre-cancerous condition that forms and is contained in the lobules. Invasive lobular carcinoma is a type of cancer that develops and break through the lobules with the potential to spread other areas of the body. DCIS is a type of cancer that forms in the lactiferous duct and considered non-invasive because it has not spread to any surrounding tissues. IBC is another uncommon but aggressive form of cancer, in which abnormal cells infiltrate the skin and lymph vessels of the breast. This type of cancer does not produce a distinct tumor or lump that can be felt and isolated within the breast. However, IBC grows rapidly and requires aggressive treatment.

Tumor is a mass of abnormal tissue and have usually two types: benign and malignant [5]. Benign is a non-cancerous tumor and not aggressive to surrounding tissues. Malignant tumors are cancerous tumor and aggressive because they invade and damage surrounding tissues. Metastatic cancer is when cancerous cells of malignant tumor spread to other parts of the body, usually through the lymph system, and form a secondary tumor [6]. The cancerous cells produce different signals that activate our immune system to produce antitumor response. The immune system consists of two types of immunity: innate immunity and acquired immunity or adaptive immunity. Innate immunity is the first and rapid response of our immune system to any pathogen. It is composed of dendritic cells, macrophages, natural killer cells, eosinophils, basophils, neutrophils, and mast cells. Adaptive immunity is activated with the stimulation of innate immune response and it is composed of T helper cells (CD4⁺), cytotoxic T cells (CD8⁺), and B cells. Cytotoxic T lymphocytes (CTLs) are preferred immune cells for targeting tumors. These cells express CD8 coreceptor and are recognised as frontline defensive cell types to fight against cancer progression. The function of CTLs is more crucial than the function of CD4⁺ T cells for the significant killing of cancer cells. The number of cytotoxic cells within the tumor microenvironment is an important foretelling marker for cancer and tumors can be classified on the base of the high and low number of CTLs [7].

Tumor cells are killed by CTLs. The majority of T cells circulate in the blood and lymph. Most are not cytotoxic and called naive T cells. The high concentration of naive T cells can be found in the lymph nodes which are scattered throughout the body. These naive T cells have the potential to become cytotoxic when activated [8]. The human body creates millions of variant T cells each with a unique ability to recognize different pathogens. But T cells cannot recognize pathogens on their own. They first encounter the pathogen by other immune cells called antigen-presenting cells (APCs) [9].

APCs are one of the first cells to detect the tumor in our body. Their role is surveillance, inspecting the environment, and taking back a sample to the nearby lymph node. Once in a lymph node, APCs search for T cells that will recognize their antigen complexed. T cells recognize the Golgi apparatus carried by APCs [10]. The reorganization is controlled by molecular interactions of these cells. APCs produce a protein called major histocompatibility complex (MHC) on its cell surface. The MHC contains virus fragment each T cell produces a unique receptor that recognizes the virus peptides MHC combinations [11]. The presence of CD8 molecules defines this class of T cells. The receptors on this particular T cells bind to virus peptides MHC complex, co-receptor CD8 binds the complex, eventually leading to T cell activation. The activated T cells divide and produce a clone of itself. The cells gain cytotoxic capability as they grow. Chemical cues attract the CTLs to the tumor site. CTL specifically recognize infected cells. The reorganization is controlled by the molecular interaction between the cells. The infected cells displayed the same virus peptide MHC complex. T cells receptor binding and clustering trigger the killing process. Cytotoxic granules are transferred to kill infected cells [12].

Mathematical modeling is commonly used in physical and engineering science to understand and explain industrial processes and complex systems [13], [14], [15], [16]. It has lots of applications in biological processes, tissue growth and development, and cancer modeling. Mathematical models provide an analytical framework to describe the interaction of immune cells and tumor. Mathematical models can be distinguished into deterministic models and stochastic models. In deterministic models, the values of dependent variables are completely determined by using values of parameters of the model, however, stochastic models show randomness such that their outcome can be viewed as probability distribution instead of unique value. Different mathematical models are reported in the literature to describe the interaction of tumor and immune cells [17], [18], [19], [20], [21], [22]. However, some of them have not parametric values verified from experimental data and some models are very complex. Complex mathematical models have more parameters and thus less well-parameterized [23]. In this article, a novel deterministic model is presented to describe the interaction of cancerous tumor and CTLs analytically. The main objective of this article is to design a less-complicated mathematical model that describe the dynamics of tumor and CTLs in breast cancer with parametric values verified from the experimental data. The values are parameters are estimated from experimental data of two research studies [24], [25] and Monte Carlo global sensitivity analysis [26] is performed to check uncertainty of the parametric values. The results produced by the deterministic model are further compared with clinical phenomena. To provide detail, we look at mathematical modelling and justification if its term in Section 2. It also includes the process of non-dimensionalization of model parameters, global sensitivity analysis of parameters, and numerical simulations of the model. In Section 3, the results are discussed and the article is concluded.

Material and methods

Mathematical modeling and parameter estimation

We design deterministic model under following assumptions: tumor cells are growing logistically in the absence of immune response [27]; cytotoxic T lymphocytes can kill tumor cells [28]; tumor cells can activate naive and noncytotoxic cells [29]; after activation of cytotoxic T cells, they will grow logistically; cytotoxic T cells will become inactive after some number of interactions with tumor cells [30]. The model describing the kinetics of tumor cell and the cytotoxic T cell is represented in Figure 1 and modelled by a system of ordinary differential equations as

(1)dX1dt=k1X1(1−X1k2)⏟Logistic growth of tumor cells−k3X1(X2k4+X2)⏟Tumor cells killed by cytolysis,

(2)dX2dt=k5X2(1−X2k6)(X1k7+X1)⏟Activation of CTLs by tumor and its logistic growth−k8X1X2⏟Inactivation of CTL after interaction with tumor−k9X2⏟Natural degradation,

where k₁, k₂, k₃, k₄, k₅, k₆, k₇, k₈, and k₉ are positive constants. It has been observed from numerous study in vitro and in vivo that tumor cells are growing exponentially when their population size is small and their growth is slowed at large population size of tumor. Under this observation, tumor growth is considered to follow logistic curve with intrinsic growth rate of k₁ and maximum carrying capacity of k₂. We used experimental literature, where available, to quantify parameters of the proposed model. However, as such data is sparse, we also rely on previous estimates of some rates and values from other modeling articles. The value of growth rate k₁ is estimated from the data [25] while the value of carrying capacity is taken from the studies [31], [32], [33]. The values of parameters k₃ and k₄ are estimated from the experimental data [34]. The parameters of the CTL equation are estimated from a clinical study of tumor cells and T lymphocyte in breast cancer patients [24]. The model equations are non-dimensionalized using

x1=X1c1, x2=X2c2, τ=c3t,

where c₃ = c₁k₃. The values of c₁, c₂ and c₃ are taken as suggested in Ref. [35]. The dimensionless model of system (1)–(2) after replacing τ by t is given below:

(3)dx1dt=f(x1,x2)=α1x1(1−x1α2)−α3x1(x2α4+x2),

(4)dx2dt=g(x1,x2)=α5x2(1−x2α6)(x1α7+x1)−α8x1x2−α9x2,

where

α1=k1c3,α2=k2c1,α3=k3c3,α4=k4c2,α5=k5c3,α6=k6c2,α7=k7c1,α8=k8c1c3,α9=k9c3

Figure 1:

A schematic diagram showing interaction between tumor cell and cytotoxic T lymphocyte.

The values of dimensionless parameters are taken from dimensional parameters and given in Table 1 (Appendix).

Steady states and nullclines

Let us consider steady states of the designed model. All the parameters of the model are supposed to be non-negative and our solution of practical interest is also non-negative values of tumor and CTLs. There always exists a trivial equilibrium point E₀ = (0, 0), which is a saddle. To examine other steady states of the system, we consider nullclines that is the curves along which dx1dt=dx2dt=0. The nullclines of the designed system (3–4) are as follows:

(5)x2=α1α4α3(1−x1α2)1−α1α3(1−x1α2)=F1(x1),

(6)x2=α6(α7+x1)α5x1(α5x1a7+x1−α8x1−α9)=F2(x1).

We consider only positive quadrant and there exist a steady state E₁ at the intersection of nullclines as shown in Figure 2. By using the dimensionless values of parameters extracted from data and listed in Table 1 (Appendix), the steady state E₁ is asymptotically stable as shown in Figure 2. According to Dulac–Bendixson criterion of closed orbit [36], there exist no closed orbit for positive values of x₁ and x₂. To illustrate, consider the function B(x1,x2)=1x1x2 and calculate:

∂(Bf)∂x1+∂(Bg)∂x2=−(α1α2x2+α5α5(α7+x1))=−L(x1,x2).

Figure 2:

Graphs showing the nullclines of designed model.

Since parameters are positive, L > 0 for positive values of x₁ and x₂. Hence Dulac–Bendixson criterion is satisfied. It follows that there is no limit cycles and no Hopf bifurcation giving rise to limit cycles occur over the domain of interest.

Global sensitivity analysis

Sensitivity analysis is the study about how uncertainty and variability of model output are affected by the uncertainty and variability of model input. It helps us in determining which parameter impact on the model output most. We can monitor that parameter to control model output. We can improve model reliability and performance. Here, we perform global sensitivity analysis to check the influence of parameters on the tumor population. We are interested in how the population of tumor cells is increasing with the change of parametric values. We consider four parameters of the model (α1, α4, α5, and α7) and examine their impact on tumor population. We take 100 sample values of each parameter for design space. All sample values are randomly generated with a uniform probability distribution. All sample values of each parameter are shown through a scatter plot in Figure 5 (Appendix).

From the scatter plot shown in Figure 3, we can see that there is no trend showing that tumor cells are increasing with the increase of parametric values of α1 and α4 (see Figure 3). However, tumor cells are increasing with an increase of α7 and decreasing with the increase of α5 as shown in Figure 3. This indicates that tumor population is sensitive to parametric values of α1 and α2. We also perform some statistics such as correlation, regression and partial correlation on the given raw data. It is shown on tornado plot through parameters and their rank (see Figure 6 in Appendix). From tornado plot, we can see that tumor population is inversely proportional to α5 and directly proportional to α7. The parameters α1 and α4 also have a direct impact with less influence.

Figure 3:

Scatter plots showing the results of global sensitivity analysis of parameters α1, α4, α7, and α5 respectively.

From the results of sensitivity analysis, it can be depicted that tumor population is sensitive to parameters α5 and α7. Therefore, we analyze the efficacy of these parameters by the contour plot to increase tumor cells (Figure 7). The parameter α5 is taken along the x-axis while the parameter α7 is taken along the y-axis. The colour bar shows the population of tumor cells. Contour plot indicates that tumor cells are at their peak with the smallest value of α5 largest value of α7. But with the larger value of α5, the influence of α7 is becoming futile.

Numerical simulations

The designed model is simulated in MATLAB Simulink environment with dimensionless parametric values listed in Table 1. The dimensionless value of initial tumor cells is 10² and the dimensionless initial number of CTLs is 10¹. Figure 4 demonstrates the dynamics of tumor with different values of α5. The arrow symbol indicates that tumor cells are decreasing with the increase of α5. Similarly, Figure 4 illustrates the tumor population with different values of α7 and it can be seen that tumor cells are increasing with an increase of α7. Hence, numerical simulations strengthen the results of sensitivity analysis. State profiles of both variables are shown in Figure 8 (Appendix). It can be seen that transients in the vicinity of E₁ exhibit decaying oscillations which are often in breast cancer patients.

Figure 4:

State profiles of tumor cells for different values of α5 and α7.

Results and discussion

By using the designed model, we can make several biological predictions about the relation of breast tumor and CTLs. The designed model has two equilibrium points namely tumor free equilibrium point (E₀) and non-zero equilibrium (E₁). The equilibrium point E₀ is a saddle while E₁ is asymptotically stable. A phase portrait of the designed, model with parametric values estimated from experimental data, is presented in Figure 2 to demonstrate asymptotic stability of E₁. Figure 8 shows decay oscillations in the population of CTLs. This cyclic fluctuation has a close agreement with experimental studies of breast cancer. In literature, it is reported that immune cells can eliminate tumors at early stages and most of the tumor cells are destroyed within two days [37], [38]. From the numerical simulations of the model shown in Figure 8, we can see that the population of CTLs is growing rapidly at early stages and showing oscillations in the neighboring of E₁. The tumor cells going to decrease with this rapid response of CTLs. From the same figure, it can be seen that tumor is not removed completely. If a tumor has a size less than 2 mm in the diameter and population-level remains less than 6 × 10⁵, then tumor remains small and stable [39]. This phenomenon is termed as “cancer without disease” and many studies have suggested that microscopic tumor never progresses to invasive [40], [41], [42]. However, our sensitivity analysis reveals that the tumor population is very sensitive to parameter α7. If we shall take a large value of α7, the tumor cells will decrease and go to zero as shown in Figures 3 and 4.

Tumor cells have the ability to escape the immune cells. Suppose that the innate immune system delays to activate CTLs and tumor cells are growing. In this situation, the model suggested that tumor cells will grow exponentially at an earlier stage, but with the activation of CTLs, tumor cells will decrease as shown in Figure 9. Also, if the innate immune system fails to activate CTLs, then model suggested that tumor cells will achieve its maximum carrying capacity at dimensionless time 30 units as shown in Figure 9. Nullclines analysis of model suggested that tumor nullcline is much dependent on the parameter α4. For the larger value of this parameter, the steady-state E₁ vanish and tumor cells approach to their maximum carrying capacity. Figure 9 shows that CTLs are not able to control tumor when there are more than approximately 1.2 × 10⁸ (x₁ = 120) tumor cells. This follows that tumor is not controlled by immune cells and there is a need for chemotherapy or immunotherapy.

Although the designed mathematical model is working perfectly to determine the interaction of tumor cells and CTLs, however, there are also some flaws. Since our objective is to describe tumor and CTLs dynamics in a less complicated way, hence designed model does not address all process of immune cells. The current model is not dealing with self-regulatory of CTLs. Moreover, most of the biological systems show noisy behaviour due to the fluctuation of its components, this phenomenon is known as stochasticity. We have yet to resolve this issue in our model. Despite all, the designed model fits the empirical data.

Conclusions

Breast cancer is the leading cause of death in women worldwide. Early and effective response of immune cells is a promising avenue for increasing survival times as well as improving the quality of the life for the breast cancer women. Finding algorithms to predict the growth of breast tumor has piqued the interest of researchers ever since the early days of cancer. The present manuscript is dealing with mathematical modeling of tumor growth and function of CTLs in breast cancer. The proposed mathematical model incorporates tumor-CTLs interaction in deterministic setting and provides a good fit with experimental data resulting from different clinical studies. We determine the equilibrium points of the model and nullclines analysis is performed to discuss stability of steady states. We have found two equilibria. One is tumor-free and shown to be saddle point. The other equilibrium point is shown to be asymptotically stable. The sensitivity analysis of parameters yields the results that are intuitively reasonable. This analysis highlights the important parameters that can be target to reduce tumor size. The results of numerical simulations of the proposed model have been validated by comparing the outcomes of clinical studies. The experimental and clinical studies have suggested that CTLs can eliminate small tumor and also larger tumor by stimulating their cytotoxicity using tumor vaccines, adoptive cell transfer therapy, or checkpoint antibodies [37], [43], [44], [45]. Numerical simulations show that CTLs are enable to remove small tumor but fail to remove larger tumor. However, in future, we are interested to extend our model by adding the effects of chemotherapy and immunotherapy.