Parameter estimation and sensitivity analysis for a model of tumor–immune interaction in the presence of immunotherapy and chemotherapy
Hesham A. Elkaranshawy1 &
Ahmed M. Makhlouf ORCID: orcid.org/0000-0002-5538-37371
Journal of the Egyptian Mathematical Society volume 30, Article number: 8 (2022) Cite this article
A mathematical model has been utilized to examine the interaction between tumor cells and immune cells. In this model, the immune cells include natural killer cells, circulating lymphocytes, CD8+T cells, CD4+T cells, and cytokines. The model not only represents the traditional role of CD4+T cells in activating CD8+T cells but also illustrates its role in killing the tumor via the secretion of cytokines. Besides, treatments with both chemotherapy and immunotherapy are considered. However, since this model was not fitted to experimental data before, parameter estimation is performed to fit the model with experimental data, first. The estimation is validated to verify the correctness of the model using the experimental data for the tumor growth. Second, numerical experiments are performed using a set of human data. Results show the mutual relations between tumor cells, and body immune cells in the absence and in the presence of therapy. Results also show that CD4+T cells could play a crucial role in immunotherapy. Third, sensitivity analysis is performed by calculating the normalized sensitivity coefficients to identify the relative influence of body parameters on the tumor cell population. The obtained results provide a tool to identify which parameters should be increased or decreased before treatment to get the optimal immune response.
According to the World Health Organization (WHO), cancer is one of the leading causes of death worldwide [1]. In the field of cancer treatment, there are many methods like surgery, radiation, hormonal therapy [2], virotherapy [3], chemotherapy, and immunotherapy. Immunotherapy is based on enhancing the effectiveness of the patient's immune system to identify and kill cancer cells. Recently, immunotherapy has been approved for the treatment of many types of cancers [4, 5]. The use of immunotherapy is safer than using chemotherapy because chemotherapy has a drawback in killing immune cells in addition to killing the tumor cells [6, 7].
Mathematical modeling can be applied to many systems, physical, chemical, biological, or else to explain them, to examine the effect of their different components, and to predict their behavior [8,9,10,11,12,13,14,15,16,17]. An increasing number of research works focused on developing mathematical models analyzing tumor growth. A model, describing tumor growth due to its internal pressure in the form of a system of partial differential equations, was proposed by Jones et al. [18]. Similar models were proposed by Tao et al. [19], Wei and Cui [20], Frieboes et al. [21], Lee et al. [22], Zhang et al. [23], Knopoff et al. [24], and Kuznetsov et al. [25]. These models are in the form of a system of differential equations also. Some proposed models involved time delayed differential equations [26] like, Rihan et al. [27], Villasana and Radunskaya [28], Xu and Bai [29] and Rihan et al. [30]. The last one includes performing parameter estimation and sensitivity analysis of the time-delayed model. Some research works proposed models including stochastic differential equations like Rihan and Rajivganthi [31] and Nutini and Sohail [32]. Other research works focused on finding suitable mathematical models that describe the interaction between the immune system and tumor cells to develop a suitable treatment strategy for each case [6, 33,34,35]. Arabameri et al. [5] proposed a model for immunotherapy using mature dendritic cells. A model was proposed by Eftimie et al. [36] and was modified by Anderson et al. [37] and by Hu and Jang [38] to consider the effect of CD4+T cells on the system. Many research works are introduced in this field; however, most of them had considered the role of some immune cells and neglected the role of the others. Hence, a mathematical model that covered the role of most of the immune system cells and incorporated immunotherapy and chemotherapy has been developed [7, 39].
In many research works, sensitivity analysis is performed to investigate the effect of changing different parameters on the output of the mathematical model [40]. Changing system parameters can drive the system from an unstable equilibrium state to a stable equilibrium state [7, 35]. Sensitivity analysis had been performed by changing each parameter by a certain small value at a specific point in time [6, 38, 41, 42]. The sensitivity of the output to each parameter is determined, and hence, the relative effectiveness of parameters can be specified. Alternatively, a normalized sensitivity coefficient can be calculated to obtain a dimensionless value for the sensitivity coefficients [40,41,42,43,44]. While the numerical method identifies the most effective parameters at a specific instant of time, the normalized sensitivity coefficient method recognizes these parameters over a wide time interval.
In this paper, a mathematical model depending on the model in [7, 39] is utilized. This model was not fitted to experimental data previously. Hence, in this study, parameter estimation [45, 46] is performed to fit the model to experimental data, numerical simulations are executed, and sensitivity analysis is performed. The model represents the interaction between immune system cells and tumor cells. The immune system cells considered are natural killers, circulating lymphocytes, CD8+T cells, CD4+T cells, and cytokines. The model not only incorporates immunotherapy and chemotherapy but also the mutual relations between these therapies, tumor cells, and body immune cells. The traditional role of CD4+T cells in activating CD8+T cells is considered as well as its role in killing the tumor via secretion of cytokines. A sensitivity analysis is implemented to the model to investigate the most effective parameters on the tumor cell population. Sensitivity analysis is performed by calculating the normalized sensitivity coefficients and it introduces the order of the most effective parameters, so it answered the question: which parameter can be changed to get the largest effect on the tumor size?
This study aims to verify that the proposed mathematical model can represent the interaction between the tumor cells and the immune system for a practical case. Hence, the model is used with confidence to perform proposed simulations.
Experimental data for the tumor growth are used in this study to verify that the model can fit experimental data. The experimental data, used in this paper, are obtained from a published work.
Parameter estimation is performed to fit the model to the given experimental data using the built-in functions "ParametricNDSolveValue" and "FindFit" in Mathematica commercial program. Consequently, numerical simulations are performed using a set of human data to show the mutual relations between therapies, tumor cells, and body immune cells. The solutions are analyzed for three cases: without any treatment, with continuous treatment, and with pulsed treatment. Finally, sensitivity analysis is performed to identify the relative influence of body parameters on the tumor cell population. Normalized sensitivity coefficients are estimated.
Mathematical model
A mathematical model for the interaction between the body's immune system and the tumor is considered. The model is proposed by Makhlouf et al. [7], and it depends upon the work of de Pillis et al. [6, 33,34,35], Eftimie et al. [36], Anderson et al. [37], and Hu and Jang [38]. The model represents the status in the absence of therapies as follows:
$$\frac{\mathrm{d}T}{\mathrm{d}t}=a T \left(1-b T\right)-c N T-D T-\frac{{c}_{1} T}{{a}_{1}+ T} I$$
$$\frac{\mathrm{d}N}{\mathrm{d}t}=e C-f N+\frac{g {T}^{2}}{h+ {T}^{2}} N-p N T$$
$$\frac{\mathrm{d}L}{\mathrm{d}t}=-m L+\frac{j {D}^{2} {T}^{2}}{k+ {D}^{2} {T}^{2}} L-q L T+\left({r}_{1} N+{r}_{2} C\right) T-u N {L}^{2}+\frac{{p}_{i} I}{{g}_{i}+I}L$$
$$\frac{\mathrm{d}Y}{\mathrm{d}t}=\frac{{\beta }_{1} T}{{\alpha }_{1}+ T} I-{\mu }_{1} Y-{\delta }_{2} T Y$$
$$\frac{\mathrm{d}C}{\mathrm{d}t}=\alpha -\beta C$$
$$\frac{\mathrm{d}I}{\mathrm{d}t}=-{\mu }_{i}I+\frac{{\beta }_{2} T}{{\alpha }_{2}+ T} Y$$
where \(T\left(t\right)\) is the tumor cells population, \(N\left(t\right)\) is the natural killer cells population, \(L\left(t\right)\) is the CD8+T cells population, \(Y\left(t\right)\) is the CD4+T cells population, \(C\left(t\right)\) is the circulating lymphocytes cell population not including natural killer cells and active CD8+T and CD4+T cells and \(I\left(t\right)\) is the concentration of the cytokine. \(D\) is the fractional kill rate, which defined as follows:
$$D=d\frac{{\left(L/T\right)}^{l}}{s+ {\left(L/T\right)}^{l}}$$
In Eq. (1), \(a T \left(1-b T\right)\) represents the tumor growth. The term \(-c N T\) represents the natural killer-induced tumor death. The term \(-D T\) represents the tumor lysis by CD8+T cells. The term \(- \frac{{c_{1} T}}{{a_{1} + T}}I\) represents the role of cytokine in killing the tumor. This last term reflects the indirect role of CD4+T cells in killing the tumor. CD4+T cells cannot kill tumor cells directly, like CD8+T cells, but it does this through the cytokines that they produce. In Eq. (2), the term \({e}_{1} C\) represents the production of natural killer cells from circulating lymphocytes. The term \(-f N\) represents the turnover of natural killer cells. The term \(\frac{g {T}^{2}}{h+ {T}^{2}}N\) represents the recruitment of natural killer cells induced due to the existence of the tumor. The term \(-p N T\) represents the death of natural killer cells by the exhaustion of tumor-killing resources. In Eq. (3), the term \(-m L\) represents the inactivation of CD8+T cells. The term \(\frac{j {D}^{2} {T}^{2}}{k+ {D}^{2} {T}^{2}}L\) represents the recruitment in CD8+T cell population due to the presence of the tumor. The term \(-qLT\) represents the death of CD8+T cells by the exhaustion of tumor-killing resources. The term \({r}_{1}NT\) represents the stimulation of CD8+T cells by natural killers-lysed tumor cell debris. The term \({r}_{2}CT\) represents the native CD8+T cells activation in the general lymphocyte population. The term \(-u N {L}^{2}\) represents the regulation and suppression of CD8+T cell activity by the natural killer cells. The term \(\frac{{p}_{i} I}{{g}_{i}+I}L\) represents the activating effect of cytokines on CD8+T cells. In Eq. (4), the term \(\frac{{\beta }_{1} T}{{\alpha }_{1}+ T}\) represents the activating effect of cytokines on CD4+T cells due to the existence of the tumor. The term \(-{\mu }_{1} Y\) represents the decaying of CD4+T cells which has a constant rate. The term \(-{\delta }_{2} T Y\) is a decaying term also which exists due to the interaction of CD4+T cells with the tumor. Equation (5) represents the circulating lymphocyte cell population. Circulating lymphocytes are assumed to be generated at a constant rate \(\alpha\) and decayed with a constant natural death rate \(\beta\). In Eq. (6), the term, \(-{\mu }_{i}I\) represents the decaying of cytokine with a constant rate \({\mu }_{i}\). The term \(\frac{{\beta }_{2} T}{{\alpha }_{2}+ T} Y\) represents the activating role of the CD4+T cells on cytokines due to the existence of the tumor [7, 39, 44].
For a set of human data, Table 1 includes the values of all the parameters included in the model [7]. Parameters \(a\) and \(b\) will be estimated in the next section to fit the model with experimental data. The definitions of all the parameters and a detailed description of each term in the model are included in [7]. Equation (6) shows that CD4+T cells activate the cytokine which in term activates CD8+T cells as can be noticed in Eq. (3) and both cytokine and CD8+T cells assist in killing tumors as presented in Eqs. (1) and (7).
Table 1 Values of patient's data
To include the effect of immunotherapy and chemotherapy, the model is modified to be:
$$\frac{\mathrm{d}T}{\mathrm{d}t}=a T \left(1-b T\right)-c N T-D T-\frac{{c}_{1} T}{{a}_{1}+ T} I-{K}_{T}\left(1-{e}^{-M}\right)T$$
$$\frac{\mathrm{d}N}{\mathrm{d}t}=e C-f N+\frac{g {T}^{2}}{h+ {T}^{2}} N-p N T-{K}_{N}\left(1-{e}^{-M}\right) N$$
$$\frac{\mathrm{d}L}{\mathrm{d}t}=-m L+\frac{j {D}^{2} {T}^{2}}{k+ {D}^{2} {T}^{2}} L-q L T+\left({r}_{1} N+{r}_{2} C\right) T-u N {L}^{2}+\frac{{p}_{i} I}{{g}_{i}+I}L-{K}_{L}\left(1-{e}^{-M}\right) L+{v}_{L}\left(t\right)$$
$$\frac{\mathrm{d}Y}{\mathrm{d}t}=\frac{{\beta }_{1} T}{{\alpha }_{1}+ T} I-{\mu }_{1} Y-{\delta }_{2} T Y+{v}_{Y}\left(t\right)$$
$$\frac{\mathrm{d}C}{\mathrm{d}t}=\alpha -\beta C-{K}_{C}\left(1-{e}^{-M}\right) C$$
$$\frac{\mathrm{d}I}{\mathrm{d}t}=-{\mu }_{i}I+\frac{{\beta }_{2} T}{{\alpha }_{2}+ T} Y+{v}_{I}\left(t\right)$$
$$\frac{\mathrm{d}M}{\mathrm{d}t}=-\gamma M+{v}_{M}\left(t\right)$$
In Eq. (14), the chemotherapy dose is represented by the function \({v}_{M}\left(t\right)\). The other term, \(-\gamma M\), represents the decaying in chemotherapy concentration which occurs with a constant rate \(\gamma\). The effect of chemotherapy on killing the tumor cells is represented in Eq. (8) by the term \(-{K}_{T}\left(1-{e}^{-M}\right) T\). The effect of chemotherapy on killing immune cells is represented in Eqs. (9), (10), and (12) by the terms \(-{K}_{N}\left(1-{e}^{-M}\right) N\), \(-{K}_{L}\left(1-{e}^{-M}\right) L\) and \(-{K}_{C}\left(1-{e}^{-M}\right) C\). Immunotherapy is represented in the model by the terms \({v}_{I}\left(t\right)\), \({v}_{L}\left(t\right)\) and \({v}_{Y}\left(t\right)\). Respectively, they are the concentrations of IL-2 cytokine therapy, CD8+T, and CD4+T adoptive immunotherapy [7].
For a set of human data, the value of each parameter is included in Table 1. The definitions of all the parameters and a detailed description of each term in the model are included in [7]. In this model, Eq. (13) shows that CD4+T cells activate the cytokine which in term activates CD8+T cells as can be noticed in Eq. (10) and both the cytokine and CD8+T cells assist in killing tumors as presented in Eq. (7) and Eq. (8).
We would like to mention here, that the cytokine used in the immunotherapy is the IL-2 cytokine because it and interferon-alpha (INF-\(\alpha\)) are the only approved cytokines from the Food and Drug Administration (FDA) in the USA. Additionally, IL-2 achieved long-term patient survival compared to the INF-\(\alpha\) [47].
Parameter estimation
In this section, two parameters are estimated to fit the basic model with experimental human data. The experimental data utilized in this section were acquired by Worschech et al. [48] and have been used by many researchers like Murphy et al. [49]. This set of human data consists of 14 points. The first 9 points are employed in the fitting process and the last 5 points to verify the values of the parameters estimated. The values of the parameters in Table 1 are employed while the tumor growth rate \(a\) and the reciprocal of the tumor carrying capacity \(b\) are estimated. The initial conditions assumed in this simulation are \(N\left(0\right)={10}^{4}\) cells, \(L\left(0\right)={10}^{2}\) cells, \(Y\left(0\right)={10}^{5}\) cells and \(C\left(0\right)={10}^{8}\) cells. The tumor size was measured in units of \({\mathrm{mm}}^{3}\) in the data set of Worschech et al. [48]. Therefore, any cell unit is converted to \({\mathrm{mm}}^{3}\). The first 9 points of the experimental data are employed in the fitting process. The intent is to determine the values of the two parameters \(a\) and \(b\) that give the best fit of the model to the given experimental data. To estimate the required parameters, the Mathematica commercial program is used to generate a parametric solution for the model in terms of the parameters a and b. This has been performed using the built-in function "ParametricNDSolveValue." Using this parametric function, the method of Levenberg–Marquardt [50, 51] is applied using the built-in Mathematica function "FindFit," to estimate these two parameters. The obtained value for \(b\) is converted to the units of \({\mathrm{cell}}^{-1}\). Hence, the estimated values for the two parameters are \(a=2.68208\times {10}^{-2}{\mathrm{ day}}^{-1}\) and \(b=1.05927\times {10}^{-10}{\mathrm{ cell}}^{-1}\). Figure 1 shows the curve fitting between the model and the first 9 points of the experimental data. To verify the results, the estimated parameters are incorporated in the model and the simulation is performed up to \(t=120\,\mathrm{ days}\). As shown in Fig. 1, the last 5 points, which were not used in the parameter estimation, are in good agreement with the model simulation.
Model fitting to experimental data
Simulation results
In this section, the mathematical model is solved numerically for proposed illustrations, and the solutions are interpolated. The solutions are analyzed for three cases: without any treatment, with continuous treatment, and with pulsed treatment. The interaction between the chemotherapies, the tumor cells, and the body immune cells is investigated. In Fig. 2a, the system (1) to (6) is solved for patient's data with initial tumor size \(T\left(0\right)=2\times {10}^{5}\), initial natural killer cells \(N\left(0\right)={10}^{4}\), initial CD8+T cells \(L\left(0\right)={10}^{2}\), initial CD4+T cells \(Y\left(0\right)={10}^{5}\), and initial circulating lymphocytes \(C\left(0\right)={10}^{8}\). With these conditions, the body's immune system overcomes the tumor, as can be seen in Fig. 2a. Using the same initial conditions but with a larger tumor, \(T\left(0\right)=2\times {10}^{7}\), the tumor is not eradicated, and treatment is needed as shown in Fig. 2b.
Simulation of the no-treatment case for patient's data
In the second case, the model (8)–(14) is simulated with continuous immunotherapy using the patient's data. The initial conditions, in this case, are: \(T\left(0\right)=2\times {10}^{7}\), \(N\left(0\right)={10}^{4}\), \(L\left(0\right)={10}^{2}\), \(Y\left(0\right)={10}^{5}\) and \(C\left(0\right)={10}^{8}\). In Fig. 3a, a continuous CD8+T cell therapy \({v}_{L}={10}^{4}\) together with cytokine IL-2 treatment \({v}_{I}={10}^{4}\) are applied. However, the tumor could not be eradicated. When this combined treatment is replaced by CD4+T cell therapy \({v}_{\mathrm{Y}}={10}^{7}\), the tumor is eliminated as shown in Fig. 3b.
Simulation of continuous treatment case for patient's data
In the third case, the model is solved for patient's data with the experimental value of initial tumor size \(T\left(0\right)=2.41758\times {10}^{8}\). Other initial conditions are: \(N\left(0\right)={10}^{4}\), \(L\left(0\right)={10}^{2}\), \(Y\left(0\right)={10}^{5}\), \(C\left(0\right)={10}^{8}\). In Fig. 4a, a continuous CD8+T cell therapy \({v}_{L}={10}^{8}\) together with cytokine IL-2 treatment \({v}_{I}={10}^{8}\) and CD4+T cell therapy \({v}_{Y}={10}^{8}\) are applied. However, the tumor could not be eradicated using this large dose of immunotherapy. The immunotherapy is failed in this case because the initial tumor is very large. Hence, chemotherapy is needed. In Fig. 4b, three pulses of chemotherapy of intensity \({v}_{\mathrm{M}}=5\) starting from the fifth day with 5 days period could eradicate the tumor. Hence, for the existing large tumor, immunotherapy failed to overcome the tumor while chemotherapy could overcome it.
Simulation of treatment for the experimental tumor size
The numerical experiments presented show interaction between tumor cells and the immune cells. They also show the crucial rule of CD4+T cells in immunotherapy. CD4+T cell therapy can be used as an alternative to both CD8+T cells and IL-2 cytokine therapies with more efficiency in tumor eradication. In the last case, we show that, for the existing large tumor, immunotherapy failed to overcome the tumor while chemotherapy could overcome it.
The sensitivity coefficient \({S}_{ij}\) for the system output \({y}_{i}\) w.r.t. the parameter \({\lambda }_{j}\) is given by [41,42,43,44]:
$${S}_{ij}=\frac{\partial {y}_{i}}{\partial {\lambda }_{j}}, \quad i=\mathrm{1,2},\dots,n,\quad j=\mathrm{1,2},\dots ,m$$
In the model (1) to (6), we have 31 parameters, they are \(a\), \({a}_{1}\), \(b\), \(c\), \({c}_{1}\), \(d\), \(e\), \(f\), \(g\), \({g}_{i}\), \(h\), \(j\), \(k\), \(l\), \(m\), \(p\), \({p}_{i}\), \(q\), \({r}_{1}\), \({r}_{2}\), \(s\), \(u\), \(\alpha\), \({\alpha }_{1}\), \({\alpha }_{2}\), \(\beta\), \({\beta }_{1}\), \({\beta }_{2}\), \({\mu }_{1}\), \({\mu }_{i}\) and \({\delta }_{2}\). We have 6 system outputs. They are \(T\), \(N\), \(L\), \(Y\), \(C\) and \(I\). We focused in this study on the effect of each parameter on the tumor size; hence, the sensitivity coefficients are calculated for the tumor size \(T\) only. The sensitivity coefficients can be calculated by solving the following system of ordinary differential equations [41,42,43,44]:
$$\frac{\mathrm{d}{\varvec{S}}}{\mathrm{d}t}={{\varvec{f}}}_{\lambda }+{\varvec{J}} {\varvec{S}}$$
where \({\varvec{S}}=\left[\begin{array}{cccc}{S}_{11}& {S}_{12}& \cdots & {S}_{1m}\\ {S}_{21}& {S}_{22}& \cdots & {S}_{1m}\\ \vdots & \vdots & \ddots & \vdots \\ {S}_{n1}& {S}_{n2}& \cdots & {S}_{nm}\end{array}\right]\), \({{\varvec{f}}}_{\lambda }=\left[\begin{array}{cccc}\frac{\partial {f}_{1}}{\partial {\lambda }_{1}}& \frac{\partial {f}_{1}}{\partial {\lambda }_{2}}& \cdots & \frac{\partial {f}_{1}}{\partial {\lambda }_{m}}\\ \frac{\partial {f}_{2}}{\partial {\lambda }_{1}}& \frac{\partial {f}_{2}}{\partial {\lambda }_{2}}& \cdots & \frac{\partial {f}_{2}}{\partial {\lambda }_{m}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial {f}_{n}}{\partial {\lambda }_{1}}& \frac{\partial {f}_{n}}{\partial {\lambda }_{2}}& \cdots & \frac{\partial {f}_{n}}{\partial {\lambda }_{m}}\end{array}\right]\), \({\varvec{J}}=\left[\begin{array}{cccc}\frac{\partial {f}_{1}}{\partial {y}_{1}}& \frac{\partial {f}_{1}}{\partial {y}_{2}}& \cdots & \frac{\partial {f}_{1}}{\partial {y}_{n}}\\ \frac{\partial {f}_{2}}{\partial {y}_{1}}& \frac{\partial {f}_{2}}{\partial {y}_{2}}& \cdots & \frac{\partial {f}_{2}}{\partial {y}_{n}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial {f}_{n}}{\partial {y}_{1}}& \frac{\partial {f}_{n}}{\partial {y}_{2}}& \cdots & \frac{\partial {f}_{n}}{\partial {y}_{n}}\end{array}\right]\) and the functions \({f}_{1}\), \({f}_{2}\), \({f}_{3}\), \({f}_{4}\), \({f}_{5}\) and \({f}_{6}\) are the right-hand sides of Eqs. (1)–(6).
The system of differential equations defined by (16) can be solved using a suitable numerical method with the initial conditions \({S}_{ij}\left(0\right)=0\). After that, the solution will be normalized to bound the range of the sensitivity functions using the following formula:
$${\left({S}_{ij}\right)}_{n}=\frac{{S}_{ij}}{\left(\frac{{y}_{i}}{{\lambda }_{j}}\right)}$$
Sensitivity analysis is performed for the patient's data with the value of initial tumor size given in the experiment [48], \(T\left(0\right)=2.41758\times {10}^{8}\), the initial natural killer cells \(N\left(0\right)={10}^{4}\), initial CD8+T cells \(L\left(0\right)={10}^{2}\), initial CD4+T cells \(Y\left(0\right)={10}^{5}\) and initial circulating lymphocytes \(C\left(0\right)={10}^{8}\). The normalized sensitivity function is plotted in Fig. 5. It shows that the most effective parameters are \(a\), \(b\) and \(l\), respectively.
Normalized sensitivity functions for the experimental data
In this paper, a mathematical model for the interaction between tumor cells and immune cells has been employed. The model has involved the role of natural killers, circulating lymphocytes, CD4+T cells, CD8+T cells, and cytokine. The model has included both chemotherapy and immunotherapy. One of the main points that have been addressed in this research work is to investigate the role of CD4+T cells and their interaction with tumors and with the treatment process. All work performed in this paper is summarized in the flowchart shown in Fig. 6.
Flowchart for the results
Experimental data have been used to estimate the tumor growth rate \(a\) and the reciprocal of the tumor carrying capacity \(b\). The estimated parameters have been verified against additional points in the experimental data set. With the estimated parameters, three numerical experiments have been performed. The first case is a no-treatment case, the second case is a continuous immunotherapy case, and the third case is a pulsed treatment case. In the first case, the immune system could eradicate a small tumor; however, the immune system could not eradicate a large tumor, and treatment is needed. In the second case, presumed CD8+T cells and IL-2 continuous therapy could not eliminate the tumor; however, a proposed CD4+T cells continuous therapy has succeeded in eliminating the tumor. In the third case, the tumor size in the experiment is considered, the immunotherapy failed in eradicating the tumor. Pulsed chemotherapy effectively eliminates the tumor in this case. These numerical experiments have shown the crucial role of CD4+T cells in immunotherapy. They also have shown that, for the existing large tumor, immunotherapy failed to overcome the tumor while chemotherapy could overcome it.
A sensitivity analysis is performed to the basic model, for the set of a patient's data with the value of initial tumor size given in the experiment, to specify the most effective parameters on the tumor cell population. It has been performed by calculating the normalized sensitivity coefficient. The list of the most effective parameters can be used in designing the most effective treatment strategy to eliminate the tumor and to avoid its recurrence after the treatment. Results revealed that the most effective parameters are the tumor growth rate \(a\), the reciprocal of the tumor carrying capacity \(b\), and the exponent of fractional tumor cell kill by CD8+T cells \(l\). In general, it is recommended to perform the sensitivity analysis for the set of patient data before deciding the treatment strategy to get the optimal effect of the treatment and to avoid tumor recurrence.
In this paper, the mathematical model for the tumor–immune interaction has been fitted to an experimental data. The prediction of the model for tumor–immune interactions fits well with the experimental data. Treatment strategies for practical cases have been considered, and the most effective body parameters have been determined using sensitivity analysis. Hence, the model can be used with confidence in making predictions of response to therapy, assessment of treatment intervention strategies, and generally in perception of interaction between tumor cells and immune cells.
All relevant data sources are included in the article.
IFN-α:
Interferon-alpha
IL-2:
WHO. "WHO". 2021. [Online]. http://www.who.int/en/news-room/fact-sheets/detail/cancer. Accessed 12 Jan 2021.
Tao, Y., Guo, Q., Aihara, K.: A partial differential equation model and its reduction to an ordinary differential equation model for prostate tumor growth under intermittent hormone therapy. J. Math. Biol. 69(4), 817–838 (2014)
Article MathSciNet MATH Google Scholar
Malinzi, J.: Mathematical analysis of a mathematical model of chemovirotherapy: effect of drug infusion method. Comput Math Methods Med 2019, 7576591 (2019)
Cesano, A., Warren, S.: Bringing the next generation of immuno-oncology biomarkers to the clinic. Biomedicines 6(14), 14–24 (2018)
Arabameri, A., Asemani, D., Hajati, J.: Mathematical modeling of in-vivo tumor-immune interactions for the cancer immunotherapy using matured dendritic cells. J. Biol. Syst. 26(1), 1–22 (2018)
de Pillis LG, Radunskaya AE. Modeling tumor-immune dynamics. In: Mathematical models of tumor-immune system dynamics. New York; 2014.
Makhlouf, A.M., El-Shennawy, L., Elkaranshawy, H.A.: Mathematical modelling for the role of CD4+T cells in tumor-immune interactions. Comput Math Methods Med 2020, 1–16 (2020)
Elkaranshawy, H.A., Ezzat, H.M., Abouelseoud, Y., Ibrahim, N.N.: Innovative approximate analytical solution for standard model of viral dynamics: hepatitis C with direct-acting agents as an implemented case. Math Probl Eng 2019, 1454739 (2019)
Elkaranshawy HA, Ezzat HA, Ibrahim NN. Dynamical analysis of a multiscale model of hepatitis C virus infection using a transformed ODEs model. In: 42nd annual international conference of the IEEE engineering in medicine & biology society (EMBC), Montreal, QC, Canada; 2020.
Keya, K.N., Kamrujjaman, M., Islam, M.S.: The influence of density in population dynamics with strong and weak Allee effect. J Egypt Math Soc 29, 4 (2021)
Kamara, A.A., Mouanguissa, L.N., Barasa, G.O.: Mathematical modelling of the COVID-19 pandemic with demographic effects. J Egypt Math Soc 29, 8 (2021)
Deif, A.S., El-Naggar, S.A.: Modeling the COVID-19 spread, a case study of Egypt. J Egypt Math Soc 29, 13 (2021)
Zafar, Z.U.A., Ali, N., Baleanu, D.: Dynamics and numerical investigations of a fractional-order model of toxoplasmosis in the population of human and cats. Chaos Solitons Fract 151, 111261 (2021)
Zafar, Z.U.A., Ali, N., Younas, S., Abdelwahab, S.F., Nisar, K.S.: Numerical investigations of stochastic HIV/AIDS infection model. Alex Eng J 60(6), 5341–5363 (2021)
Zafar, Z.U.A., Rezazadeh, H., Inc, M., Nisar, K.S., Sulaiman, T.A., Yusuf, A.: Fractional order heroin epidemic dynamics. Alex Eng J 60(6), 5157–5165 (2021)
Elkaranshawy, H.A., Ezzat, H.M., Ibrahim, N.N.: Lyapunov function and global asymptotic stability for a new multiscale viral dynamics model incorporating the immune system response: Implemented upon HCV. PLoS ONE 16(10), e0257975 (2021)
Elkaranshawy, H.A., Ali, A.M.E., Abdelrazik, I.M.: An effective heterogeneous whole-heart mathematical model of cardiac induction system with heart rate variability. Proc Inst Mech Eng [H] 235(3), 323–335 (2021)
Jones, D.S., Plank, M.J., Sleeman, B.D.: Differential equations and mathematical biology, 2nd edn. CRC Press, Boca Raton (2009)
Book MATH Google Scholar
Tao, Y., Guo, Q., Aihara, K.: A mathematical model of prostate tumor growth under hormone therapy with mutation inhibitor. J Nonlinear Sci 20(2), 219–240 (2010)
Wei, X., Cui, S.: Existence and uniqueness of global solutions for a mathematical model of antiangiogenesis in tumor growth. Nonlinear Anal Real World Appl 9(5), 1827–1836 (2008)
Frieboes, H.B., Jin, F., Chuang, Y.L., Wise, S.M., Lowengrub, J.S., Cristini, V.: Three-dimensional multispecies nonlinear tumor growth-II : tumor invasion and angiogenesis. J Theor Biol 264(4), 1254–1278 (2010)
Lee, H.G., Kim, Y., Kim, J.: Mathematical model and its fast numerical method for the tumor growth. Math Biosci Eng 12(6), 1173–1187 (2015)
Zhang, Y., He, Y., Chen, H.: Boundary element method for a free boundary problem modeling three dimensional tumor growth. Comput Math Appl 73(7), 1627–1641 (2017)
Knopoff, D.A., Fernández, D.R., Torres, G.A., Turner, C.V.: adjoint method for a tumor growth PDE-constrained optimization problem. Comput Math Appl 66(6), 1104–1119 (2013)
Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A., Perelson, A.S.: Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bull Math Biol 56(2), 295–321 (1994)
Article MATH Google Scholar
F. A. Rihan, Delay Differential Equations and Applications to Biology, Springer, 2021.
Rihan, F.A., Lakshmanan, S., Maurer, H.: Optimal control of tumour-immune model with time-delay and immuno-chemotherapy. Appl Math Comput 353, 147–165 (2019)
MathSciNet MATH Google Scholar
Villasana, M., Radunskaya, A.: A delay differential equation model for tumor growth. J Math Biol 47(3), 270–294 (2003)
Xu, S., Bai, M.: Stability of solutions to a mathematical model for necrotic tumor growth with time delays in proliferation. J Math Anal Appl 421(1), 955–962 (2015)
Rihan, F.A., Abdel Rahman, D.H., Lakshmanan, S., Alkhajeh, A.S.: A time delay model of tumour-immune system interactions: Global dynamics, parameter estimation, sensitivity analysis. Appl Math Comput 232, 606–623 (2014)
Rihan, F.A., Rajivganthi, C.: Dynamics of tumor-immune system with random noise. Mathematics 9(21), 2707 (2021)
Nutini, A., Sohail, A.: Deep learning of the role of interleukin IL-17 and its action in promoting cancer. Bio-Algorithms Med-Syst 16(4), 20200052 (2020)
de Pillis, L.G., Fister, K.R., Gu, W., Head, T., Maples, K., Neal, T., Murgan, A., Kozai, K.: Optimal control of mixed immunotherapy and chemotherapy of tumors. J Biol Syst 16(1), 51–80 (2008)
de Pillis, L.G., Savage, H., Radunskaya, A.E.: Mathematical model of colorectal cancer with monoclonal antibody treatments. Br J Med Med Res 4(16), 3101–3131 (2014)
de Pillis, L.G., Gua, W., Radunskaya, A.E.: Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. J Theor Biol 238(4), 841–862 (2006)
Eftimie, R., Bramson, J.L., Earn, D.J.D.: Modeling anti-tumor Th1 and Th2 immunity in the rejection of melanoma. J Theor Biol 265(3), 467–480 (2010)
Anderson, L., Jang, S., Yu, J.L.: Qualitative behavior of systems of tumor-CD4+-cytokine interactions with treatments. Math Methods Appl Sci 38(17), 4330–4344 (2015)
Hu, X., Jang, S.R.J.: Dynamics of tumor-CD4+-cytokine-host cells interactions with treatments. Appl Math Comput 321, 700–720 (2018)
Elkaranshawy HA, Makhlouf AM, Abouelseoud Y. Using Padé approximant method to solve the mathematical model of tumor-immune interactions. In:2020 42nd annual international conference of the ieee engineering in medicine & biology society (EMBC), Montreal, QC, Canada; 2020.
Zi, Z.: Sensitivity analysis approaches applied to systems biology models. IET Syst Biol 5(6), 336–346 (2011)
Gul, R., Schütte, C., Bernhard, S.: Mathematical modeling and sensitivity analysis of arterial anastomosis in the arm. Appl Math Model 40(17–18), 7724–7738 (2016)
Kirch, J., Thomaseth, C., Jensch, A., Radde, N.E.: The effect of model rescaling and normalization on sensitivity analysis on an example of a MAPK pathway model. EPJ Nonlinear Biomed Phys 4, 3 (2016)
Bighamian, R., Parvinian, B., Scully, C.G., Kramer, G., Hahn, J.O.: Control-oriented physiological modeling of hemodynamic responses to blood volume perturbation. Control Eng Pract 73, 149–160 (2018)
Makhlouf, A.M., Elkaranshawy, H.A.: Sensitivity analysis for a mathematical model of tumor-immune interactions. UPB Sci Bull Ser A 83(2), 317–326 (2021)
MathSciNet Google Scholar
Wang, F., Idrees, M., Sohail, A.: "AI-MCMC"for the parametric analysis of the hormonal therapy of cancer. Chaos Solitons Fract 154, 111618 (2022)
Yu, Z., Arif, R., Fahmy, M.A., Sohail, A.: Self organizing maps for the parametric analysis of COVID-19 SEIRS delayed model. Chaos Solitons Fract 150, 111202 (2021)
Berraondo, P., Sanmamed, M.F., Ochoa, M.C., Etxeberria, I., Aznar, M.A., Pérez-Gracia, J.L., Rodríguez-Ruiz, M.E., Ponz-Sarvise, M., Castañón, E., Melero, I.: Cytokines in clinical cancer immunotherapy. Br J Cancer 120, 6–15 (2019)
Worschech, A., Chen, N., Yu, Y.A., Zhang, Q., Pos, Z., Weibel, S., Raab, V., Sabatino, M., Monaco, A., Liu, H., Monsurró, V., Buller, R.M., Stroncek, D.F., Wang, E., Szalay, A.A., Marincola, F.M.: Systemic treatment of xenografts with vaccinia virus GLV-1h68 reveals the immunologic facet of oncolytic therapy. BMC Genomics 10, 301 (2009)
Murphy, H., Jaafari, H., Dobrovolny, H.M.: Differences in predictions of ODE models of tumor growth: a cautionary example. BMC Cancer 16, 163 (2016)
Bates, D.M., Watts, D.G.: Nonlinear regression analysis and its applications. Wiley, New York (1988)
Gill, P.E., Murray, W., Wright, M.H.: Practical optimization, California. Emerald Group Publishing Limited, Boston (1982)
The authors would like to thank Dr. Lamiaa El-Shennawy for her valuable advice and suggestions.
Department of Engineering Mathematics and Physics, Faculty of Engineering, Alexandria University, Alexandria, Egypt
Hesham A. Elkaranshawy & Ahmed M. Makhlouf
Hesham A. Elkaranshawy
Ahmed M. Makhlouf
H. A. Elkaranshawy proposed the model and the methodology of this article. A. M. Makhlouf investigate the model and analyzed and interpreted the results. All authors read and approved the final manuscript.
Correspondence to Ahmed M. Makhlouf.
Elkaranshawy, H.A., Makhlouf, A.M. Parameter estimation and sensitivity analysis for a model of tumor–immune interaction in the presence of immunotherapy and chemotherapy. J Egypt Math Soc 30, 8 (2022). https://doi.org/10.1186/s42787-022-00143-0
Tumor–immune interactions