Mathematical modeling and dynamical analysis of anti-tumor drug dose-response

1.   Introduction

Cancer is a major public health concern globally. Using anti-tumor drugs is a way to control tumor growth. To achieve a better therapeutic effect, it is essential to choose the appropriate drug dose. A lower dose may not achieve the treatment effect but stimulate tumor growth and a higher dose is harmful to the patient's health. Therefore, it is necessary to study the quantitative relationship of drug dose-response. Hormesis is a biphasic dose-response phenomenon, which is generally defined as the negative effects of chemicals on organisms that occur at high doses (such as growth and development inhibition), but they are beneficial at low doses (such as stimulating growth and development). A large number of experiments have shown that the hormesis appeared in the dose-response relationship of various anti-tumor drugs ^[1]. Modeling this phenomenon will contribute to finding the appropriate dose ^[2].

The hormesis phenomenon in the dose-response relationship was first proposed by Edward J. Calabrese. In 2005, he and his collaborators evaluated the dose-response relationship between human tumor cell lines and various drugs including anti-tumor drugs. The biphasic dose-response similar to hormesis has been universally confirmed in 136 tumor cell lines with more than 30 tissue types and 120 different agents ^[1]. In 2009, Nascarella et al. proposed several observable indicators to quantitatively analyze whether hormesis occurs in the dose-response of a certain drug ^[3]. In 2014, Pearce et al. found that a particular class of tumor-directed immune reactants, anticancer antibodies, stimulated tumor growth at low doses and inhibited growth at higher doses ^[4]. Their experiments on multiple murine models (normal mice, mice with high immunocompetence, mice with poor immunogenicity, etc.) all showed the hormesis, indicating that the effect is widespread ^[5]. In 2015, Yoshimasu et al. proposed a formula to fit the dose-response relationship of anti-tumor drugs in vitro. Their theoretical model helps predict how hormesis affects patients with malignant tumors ^[6]. In 2019, Qian Li et al. proposed and expanded a discrete-time tumor growth model that combines radiotherapy and CTLA4 pathway inhibitors to describe the development of the tumor microenvironment ^[7]. They found that the intensity of radiotherapy had a significant impact on the kinetics of tumor growth. However, few studies have used dynamical models to study the hormesis in anti-tumor drug dose-response.

In this study, we combine experimental data with a dynamical model to quantitatively explore the hormesis of anti-tumor drugs. This study is structured as follows. Section 2 presents a mathematical model and a nondimensional model of the interactions among tumor cells, natural killer cells and M2-polarized tumor-associated macrophages. Section 3 provides the theoretical analysis. In Section 4, the data and parameter fitting are presented. Bifurcation analysis and quantitative simulations are conducted in Section 5. Finally, some conclusions and future work are addressed in Section 6.

2.   A mathematical model and its nondimensionalization

Nonhuman sialic acid N-glycolyl-neuraminic acid (Neu5Gc) is one of the three major mammalian sialic acids. It is found to have a significant role in human cancer immunity and tumorigenesis ^[8]. Studies have demonstrated that anti-Neu5Gc antibodies can stimulate or inhibit progression of tumors depending on the dosage used ^[5,9,10]. Experiments analyse the growth of Neu5Gc-positive tumors in Neu5Gc-deficient mice, following administration of increasing concentrations of anti-Neu5Gc antibodies ^[4]. They find the phenomenon that lower antibody concentrations stimulate tumor growth and higher concentrations inhibit growth (Figure 1).

In the experimental mice, the anti-Neu5Gc antibodies are able to affect tumor growth by interacting with natural killer (NK) cells ^[5]. When the Fab portion of the anti-Neu5Gc antibody binds to an antigen, its Fc portion is recognized and bound by Fc$\gamma$ receptors on NK cells, leading to activation of NK cells ^[11]. This antibody-dependent cellular cytotoxicity (ADCC) mediates NK cells to kill tumor cells. Besides, higher doses of antibodies promote NK cell-mediated ADCC. The experiments on mice show that when the antibody dose is higher, the number of NK cells is larger than the control group. When the antibody dose is lower, the number of NK cells is smaller than the control group. This indicates that antibody dose can affect NK cell growth and higher doses favor NK cell expansion.

