Nonlocal Reaction–Diffusion Model of Viral Evolution: Emergence of Virus Strains

Abstract

This work is devoted to the investigation of virus quasi-species evolution and diversification due to mutations, competition for host cells, and cross-reactive immune responses. The model consists of a nonlocal reaction–diffusion equation for the virus density depending on the genotype considered to be a continuous variable and on time. This equation contains two integral terms corresponding to the nonlocal effects of virus interaction with host cells and with immune cells. In the model, a virus strain is represented by a localized solution concentrated around some given genotype. Emergence of new strains corresponds to a periodic wave propagating in the space of genotypes. The conditions of appearance of such waves and their dynamics are described.

1. Introduction

Human infections with rapidly evolving viruses such as the human immunodeficiency virus (HIV) or the hepatitis C virus (HCV) remain a challenge for health-care systems. Infections are usually initiated by one or few virions that then replicate and generate a swarm of progeny viruses with distinct but related genomes [1,2]. Collectively these swarms of viruses are called a virus quasi-species [3,4,5,6,7]. This quasi-species nature enables viruses to rapidly evolve within an infected host organism and adapt to constraints mediated by immune responses or antiviral drugs [8,9]. It also allows viruses to broaden their host cell tropism and to spread to diverse tissues [10]. Well studied examples for virus adaptation are the development of drug resistance or the generation of variants within virus-specific cytotoxic T lymphocyte (CTL) epitopes that diminish immune recognition and destruction of infected cells [11,12,13]. Since the immune system can also adapt to respective virus changes [14], an increase in the number of CTL target regions over time of infection as well as successive shifts in the hierarchy of immunodominance have been observed [15,16].

Gaining a mechanistic understanding of the dynamic interplay between the processes of virus replication, mutation, and elimination by immune responses and drug-based treatment requires the development of mathematical models which could be used to predict the generation of viral variants that escape the immune recognition and confer resistive to antiviral drugs. The existing models of virus evolution are based on the concept of quasi-species. i.e., an ensemble of related genomes [4]. The models can be formulated either as deterministic high-dimensional systems of ODEs, describing the densities of individual strain [15,17] or stochastic models with genetic algorithms [18]. The cooperative interactions in viral populations are considered to be key for linking the quasi-species dynamics in a changing virus-host environment with the genetic markers of viral evolution and the disease pathogenesis [10,19]. This implies that nonlocal interactions between the quasi-species in the genotype space need to be considered to predict the evolution of viruses to form distinct phenotypes.

Nonlocal reaction–diffusion equations represent an appropriate framework to describe evolution of biological species [20,21,22]. These equations take into account nonlocal consumption of resources characterizing intraspecific competition and possibly leading to the emergence of multi-modal population density distributions. Considered to be depending on a morphological characteristic and on time, localized in-space distributions can be interpreted as biological species, and the emergence of multi-modal distributions corresponds to the appearance of new species. In this work we will study virus quasi-species and will analyze the emergence of new strains in the space of genotypes. We consider the nonlocal reaction–diffusion equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qJ\left(u\right))-uf\left(S\left(u_{\tau }\right)\right)-\sigma \left(x\right)u$$

(1)

introduced in the previous work [23] devoted to the existence and dynamics of virus strains, but not to the emergence of new strains since this question is different both from the biological and modelling points of view. Here $u(x,t)$ is a dimensionless virus density distribution depending on its genotype x considered to be a continuous variable and on time t. The diffusion term in the right-hand side of this equation characterizes virus mutations, and the other terms describe virus reproduction, its elimination by immune response and by genotype-dependent mortality, either natural or caused by an antiviral treatment.

We describe the virus reproduction and immune response terms in more detail. We begin with the diffusion term. Assuming that there is a sequence of reversible mutations with consecutive genotypes $x_i$, we can write the equation for the density $u_i$ of virus with genotype $x_i$:

$$\frac{d{u}_i}{dt}=\mu (u_{i-1}-u_i)+\mu (u_{i+1}-u_i),$$

where $\mu$ is the frequency of mutations. This equation represents a discretization of the diffusion equation with the diffusion coefficient proportional to $\mu$.

The virus reproduction rate is conventionally considered either proportionally to its density u or, if we take into account the limitation on the quantity of the host cells where virus can multiplicate, as a logistic term $ru(1-qu)$. Here r is a proportionality coefficient, and q is a positive constant corresponding to the inverse carrying capacity in population dynamics. The latter case implicitly signifies that there is one-to-one correspondence between the virus genotype and the type of infected cells. In a more general and biologically realistic setting we should accept that viruses with different genotypes can infect the same cells. In this case, we replace the conventional logistic term by the term $ru(1-qJ(u\left)\right)$, where $J\left(u\right)={\int }_{-\infty }^{\infty }\varphi (x-y)u(y,t)dy$. The kernel $\varphi (x-y)$ in this integral characterizes the efficacy of host cell infection depending on the difference in genotypes. In general, it is a decreasing function of the modulus of its argument. Its exact form in the applications is not known, and different examples will be considered below. The integral is taken in the infinite limits for convenience of presentation. It implies that the genotype space is sufficiently large, and it can be mathematically approximated as a real line.

