Dynamic Analysis of a Model on Tumor-Immune System with Regulation of PD-1/PD-L1 and Stimulation Delay of Tumor Antigen

Abstract

We propose and investigate a mathematical model on interaction between tumor and the immune system, where the regulation of PD-1/PD-L1 and the stimulation delay of tumor antigen for the immune system are considered. Though delay will not change the structure of equilibria, the global dynamics in the case without delay is simple compared with that in the case with delay. Theoretic analysis and numerical simulations show that the incorporation of delay leads to complex dynamics, including the appearance of oscillating solutions, periodic solutions from Hopf bifurcation, and homoclinic orbits, etc. The effect of the immunotherapy including anti-PD-1/PD-L1 inhibitor and tumor vaccine is also discussed.

Appendices

High Order Equilibrium of System (4)

1.1 Local Dynamics of System (4) Near $P_0$

To investigate the dynamical behavior of system (4) near $P_0$ when $\beta _2=1$ and $c\ne \beta _1$, we let $u=x-1$ and $v=y$ to transform (4) into

$$\begin{aligned} \begin{aligned}&u'=-u+(c-\beta _1)v-\beta _1uv, \\&v'=-rv[(u+v)+\beta _0(1+u)v]. \end{aligned} \end{aligned}$$

(14)

Further, let $u=w+(c-\beta _1)v$. Then (14) becomes

$$\begin{aligned} \begin{aligned}&w'=-w-\beta _1v(w+cv-\beta _1v)\triangleq F(w, v), \\&v'=-rv[(1+\beta _0v)w+(c+\beta _0+1-\beta _1)v +\beta _0(c-\beta _1)v^2] \triangleq G(w, v). \end{aligned} \end{aligned}$$

(15)

By the center manifold theorem [8], we assume that the local center manifold of (15) at the origin is expressed by the function

$$\begin{aligned} w=h(v)=av^2+o(v^2), \end{aligned}$$

where a is a constant to be determined. From $\frac{dw}{dt}=[2av+o(v)]\frac{dv}{dt}$ and (15) we have

$$\begin{aligned}{}[2av+o(v)]vG(h(v),v)-h(v)-vF(h(v), v)=0. \end{aligned}$$

(16)

By letting the coefficient of $v^2$ in the left-hand side of (16) be zero, we get $a=\beta _1(\beta _1-c)$.

Substituting the obtained local center manifold $w=\beta _1(\beta _1-c)v^2+o(v^2)$ into the second equation of (15) yields

$$\begin{aligned} v'=r(\beta _1-1-c-\beta _0)v^2+r(\beta _1-\beta _0)(c-\beta _1)v^3+o(v^3). \end{aligned}$$

(17)

Then when $\beta _1\ne 1+c+\beta _0$, $P_0$ is a saddle-node. Further, according to the first equation of (15), when $\beta _1< 1+c+\beta _0$, $P_0$ is attracting in the part corresponding to the node while when $\beta _1> 1+c+\beta _0$, $P_0$ is repelling in the part corresponding to the node. When $\beta _1=1+c+\beta _0$, (17) becomes

$$\begin{aligned} v'=-r(1+\beta _0)(1+c)v^3+o(v^3). \end{aligned}$$

It follows that, when $\beta _2=1$ and $\beta _1= 1+c+\beta _0$, $P_0$ is a locally asymptotically node by the first equation of (15).

In summary, this gives the local stability of $P_0$ when $\beta _2=1$ and $c\ne \beta _1$. $\square$

1.2 Local Dynamics of System (4) Near $P^*_*$

In order to study the dynamical behavior of (4) near $P^*_*$, we translate $P^*_*$ to the origin by the transformation $u=x-x^*_*$ and $v=y-y^*_*$. Then (4) becomes

$$\begin{aligned} \begin{aligned}&u'=-(1+\beta _1y^*_*)u+(c-\beta _1x^*_*)v-\beta _1uv, \\&v'=-r(v+y^*_*)[(\beta _0y^*_*+\beta _2^*)u +(1+\beta _0x^*_*)v+\beta _0uv]. \end{aligned} \end{aligned}$$

(18)