This hormesis phenomenon has been reproduced in multiple mouse models^[5]. In this study, we construct a dynamic model of the interactions among tumor cells, NK cells and M2-polarized tumor-associated macrophages (TAMs) to describe the hormesis of anti-Neu5Gc antibodies. Based on both the knowledge of immune system function and experiments reported in references ^[4] and ^[5], we have the following model assumptions: 1) The tumor cells grow logistically in the absence of an immune response. 2) NK cells kill tumor cells when they interact via the law of mass action. TAMs stimulate tumor cells to proliferate when they interact via the law of mass action. 3) NK cells and TAMs may be stimulated to proliferate upon contact with tumor cells. 4) The apoptotic rates of NK cells and TAMs are independent of tumor cell density. 5) The effect of stimulation or inhibition could occur independently of any adaptive immunity. 6) Antibody dosage affects NK cell growth and higher doses favor NK cell expansion.

Using above assumptions, the graphical representation of the mathematical model, including three different cells and various interactions, is depicted in Figure 2.

Then, we establish the following ordinary differential equations (ODEs) containing three variables:

where $T$ represents tumor cell population, $N$ represents NK cell population and $M$ represents TAM population. $r$ is the intrinsic growth rate of the tumor, $K$ is the carrying capacity of the tumor, $a$ is the killing rate of NK cells on tumor cells and $b$ is the rate of proliferation of tumor when stimulated by TAMs. $k_1$ and $k_2$ are the rates of proliferation of NK cells and TAMs when stimulated by tumor, respectively. $c_1$ and $c_2$ are the apoptotic rates of NK cells and TAMs, respectively. The parameters $r, K, a, b, k_1, k_2, c_1$ and $c_2$ are constants and assumed to be non-negative. In our model, antibody dosage affects proliferation and apoptosis rates of NK cells (parameters $k_1$ and $c_1$).

For convenience of theoretical analysis, we nondimensionalize the model (2.1). Time is scaled relative to the apoptotic rate of NK cells ($c_1$), and we make the following substitutions:

where $N_0$ represents initial population of NK cells and $M_0$ represents initial population of TAMs. Then we obtain the nondimensional system of equations:

where

In this nondimensional system, antibody dose affects the parameter $D$. When the antibody dose is low, the parameter $D$ would be small and when the antibody dose is high, the parameter $D$ would be large.

3.   Theoretical analysis of the model system

3.1.   Positivity and boundedness

According to the comparison theorem ^[12], we can easily obtain the following theorem.

Theorem 3.1. Let the initial conditions of system (2.3) be positive, then the solutions $(x(\tau), {\rm{ y}}(\tau), {\rm{ z}}(\tau))$ of the system are nonnegative for all $\tau \ge 0$.

Theorem 3.2. System (2.3) has a nonnegative invariant set $\Omega = \{ (x(\tau), y(\tau), z(\tau)) \in {R^3}\left| 0 \le x(\tau) \le {\alpha}, 0 \le \right.y(\tau) \le {\beta}, 0 \le z(\tau) \le {\gamma}\}$, where ${\alpha} = \min\left\{\frac{F}{E}, \frac{1}{D}\right\}, {\beta} = \frac{A}{B}, \gamma = \frac{A}{C}\left(\alpha-1\right).$

Proof. Let $(x(\tau), y(\tau), z(\tau))$ be a solution of system (2.3) with $(x(0), y(0), z(0)) \in \Omega$. When $0 \le x(\tau) \le \alpha \le \frac{1}{D}$, from system (2.3) we obtain

which implies that $y$ is monotonic decreasing, so $0 \le y(\tau) \le y(0) \le {\beta}, {\rm{ }}\forall {\rm{ \mathsf{ τ}}} \ge {\rm{0}}{\rm{.}}$

Similarly, when $0 \le x(\tau) \le \alpha \le \frac{F}{E}$, from system (2.3) we obtain

which implies that $z$ is monotonic decreasing, so $0 \le z(\tau) \le {\gamma}, {\rm{ }}\forall {\rm{ \mathsf{ τ}}} \ge {\rm{0}}{\rm{.}}$

From system (2.3) we obtain

