Modeling free tumor growth: Discrete, continuum, and hybrid approaches to interpreting cancer development

1.   Introduction

Cancer, as the second-most common fatal disease, remains a serious threat to human health, and millions of people suffer from cancer around the world. Cancer cells proliferate and spread rapidly, which causes many difficulties in treatment. The function of the immune system is to prevent tumor cells from invading the body, and the body restricts the growth of tumor cells. In view of molecular cell biology, the immune system has a natural defense for the body, since it can recognize and destroy tumor cells [1,2]. It is obvious that immune system treatment is effective for the treatment of cancer and has more advantages than traditional treatment methods, such as chemotherapy, operation, and so on. Just owing to the character of the immune system, immunotherapy arises gradually, and immunotherapy has become the primary treatment instead of other methods in recent decades. Due to the complexity associated with the function of the immune system in the process of tumor proliferation, pathologists and clinicians think that treating malignant tumors is still very hard, such that human beings cannot overcome it even today.

It is necessary to construct a mathematical model for understanding the mechanism of the immune system, numerous linear and nonlinear mathematical models endowed with functional response have been considered by mathematicians, which play a crucial role in predicting or controlling tumor growth. It is also necessary to study the mechanisms underlying the interaction between the tumor-immune system to understand immunotherapy. Further development of immunotherapy depends on a better comprehension of the interaction between the tumor system and immune system. Tumor cells can escape the monitoring of the immune system to proliferate infinitely, which is called immune suppression. Many researchers are devoted to removing immune suppression so that cancer cells can be controlled by the immune system. Mathematical models, for example, ordinary differential equation (ODE), partial differential equation (PDE) and partial functional differential equation (PFDE), are extensively used to investigate the dynamics of the immune system. Many researches study how the immune system influences tumor cells by modeling and analyzing mathematical models. Experts attempt to solve the problem based on theoretical analysis and experimental data. There is no doubt that constructing mathematical models provides us with a new view to recognize the tumor-immune system.

In recent years, immunotherapy models have been extensively utilized to predict the behaviors of tumor cells. Ref. [3] constructed a tumor-immune system with two distributed delays to describe the interactions between tumor cells and immune cells. Ref. [4] proposed a tumor-immune model with time delay, and explained that the time delay of the immune reaction could change the stability of equilibrium. Ref. [5] analyzed a two-dimensional tumor-immune system with two delays. By using bifurcation theory, the stability of equilibrium and the existence of Hopf bifurcation were given when the time delays were regarded as the bifurcation parameters. Relevant scholars proposed all kinds of differential equations [6-8]. Refs. [9,10] showed that different sizes of delays could cause a change in equilibrium stability. Refs. [11-13] studied the various dynamic properties associated with the changing time delay between tumor cells and immune cells. Ref. [14] studied a tumor-immune system in view of optimal control, and showed the method associated with minimizing the cost of immune-chemotherapy and reducing the load of tumor cells. Ref. [15] investigated a tumour-immune system with delays and fractional-order, and provided the necessary and sufficient conditions for stability of the steady states and Hopf bifurcating periodic solutions. Ref. [16] presented a delay differential equation to analyze the dynamics between effector cells and tumor cells. The existence and stability of possible steady states and the local stability of Hopf bifurcating periodic solutions are discussed. Das et al. studied a series of studies associated with tumor-immune interaction. For example, Ref. [17] studied the deterministic and stochastic dynamics of tumor-immune interactions. Refs. [18,19] considered the dynamics of time-delayed tumor-immune systems. Refs. [20,21] compared the deterministic and noise-induced chaotic dynamics, and described the growth and proliferation process of tumor cells. The analysis showed that the onset of chaos in the system can be predicted. Refs.[22,23] established the optimal treatment strategies that maximize the number of immune-effector cells, minimized the number of cancer cells, and detrimental effects caused by the amount of drugs.