It follows from $\det J(P^*_*)=0$ that $c-\beta _1x^*_*=-\frac{(1+\beta _1y^*_*)(1+\beta _0x^*_*)}{\beta _0y^*_*+\beta _2^*}$. Then (18) is rewritten as

$$\begin{aligned} \begin{aligned}&u'=-au-\frac{am}{b}v-\beta _1uv, \\&v'=-r(v+y^*_*)(bu+mv+\beta _0uv], \end{aligned} \end{aligned}$$

(19)

where $a=1+\beta _1y^*_*$, $b=\beta _0y^*_*+\beta _2^*$, and $m=1+\beta _0x^*_*$.

Furthermore, using the reversible transformation of variables

$$\begin{aligned} \begin{pmatrix} u\\ v \end{pmatrix} =\begin{pmatrix} a &{} -\frac{am}{b} \\ r y^*_*b &{} a \end{pmatrix} \begin{pmatrix} w\\ z \end{pmatrix}, \end{aligned}$$

that is,

$$\begin{aligned} \begin{pmatrix} w\\ z \end{pmatrix}= \frac{1}{a(a+ry^*_*m)} \begin{pmatrix} a &{} \frac{am}{b} \\ -r y^*_*b &{} a \end{pmatrix} \begin{pmatrix} u\\ v \end{pmatrix}, \end{aligned}$$

we can transform (19) into

$$\begin{aligned}&w'= -(a+ry^*_*m)w+p_{20}w^2+p_{11}wz+p_{02}z^2 +p_{30}w^3+p_{21}w^2z+p_{12}wz^2 +p_{03}z^3, \nonumber \\&z'= q_{20}w^2+q_{11}wz+q_{02}z^2 +q_{30}w^3+q_{21}w^2z+q_{12}wz^2 +q_{03}z^3, \end{aligned}$$

(20)

where

$$\begin{aligned} p_{20}= & {} -\frac{ry^*_*[ab(\beta _1+mr) +rmy^*_*(bmr+a\beta _0)]}{a+ry^*_*m}, \\ p_{11}= & {} -\frac{a [ab(\beta _1+mr)+ry^*_*m(mrb-mr\beta _0y^*_* +a\beta _0-b\beta _1) ]}{b(a+ry^*_*m)}, \\ p_{02}= & {} \frac{a^2m(b\beta _1+mr\beta _0y^*_*)}{(a+ry_*^*m)b^2}, \\ p_{30}= & {} -\frac{mr^3{y^*_*}^2b\beta _0a}{a+ry^*_*m}, \\ p_{21}= & {} -\frac{mr^2y^*_*\beta _0a(2a-ry^*_*m)}{a+ry^*_*m}, \\ p_{12}= & {} -\frac{mr\beta _0a^2(1-2ry^*_*m)}{b(a+ry^*_*m)}, \\ p_{03}= & {} \frac{m^2ra^3\beta _0}{(a+ry^*_*m)b^2}, \\ q_{20}= & {} -\frac{r^2y_*^*b((ay^*_*\beta _0+ba -by^*_*\beta _1+bmry^*_*))}{a+ry^*_*m}, \\ q_{11}= & {} -\frac{r[(b\beta _1-a\beta _0)(rmy^*_*-a)y^*_* -ab(a+rmy^*_*)]}{a+ry^*_*m}, \\ q_{02}= & {} \frac{ray^*_*m(a\beta _0-b\beta _1)}{b(a+ry_*^*m)}, \\ q_{30}= & {} -\frac{r^3{y^*_*}^ 2ab^2\beta _0}{a+ry^*_*m}, \\ q_{21}= & {} -\frac{r^2y^*_*a\beta _0b(2a-ry^*_*m)}{a+ry^*_*m}, \\ q_{12}= & {} -\frac{ra^2\beta _0(a-2ry^*_*m)}{a+ry^*_*m}, \\ q_{03}= & {} \frac{ra^3\beta _0m}{b(a+ry*_*m)}. \end{aligned}$$

Then the local center manifold of (20) at the origin can be obtained similarly as that of (15), which is

$$\begin{aligned} w=\frac{a^2m(b\beta _1+mr\beta _0y^*_*)}{b^2(a+ry^*_**m)^2}z^2+o(z^2). \end{aligned}$$

Substituting it into the second equation of (20) yields