According to the comparison theorem ^[12], the lower bound of $x(\tau)$ is given by the lower bound of the solution of equation

Thus, we obtain $x(\tau) \ge 0, {\rm{ }}\forall {\rm{ \mathsf{ τ}}} \ge {\rm{0}}.$ The upper bound of $x(\tau)$ is given by the upper bound of the solution of equation

Thus, we obtain $x(\tau) \le \frac{A + C{\gamma}}{A} = \alpha, {\rm{ }}\forall {\rm{ \mathsf{ τ}}} \ge {\rm{0}}.$ In summary, $0 \le x(\tau) \le \alpha, {\rm{ }}\forall {\rm{ \mathsf{ τ}}} \ge {\rm{0}}.$

Therefore, the solutions of Eq (2.3) are bounded and the theorem is proven.

3.2.   Stability of equilibrium points

Considering the biological significance, the system (2.3) has four equilibrium points.

1) $E_0(0, 0, 0)$: Tumor cells, NK cells and TAMs are eliminated.

2) $E_1(1, 0, 0)$: The tumor cell population reaches the carrying capacity while NK cells and TAMs are eliminated.

3) ${E_2}(\frac{1}{D}, \frac{{A(D - 1)}}{{BD}}, 0)$: The TAMs are eliminated. Tumor cells and NK cells coexist.

4) When$\frac{F}{E} = \frac{1}{D}$, there exists a line of equilibrium points ${E_3}(\frac{1}{D}, {y^*}, \frac{{A - AD + BD{y^*}}}{{CD}}).$ For biological and dynamical relevance, we only take the subset of this line in the first octant. Thus, when $D > 1$, it requires that $y^* \geq \frac{A(D-1)}{BD}$; when $0 < D \leq 1$, we have that $y^* \geq 0$. The tumor cells, NK cells and TAMs coexist.

Remark: The system (2.3) has another equilibrium point $(\frac{F}{E}, 0, \frac{{A(F - E)}}{{CE}})$. However, considering the biological significance, tumor cells should not exceed the maximum carrying capacity of the environment. Then we have $F/E \le 1$. When $F = E$, this equilibrium point becomes $E_1(1, 0, 0)$. When $F < E, \frac{{A(F - E)}}{{CE}} < 0$, so this equilibrium point does not exist.

According to the stability theory of ODEs ^[13], we can easily obtain the Theorems 3.3 and 3.4.

Theorem 3.3. The tumor-free equilibrium point $E_0(0, 0, 0)$ is always unstable.

Theorem 3.4. When $D < 1, E < F$, the equilibrium point $E_1(1, 0, 0)$ is locally asymptotically stable.

Theorem 3.5. When $D > max\{1, \frac{E}{F}\}$, the equilibrium point ${E_2}(\frac{1}{D}, \frac{{A(D - 1)}}{{BD}}, 0)$ is locally asymptotically stable.

Proof. The Jacobian matrix $J({E_2})$ at equilibrium point $E_2$ is given by

and the corresponding characteristic equation is

where

According to Routh-Hurwitz criterion ^[14], the roots of Eq (3.7) have negative real parts equivalent to

Therefore, the equilibrium point $E_2$ is locally asymptotically stable if

When $D > max\{1, \frac{E}{F}\}$, the Eq (3.10) holds, so the equilibrium point $E_2$ is locally asymptotically stable.

Theorem 3.6. When $Q = (F-1)By^*-AF(1-\frac{1}{D}) < 0$, the equilibrium point ${E_3}(\frac{1}{D}, {y^*}, \frac{{A - AD + BD{y^*}}}{{CD}})$ is stable. When $Q > 0$, the equilibrium point ${E_3}$ is unstable. When $Q = 0, B = Ck, D = E = F = 1$, the equilibrium point ${E_3}$ is stable on the surface $z = ky^{E/D}, \forall k \ge 0$.

Proof. The Jacobian matrix $J({E_3})$ at equilibrium point $E_3$ is given by

and the corresponding characteristic equation is

When $Q > 0$, Eq (3.12) has a positive root, so $E_3$ is unstable.

When $Q < 0$, Eq (3.12) has one zero root and two negative roots, and we use the Center Manifold Theorem ^[13] to solve the stability of $E_3$. Let

