Parameter estimation and sensitivity analysis for a model of tumor–immune interaction in the presence of immunotherapy and chemotherapy

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.

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 timedelayed 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

Method
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: where T (t) is the tumor cells population, N (t) is the natural killer cells population, L(t) is the CD8 + T cells population, Y (t) is the CD4 + T cells population, C(t) is the circulating lymphocytes cell population not including natural killer cells and active CD8 + T and CD4 + T cells and I(t) is the concentration of the cytokine. D is the fractional kill rate, which defined as follows: In Eq. (1), aT (1 − bT ) represents the tumor growth. The term −cNT represents the natural killer-induced tumor death. The term −DT represents the tumor lysis by CD8 + T cells. The term − 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 −fN represents the turnover of natural killer cells. The term gT 2 h+T 2 N represents the recruitment of natural killer cells induced due to the existence of the tumor. The term −pNT represents the death of natural killer cells by the exhaustion of tumor-killing resources. In Eq. (3), the term −mL represents the inactivation of CD8 + T cells. The term jD 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 −uNL 2 represents the regulation and suppression of CD8 + T cell activity by the natural killer cells. The term p i I g i +I L represents the activating effect of cytokines on CD8 + T cells. In Eq. (4), the term β 1 T α 1 +T represents the activating effect of cytokines on CD4 + T cells due to the existence of the tumor. The term −µ 1 Y represents the decaying of CD4 + T cells which has a constant rate. The term −δ 2 TY 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 α and decayed with a constant natural death rate β . In Eq. (6), the term, −µ i I represents the decaying of cytokine with a constant rate µ i . The term β 2 T α 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).
To include the effect of immunotherapy and chemotherapy, the model is modified to be: , v L (t) and v Y (t) . 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-α ) 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-α [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 (0) = 10 4 cells, L(0) = 10 2 cells, Y (0) = 10 5 cells and C(0) = 10 8 cells. The tumor size was measured in units of mm 3 in the data set of Worschech et al. [48]. Therefore, any cell unit is converted to 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 cell −1 . Hence, the estimated values for the two parameters are a = 2.68208 × 10 −2 day −1 and b = 1.05927 × 10 −10 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 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.

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 (0) = 2 × 10 5 , initial natural killer cells N (0) = 10 4 , initial CD8 + T cells L(0) = 10 2 , initial CD4 + T cells Y (0) = 10 5 , and initial circulating lymphocytes C(0) = 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 (0) = 2 × 10 7 , the tumor is not eradicated, and treatment is needed as shown in Fig. 2b.
In the third case, the model is solved for patient's data with the experimental value of initial tumor size T (0) = 2.41758 × 10 8 . Other initial conditions are: N (0) = 10 4 , L(0) = 10 2 , Y (0) = 10 5 , C(0) = 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 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.
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.

Sensitivity analysis
The sensitivity coefficient S ij for the system output y i w.r.t. the parameter j is given by [41][42][43][44]: 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 , α, α 1 , α 2 , β, β 1 , β 2 , µ 1 , µ i and δ 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]:   (1)-(6). The system of differential equations defined by (16) can be solved using a suitable numerical method with the initial conditions S ij (0) = 0 . After that, the solution will be normalized to bound the range of the sensitivity functions using the following formula: Sensitivity analysis is performed for the patient's data with the value of initial tumor size given in the experiment [48], T (0) = 2.41758 × 10 8 , the initial natural killer cells N (0) = 10 4 , initial CD8 + T cells L(0) = 10 2 , initial CD4 + T cells Y (0) = 10 5 and initial circulating lymphocytes C(0) = 10 8 . The normalized sensitivity function is plotted in Fig. 5. It shows that the most effective parameters are a , b and l , respectively.

Discussion
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.
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 (17)  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.

Conclusions
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.