A nonlinear delayed model for the immune response in the presence of viral mutation

TL;DR

This paper introduces a delayed nonlinear model of immune response to viral infection, accounting for wild-type and mutant viruses, revealing oscillatory and chaotic dynamics caused by immune response delays.

Contribution

It presents a novel delayed nonlinear model capturing complex immune-virus interactions and analyzes the effects of immune response delays on system dynamics.

Findings

Delayed responses induce oscillations and chaos in immune-virus dynamics

Numerical analysis confirms complex behaviors due to delays

Analytical results support numerical findings

Abstract

In this paper we consider a delayed nonlinear model of the dynamics of the immune system against a viral infection that contains wild-type virus and one mutant. A finite response time of the immune system was considered in order which leads to sustained oscillatory behavior as well as chaotic behavior, triggered by the presence of delays. We present a thoroughly numeric analysis and some analytical results.

Equations

$$\dot{y}$$

$$\dot{y_m}$$

$$\dot{z}$$

$$y$$

$${\dot{x}}_i=\sum _{j=1}^k{F}_j^i\left(x_1(t-{\tau }_j),x_2(t-{\tau }_j),\cdots \right).$$

$${\dot{x}}_i\approx \sum _{j=1}^k\left(F_j^1(x_1,\cdots ){|}_{x^{\ast }}+\frac{\partial F_j^i}{\partial x_1}{|}_{x^{\ast }}(x_1(t-{\tau }_j)-x_1^{\ast })+\frac{\partial F_j^i}{\partial x_2}{|}_{x^{\ast }}(x_2(t-{\tau }_j)-x_2^{\ast })+\cdots \right).$$

$$J_{j}=\left[\begin{array}[]{ccc}\frac{\partial F_{j}^{1}}{\partial x_{1}}&\frac{\partial F_{j}^{1}}{\partial x_{2}}&\cdots\\ \vdots&\vdots&\vdots\\ \frac{\partial F_{j}^{n}}{\partial x_{1}}&\frac{\partial F_{j}^{n}}{\partial x_{2}}&\cdots\\ \end{array}\right];$$

$$J=\left[\begin{array}{ccc}r(1-\alpha )(1-\frac{2y^{\ast }}{K})-a-p{z}^{\ast }& 0& -p{y}^{\ast }\\ 0& {\alpha }_m{r}_m(1-\frac{2y_m}{K_m})-a_m-p_m{z}^{\ast }& -p_m{y}_m^{\ast }\\ \left(\frac{c{z}^{\ast }}{(1+d{y}^{\ast }{)}^2}-q{z}^{\ast }\right)e^{-\lambda {\tau }_1}& \left(\frac{c_m{z}^{\ast }}{(1+d_m{y}_m{)}^2}-q_m{z}^{\ast }\right)e^{-\lambda {\tau }_2}& \frac{c{y}^{\ast }}{(1+d{y}^{\ast })}+\frac{c_m{y}_m^{\ast }}{(d_m{y}_m^{\ast }+1)}-q{y}^{\ast }-q_m{y}_m^{\ast }-b\end{array}\right]$$

$$\tilde{J}=\left[\begin{array}[]{ccc}r(1-\alpha)-a&0&0\\ 0&r_{m}-a_{m}&0\\ 0&0&-b\end{array}\right].,$$

$$J_{SF}=(-0.4825880248-1.960851071\lambda )\exp (-\lambda \tau )-0.7961892098{\lambda }^2-{\lambda }^3$$

$$-0.08061172163\lambda +8.061172243\dot{1}{0}^{-11}.$$

$$F_{SF}=0.4726938145y^4+y^6-3.838438673y^2-\mathrm{0.2328912017.}$$

Full text

1 Introduction

We consider a nonlinear set of delay differential equations (DDEs) to model the interaction of the immune system with an external pathogen, e.g., a viral infection. Our model follows one presented in Ref. [1] in which a time delay takes into account the non-instantaneous immune response caused by a sequence of events (e.g., activation of antigenic response or production of immune cells) that occurs within a finite time period. In addition, the presence of sustained aperiodic oscillations and chaotic trajectories observed in real data [2, 3, 4] indicates that time delays are needed to allow bifurcations that cause chaotic behavior even in models that are one- and two-dimensional [5]. In ordinary differential equations (ODEs) a minimum set of three coupled equations is required.