The motivation of this paper is as follows. First, Ref. [24] showed that there existed a time delay during the immune system recognition of tumor cells. Ref. [25] pointed out that delay played an important role in the interaction between the immune system and tumor cells. The different delays could cause different phenomena of the tumor-immune system, which is useful for understanding the tumor-immune system and providing treatment for tumors. Second, Ref. [26] indicated that tumor cells could spread in the body and revealed some factors causing the diffusion of tumor cells. Therefore, it is necessary and reasonable to introduce delay and diffusion into the tumor-immune system. Third, Ref. [27] considered some results on the nonlinear dynamics of delayed differential equation models describing the interaction between tumor cells and effector cells of the immune system, however, it did not consider the influence of diffusion of tumor cells and immune cells. For the immune reaction of organisms, it must be enough time to identify nonself cells, and the growth of tumor cells also needs some time. Based on the model of Ref. [27], we introduce two delays and diffusion, and propose a tumor-immune system with two delays and diffusion.

The paper is organized as follows. In Section 2, we establish a partial functional differential equation (PFDE) associated with tumor-immune system, and study the dynamics of this system. In Section 3, the stability of nonnegative equilibrium and the existence of Hopf bifurcation are presented. The normal form of Hopf bifurcation is deduced in Section 4. Numerical simulations to test our findings are shown in Section 5. In Section 6, we give a conclusion for our work.

2.   Mathematical modeling

The essential mechanisms of interaction between tumor cells and immune cells is shown as follows (see Figure 1).

Immune cells assault tumor cells, and the multiplication of immune cells is stimulated by the appearance of tumor cells. The body can not recognize tumor cells immediately, and there exists a time delay $\tau_1$ for responding to the entrance of nonself cells, and the response time may be short, however, it cannot be ignored. Meanwhile, there also exists a time delay $\tau_2$ when tumor cells are recognized by immune cells. Furthermore, tumor cells and immune cells diffuse in the body under conditions of limited resources and space, and the growth and interaction between tumor cells and immune cells not only depend on time, but are also affected and restricted by the living spatial environment. Hence, it is also indispensable to investigate the complex interaction caused by spatial factors. Considering the above factors, we construct the following model to analyze the dynamics of the tumor-immune system:

where $\Omega = [0, \pi]$ is a bounded open domain in $\mathrm{R}^N\; (N\geq1)$ with a smooth boundary $\partial\Omega$, and $\Delta$ denotes the Laplacian operator in $\mathrm{R}^N\; (N\geq1)$. $u$ and $v$ are the numbers of tumor cells and immune cells, respectively. The description of the parameters is given in Table 1, and these parameters are all positive. We also assume that $u_0, v_0\in \mathfrak{C} = C([-\tau, 0], X)$ and $X$ is defined by $X = \{u, v\in W^{2, 2}(\Omega):\frac{\partial u(x, t)}{\partial x} = \frac{\partial v(x, t)}{\partial x} = 0, \; x\in\partial \Omega\}$. There exists a time delay for producing a proper immune response in host cells, and considering the time delay of tumor growth $\tau_1$, the growth of tumor cells follows $ru(1-\frac{u(x, t-\tau_1)}{\alpha})$. $\beta u(x, t-\tau_2) v(x, t-\tau_2)$ means that immune cells are activated due to the stimulus of tumor cells.

3.   Bifurcation Analysis

System (2.1) has two boundary equilibria $E_0 = (0, 0), E_1 = (\alpha, 0)$, and it has a unique positive constant steady equilibrium $E^* = (u^*, v^*) = (\frac{d}{\beta-p}, \frac{r}{m} (1-\frac{d}{\alpha (\beta-p)}))$ when $d < \alpha (\beta-p)$ and $\beta > p$. For convenience, we show this condition associated with the existence of positive equilibrium $E^*$ by using the following assumption.

Denote $U(x, t)$ = $(u(x, t), v(x, t))^\mathrm{T}$, then linearized system for Eq. (2.1) can be rewritten as a differential equation at equilibrium $E = (u_0, v_0)$ (here $(u_0, v_0) = (0, 0), \; (\alpha, 0), \; (u^*, v^*))$:

where

Hence, the characteristic equation of (2.1) at $(u_0, v_0)$ is given as follows:

where

3.1.   The case for $\tau_{1}=0$, $\tau_{2}=0$