To facilitate subsequent calculations, we make the following variable substitutions to move the equilibrium point to the origin:

Then we obtain the new model:

The Jacobian matrix $J_0$ at $(0, 0, 0)$ is given by

and the corresponding eigenvalues are $\lambda_1 = 0, \lambda_2 = \frac{1}{2}(-\frac{A}{D}-P), \lambda_3 = \frac{1}{2}(-\frac{A}{D}+P)$. Construct a matrix $\overline{M}$ whose columns are the eigenvectors of $J_0$:

Let $\overline{T} = {\overline{M}}^{-1}$, then

and

Make variable substitutions:

Then the system becomes:

The model (3.21) satisfies the conditions of the Center Manifold Theorem ^[13], and there must exist a central manifold

Let

Then $h(y_1)$ satisfies the following partial differential equation and boundary conditions

$h(y_1) = 0$ is a solution of the above equation. Substituting $h(y_1) = 0$ into the first equation of model (3.21), we obtain the equation satisfied by the solution on the central manifold is

Obviously the zero solution of this equation is stable. Thus the zero solution of model (3.21) is stable.

When $Q = 0$, the Center Manifold Theorem ^[13] used in previous proof is inapplicable. Instead, to analyze the stability of $E_3$, we reduce model (2.3) to a two-dimensional system (3.26) along the surface $z = ky^{E/D}$, for any $k \ge 0$, which is invariant under (2.3) when $F/E = 1/D$.

When $Q = 0, B = Ck, D = E = F = 1$, the Jacobian matrix $J({E_3})$ at equilibrium point $E_3$ is given by

and the corresponding characteristic equation is

The eigenvalue $\lambda_{1, 2} = \frac{1}{2}(-A\pm \sqrt{A^2-4y}) < 0$, so $E_3$ is stable.

In summary, when $Q < 0$, the equilibrium point ${E_3}$ is stable. When $Q > 0$, the equilibrium point ${E_3}$ is unstable. When $Q = 0, B = Ck, D = E = F = 1$, the equilibrium point ${E_3}$ is stable along the surface $z = ky^{E/D}, \forall k \ge 0$.

Corollary: When $D > 1, F \le 1$, the equilibrium point ${E_3}(\frac{1}{D}, {y^*}, \frac{{A - AD + BD{y^*}}}{{CD}})$ is stable. When $D < 1, F \ge 1$, the equilibrium point ${E_3}$ is unstable.

Proof. It's easy to verify that when $D > 1, F \le 1$, $Q = (F-1)By^*-AF(1-\frac{1}{D}) < 0$, the equilibrium point ${E_3}$ is stable. When $D < 1, F \ge 1$, $Q > 0$, the equilibrium point ${E_3}$ is unstable.

The above theoretical results are summarized in Table 1. Under the condition of $E < F$, the equilibrium point $E_1$ is stable when $D < 1$ and the equilibrium point $E_2$ is stable when $D > 1$. This indicates that parameter $D$ is the key parameter to change the steady state of the system. Increasing the value of parameter $D$ will change the steady state of the system from $E_1$ to $E_2$, which corresponds to the process of increasing the antibody dose. Therefore, we suppose that equilibrium point $E_1$ corresponds to the steady state under the low dose of antibody and equilibrium point $E_2$ corresponds to the steady state under the high dose of antibody.

4.   Data and parameter fitting

4.1.   Collection of data

We use the experimental data obtained from references ^[4] and ^[5]. These experiments analyzed the growth of Neu5Gc-positive tumors in Neu5Gc-deficient mice, following administration of increasing concentrations of anti-Neu5Gc antibodies.

The experimental data ^[5] includes the growth process of tumor in the first 10 days of certain mice (called LLC mice) under four different antibody doses (other experimental conditions remain the same), as shown in Figure 3. We also use the Siglec-E-deficient Siglece-null (SigE-/-) mice which have an enhanced innate immune response ^[5]. The conversion of tumor volume to the number of tumor cells is as follows: the volume of a single mammalian cell is about 100–10000 $\mu {m^3}$ ^[15]. According to reference ^[16], mean cell density in solid (tumor) tissue is $10^6/mm^3$. Therefore, we assume that the volume of a single tumor cell is 1000 $\mu {m^3}$ and there are about ${10^9}$ tumor cells in $1c{m^3}$ tumor.