$$\begin{aligned} z'=-\frac{ray^*_*m(\beta _1\beta _2^*-\beta _0)}{b(a+ry^*_*m)}z^2+o(z^2), \end{aligned}$$

where $b\beta _1-a\beta _0=\beta _1\beta _2^*-\beta _0$ has been used. From $\varDelta (\frac{\beta _0}{\beta _1}) =\frac{(\beta _1^2+\beta _1-\beta _0\beta _1+c\beta _0)^2}{\beta _1^2}>0$ for $\beta _1>1+c+\beta _0$, we know that $\beta _1\beta _2^*-\beta _0\ne 0$ since $\varDelta (\beta _2^*)=0$. Therefore, $P^*_*$ is a saddle-node (see, for example, [8]). $\square$

Proof of Proposition 1

During the discussion, Maple has been used to do some symbolic calculations.

• (i) A straightforward calculation gives

  $$\begin{aligned} \beta _2^{**}(1+\beta _0)=\frac{2(1+\beta _0) {[}\sqrt{(1+\beta _0)^2+c^2(c\beta _0+3+3\beta _0)} -(1+\beta _0)]}{c^2(c\beta _0+3+3\beta _0)}>0 \end{aligned}$$

  and

  $$\begin{aligned}&\beta _2^{**}(1+\beta _0)-1 \\= & {} \frac{2(1+\beta _0) \sqrt{(1+\beta _0)^2+c^2(c\beta _0+3+3\beta _0)} -[2(1+\beta _0)^2 +c^2(c\beta _0+3+3\beta _0)]}{c^2(c\beta _0+3+3\beta _0)}. \end{aligned}$$

  Since

  $$\begin{aligned}&[2(1+\beta _0)\sqrt{(1+\beta _0)^2+c^2(c\beta _0 +3+3\beta _0)}]^2 -[2(1+\beta _0)^2 +c^2(c\beta _0+3+3\beta _0)]^2 \\= & {} - c^4(c\beta _0+3+3\beta _0)^2 \\< & {} 0, \end{aligned}$$

  it follows that $\beta _2^{**}(1+\beta _0)<1$.

• (ii) On the one hand, we have

  $$\begin{aligned} \beta _2^{**}(1+\beta _0+c)=\frac{2(1+c+\beta _0) \sqrt{\varDelta _0}-Q}{c^2[3(1+c)+\beta _0(3+c)]}, \end{aligned}$$

  where

  $$\begin{aligned} \varDelta _0= & {} (1+c+\beta _0)^2+c^2(c+2)(c+2\beta _0+2), \\ Q= & {} 2(1+c+\beta _0)^2+c^2(1+c)(\beta _0+1). \end{aligned}$$

  Since

  $$\begin{aligned}&4(1+c+\beta _0)^2\varDelta _0^2-Q^2 \\= & {} c^2[3(1+c+\beta _0)+c\beta _0] [4(1+c+\beta _0)^2+c^2(1+\beta _0)-(\beta _0-1)c^3], \end{aligned}$$

  we can easily see that $\beta _2^{**}(1+\beta _0+c)$ can be positive or negative. On the other hand,

  $$\begin{aligned}&\beta _2^{**}(1+\beta _0+c)-1 \\= & {} \frac{\begin{array}{l} 2 \{(1+c+\beta _0)\sqrt{(c^2+c+1)^2+2 c\beta _0(c+1)^2+2c^2(c+1)+\beta _0(\beta _0+2)} \\ - (1+c+\beta _0)^2-c^2[2(c+1)+\beta _0(2+c)] \}\end{array}}{c^2(c\beta _0+1+c+\beta _0)} \end{aligned}$$

  and

  $$\begin{aligned}&(1+c+\beta _0)^2[(c^2+c+1)^2+2 c\beta _0(c+1)^2+2c^2(c+1)+\beta _0(\beta _0+2) ] \\&- \{ (1+c+\beta _0)^2+c^2[2(c+1)+\beta _0(2+c)]\}^2 \\= & {} -c^4(1+\beta _0)(c+1)(3+3\beta _0+3c+\beta _0c) \\< & {} 0 \end{aligned}$$

  thus $\beta _2^{**}(1+\beta _0+c)<1$.