The reproduction rate with the integral $J\left(u\right)$ corresponds to nonlocal consumption of resources in population dynamics [22,24]. If we replace the kernel $\varphi \left(x\right)$ by the $\delta$-function, then we obtain the previous “local” case. Finally, if cell contamination is independent of virus genotype, then we have the integral $I\left(u\right)={\int }_{-\infty }^{\infty }u(y,t)dy$ corresponding to the global consumption of resources. Behavior of solutions of Equation (1) can be essentially different in these three cases.

Virus elimination by immune cells is proportional to the virus density u and to the concentration of immune cells C. Since immune response is stimulated by the antigen (virus), then the concentration of immune cells can be considered to be a function of virus density $C=f\left(u\right)$. The function $f\left(u\right)$ characterizes the intensity of immune response. It is positive and growing for some limited interval of values of u where the antigen stimulates the immune response, and it is decreasing for u sufficiently large due to the exhaustion and death of immune cells provoked by large virus concentrations [25,26,27]. The qualitative form of this function is well described by the dependence $f\left(u\right)=(k_{1}u+k_2)e^{-k_{3}u}$ with some positive constants $k_1,k_2,$ and $k_3$. The approximation of the concentration of immune cells as a function of virus density can be derived from a more complete model under the assumption of large reaction rate constants [28].

Clonal expansion of immune cells requires several cell proliferations and differentiations, and it usually takes 3–4 days. Therefore, the rate of virus elimination by immune cells can take into account this time delay if it is not negligible in the time scale of virus evolution. In this case, the corresponding term becomes $uf\left(u_{\tau }\right)$, where $u_{\tau }=u(x,t-\tau )$. Furthermore, similar to virus reproduction, virus elimination is also nonlocal in the genotype space. In [23] it was described by the term ${\int }_{-\infty }^{\infty }\theta (x-y)f\left(S\left(u_{\tau }\right)(y,t)\right)dy$, where $S\left(u_{\tau }\right)(y,t)={\int }_{-\infty }^{\infty }\psi (y-z)u(z,t-\tau )dz$. The inner integral S characterizes a cross-reactive stimulation of immune response by different antigens, while the outer integral describes a cross-reactive virus elimination by different immune cells. Both assumptions are biologically justified. However, this model becomes excessively complex, and we will restrict ourselves here only to the inner integral assuming the outer term is local, i.e., that the kernel $\theta \left(x\right)$ is replaced by the $\delta$-function.

The last term in the right-hand side of Equation (1) describes virus mortality with the rate depending on its genotype. The viability interval, i.e., the rate of genotypes where virus multiplication rate exceeds its mortality can depend on its intrinsic features and on an antiviral treatment.

Some particular cases of Equation (1) are studied in the literature. Considered without immune response and genotype-dependent mortality, this nonlocal reaction–diffusion equation and some its variations were widely studied in relation to various applications [29,30,31,32] and from the point of view of their mathematical properties [33,34,35,36]. One of the main features of this equation is that its homogeneous in-space stationary solution can become unstable leading to the emergence of periodic in-space solutions. We will return to this question below. The local equation ($J\left(u\right)\to u,S\left(u\right)\to u$) with time delay in the immune response term was suggested in [23,26,27,37] as a model of virus spread in tissues. The presence of time delay can lead to complex patterns of wave propagation.

In our previous work [23] we studied the existence and dynamics of virus strains considered to be localized solutions (pulses) in the space of genotypes, with the understanding that a virus strain can be characterized by its most frequent genotype with a narrow density distribution around it. The existence and stability of such solutions is not a priori given. In the local bistable equation such solutions exist but they are not stable. In the local monostable equation such solutions do not exist. It was previously shown that stable pulses exist for the nonlocal bistable equation [38]. In [23], it was revealed that persistent virus strains can exist due to the interaction of nonlocal (global) virus reproduction with immune response or with genotype-dependent mortality rate. This modelling approach allows us to investigate the competition of different strains and the emergence of resistant strains due to treatment.

In this work we will study the question about the emergence of new strains. From the modelling point of view, these two cases, existence of stable strains and emergence of new strains are complementary, they are not observed at the same time. The former corresponds to stable pulses while the latter to periodic travelling waves. A nonlocal reaction–diffusion equation describing the emergence of biological species was suggested in [22]. It represents a particular case of Equation (1) without immune response and genotype-dependent mortality. We will show that immune response plays an important role in the dynamics of virus quasi-species.

2. Bifurcations of Periodic Structures

An important property of nonlocal reaction–diffusion equations is that the homogeneous in-space stationary solution can lose its stability with respect to spatial perturbations leading to the emergence of periodic solutions. This property was revealed and studied in some models of population dynamics [22,39]. Here we will study it for the models of infection development described by the equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qJ\left(u\right))-uf\left(S\left(u_{\tau }\right)\right)$$

(2)

similar to Equation (1) but without the genotype-dependent mortality term. We will analyze various particular cases of this equation.

2.1. Single Nonlocal Term

We consider the reaction–diffusion equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qJ\left(u\right))-f\left(u\right)u$$

(3)

for all real x, where

$$J\left(u\right)={\int }_{\infty }^{\infty }\varphi (x-y)u(y,t)dy,$$

the function $\varphi \left(x\right)$ is bounded and non-negative. We will suppose that ${\int }_{-\infty }^{\infty }\varphi \left(x\right)dx=1$. In what follows, we set $r=q=1$. Let $u_0>0$ be a solution of the equation

$$f\left(u\right)=1-u.$$

(4)