4.2.   Parameter estimation

The model includes three variables and eight parameters. We obtain two parameters from the literature and identify the other six parameters by fitting experimental data. In this study, we use experimental data which is obtained under four different antibody doses. In our model, antibody dosage affects proliferation and apoptosis rates of NK cells (parameters $k_1$ and $c_1$). Therefore, the values of the parameters $c_1$ and $k_1$(corresponding to the parameter $D$) estimated by four sets of data are different.

To identify the parameters in the model, we convert the parameter identification into the problem of optimization by minimizing the following objective function. Mathematically, the objective function is defined as the error between the simulation results and the time series experimental data. The formulation can be expressed as

$y_i^D({t_j})$ represents the measured data of component $i$ at time-point $t_j$, which in our case, is the time series data of mouse experiments obtained from the literature. ${y_i}({t_j}, K)$ represents the $i$th component of the solutions of ODEs at time $t_j$ with parameter set $K$. The numerical solution of the ODEs is solved by MATLAB, with the initial value set to be the approximation of the first data point of the mouse data. The particle swarm optimization (PSO) ^[17] algorithm is used to optimize this object function. We first estimate six parameters of the model using data from the control group (without antibody). Then we fix the parameters $r, K, b$ and $k_2$ and only change $k_1$ and $c_1$ to fit the experimental data under three different drug doses. The obtained optimization parameters of the control group are listed in Table 2. The optimal values of parameters $c_1$ and $k_1$ of the LLC mice model when the dose is 14 ug, 28 ug and 56 ug are listed in Table 3 and the estimated values of the other four parameters are the same as in Table 2.

We use the estimated parameter values to simulate the growth process of the tumor in the first 12 days under different antibody doses and the results are depicted in Figure 3.

From Figure 3 we can observe that the numerical simulations of the model fit the experimental data well. When the antibody dose is 14 ug and 28 ug, the tumor volume in mice is always larger than the control group, indicating that these two lower antibody doses would promote tumor growth. When the antibody dose is 56 ug, the tumor volume in mice is always smaller than the control group, pointing that the higher antibody dose can inhibit tumor growth. Besides, 28 ug of antibody has a more significant effect on promoting tumor growth than 14 ug. The above simulation results are consistent with the hormesis phenomenon.

4.3.   Parameter uncertainty and identifiability

In our study, we use 24 data points under four different antibody doses to fit six parameters in our model, so the parameter uncertainty needs to be translated into confidence intervals for model predictions ^[20]. After estimating the optimal parameters, we exploit the profile likelihood to analyze the identifiability of the parameters ^[21] and use the software package ^[22] to estimate the parameters' confidence intervals. These finite sample confidence intervals are listed in Tables 2 and 3 and the confidence level is 95%.

The profile likelihood of six parameters in the LLC mice model are shown in Figure 4 and the profile likelihood of six parameters in the SigE-/- mice model are shown in Figure 5.

To be specific, the profile likelihood for $r$, $b$ and $c_1$ in LLC mice model and $r$, $k_2$ and $c_1$ in SigE-/- mice model show a steep concave shape, indicating that the optimization routines would reach their minimum rapidly. The profile likelihood for $K$, $k_1$ and $k_2$ in LLC mice model and $K$, $k_1$ and $b$ in SigE-/- mice model also show a concave shape, however, the curves on one side of the vertical dashed lines decrease slowly, indicating that their optimization routines would reach their minimum slowly.

4.4.   Sensitivity analysis

The global sensitivity of parameters reflects how the system responds to the perturbation of parameters in the model. To obtain the sensitivity of the input parameters to all variables in the model, the sensitivity function ${s_j}(t)$ of the parameter $P_j$ at time $t$ is defined as follows ^[23,24]:

where $O(t)$ is the model output at time $t$ and $\Delta {P_j}$ is a small perturbation which is $10\%$ in our situation. $\frac{S_j}{T_0} = \frac{\int_0^{T_0} {{s_j}(t)} dt}{T_0}$ is the sensitivity value of parameter $P_j$ ($T_0$ is the total time length of simulation time). We choose different $T_0$ and display the sensitivity analysis to the perturbation of parameters in the model in Figure 6.