Because the fundamental underlying mechanisms are non-instantaneous, several biological models have recently been modeled using delay differential equations. Among these are a predator prey model with delays [6], a model for the dynamics of the hormonal control of the menstrual cycle [7], a model for human respiration [8], a model for dioxide carbon levels in the blood [10, 9], and a number of models for viral dynamics [2, 5, 11, 1, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

In previous research [2] we analyzed the cellular immune response and found that stationary solutions bifurcate to an unstable fixed point when delays are longer than a critical immune response time $\tau_{c}$. We found that increasing the time delay causes the system to suffer a series of bifurcations that can evolve into a chaotic regime. We used two coupled delayed equations to model the interaction of the immune system with a target population [5]. We used some analytical tools to analyze delayed systems [11], and we published new results for the model originally presented in Ref. [2]. Here we consider a three-dimensional version of a model that previously appeared in the literature [1, 22] for the dynamics of the population of virus $y(t)$ and of immune cells $z(t)$, and also a mutant population of virus $y_{m}(t)$.

Delay differential equations require both the initial conditions and the history of the dynamic variable values of $t<\tau$. Because we are using models with discrete delays, $\tau$ is constant. This is in contrast to a system with distributed delays in which $\int_{t-r}^{t}k(t-s)x(s)ds=\int_{0}^{r}k(z)x(t-z)dz$, where $0\leq r\leq\infty$ is the distributed delay and the kernel $k$ is normalized, and thus $\int_{-\infty}^{\infty}k(y)dy=1$. For an identically null $k(u),\forall u>u_{\rm max}$ the delay can be represented by integrals of type $\int_{-\infty}^{t}M_{1}(s)k(t-s)ds=\int_{0}^{\infty}M_{1}(t-u)k(u)du$. These are “bounded delays” because they represent the values of $M_{1}$ at a past time $(t-u_{\rm max},t)$. A discrete delay is a particular kind of bounded delay. More complicated forms are also possible, e.g., delays of type $x(t-r[x(t)])$ distributed over space.

Introducing delays allows us to model richer behavior, e.g., the well-known logistic equation governing the dynamics of a population density $N(t)$: $\dot{N}(t)=N(t)\left(1-\frac{N(t)}{K}\right)$, with $r$ the growth rate and $K$ the carrying capacity. Note that for every initial condition $N(0)$ the system ultimately reaches the stable equilibrium $N(t)\rightarrow K$. A delayed version of this model can be used for a species population that gathers and stores food, i.e., when resources vanish, the species population starves within finite time $\tau$. Reference [23] assumes this and analyzes the delayed system $\dot{N}(t)=N(t)r\left(1-\frac{N(t-\tau)}{K}\right)$. This delayed version of the logistic equation can model chaotic behavior that instantaneous one dimensional models cannot because ODE systems need at least a three-dimensional state space to model chaos, as demonstrated in Lorenz’s seminal work [24]. Here the number of initial conditions is equal to the number of degrees of freedom. In delayed systems the number of degrees of freedom is infinity and chaos occurs in even one dimensional systems, as in the case for one-dimensional non-invertible maps.

We present the model in the next section. In section 3 we present some analytical and numeric results, and in section 4 we present our conclusions.

2 Model

Our model is based on research described in Refs. [1, 22] that uses a two-dimensional model for the dynamics of the population of virus $y(t)$ and of immune cells $z(t)$. We use time-lagged response for the immune system, following previous research demonstrating its importance in the appearance of the Hopf bifurcations [3], chaotic trajectories [2], and sustained oscillatory behavior rarely seen in the instantaneous version of the model [4]. Here we extend the model to a spreading population of mutant virus $y_{m}(t)$,

[TABLE]

where $r(1-\alpha)$ is the growth rate of the viral population for $y\approx 0$. This rate decreases and reaches zero when the virus population equals $K$. The virus population decays with $a$. We then have a net rate of $r(1-\alpha)-a$ and a carrying capacity of $\frac{K(r(1-\alpha)-a)}{r(1-\alpha)}$. The viruses are eliminated by cells of the immune system according a rate $p$. The term $y_{m}$ represents the concentration of the mutant viruses. Its net growth rate and carrying capacity are, respectively, $r_{m}\alpha_{m}-a_{m}$ and $\frac{K_{m}(r\alpha_{m}-a)}{r\alpha_{m}}$. They are eliminated at a rate $p_{m}$. The immune cell concentration $z$ grows proportionally to the virus population according to a saturation term. The $\tau_{2}$ value is the delay in the immune response to the viral infection. The delay $\tau_{1}$ refers to the processes used by the organism to prepare the cells to fight the virus. Immune cells are attacked and destroyed by the original viruses and their mutant version with rates $q$ and $q_{m}$, respectively. The terms $1/(1+dy(t-\tau_{1}))$ and $1/[1+d_{m}y_{m}(t-\tau_{2})]$ shows that the immune response is proportional to the product of the virus (either $y$ or $y_{m}$) and the population of immune cells $z$, but saturates when the virus population is higher. Numerical estimations of the parameters are provided in Ref. [22].

[TABLE]

Identical numeric values are assumed for $K_{m}$, $r_{m}$, $a_{m}$, $c_{m}$, $d_{m}$, and $q_{m}$. The $p_{m}=0.9<p$ value is an exception because here it is more difficult for the immune system to eliminate cells infected by the mutant virus. We also assume $\alpha=1$ and $\alpha_{m}=0.05$, which indicates that the mutation is a residual portion of the replication mechanisms.

3 Results

Reference [1] presents several analytical results for the two-dimensional version presented in (1), which does not take into account the mutant population $y_{m}$. Because our model is three-dimensional it is cumbersome to analyze, and we focus on numeric results. Similar to the procedure used in the logistic map, we focus on the emergence of bifurcations and chaos as time-delay values increase. The system in (1) has a total of 11 equilibrium points. Six are facial points (with at least one null component). Those with simple algebraic expressions are

[TABLE]

The others present cumbersome algebraic expressions, which we omit here for sake of simplicity.

Note that the stability of the fixed points of a $n$-dimensional system with $k$ delays can be analyzed using the usual Jacobian evaluated at the equilibrium point [11]. Each $\dot{x_{i}},i=1,\cdots,n$ be written

[TABLE]

Performing a series expansion around the equilibrium point $x^{*}=(x_{1}^{*},\cdots,x_{n}^{*})$, we obtain for each $x_{i},i=1,\cdots,n$

[TABLE]

We then have a linear system of variables $\tilde{x}_{i}\equiv x_{i}-x_{i}^{*}$ with $k$ Jacobian matrices that take the form

[TABLE]

evaluated at the fixed points. The stability of a particular fixed point is determined by the eigenvalues of its corresponding Jacobian. Bifurcations occur whenever one eigenvalue crosses the imaginary axis as one or more parameters, including the delays, change. Typical bifurcations involve a turning point when the eigenvalue is initially null, and a Hopf bifurcation when a pair of complex eigenvalues crosses the imaginary axis [11].

The general expression for the Jacobian is

[TABLE]

The Jacobian for the origin is thus

[TABLE]

which holds for all values of $\tau_{1},\tau_{2}$. The eigenvalues are $-\alpha r-a+r$, $r_{m}-a_{m}$, and $-b$. Stability (with only negative eigenvalues) can be achieved for smaller $r$ and $r_{m}$ (the viral growth rate), and larger $a,a_{m}$ (the natural population decay of the virus). Here ultimately the system loses all of its viruses and has no immune cells irrespective of the delay. For the set of parameters chosen here, however, the origin is unstable $\forall\tau_{1},\tau_{2}$. Note that for the three fixed points in (3) the stability is unchanged when there are non-null delays. This can be seen from (4) by substituting $z=0$. Note that $r-\alpha$, $r-a$, and $-r_{m}+a_{m}$ are common eigenvalues, a condition that renders the origin unstable for all three.

For null delays and the chosen set of parameters, only two of the 11 equilibria are stable. Because $y,y_{m}$ and $z$ are densities and therefore positive quantities, one stable equilibrium is physically irrelevant: $y=5.568989996$, $y_{m}=5.046486447$, and $z=-7.88108099$. The other stable equilibrium is the spiral focus (SF): $y=0.06265629108$, $y_{m}=0.3385711289$, and $z=2.580953047$. For the parameters used, the remaining equilibria are all unstable and comprise six facial equilibria and two physically-irrelevant equilibria that have at least $y<0$ or $y_{m}<0$ or $z<0$. Here we focus on how increasing the value of the time delay alters the stability of the stable SF solution.

A theorem presented in Ref. [25] describes the conditions for switches in stability when there are delays and finds a critical $\tau^{*}>0$ above which the equilibrium point is always unstable. The theorem states:

Theorem 1

Let a characteristic equation of a given fixed point be written $R(\lambda)+S(\lambda)\exp(-\lambda\tau)=0$. $R(\lambda)$ and $S(\lambda)$ are analytical in the right half plane and $\Re{\lambda}>-\delta,\delta>0$. When the following properties hold:

• (i)

$R(\lambda)$* and $S(\lambda)$ have no common zero;*

• (ii)

$\overline{R(-Is)}=R(Is)$,$\overline{S(-Is)}=R(Is)$, where the bar indicates the conjugate and $I=\sqrt{(}-1)$;

• (iii)

$R(0)+S(0)=0$;

• (iv)

The half right plane possesses at most a finite number of roots of $R(\lambda)+S(\lambda)\exp(-\lambda\tau)=0$ when $\tau=0$; and

• (v)

$F(y)=\left|R(Iy)\right|^{2}-\left|S(Iy)\right|^{2}$* when real $y$ has at most a finite number of zeros.*

Then the following statements are true:

• (a)

If $F(y)=0$ has no positive real roots, and if the associated fixed point is stable (unstable) for null delays, it will remain stable (unstable) for all delays.

• (b)

If $F(y)=0$ has at least one positive root and all roots are simple, stability switches can occur with increasing $\tau$. There exists a $\tau^{*}>0$ above which the fixed point is unstable for all $\tau>\tau^{*}$. As $\tau$ varies from zero to $\tau^{*}$ at most a finite number of stability switches may occur.

Reference [5] used this theorem to analyze their model. Here we consider equal delays $\tau_{1}=\tau_{2}=\tau$. The Jacobian of the spiral focus stable equilibrium ($y=0.06265629108$, $y_{m}=0.3385711289$, $z=2.580953047$ is

[TABLE]

This equation is clearly of type $R(\lambda)+S(\lambda)\exp(-\lambda\tau)$. Thus $F(Iy)=\left|R(Iy)\right|^{2}-\left|S(Iy)\right|^{2}$ yields

[TABLE]

The roots of this equation are $\pm 1.330454624,\pm 0.2455259061I,\pm 1.477335558I$, and thus the condition $b$ of item $v$ of the theorem holds. Figure 1 shows the expected stability switches. We plot $z(t)$ versus $y(t)$ for $\tau_{1}=\tau_{2}=0$ (the stable case), $\tau_{1}=\tau_{2}=0.2$, $\tau_{1}=\tau_{2}=5$, and $\tau_{1}=\tau_{2}=15$. There is still stability for delay $\tau=\tau_{1}=\tau_{2}=0.2$, but this is lost in $\tau=5$ and $\tau=15$. These results demonstrate how the introduction of delays can change the stability of a stable solution and promote a richer dynamics for the system.

For the sake of comparison, we use the two-dimensional model proposed in Ref. [1] (which has no mutant virus) and plug $y_{m}=0$ and $K_{m},r_{m}a_{m},c_{m},d_{m},q_{m},\tau_{2}=0$ into (1). Figure 2 shows the maxima values of $z(t)$ versus $\tau_{1}$. Note that there is a series of bifurcations that switches between sustained oscillations and chaotic behavior, with windows of periodic behavior (e.g., around $\tau_{1}=14$).

When we use the term $\frac{c_{m}y_{m}(t)z(t)}{1+d_{m}y_{m}(t)}$ with $\tau_{2}=0$ to introduce the mutant component into our model, it changes the dynamics of (1). Figure 3 shows that periodic orbits are present but not chaotic behavior. Although merely inserting a new equation into the system does not enrich the dynamics, the situation changes completely when $\tau_{2}\neq 0$. Figure 4 shows that this time delay causes more complex patterns to emerge, including regions of chaotic behavior.

4 Conclusion

We have considered a nonlinear set of delay differential equations to model the interaction between an immune system and an external pathogen, e.g., a viral infection. We extend the previous model considered in [1] by introducing a new variable that takes into account mutant viruses. We find a series of bifurcations that lead to chaotic behavior, an outcome that agrees with the results observed in real data [2, 3, 4] and that corroborates previous work indicating the need for the time delays that generate richer behavior [1].

Acknowledgments

AC thanks the Alagoas State Research Agency FAPEAL for support through major projects (PPP - 20110902-011-0025-0069 / 60030-733/2011), also CNPq for PDE (207360/2014-6) and Universal (423713/2016-7) grants. DM acknowledges a scholarship by the Brazilian funding agency CAPES. The Boston University work was supported by DTRA Grant HDTRA1-14-1-0017, by DOE Contract DE-AC07-05Id14517, and by NSF Grants CMMI 1125290, PHY 1505000, and CHE-1213217.

5 References