Then $u_0$ is a stationary solution of Equation (3). We will study stability of this stationary solution.

To linearize Equation (3) about $u=u_0$, we look for its solution in the form $u=u_0+v{e}^{\lambda t}$, where v is a small perturbation and we obtain the eigenvalue problem

$$D{v}^{″}-u_0(J\left(v\right)+f^{\prime }\left(u_0\right)v)=\lambda v.$$

(5)

Applying the Fourier transform, we get

$$\lambda =-D{\xi }^2-u_0(\tilde{\varphi }\left(\xi \right)+f^{\prime }\left(u_0\right)),$$

(6)

where $\tilde{\varphi }\left(x\right)$ is the Fourier transform of the function $\varphi \left(x\right)$. If we replace $\varphi \left(x\right)$ by the $\delta$-function, instead of (6) we have

$$\lambda =-D{\xi }^2-u_0(1+f^{\prime }\left(u_0\right)).$$

(7)

Assuming that

$$1+f^{\prime }\left(u_0\right)>0,$$

(8)

we conclude from (7) that $\lambda <0$ for all real $\xi$.

Let us now analyze equality (6). Since $\tilde{\varphi }\left(0\right)=1$, then $\lambda \left(0\right)<0$. Suppose that $\varphi \left(x\right)$ is an even function, $\varphi \left(x\right)=\varphi (-x)$ for all x. Then

$$\tilde{\varphi }\left(\xi \right)={\int }_{-\infty }^{\infty }\varphi \left(x\right)cos\left(\xi x\right)dx,$$

$\tilde{\varphi }\left(0\right)=1$, and $\tilde{\varphi }\left(\xi \right)<1$ for all $\xi \ne 0$. Assuming that $\lambda =0$ in (6) for some $\xi$, we obtain the stability boundary:

$$\tilde{\varphi }\left(\xi \right)=-f^{\prime }\left(u_0\right)-D{\xi }^2/u_0.$$

(9)

This equality should be satisfied for some real $\xi$. Its value is related to the wave number of the corresponding eigenfunction. If we consider a bounded interval with periodic boundary conditions, then $\xi =2\pi k/L$, where L is the length of the interval and k = 1, 2, 3 …

2.2. Examples

Consider the functions