From Figure 6, we can see that the variable $T$ (Tumor) in the model is most sensitive to the perturbation of $r$ (tumor intrinsic growth rate), and is sensitive to the perturbation of $a$ (NK cell killing rate to tumor) and $c_1$(NK cell apoptotic rate) too. This indicates that inhibiting tumor proliferation is effective for tumor control and changing the killing rate of NK cells or the NK cell apoptotic rate will also have a certain impact on tumor proliferation. The variable $N$ (NK cells) in the model is most sensitive to the perturbation of $c_1$. Besides, $c_1$ is the only parameter that has a significant impact on NK cells among the eight parameters of the model.

5.   Numerical results

5.1.   Simulating the influence of the key parameter D on the dynamics of model

The theoretical analysis presented in Section 3 indicates that parameter $D = \frac{Kk_1}{c_1}$ is a key parameter to change the steady state of the system. This suggests that the ratio of NK cell proliferation rate ($k_1$) to apoptotic rate ($c_1$) is a capital factor affecting the steady state of the system. In this section, we verify this result from the perspective of dynamics. We simulate the growth process of the tumor in first 100 days under four different antibody doses, as shown in Figure 7(a).

The cases where the antibody dose is 14 ug, 28 ug and control group correspond to $D < 1$ in the dimensionless model. The system will eventually reach the equilibrium point $E_1$. The tumor growth rate at the early stage is different under three antibody doses and the tumor cells proliferate fastest when the dose is 28 ug. The case of 56 ug belongs to the situation of $D > 1$ in the dimensionless model. Tumor growth will first go through an oscillating process and finally reach the steady state $E_2$.

We also simulate the dynamic of NK cells in first 100 days, as shown in Figure 7(b). When the antibody dose is low (14 ug, 28 ug), the number of NK cells is less than the control group and the promoting effect of M2-like TAMs is greater than the killing effect of NK cells, which generally promotes tumor growth. When the antibody dose is high (56 ug), the number of NK cells is more than the control group and the killing effect of NK cells is greater than the promoting effect of M2-like TAMs, which generally inhibits tumor growth. These simulation results are consistent with the hormesis phenomenon.

5.2.   Bifurcation analysis of parameters

To further understand the relationship between the steady state of model and parameter $D$, we use one-parameter bifurcation analysis and find that the stability of the equilibrium point changes as the parameter $D$ increases, as shown in Figure 8.

When $D < 1, E_1$ is stable; when $D > 1, E_1$ is unstable and $E_2$ is stable. By increasing the parameter $D$, the system undergoes a transition from monostability in a higher state to monostability in a lower state. In addition, we draw trajectories in a 3-dimensional phase-space with different initial values to represent that the value of the dimensionless parameter $D$ determines the dynamics of the system (Figure 9).

5.3.   Nonlinear relationship between antibody dose and model parameter D

After simulating the dynamics of tumor and immune cells, we would like to obtain a quantitative drug dose-response relationship curve, which reflects the hormesis phenomenon, to find the threshold point of drug dose. As mentioned above, the parameter $D$ is positively correlated with the ratio of NK cell proliferation to apoptosis and is a key parameter to change the steady state of the system. Therefore, we look for the quantitative relationship between $log (1/D)$ and drug (antibody) dose.

There are the following reasons for choosing $log (1/D)$ to represent the biological response. When $D > 1$, the system will eventually reach the steady state ${E_2}(\frac{1}{D}, \frac{{A(D - 1)}}{{BD}}, 0)$ and the steady state of tumor cell is $1/D$.

When $D < 1$, the system will ultimately reach the equilibrium point ${E_1}(1, 0, 0)$, which is not related to $D$. However, different antibody dosage will influence the early growth rate of the tumor. The parameter $D$ is positively correlated with the ratio of NK cell proliferation rate to apoptosis rate. The smaller $D$ or the slower $k_1$ will lead to the acceleration of tumor growth. Therefore, $log (1/D)$ is positively correlated with tumor growth rate.