When $\tau_{1} = 0$, $\tau_{2} = 0$, Eq. (3.1) for equilibrium $E_0 = (0, 0)$ becomes

For Eq. (3.1) with $n = 0$, the product of two eigenvalues $-rd < 0$, obviously, equilibrium $E_0 = (0, 0)$ is always unstable.

When $\tau_{1} = 0$, $\tau_{2} = 0$, Eq. (3.1) for equilibrium $E_1 = (\alpha, 0)$ becomes

When $(H_1)$ does not hold,

thus, equilibrium $E_1 = (\alpha, 0)$ is always locally asymptotically stable when $(H_1)$ does not hold.

When $(H_1)$ holds, there also exists positive equilibrium $E_2$. The stability analysis for $E_1 = (\alpha, 0)$ is similar to $E^* = (u^*, v^*)$, and we only show the general stability results for $E = (u_0, v_0)$ (here $(u_0, v_0) = (\alpha, 0), (u^*, v^*)$). Thus, we show following assumption.

When $\tau_{1} = 0$, $\tau_{2} = 0$, Eq. (3.1) becomes

note that under $(H_2)$, $A_{n} +B_{n}+ D_{n}\geq A_{0} +B_{0}+ D_{0} = d+pu_0 +\frac{2r u_0}{\alpha}-r+mv_0-\beta u_0 > 0$. Denote

where

Thus, we obtain following theorem.

Theorem 3.1. For system (2.1) with $\tau_{1} = 0$, $\tau_{2} = 0$, the stability results for equilibria are given as follows.

(1) Equilibrium $E_0 = (0, 0)$ is always unstable.

(2) Equilibrium $E_1 = (\alpha, 0)$ is always locally asymptotically stable when $(H_1)$ does not hold.

(3) When $(H_1)$ and $(H_2)$ hold, the stability results for equilibrium $E_1$ (or $E^*$) are given as follows:

(3-1) If $\xi_0 < 0$, equilibrium $E_1$ (or $E^*$) is unstable.

(3-2) If $\xi_0 > 0$ and $-\frac{\xi_1}{2\xi_2} < 0$, equilibrium $E_1$ (or $E^*$) is locally asymptotically stable.

(3-3) If $\xi_0 > 0$, $-\frac{\xi_1}{2\xi_2} > 0$, $\mathfrak{F}(\lfloor-\frac{\xi_1}{2\xi_2}\rfloor) > 0$ and $\mathfrak{F}(\lfloor-\frac{\xi_1}{2\xi_2}\rfloor+1) > 0$, equilibrium $E_1$ (or $E^*$) is locally asymptotically stable.

(3-4) If $\xi_0 > 0$, $-\frac{\xi_1}{2\xi_2} > 0$, and $\mathfrak{F}(\lfloor-\frac{\xi_1}{2\xi_2}\rfloor) < 0$ or $\mathfrak{F}(\lfloor-\frac{\xi_1}{2\xi_2}\rfloor+1) < 0$, equilibrium $E_1$ (or $E^*$) is unstable.

3.2.   The case for $\tau_{2}=0$, $\tau_{1}\neq0$

The stability analysis for $E_1 = (\alpha, 0)$ is similar to $E^* = (u^*, v^*)$, thus, we only show the general stability results for $E = (u_0, v_0)$ (here $(u_0, v_0) = (\alpha, 0), (u^*, v^*)$).

When $\tau_{2} = 0$, $\tau_{1}\neq0$, the characteristic equation of system (2.1) is given as follows:

where $L_{1} = A_n+D_n$, $L_{2} = E_n+G_n$, $L_{3} = B_n, \; L_{4} = C_n+F_n$, with

We might suppose that $\pm\mathrm{i} \omega\; (\omega > 0)$ are a pair of purely imaginary roots of Eq. (3.5). Substituting them into Eq. (3.5) and separating the real and imaginary parts, we obtain