• (iii) It is easy to see that the equation $\beta _2^{**}(\beta _1)=1$ has two distinct roots with different signs. A straightforward calculation finds that the positive root is $\beta _1={\hat{\beta }}_1$, where

  $$\begin{aligned} {\hat{\beta }}_1=\frac{2c+(\beta _0+1)(c+1) +\sqrt{(\beta _0+1)[\beta _0(c-1)^2+(3c+1)^2]}}{2}. \end{aligned}$$

  Then

  $$\begin{aligned} {\hat{\beta }}_1-(1+c+\beta _1)=\frac{ (\beta _0+1)(c-1)+\sqrt{(\beta _0+1)[\beta _0(c-1)^2+(3c+1)^2]}}{2}>0 \end{aligned}$$

  since $(\beta _0+1)[\beta _0(c-1)^2+(3c+1)^2]>(\beta _0+1)^2(c-1)^2$.

• (iv) For simplicity of notation, denote $r_1=2\beta _1+(\beta _0-1+\beta _1)c$, $r_2=c(\beta _1-1-\beta _0)(c\beta _0+\beta _1)+2\beta _1^2$, and $\varDelta _0=(1+c) (\beta _1-c)(c\beta _0+\beta _1)$. We can rewrite $\beta _2^*(\beta _1)$ and $\beta _2^{**}(\beta _1)$ as

  $$\begin{aligned} \beta _2^*(\beta _1)=\frac{r_1-2\sqrt{\varDelta _0}}{c^2}, \qquad \beta _2^{**}(\beta _1)=\frac{2 \beta _1\sqrt{\varDelta _1}-r_2}{c^2(c\beta _0+3\beta _1)}. \end{aligned}$$

  Then

  $$\begin{aligned} \beta _2^*(\beta _1)-\beta _2^{**}(\beta _1) =\frac{ [r_1(c\beta _0+3\beta _1)+r_2 ] -2[ (c\beta _0+3\beta _1)\sqrt{\varDelta _0}+ \beta _1\sqrt{\varDelta _1} ]}{c^2(c\beta _0+3\beta _1)}. \end{aligned}$$

  A direct calculation gives

  $$\begin{aligned}&{[}r_1(c\beta _0+3\beta _1)+r_2]^2 -4\left[ (c\beta _0+3\beta _1)\sqrt{\varDelta _0}+ \beta _1\sqrt{\varDelta _1} \right] ^2\\&\quad = 4(c\beta _0+3\beta _1) ( \varDelta _2- \beta _1\sqrt{\varDelta _0\varDelta _1}), \end{aligned}$$

  where

  $$\begin{aligned} \varDelta _2= & {} c^3\beta _0(\beta _0+c\beta _0+2\beta _1+1) +\beta _1(c^2+2\beta _1^2) \\&+c\beta _1 [c(\beta _1-\beta _0)(\beta _0+c\beta _0+\beta _1) +2\beta _1(\beta _1-1) ] \\> & {} 0 \end{aligned}$$

  for $\beta _1>1+c+\beta _0$. Furthermore,

  $$\begin{aligned} \varDelta _2^2- \beta _1^2{\varDelta _0\varDelta _1} = c^4[\beta _0(c+1)+(1+\beta _1) ]^2 [\beta _1(\beta _1+1)+\beta _0(c-\beta _1) ]^2>0 \end{aligned}$$

  for $\beta _1>1+c+\beta _0$. Therefore, we have $\beta _2^*(\beta _1)-\beta _2^{**}(\beta _1)>0$ for $\beta _1>1+c+\beta _0$.

• (v) For the quadratic function $\varPhi _0 (\beta _2 )$ of $\beta _2$, it is easy to verify that $\varPhi _0(\beta _2^{(0)}(\beta _1))<0$. It follows easily that when $\varDelta _1>0$, the two roots of $\varPhi _0 (\beta _2)$, ${\hat{\beta }}_2^{**}(\beta _1)$ and $\beta _2^{**}(\beta _1)$, satisfy ${\hat{\beta }}_2^{**}(\beta _1)<\beta _2^{(0)}(\beta _1)<\beta _2^{**}(\beta _1)$. $\square$ This completes the proof of Proposition 1.