In summary, using $log (1/D)$ to indicate the biological response can better reflect the changes in the tumor growth process caused by different antibody doses. According to reference ^[5], we infer that the quantitative relationship between the biological response $log(1/D)$ and the antibody dose (denoted as $\theta$) is as follows:

We use 14 ug, 28 ug and 56 ug experimental data for LLC mice to fit the quantitative relationship between $log(1/D)$ and drug (antibody) dose in LLC mice, as shown in Figure 10. The coefficients of the drug dose-response curve are listed in Table 4.

The blue dashed line intersects with the dose-response curve at the point (47.16, 1.55), suggesting that the antibody dose of 47.16 ug has no effect on tumor growth. This point is the threshold that we obtain from promoting tumor growth to inhibiting tumor growth. In clinical treatment, it will have the adverse effect of stimulating tumor growth if the antibody dose is less than this threshold. Therefore, the antibody dose used in clinical should be larger than this threshold to inhibit tumor growth.

The grey dashed line represents the case of $D = 1$. This line corresponds to the critical point where the steady state of the system changes from $E_1$ to $E_2$. The upper part of the line corresponds to the steady state $E_1$ and the lower part corresponds to the steady state $E_2$. The two dashed lines together divide the dose-response relationship curve into three parts (Figure 10).

1) In first part above the blue dashed line, the value of $log(1/D)$ is greater than that of the control group and the early growth rate of tumor is larger than that of the control group. The antibody dose in this interval stimulates tumor growth.

2) For second part between the two dashed lines, the value of $log(1/D)$ is less than the control group but larger than zero. The early growth rate of tumor is smaller than that of the control group, but the system steady state has not changed. The antibody dose in this interval can slow down the early growth of tumor, but still cannot inhibit tumor growth at last, and the tumor will reach the maximum carrying capacity of the environment.

3) In third part below the grey dashed line, the value of $log(1/D)$ is less than zero. The early growth rate of tumor is smaller than that of the control group and the system steady state changes from $E_1$ to $E_2$. The antibody dose in this interval can slow down the early proliferation of the tumor. After the oscillating process, it can finally achieve the effect of inhibiting tumor growth and the tumor will be eliminated eventually as the parameter $D$ gradually increases.

The antibody dose-response curve in Figure 10 well displays the characteristics of hormesis effect that low-dose antibodies stimulate tumor growth and high-dose antibodies inhibit tumor growth. It also clearly reflects the threshold point of the dose from promoting tumor growth to inhibiting tumor growth.

5.4.   The influence of the innate immune response on the dose-response

To investigate the effect of innate immune response on the drug dose-response, we further simulate the data of the Siglec-E-deficient Siglece-null (SigE-/-) mice. Siglec-E in the mouse is the functionally equivalent homolog of the inhibitory human Siglec-9, being expressed prominently on myeloid cells ^[5]. Siglecs are a family of sialic-acid-binding immunoglobulin-like lectins. They could regulate the growth and survival of cells in the innate and adaptive immune systems, either by inhibition of proliferation or induction of apoptosis ^[25]. In SigE-/- mice, the deficiency of Siglecs affects the regulation of innate immune system cells' growth and leads to an enhanced innate immune response. We identify the parameters in model (2.1) using the data of SigE-/- mice. The obtained optimization parameters of the control group are listed in Table 5. The optimal values of parameters $c_1$ and $k_1$ of the SigE-/- mice model when the dose is 11 ug and 28 ug are listed in Table 6 and the estimated values of the other four parameters are the same as in Table 5.

We use these parameters to simulate the growth process of the tumor in the first 16 days under different antibody doses and the results are depicted in Figure 11.

As shown in Figure 11, when the antibody dose is 11 ug, the tumor volume in mice is larger than that of the control group. The number of NK cells is less than that of the control group and the promoting effect of M2-like TAMs is greater than the killing effect of NK cells, which generally promotes tumor growth. When the dose is 28 ug, the tumor volume in the mice is smaller than that of the control group. The number of NK cells is more than that of the control group, and the killing effect of NK cells is dominant, which generally inhibits tumor growth. The above results show that the hormesis also appears in mice with enhanced innate immune response.