thus, $R: = \mathrm{sin}(\omega \tau_1) = \frac{(\omega^{2}-L_{2}) L_{3} \omega+L_{1} L_{4} \omega}{L_{3}^{2} \omega^{2}+L_{4}^{2}}$, $Q: = \mathrm{cos}(\omega \tau_1) = \frac{(\omega^{2}-L_{2}) L_{4}- L_{1} L_{3}\omega^{2}}{L_{3}^{2} \omega^{2}+L_{4}^{2}},$ which implies that

Let $z = \omega^{2},$ then the above equation can be rewritten in the following form

Under assumption

Eq. (3.6) has no positive roots.

Under assumption

Eq. (3.6) has one unique positive root $\omega_{0}$.

Under assumption

Eq. (3.6) has two positive roots $\omega_{1, 2}$, and the critical time delay can be determined as,

where $R: = \mathrm{sin}(\omega \tau_1) = \frac{(\omega^{2}-L_{2}) L_{3} \omega+L_{1} L_{4} \omega}{L_{3}^{2} \omega^{2}+L_{4}^{2}}$, $Q: = \mathrm{cos}(\omega \tau_1) = \frac{(\omega^{2}-L_{2}) L_{4}- L_{1} L_{3}\omega^{2}}{L_{3}^{2} \omega^{2}+L_{4}^{2}},$ $l = 0, 1, 2;\; j = 0, 1, 2, \cdots.$ Then the following transversality conditions yield,