$${\varphi }_1\left(x\right)=\sqrt{\frac{a}{\pi }}{e}^{-a{x}^2},{\varphi }_2\left(x\right)=\frac{a}{2}{e}^{-a\left|x\right|},{\varphi }_3\left(x\right)=\left\{\begin{array}{ccc}1/\left(2N\right)& ,& \left|x\right|\le N\\ 0& ,& \left|x\right|>N\end{array}\right..$$

Then

$${\tilde{\varphi }}_1\left(\xi \right)=e^{-{\xi }^2/\left(4a\right)},{\tilde{\varphi }}_2\left(\xi \right)=\frac{a^2}{a^2+{\xi }^2},{\tilde{\varphi }}_3\left(\xi \right)=\frac{1}{\xi N}sin\left(\xi N\right).$$

Figure 1 (left) shows a graphical solution of Equation (9) for the function ${\varphi }_3\left(x\right)$. The curves corresponding to the functions in the left-hand side and in the right-hand side of this equation touch each other. The corresponding values of parameters belong to the stability boundary. For lesser values of the diffusion coefficient, the homogeneous in-space stationary solution $u_{-}$ loses its stability resulting in the emergence of a periodic in-space solution.

Let us note that Fourier transforms of the functions ${\varphi }_1\left(x\right)$ and ${\varphi }_2\left(x\right)$ are positive. Therefore, if $f^{\prime }\left(u_{-}\right)\ge 0$, then Equation (9) does not have solution, and solution $u_{-}$ is stable. If $f^{\prime }\left(u_{-}\right)<0$, it has a solution for sufficiently small values of the diffusion coefficient. Thus, emergence of periodic solutions is determined by the interaction of virus mutations, nonlocal competition for host cells and immune response. In terms of virus population distribution in the space of genotypes, these periodic solutions correspond to different virus strains.

2.3. Double Nonlocal Equation

Consider Equation (2) with two nonlocal terms $J\left(u\right)$ and $S\left(u\right)$ and without time delay ($\tau =0$). For simplicity of presentation, we set $f\left(u\right)=bu$. In the stationary case, we obtain the equation

$$D{u}^{″}+u(1-J\left(u\right)-bS\left(u\right))=0.$$

Let us recall that $J\left(u\right)={\int }_{-\infty }^{\infty }\varphi (x-y)u\left(y\right)dy$, $S\left(u\right)={\int }_{-\infty }^{\infty }\psi (x-y)u\left(y\right)dy$. Assuming that ${\int }_{-\infty }^{\infty }\varphi \left(y\right)dy={\int }_{-\infty }^{\infty }\psi \left(y\right)dy=1$, we get a homogeneous in-space stationary solution $u_{-}=1/(1+b)$ of this equation. Linearizing the equation about this stationary solution, we obtain the eigenvalue problem:

$$D{v}^{″}-u_{-}(J\left(v\right)+bS\left(v\right))=\lambda v.$$

(10)

Applying the Fourier transform to (10), we get

$$\lambda =-D{\xi }^2-u_{-}\left(\tilde{\varphi }\left(\xi \right)+b\tilde{\psi }\left(\xi \right)\right).$$

(11)

Consider the functions

$$\varphi \left(x\right)=\left\{\begin{array}{ccc}1/\left(2N_1\right)& ,& \left|x\right|\le N_1\\ 0& ,& \left|x\right|>N_1\end{array}\right.,\psi \left(x\right)=\left\{\begin{array}{ccc}1/\left(2N_2\right)& ,& \left|x\right|\le N_2\\ 0& ,& \left|x\right|>N_2\end{array}\right.$$

Then

$$\tilde{\varphi }\left(\xi \right)=\frac{1}{\xi N_1}sin\left(\xi N_1\right),\tilde{\psi }\left(\xi \right)=\frac{1}{\xi N_2}sin\left(\xi N_2\right).$$

Stability boundary is determined by the equation

$$\tilde{\varphi }\left(\xi \right)+b\tilde{\psi }\left(\xi \right)=-D{\xi }^2/u_{-}$$

(12)

obtained from Equation (11) with $\lambda =0$. An example of graphical solution of this equation is shown in Figure 1. The stability boundary corresponds to the case where the functions in the left-hand side and in the right-hand side of this equation touch each other (solid and point lines). If they intersect (dashed and point lines), then the stationary solution is unstable. For the fixed values of $N_1$, $N_2$ and b, stability conditions are determined by the diffusion coefficient.

Proposition 1.

For any positive values $N_1$, $N_2$ and b there exists a critical value $D_c$ of the diffusion coefficient such that the stationary solution $u_{-}=1/(1+b)$ is stable for $D>D_c$ and unstable for $D<D_c$.

The proof of this proposition is straightforward. It is sufficient to note that the function $\varphi \left(\xi \right)+b\tilde{\psi }\left(\xi \right)$ has negative values for any values of parameters.

Dependence of the stability conditions on $N_1$ and $N_2$ is more complex. Both can have stabilizing or destabilizing effect on the solution. In the example in Figure 1, decreasing the value $N_2$ leads to the instability of the solution (solid and dashed lines).

2.4. Delay Equation

Stability of solutions of the delay equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qu)-f\left(u_{\tau }\right)u$$

(13)

without nonlocal terms was studied in [37]. Spatial perturbations of the homogeneous in-space solution can lead to a complex spatiotemporal behavior with standing wave, travelling waves and aperiodic dynamics.

2.5. Nonlocal Delay Equation

Consider now Equation (2) with a single nonlocal term and with time delay:

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qJ\left(u\right))-uf\left(u_{\tau }\right).$$

(14)

In what follows, we set $r=q=1$. Linearizing this equation about a stationary solution $u=u_0$, we obtain the eigenvalue problem

$$D{v}^{″}-u_0\left(J\left(v\right)+f^{\prime }\left(u_0\right)e^{-\lambda \tau }v\right)=\lambda v.$$

Applying the Fourier transform, we get

$$\lambda =-D{\xi }^2-u_0\left(\tilde{\varphi }\left(\xi \right)+f^{\prime }\left(u_0\right)e^{-\lambda \tau }\right).$$

For $\lambda =0$ we obtain Equation (9) for the stability boundary. Therefore, as it can be expected, time delay does not influence the bifurcation of stationary solution. Consider now the bifurcation of time periodic solution. We set $\lambda =i\nu$, where $\nu$ is a real number. Then separating the real and imaginary part in the last equation we obtain:

$$\nu =u_0{f}^{\prime }\left(u_0\right)sin\left(\nu \tau \right),D{\xi }^2+u_0\tilde{\varphi }\left(\xi \right)+u_0{f}^{\prime }\left(u_0\right)cos\left(\nu \tau \right)=0.$$

Set $z=\nu \tau$. Then

$$cosz=-\frac{D{\xi }^2/u_0+\tilde{\varphi }\left(\xi \right)}{f^{\prime }\left(u_0\right)},\tau =\frac{z}{u_0{f}^{\prime }\left(u_0\right)sinz}.$$

(15)

We find the value of z from the first equation and the value of $\tau$ from the second equation. The first equation has a solution if and only if

$$|\frac{D{\xi }^2/u_0+\tilde{\varphi }\left(\xi \right)}{f^{\prime }\left(u_0\right)}|\le 1.$$

(16)

In the case of the local equation where $\tilde{\varphi }\left(\xi \right)=1$, the minimum of the numerator is reached at $\xi =0$. Therefore, the loss of stability occurs with the space independent perturbations assuming that $|f^{\prime }\left(u_0\right)|<1$. For the nonlocal equation the loss of stability can occur for $\xi \ne 0$ and for $|f^{\prime }\left(u_0\right)|\ge 1$.

Consider the case where the function $f\left(u\right)$ is linear, $f\left(u\right)=k_{1}u$. If $k_1>1$ and $N_1$ is sufficiently small, then the homogeneous in-space solution $u_0$ can lose its stability with respect to temporal perturbation for the time delay $\tau$ large enough. If $\tau$ is less than a critical value then temporal and spatial perturbations decay (Figure 2, left). If we increase N for the other parameters fixed, then the constant solution $u_0$ can lose its stability with respect to spatiotemporal perturbations. Figure 2 (right) shows the emergence of a spatiotemporal pattern at the center of the interval. It propagates and gradually fills the whole spatial domain. Let us now take $k_1<1$. Then time oscillations in the local problem ($N_1=0$) decay for any time delay. In the nonlocal problem and $N_1$ sufficiently large various spatiotemporal patterns can be observed (Figure A1 in Appendix A).

3. Emergence of Strains as Periodic Wave Propagation

3.1. Propagation Of Waves

We study in this section propagation of described by Equation (1) assuming for simplicity that the second integral term $S\left(u\right)$ becomes local, i.e., the kernel $\psi \left(x\right)$ is the replaced by the $\delta$-function:

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qJ\left(u\right))-uf\left(u_{\tau }\right).$$