We fit both the quantitative relationship between $log(1/D)$ and antibody dose in LLC mice and SigE-/- mice, as shown in Figure 12. To fit the dose-response curve, we use 14 ug, 28 ug and 56 ug data for LLC mice and use control group, 11 ug and 28 ug data for SigE-/- mice. The dose-response curve of mice with an enhanced innate immune response (SigE-/-) is on the left of the dose-response curve of ordinary mice (LLC), indicating that enhancing the innate immune response will shift the dose-response curve to the left.

Additionally, the antibody dose of 28 ug is used in two mouse experiments. In LLC mice, the value of $log(1/D)$ corresponding to 28 ug is greater than that of the control group, indicating that 28 ug antibody can stimulate tumor growth. In the experiment of mice with enhanced innate immune response, the value of $log(1/D)$ corresponding to 28 ug is smaller than that of the control group, showing that 28 ug antibody can inhibit tumor growth. So with enhanced innate immune response, a lower dose of antibody is required to achieve the same stimulation and inhibition effect.

5.5.   The great effect of intrinsic growth rate on tumor growth process

In this section, we investigate the influence of the intrinsic growth rate of tumor ($r$) on tumor growth. We only change the parameter $r$ and keep other parameters and initial values of the model (2.1) unchanged.

First, when the drug dose is 14 ug and 28 ug (the corresponding value of parameter $D$ at the two doses can be obtained using (5.1)), we change the parameter $r$ to obtain the dynamics of tumor cells in four cases, as shown in Figure 13(a), (b). We find that changing the intrinsic growth rate of tumor cells will affect the proliferation of tumor. The greater intrinsic growth rate of tumor cells will lead to the quicker proliferation of tumor in the early stage, and the number of tumor cells will reach a steady state faster.

Second, with the drug dose of 56 ug, the parameter $r$ is changed to obtain the dynamics of tumor cells in three cases, as shown in Figure 13(c). We find that reducing the inherent growth rate of tumor cells can reduce the peak value of tumor cell number during the growth process. Meanwhile, reducing the inherent growth rate of tumor cells will prolong the oscillation process of tumor cell growth, which causes the tumor to take a longer time to reach a steady state.

Besides, changing the intrinsic growth rate of tumor cells will not change the steady state that the number of tumor cells finally reaches. This result is consistent with the stability analysis, that is, changing the parameter $r$ will not change the steady state of the system. This indicates that only reducing the intrinsic growth rate of tumor cells cannot effectively inhibit tumor and other treatments are needed, such as increasing the dosage of drugs, enhancing the innate immune response and so on.

6.   Conclusion and discussion

In this study, we construct a dynamic model to study the hormesis in the anti-tumor drugs dose-response and investigate the mechanism of tumor inhibition or promotion by the immune system. We fit out the quantitative dose-response curve of drugs (antibodies), which can well reflect the hormesis, and obtain the critical threshold point of antibody dose changing from promoting tumor growth to inhibiting tumor growth. In clinical applications, if the antibody dose is lower than this threshold, it will have the adverse effect of stimulating tumor growth. Therefore, to inhibit tumor growth, the antibody dose should be larger than this threshold. Comparing the experiments of mice with enhanced innate immune response and ordinary mice, we find that improving the innate immune response can shift the dose-response curve to the left. Besides, we find that changing the intrinsic growth rate of tumor will affect the proliferation of tumor and prolong the oscillation process of tumor growth. However, only reducing the intrinsic growth rate of tumor cannot effectively inhibit tumor growth and other treatments are needed, such as increasing the dosage of drugs, enhancing the innate immune response and so on.

In this work, we only model tumor cells, NK cells and M2-polarized macrophages, and do not consider other immune cells, cytokines and other factors. Therefore, other immune cells and cytokines can be added to the existing model to further study the dynamics of their interactions. In addition, this study only considers the impact of enhancing innate immune response on the dose-response curve of anti-tumor drugs, and the impact of other immunotherapy can be considered in subsequent studies.

Acknowledgments

This work is supported by the Key Program of the National Nature Science Foundation of China (No.11831015), the Chinese National Nature Science Foundation (No.61672388) and the National Key Research and Development Program of China (No.2018YFC1314600).

Conflict of interest

All authors declare no conflicts of interest in this paper.