where $\Delta = L_2^2 z+L_3^2 z^2, \; l = 0, 1, 2;\; j = 0, 1, 2, \cdots.$ Obviously, under $(H_1)\; and\; (H_4)$, $h(z) = 0$ has one positive real root $z_0$ and $h^{'}(z_0) > 0;$ under $(H_1)\; and\; (H_5)$, $h(z) = 0$ has two positive real roots $z_1$ and $z_2$. Suppose $z_1 < z_2$, then $h^{'}(z_1) < 0$ and $h^{'}(z_2) > 0$.

Theorem 3.2. For system (2.1) with $\tau_{2} = 0, \; \tau_1 > 0$, the following conclusions hold when assumptions $(H_1), (H_2)$ and the stability conditions of Theorem 3.1 (3-2) or (3-3) hold,

(1) If $(H_3)$ is satisfied, then equilibrium $E_1$ (or $E^*$) is locally asymptotically stable for all $\tau_1 > 0.$

(2) If $(H_4)$ is satisfied, then equilibrium equilibrium $E_1$ (or $E^*$) is locally asymptotically stable for all $0 < \tau_1 < \tau_{10}^{(0)},$ and unstable for $\tau_1 > \tau_{10}^{(0)},$ furthermore, system (2.1) undergoes a Hopf bifurcation at equilibrium $E_1$ (or $E^*$) when $\tau_1 = \tau_{10}^{(j)}, j = 0, 1, 2, \cdots.$

(3) If $(H_5)$ is satisfied, then there is a positive integer m such that the equilibrium $E_1$ (or $E^*$) is locally asymptotically stable when $\tau_1 \in[0, \tau_{12}^{(0)})\bigcup \cup_{k = 0}^{(m-1)} (\tau_{11}^{(k)}, \tau_{12}^{(k+1)})$, and unstable $\tau_1\in \cup_{k = 0}^{(m-1)}(\tau_{12}^{(k)}, \tau_{11}^{(k)}) \bigcup(\tau_{12}^{(m)}, +\infty),$ furthermore, system (2.1) undergoes a Hopf bifurcation at equilibrium $E_1$ (or $E^*$) when $\tau_1 = \tau_{1l}^{(j)}, \; l = 1, 2\; ; j = 0, 1, 2, \cdots$.

3.3.   The case for $\tau_{2}\neq 0, \tau_{1}\neq 0$

The first stable interval mentioned in Theorem 3.2, on which the equilibrium $E_1$ (or $E^*$) is locally asymptotically stable, is denoted by stable region $I$. That is, assumptions $(H_1), (H_2)$ and the stability conditions of Theorem 3.1 (3-2) or (3-3) hold, and if $(H_3)$ also holds, $I = \{\tau_1|\tau_1 > 0\}$. If $(H_4)$ also holds, $I = \{\tau_1|0 < \tau_1 < \tau_{10}^{(0)}\}$. If $(H_5)$ holds, $I = \{\tau_1|0 < \tau_1 < \tau_{12}^{(0)}\}$. Next, we choose $\tau_{1} = \tau_{1}^{*} \in I.$ Regarding $\tau_{2}$ as a parameter, the characteristic equation of (2.1) at $(u_0, v_0)$ (here $(u_0, v_0) = (\alpha, 0), (u^*, v^*)$) is rewritten as follows:

where

Letting $\lambda = \mathrm{i} \omega (\tau_{2}) (\omega > 0)$ be the root of the above equation, then separate the real and imaginary parts for above equation,

which leads to,

Suppose

then, we get $F(0) = G_n^{2}+C_n^{2}-E_n^{2}-F_n^{2}+2C_n G_n-2E_nF_n < 0, F(+\infty) > 0.$ Hence, $F(\omega) = 0$ has definite positive roots $\omega_{2i}, \; i = 0, 1, 2.$ For every fixed $\omega_{i}$, there is a sequence of $\tau_{2i}^{(j)}, j = 0, 1, 2, \cdots,$ defined by

where $j = 0, 1, 2, \cdots,$

with

Let $\tau_{20} = \min \tau_{2i}^{(j)}, i = 0, 1, 2, \; j = 0, 1, 2, \cdots$, when $\tau_{2} = \tau_{20}$, Eq. (3.1) has a pair of purely imaginary roots $\pm \mathrm{i} \omega_{20}$ for $\tau_{1}^{*} \in I$. Assume

then we have,

Theorem 3.3. Assuming that $(H_1), \; (H_6)\; and \; (H_7)$ are satisfied, and one of the assumptions $(H_3), (H_4)$ and $(H_5)$ holds, then for $\tau_{1}^{*} \in I,$system (2.1) undergoes a Hopf bifurcation at $E = (u_0, v_0)$ when $\tau_{2} = \tau_{20}^{(j)}$, here $(u_0, v_0) = (\alpha, 0)$ or $(u^*, v^*)$ with $u^* = \frac{d}{\beta-p}, \; v^* = \frac{r}{m} (1-\frac{d}{\alpha (\beta-p)})$, and the equilibrium $E$ is locally asymptotically stable when $\tau_{2}\in[0, \tau_{20})$.

4.   Normal formal of Hopf bifurcation

In this section, we derive the normal formal of Hopf bifurcation at $E^{*} = (u^*, v^*)$ when $\tau_{2} = 0, \; \tau_{1} \neq 0$ for the system (2.1). We denote the critical value $\tau = \tau_2 = \tau_{c}$. Define $\overline{u}(x, t) = u(x, \tau t)-u^{*}, \; \overline{v}(x, t) = v(x, \tau t)-v^{*},$ and drop the bar for convenience, the system (2.1) can be rewritten as,

Let $h = (h_{11}, \; h_{12})^{T}$ be the eigenvector of the linear operator corresponding to the eigenvalue $\mathrm{i} \omega \tau$, and let $h^{*} = (h_{21}, \; h_{22})^{T}$ be the normalized eigenvector of the adjoint operator of the linear operator corresponding to the eigenvalues $-\mathrm{i} \omega \tau$ satisfying the inner product $< \overline{h^{*}}, h > = \overline{h^{*}}^{T}\cdot h = 1$. By a simple calculation, we get

Suppose the solution of Eq. (4.1) is

where

The derivation with respect to $t$ is

where the differential operator $D_{i} = \frac{\partial}{\partial T_{i}}, \; i = 0, \; 1, \; 2, \cdots$.

Denote

We obtain

We take perturbations as $\tau = \tau_{c}+\varepsilon \mu$, to deal with the delayed terms, we expend $u(x, t-1), \; v(x, t-1)$ at $u(x, T_{0}-1, T_{1}, \; T_{2}, \cdots)$ and $v(x, T_{0}-1, T_{1}, \; T_{2}, \cdots)$, respectively, that is,

where $u_{j, 1} = u_{j}(x, T_{0}-1, T_{1}, \; T_{2}, \cdots), \; v_{j, 1} = v_{j}(x, T_{0}-1, T_{1}, \; T_{2}, \cdots), \; j = 1, 2, 3, \cdots.$

Substituting Eqs. (4.3)-(4.5) into Eq. (4.1), for the $\varepsilon$-order terms, we obtain

Since $\pm \mathrm{i} \omega \tau$ are the eigenvalues of the linear part of Eq. (4.1), the solution of Eq. (4.6) can be expressed in the following form

where c.c. means the complex conjugate of the preceding terms, and $h$ is given in Eq. (4.2).

For the $\varepsilon^{2}$-order terms, we obtain

Substituting Eq. (4.7) into the right side of Eq. (4.8), we obtain the coefficient vector of term $\mathrm{e}^{\mathrm{i}\omega \tau_{c} T_{0}},$ denoted as $m_{1}$, by solvability conditions, that is, $< h^{*}, \; m_{1} > = 0,$ we obtain

where

Suppose the solution of Eq. (4.8) is

Denoted $c_k = < \cos(n x) \cos(n x), \cos(k x) > = \int_{0}^{\pi} \cos(n x) \cos(n x) \cos (kx) \mathrm{d}x$, thus

with

For the $\varepsilon^{3}$-order terms, we obtain

Substituting solutions Eq. (4.7) and Eq. (4.10) into the right side of Eq. (4.12), we obtain the coefficient vector of term $\mathrm{e}^{\mathrm{i} \omega \tau_{c} T_{0}},$ denoted as $m_{2}$, by solvability conditions, let $< h^{*}, \; m_{2} > = 0,$ we obtain

where

with

According to the above analysis, the normal form of Hopf bifurcation for system (2.1) reduced on the center manifold is

making $G \rightarrow G/\varepsilon$, thus, Eq. (4.14) becomes:

where $M$ and $\chi$ are given by Eq. (4.9) and Eq. (4.13), respectively.

Let $G = r{\text{e}^{\mathrm{i}\theta}}$ and substitute it into Eq. (4.15), and we obtain the Hopf bifurcation normal form in polar coordinates:

According to the bifurcation normal form Eq. (4.16) in polar coordinates, we have following theorem.

Theorem 4.1. For system (4.16), if $\frac{\operatorname{Re}(M)\mu}{\operatorname{Re}(\chi)} < 0$ holds, system (2.1) exists periodic solutions near equilibrium $E = (u_0, v_0)$, here $(u_0, v_0) = (\alpha, 0)$ or $(u^*, v^*)$ with $u^* = \frac{d}{\beta-p}, \; v^* = \frac{r}{m} (1-\frac{d}{\alpha (\beta-p)})$.

(1) If $\operatorname{Re}(M)\mu < 0$, the bifurcating periodic solutions are unstable.

(2) If $\operatorname{Re}(M)\mu > 0$, the bifurcating periodic solutions are locally asymptotically stable.

5.   Numerical simulations

According to Ref. [12] we choose the following parameter (see Table 2) for numerical simulation. Obviously, $(H_1)$ holds, and system (2.1) has two boundary equilibria $E_0 = (0, 0)$ and $E_1 = (\alpha, 0) = (5\times10^6, 0)$, and one nontrivial equilibrium $E^* = (u^*, v^*) = (2.5948\times10^5, 7.8643\times10^5)$. Actually, we only care about the stability of the nontrivial equilibrium $E^*$. In this section, we plot numerical simulations by using MATLAB software.

When $\tau_{1} = 0, \; \tau_{2} = 0$, according to Theorem 3.1, $E_{0} = (0, 0)$ is always unstable. This means that if system (2.1) without delays, the immune system cannot restrain the growth of tumor cells, and tumor cells proliferate quickly in the body. For equilibrium $E^*$, assumption $(H_2)$ and the conditions in Theorem 3.1 (3-2) hold, thus equilibrium $E^* = (u^*, v^*)$ is locally asymptotically stable. We choose the initial value $(2.5947 \times 10^{6}, \; 7.8642 \times 10^{5})$, $E^* = (u^*, v^*)$ is locally asymptotically stable (see Figure 2). Although tumor cells and immune cells can coexist at this moment, the immune system can suppress the growth of tumor cells effectively.

By a simple calculation, $(H_5)$ holds, from Eq. (3.7), we obtain $\tau_{12}^{(0)} = 14.8403$, we choose $\tau_{1} = 3 \in (0, \tau_{12}^{(0)})$ and choose the initial value $(2.5947 \times 10^{6}, \; 7.8642 \times 10^{5})$, according to Theorem 3.2, $E^* = (u^*, v^*)$ is locally asymptotically stable (see Figure 3). From Figure 3, we conclude that if tumor cells need a long time to proliferate, the tumor cells cannot diffuse, and it is easy to cure the tumor at this time. Thus, the diffusion of tumor cells can be controlled by making the size of the delay small. We reduce the diffusion of tumor cells by controlling the time of the proliferation of tumor cells.

We choose $\tau_1 = 15$ and the initial value $(2.5947 \times 10^{6}, \; 7.8642 \times 10^{5})$, from (4.15), we obtain $\mathrm{R e}(M) > 0, \; \mathrm{R e}(\chi) < 0$, thus, according to Theorem 4.1, system (2.1) exhibits stable periodic solutions near $E^{*}$, and the direction of Hopf bifurcation is forward. This means that, when $\tau$ is close to $\tau_{c}$, the number of tumor cells can be controlled, and the number of tumor cells and immune cells varies periodically (see Figure 4).

When $\tau_{2}\neq 0, \tau_{1}\neq 0$, the stable interval mentioned in Theorem 3.2 $I = \{\tau_1|0 < \tau_1 < \tau_{12}^{(0)} = 14.8403\}$, on which the equilibrium $E^*$ is locally asymptotically stable for $\tau_2 = 0$, and $(H_2)$ holds, from (3.8), we obtain $\tau_{21}^{(0)} = 4.3239$. We choose $\tau_{1} = 1 \in I, \; \tau_{2} = 1$, according to Theorem 3.3, we choose the initial value $(2.5947 \times 10^{6}, \; 7.8642 \times 10^{5})$, $E^* = (u^*, v^*)$ is locally asymptotically stable (see Figure 5). The time of tumor growth is fixed within a certain range, and the tumor cells cannot escape the control of the immune system.

Remark: Comparing the above figures, we have the following results. When the time delay is less than the critical value, the equilibrium is locally asymptotically stable, that is, the immune system can restrain the growth of tumor cells effectively at this moment, and tumor cells may not proliferate insanely in the body. Tumor cells originate from the malignant growth of normal cells. It is shown that the immune system can restrain tumor growth, although immune cells and tumor cells could coexist, and the immune system has higher efficiency (see Figures 1, 2 and 5). As the time delay of immune cells identifying tumor cells increased, immune system is able to identify itself from nonself substances, and immunotherapy can work effectively. Then, the numbers of tumor cells and immune cells demonstrate stable periodic variation (see Figure 3). When the time delay of immune cells identifying tumor cells is larger, the numbers of tumor cells are larger, that is, the immune system loses function, and the body may be damaged. On the other hand, the growing speed of tumor cells depends on the mood and diet of the patient, and it is very important to maintain an optimistic mood and a regular diet.

6.   Conclusion

We constructed a tumor-immune system with two delays and diffusion, and investigated how the delays affect the dynamics of the system. We analyzed the existence and stability of equilibria, and studied the dynamic properties of Hopf bifurcation. The nontrivial equilibrium is locally asymptotically stable under suitable parameters for the system with two delays, and the system occurred stable periodical solutions when the delay of tumor cell proliferation passed through the critical value.

Conflict of interest

The authors declare that they have no competing interests.

Funding

This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No. 2572021DJ01) and the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2019A001).

Availability of data and materials

The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

The idea of this research was introduced by Yuting Ding. All authors contributed to the main results and numerical simulations.

Appendix A

In this section, we show the code of MATLAB software for simulating partial differential equations (that is, Figure 2).

Appendix B

In this section, we show the code of MATLAB software for simulating partial functional differential equations (that is, Figures 3-5).