(17)

Local and delay equations.

To study the behavior of solutions, we begin with the local case. Then we get conventional reaction–diffusion equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+F\left(u\right),$$

(18)

where the function $F\left(u\right)=ru(1-qu)-uf\left(u\right)$ can have different numbers of zeros depending on the values of parameters. Besides the zero $u_{+}=0$, there can exist up to three positive zeros, the maximal zero $u_{-}$ and possibly one or two intermediate zeros $u_1$ and $u_2$.

Monostable case. If there is only one positive zero $u_{-}$, then $F^{\prime }\left(u_{+}\right)>0$, $F^{\prime }\left(u_{-}\right)<0$. The $[u_{+},u_{-}]$-waves, i.e., the waves with the limits $u(\pm \infty )=u_{\pm }$ at infinity, exist for all values of the speed greater than or equal to some minimal speed $c_0$. These waves are stable in appropriate weighted spaces [24].

Bistable case. In the bistable case, $F^{\prime }\left(u_{\pm }\right)<0$, and there is an additional zero $u_1\in (u_{+},u_{-})$. The $[u_{+},u_{-}]$-wave exists for a single value of speed $c_1$, and this wave is globally asymptotically stable.

Monostable–bistable case. In this case, there are two intermediate zeros, $u_1,u_2$, $u_1<u_2$, and $F^{\prime }\left(u_{+}\right)>0$, $F^{\prime }\left(u_1\right)<0$, $F^{\prime }\left(u_2\right)>0$, $F^{\prime }\left(u_{-}\right)<0$. The monostable $[u_{+},u_1]$-waves, i.e., the waves with the limits $u(\pm \infty )=u_{\pm }$ at infinity, exist for all values of the speed greater than or equal to some minimal speed $c_0$. The bistable $[u_1,u_{-}]$-wave exists for a single value of speed $c_1$. If $c_1>c_0$, then there exist $[c_{+},c_{-}]$-waves for all speeds $c\in [c_0,c_1)$. If $c_1\le c_0$, then such waves do not exist, and there is a system of two waves propagating one after another with different speeds and a growing distance between them. All these properties including convergence of solutions of the Cauchy problem to waves and systems of waves can be found in [40].

Delay reaction–diffusion equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-qu)-uf\left(u_{\tau }\right)$$

(19)

is a particular case of Equation (17) where the integral $J\left(u\right)$ is replaced by u. In numerical simulations this equation is considered in a sufficiently long interval $0<x<L$ with the homogeneous Neumann boundary conditions and with some initial condition $u(x,t)=u_0\left(x\right),-\tau \le t\le 0$, where $u_0\left(x\right)=u_0$ for $0\le x\le x_0$ and $u_0\left(x\right)=0$ otherwise. Here $u_0$ and $x_0$ are some positive constants. Equation (19) was introduced in [27] to study spatial models of infection development in tissues. The influence of time delay on the wave propagation manifests itself in the most spectacular way in the monostable–bistable case where $c_0>c_1$, and there is a system of two waves propagating with different speeds. The presence of time delay can lead to the emergence of complex spatiotemporal structures between the two waves. Some examples of numerical simulations are shown in Figure 3.

Existence of waves described by this equation was proved in [41,42]. Since conventional monotonicity conditions and the maximum principle are not applicable in this case, the proof of the wave existence requires sophisticated mathematical techniques. There are only a few works where the wave existence is proved for the delay reaction–diffusion equation in the bistable case without the monotonicity condition (see also [33]).

3.1.1. Nonlocal Equation

The presence of the nonlocal term in Equation (17) can influence the regimes of wave propagation presented above for the local equation. Let us recall that the homogeneous in-space stationary solution $u_{-}$ can be stable or unstable depending on the values of parameters. In particular, for a sufficiently small diffusion coefficient or for a sufficiently large N in the definition of the kernel $\varphi \left(x\right)$:

$$\varphi \left(x\right)=\frac{1}{2N}\left\{\begin{array}{ccc}1& ,& \left|x\right|\le N\\ 0& ,& \left|x\right|>N\end{array}\right.$$

(20)

this solution loses its stability resulting in the bifurcation of a periodic in-space stationary solution. If this solution is stable, then the regimes of wave propagation are the same as before (monostable, bistable, monostable–bistable). Suppose now that it is unstable and consider, first, the monostable case. Then there are two transition: the first one is provided by the $[u_{+},u_{-}]$-wave, the second one is the transition from the constant solution $u_{-}$ to the periodic solution $u_p\left(x\right)$. If the speed $c_0$ of the former is greater than the speed $c_p$ of the latter, then they propagate one after another one with a growing distance between them (Figure 4, left). If $c_p>c_0$, then they merge and form a single periodic wave (Figure 4, right). These different regimes are also observed in the bistable and monostable–bistable cases.

Next, consider Equation (17) with time delay and, for simplicity of presentation, with a linear function $f\left(u\right)$. In this case, we have a monostable equation with two stationary points $u_{+}=0$ and $u_{-}>0$. Stability analysis of the homogeneous in-space stationary solution $u_{-}$ with respect to spatial and temporal perturbations was carried out in Section 2. If this point is stable, then we observe propagation of a usual $[u_{+},u_{-}]$-wave with a constant speed and a constant profile. However, this wave is not necessarily monotonic with respect to x, as it is the case for the local equation. Damped oscillation behind the wave occur for N sufficiently large (Figure 5, left).

Behavior of solutions becomes more complex if the solution $u_{-}$ is unstable. If the spatiotemporal perturbation of this solution propagates in space with the speed less than the speed of the $[u_{+},u_{-}]$-wave, then this wave propagates with a constant speed and a constant profile, possibly with decaying spatial oscillations behind the wave front. This wave is followed by the region of spatiotemporal oscillations (Figure 5, right). If the perturbations of the solution $u_{-}$ propagate faster than the $[u_{+},u_{-}]$-wave, then they merge and form an oscillating wave propagating with a variable speed (not shown).

Let us recall that in the monostable case the waves exist for all values of the speed greater than or equal to the minimal speed $c_0$. The value of the minimal speed is determined by the linearized problem at 0, and it does not depend on N and $\tau$ if the wave remains stable.

3.1.2. Bifurcations of Waves and Pulses

Existence and stability of pulses and waves for Equation (17) depends on the width N of the support of the kernel $\varphi \left(x\right)$.

Let us recall that the local reaction–diffusion equation with a bistable function $F\left(u\right)$ has a positive stationary solution decaying at infinity (pulse solution) if and only if $I_F={\int }_{u_{+}}^{u_{-}}F\left(u\right)du>0$. This stationary solution is unstable. It also has a stable $[u_{+},u_{-}]$-wave whose speed is positive under the same condition on the integral $I_F$. Similar properties hold for the nonlocal equation with N sufficiently small. With an increase of N, the wave becomes non-monotone as a function of x but it still has a constant speed and profile.

Periodic structures and waves appear as the width N of the support of the kernel $\varphi \left(x\right)$ exceeds a critical value $N_c^1$ for which the corresponding eigenvalue of the linearized problem crosses the origin. If we further increase the value of N, then instead of periodic waves we observe stable pulses. Thus, there are two bifurcations with a transition from simple waves to periodic waves (through an intermediate regime with two waves) and from periodic waves to pulses. The first bifurcation occurs due to the essential spectrum crossing the origin. The second one is a nonlocal bifurcation where the speed of the periodic wave decreases as N approaches a critical value $N_c^2$, and it becomes zero for N exceeding the critical value. At the same time, the spikes of the periodic wave become pulses. Let us note that multiple pulses are not stationary solutions of Equation (17), they slowly move from each other with a decaying speed.

3.2. Emergence of Strains

Virus density distribution $u(x,t)$ as a function of its genotype x and time t characterizes the existence of virus strains and their evolution. In this context, a strain is a positive localized solution of Equation (1), i.e., a solution with maximum at some $x_0$ (most frequent genotype) and rapidly decaying as the distance $|x-x_0|$ increases. Existence and stability of stationary localized solutions (pulses) of Equation (1) was studied in [23]. They correspond to persisting virus strains. In this section, we will study the emergence of new strains due to propagating of periodic waves. As it was discussed in the previous section, stable pulses and waves are mutually exclusive.

3.2.1. Initiation of Periodic Waves

Suppose that the initial virus distribution $u_0\left(x\right)$ represents a non-negative function with a narrow support at the center of the interval. For a properly chosen values of parameters the solution of Equation (1) with this initial condition develops to a periodic travelling wave. At the first stage of this dynamics, the solution rapidly growth remaining localized at the center of the interval (Figure 6a). It reaches a maximal level, and then the peak gradually decreases and becomes wider (Figure 6b,c). At some moment of time, two other peaks appear from each side of the first one. After a short transient period, they converge to approximately the same height. If the interval is sufficiently large, then other peaks will appear after some time gradually filling the whole interval (Figure 4, right). This is a typical dynamics of the initiation of periodic waves which occurs under the conditions presented in Section 2.

Let us also note that new strains (peaks) appear at some distance from the first one by a genetic jump and not as a gradual evolution of the original strain. The value of the virus density between them is close to zero. This is not the case for a large diffusion coefficient (mutation rate) where new strains appear continuously (Appendix A, Figure A2). If the diffusion coefficient is large enough, then the strains do not form, and viruses with any genotype exist. This is determined by the stability of the stationary points (Section 2).

3.2.2. the Influence of Immune Response

We consider the function of immune response in the form $f\left(u\right)=(k_{1}u+k_2)e^{-k_{3}u}$. In order to explain the influence of immune response on the emergence of strains, consider the function $F\left(u\right)=ru(1-qu-f(u)$). If $k_2=k_3=0$, then $F\left(u\right)$ has a single positive zero $u_{-}$, and this is a monostable case. For the values of $k_1$ sufficiently small, behavior of solution of Equation (1) is similar to the case where $f\left(u\right)\equiv 0$ with the propagation of a periodic wave and the emergence of new strains (Figure 6). If $k_1$ is large enough, then the equilibrium $u_{-}$ becomes stable, and there is a stationary $[u_{+},u_{-}]$-wave without emerging peaks, as it is the case of a periodic wave.

Let us recall that the growing branch of the function $f\left(u\right)$ corresponds to the antigen stimulated immune response while the decreasing characterizes death or exhaustion of immune cells due to high virus concentration. Thus, if we consider only the growing branch, then a strong immune response (large $k_1$) does not eliminate infection but prevents the formation of virus strains. Instead of the localized solutions with separated peaks, the virus density distribution converges to a constant positive solution.

In the case where the decreasing branch of the immune response function is present ($k_2=0$, $k_3\ne 0$), the function $F\left(u\right)$ can have up to three positive zeros $u_1<u_2<u_{-}$. As we discussed above, this is a monostable–bistable case where the behavior of solutions depends on the values of parameters and on the choice of initial condition. Set $u_0\left(x\right)=u_0$ for $x_1\le x\le x_2$ and $u_0\left(x\right)=0$ otherwise. If $u_0$ is small enough, then the solution represents a monostable wave (Figure 7, left) without strain formation. If $u_0$ is sufficiently large, then for the same values of parameters as before, the central peak is formed followed after some time by the appearance of two monostable waves of low amplitude (Figure 7, right).

Furthermore, the width $x_2-x_1$ of the support of the initial condition can also influence this behavior. A counterintuitive result is that increasing the support of the initial condition leads to the disappearance of the high amplitude peak and to the convergence of solution to the low amplitude monostable wave. The explanation of this effect is that two pulses (peaks) form if the support is sufficiently wide. They compete with each other, their amplitude becomes less than for a single pulse, it is not sufficient to overcome the threshold and to form a stable central pulse.

Under further increase of $k_3$, propagation of a periodic wave, as it is described above, is observed. If $k_2\ne 0$, then $f\left(0\right)>0$, i.e., immune response is nonzero even without antigen due to memory cells. Depending on the values of parameters, solution can form a stable pulse, vanish, or initiate simple or periodic waves described above.

3.2.3. Effect of the Delay of the Antiviral Immune Response

In the case of a nonlocal equation with time delay in the immune response term, spatial structures presented in the previous paragraph can become oscillating. Some examples of virus strain evolution are shown in Figure 8. The left and middle figures are obtained for the same values of parameters with different initial conditions. If the initial virus load is large enough, then there is a dominating virus strain and some other strains with low virus density and a complex spatiotemporal behavior. If the initial virus load is sufficiently small, then the dominating central strain does not exist, and there is only a variety of different genotypes with low densities. Changing the properties of the immune response (function $f\left(u\right)$), we observe stable stationary strains similar to those in Figure 4 (right image) after the initial front propagation.

3.2.4. the Influence of Genotype-Dependent Mortality

We will finish this section with the analysis of the genotype-dependent mortality on the emergence and evolution of virus quasi-species. In order to show this influence more precisely, we consider it in the case without immune response, $f\left(u\right)\equiv 0$. Set $\sigma \left(x\right)=0$ for $x_1^{*}\le x\le x_2^{*}$ and $\sigma \left(x\right)=1$ otherwise. Figure 9 (left) shows the emergence of virus strains for $x_1^{*}=0.3,x_2^{*}=0.7$ and $N=0.09$. The initial condition has a support at the center of the interval. Similar to the case of initiation of a periodic wave, there is one strain in the beginning of the simulation, and two other strain appear sometime later. These three strains fill the whole admissible interval where $\sigma \left(x\right)=0$, and new strains do not appear outside of this interval because virus mortality rate is greater there that its reproduction rate.

In the middle figure, we begin the simulation with $N=0.08$ with five emerging strains. When they become steady and do not evolve any more, we change the value of N to $N=0.09$, as in the previous simulation. However, this time we observe the regime with four strains instead of three strains observed previously. Hence, different stationary regimes can exist for the same values of parameters, and the initial condition determines the convergence to each of them.

Finally, consider the case where $N=0.2$ (Figure 9, right). As before, there is single peak of solution at the center of the interval in the beginning of the simulation. However, after some time, it disappears giving rise to two other peaks. Such behavior is determined by the size of the admissible interval: it cannot support three wide peaks, and the two of them from the sides suppress the one at the center of the interval. Thus, virus tends to fill the admissible interval in the most efficient way, that is, to maximize its total density.

4. Discussion

4.1. Virus Quasi-Species

Speciation is considered to be a general property of the living matter [43]. It manifests itself in the emergence of biological species and in a variety of other systems [44]. In the framework of mathematical modelling, speciation appears due to a non-homogeneous density distributions $u(x,t)$. In biological populations, x can be a morphological characteristic or some characterization of the genotype.

Describing virus quasi-species dynamics, we observe some similarities with the general speciation theory due to the competition for host cells but also some specific features because of the presence of the immune response and of the genotype-dependent mortality. If we consider the virus density distribution $u(x,t)$ as a function of its genotype x and of time t, then a virus quasi-species (strain or variant genome) corresponds to a localized solution with a maximal value at some genotype $x_0$ and rapidly decaying as the distance $|x-x_0|$ increases. In the case of persistent strains, the most frequent genotype $x_0$ is fixed but it can be also time-dependent for the evolving strains.

Existence of virus strains considered to be a positive stationary density distribution decaying at infinity (in the approximation of an infinite genotype space) is not a priori given. In the previous work [23] we revealed two mechanisms providing the existence and stability of such solutions. In the first case, the existence of virus strains is determined by the genotype-dependent mortality where the virus can survive only inside some viability interval. The maximum of the density distribution is achieved in the middle of the corresponding viability interval. The second mechanism is determined by the immune response under the assumption that the immune response function $f\left(u\right)$ decreases for large u. In this case, the virus can survive and form a persistent strain if its concentration is sufficiently high, and if the competition with other strains for host cells occurs in a sufficiently wide range of genotypes.

In mathematical terms, virus strains correspond to stable pulse solutions of the corresponding nonlocal reaction–diffusion equations. Existence of stable pulses does not occur for the conventional local equations.

4.2. Emergence of New Quasi-Species: Summary of the Results

It is important to note that stable pulses and periodic travelling waves are mutually exclusive, they are not observed for the same values of parameters. In this work we study periodic travelling waves. Emergence of new peaks in the virus density distribution during the wave propagation corresponds to the emergence of new virus strains.

From the mathematical point of view, the conditions of the emergence of periodic travelling waves can be determined by the linear stability analysis of the homogeneous in-space stationary solutions. In order to show the influence of different factors on the stability conditions, this analysis is carried out in Section 2 for the single nonlocal term, for both nonlocal terms and for the nonlocal delay equation. In the presence of a single nonlocal term, periodic spatial structures bifurcate from the constant solution if the diffusion coefficient D is sufficiently small and if the kernel $\varphi \left(x\right)$ of the integral satisfies certain conditions. In particular, for a piece-wise constant kernel, its support should be sufficiently large.

In the case of two nonlocal terms, the qualitative behavior of solutions is similar. The interaction of the nonlocal terms can have stabilizing or destabilizing effect depending on the values of parameters. The influence of time delay (and a single nonlocal term) on the stability conditions depends on the immune response function resulting in temporal or spatiotemporal oscillations.

The transition between a virus-free equilibrium and an infected equilibrium is provided by travelling waves (Section 3). In the case without nonlocal terms or time delay, this is a conventional wave with a constant speed and profile or two waves propagating one after another with different speeds. The nonlocal terms can result in the emergence of periodic waves, while time delay can lead to a complex spatiotemporal pattern formation.

Let us note that the qualitative behavior of solutions is quite robust, and it is not very sensitive to the particular choice of the immune response function and of the integral kernels. The range of genotypes is supposed to be sufficiently large to neglect the influence of the boundaries. Mathematically, we can consider the whole real axis. In numerical simulations, we consider a bounded sufficiently large interval. In most cases, we stop the simulations before the wave approaches the boundary of the interval. In some cases (Figure 8), we continue the simulations to reveal spatiotemporal pattern formation inside the interval. In this case, periodic boundary conditions are more convenient because they do not influence the behavior of the solution inside the interval. Dirichlet or Neumann boundary conditions can influence the behavior of solutions. Thus, periodic boundary conditions do not have here biological significance, but they are more appropriate for mathematical modelling.

4.3. Biological Interpretations

Suppose that the initial viral load is localized in a narrow interval of genotypes (e.g., a founder virus in the HIV case). Due to its multiplication and mutations, the density distribution grows and widens. Viruses with different genotype begin to compete for the host cells leading to the appearance of new strains at some distance to avoid the competition with the existent strains. This description of speciation is generic [22], it is not specific for virus quasi-species. Specific features of virus diversification and strain emergence are related to the immune response.

If we take into account the cross-reactivity in the immune response, i.e., different antigens can stimulate the same immune cells, then the immune response interferes with virus competition for the host cells. This interaction is quite complex and can act in different ways on different strains. However, if the mutation rate (diffusion in the genotype space) is sufficiently small, then the speciation of virus quasi-species will necessarily occur. The critical conditions leading to the emergence of new strains depend on the parameters of the problem, and the immune response can have both stabilizing and destabilizing effect.

The influence of immune response becomes easier to determine if we neglect cross-reactivity. Assuming that the immune response function $f\left(u\right)$ is increasing due to the stimulation of immune response by the virus antigens, the model predicts that immune response acts to suppress the formation of new strains. If the growth rate of the function $f\left(u\right)$ is sufficiently high, then the speciation of the virus density distribution completely disappears. In this case, instead of a discrete set of virus strains our study predicts that a uniform density distribution as a function of virus genotype will take place.

The situation becomes again more complex if we consider also a decreasing branch of the immune response function, which appears due virus-induced death of immune cells, and time delay in the immune response required for the clonal expansion of immune cells. In this case, along with periodic travelling waves described above, complex nonlinear dynamics of solutions can take place with various patterns of emerging and disappearing strains.

Genotype-dependent virus mortality restrains the evolution of virus species to the admissible interval. The emergence and the evolution of virus strains within the viability interval depend on the values of parameters and on the initial viral load. Let us note that different virus density distributions can be observed for the same values of parameters.

Concluding this discussion, we point out the limitations of the model considered in this work due to various simplifying assumptions. We do not take into account the presence of different immune cells and of cytokines participating in the immune response, complex intracellular regulation of cell fate and of virus multiplication. On the other hand, these and other simplifications allow us to reveal some generic properties of the evolution of virus quasi-species, which would be more difficult to identify in a more complex model. This modelling framework provides a starting basis for further investigations and for the introduction of more detailed models.