Introduction

Mathematical model

To formulate the equation of the virus dynamics, we follow the view that antigen dose, time period during which it is available and its “geographical” distribution within this host influence the duration and extent of the immune responses [38]. The underlying regulation implies a bell-shaped relationship between viral load and the magnitude of the antiviral T cell response so that high antigen amount leads to exhaustion of T cells [39, 40]. The biological scheme of the model and the relevant processes are shown in Fig 1(a). To make the model analytically tractable, we do not consider explicitly the population dynamics of the antiviral T cell response. The magnitude of the clonal expansion is assumed to be a non-linear function of the viral load with a time lag. The common practice of using time delay in the life sciences is when there are some hidden variables and processes which are not well understood but are known to take some time to react. Although the model is general enough, we develop it with a particular interest to study non-, or weakly cytopathic infections such as LCMV, SIV and HIV-1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Schematic representation of virus infection dynamics regulation (left) and qualitative forms of the function f(v) (right).

Low level infection stimulates immune response while high level infection down-regulates it. The former corresponds to the growing branch of the function f(v) while the latter to its decreasing branch.

https://doi.org/10.1371/journal.pone.0168576.g001

We consider the following equation for the concentration of virus v in lymphoid organs such as the spleen or lymph nodes or in any other permissive tissues: (1) Here v = v(x, t), v_τ = v(x, t − τ), the function f(v_τ) will be specified below. The first term in the right-hand side of this equation describes virus diffusion, the second term its production and the last term its elimination by CTL cells. The parameter D is the diffusion coefficient (or diffusivity) and k stands for the replication rate constant. The parameterized function f(v_τ) characterizes the virus induced clonal expansion of T cells, i.e. the number and function of these cells generated after some time-delay τ, i.e. it depends on the concentration of virus some time before. The qualitative form of this function is shown in Fig 1(b) [41].

Eq (1) is a quasi-linear functional reaction-diffusion equation. The fundamentals of the general theory of these equations can be found in [42, 43]. It is worthwhile to note that for t ∈ [t₀, t₀ + τ] each solution of Eq (1) satisfies the usual non-autonomous reaction-diffusion equation (2) with some boundary conditions and with the initial condition v(t₀, x) = v₀(t₀, x), where v₀(s, x), s ∈ [t₀ − τ, t₀] is the prescribed initial function for Eq (1) and a(t, x) = f(v₀(t − τ, x)). After solving Eq (2) on [t₀, t₀ + τ] we can again apply the same argument to (1.1) for t ∈ [t₀ + τ, t₀ + 2τ] and so on. This procedure is known as the method of steps [44], it is simple and efficient but it works only on finite time intervals. Consequently, other approaches should be applied in order to study global in time behaviour of solutions to Eq (1). One can expect (and this is confirmed analytically, cf. [44]) that the small delays are “harmless” in the sense that the asymptotic properties of solutions of delay differential equations varies continuously as τ → 0+. However, if τ becomes large, the delayed nonlinearity might trigger essential qualitative changes in the dynamics of the system. It also presents a considerable technical complication from the mathematical point of view.

In this respect, the delayed reaction-diffusion equations whose reaction term g is either of logistic type (i.e. g is as in (1.1), when v is not separated multiplicatively from v_τ) or of the Mackey-Glass type (when v is separated multiplicatively from v_τ, i.e. g = −kv + b(v_τ)) are between the most studied ones. It is an interesting point of discussion whether the Mackey-Glass type models reflect more adequately the biological reality than the logistic models, see e.g. [45] and especially [44] for further details. In any case, in many regards both models exhibit similar types of qualitative behaviour of solutions. In particular, this can be seen from the existence, uniqueness and stability results describing travelling waves (both monostable and bistable) for various subclasses of the logistic type and the Mackey-Glass type equations, see [1, 46–52] among many other references. We note that the investigation of delayed logistic models is more difficult and technically involved than the studies of the Mackey-Glass type systems, precisely because of the multiplicative non-separateness of v and v_τ. For example, so far no analytical results on the existence and uniqueness of monostable (or even bistable, cf. [53–58]) waves in the delayed Eq (1) with non-monotone nonlinear response f were available in the literature. This work presents the first contribution in this respect: below, we propose a novel approach to the both aforementioned problems. This approach combines the phase plane method with simple fixed points arguments.

We will also study nonlinear dynamics of reaction-diffusion waves described by Eq (1). Travelling waves in populations dynamics are extensively studied in ecology and epidemiology (see [59–61] and the references therein). In the case of reaction-diffusion systems they can be accompanied by pattern formation and spreading [59, 62–65]. Pattern formation and wave propagation in the case of delay reaction-diffusion equations are discussed in [66–68].

Virus spread without time delay

The above model is used to identify and characterize some fundamental types and patterns of the spatiotemporal dynamics of virus infections in tissues. We start the analysis of the model with a simplifying assumption of τ = 0. In this case, virus distribution in the tissues such as spleen or lymph node can be described by the reaction-diffusion equation (3) The above model governs a spatio-temporal dynamics of the virus infection resulting from a balance of the diffusion, replication, and elimination by the cytotoxic T cell response. The T cell response depends in a bell-shaped way on the viral load. In contrast to Eq 1, the time needed for the cells to expand and migrate to the site of infection is assumed to be taken into account indirectly via the characteristics of the function f(⋅) rather that explicit time delay τ. Here v = v(x, t) is the dimensionless normalized virus concentration, x ∈ R₊, t ∈ [0, +∞). By a re-scaling of variables this equation can be reduced to the same equation with D = 1, k = 1. We will consider this equation on the whole axis, i.e., in Ω = ((x, t): −∞ < x < +∞, 0 ≤ t < +∞) with a non-negative initial condition u(x, 0) = u₀(x).

For what follows it is convenient to introduce the function (4) Depending on the form of the function f(v) we will get different behavior of solutions. Let us note that F(0) = 0. This function can have other zeros for v > 0. We will also suppose that the function f(v) is non-negative, it is continuous together with its second derivatives.

Under the above assumptions, this model of virus infection implies a natural description of the spatiotemporal dynamics of the virus in terms of travelling- or standing waves. The travelling wave represents a pattern in the viral density that travels across the tissue in an uninterrupted fashion. Another possibility represented by a standing wave would be a viral density pattern that appear to be confined to some region and standing still. Note that standing waves are the partial case of traveling waves correspond to a zero velocity of the wave front. Hence, we do not consider them as a separate case.

We will study the conditions for the existence and properties of travelling wave solution of Eq (3), that is solutions of the form v(x, t) = w(x − ct) that moves at a constant speed c in the positive x-direction. The function w(x) (we keep the same notation x for the argument of w(⋅)) satisfies the equation (5) with appropriate conditions at infinity equal one of the steady state, i.e., lim_{x → −∞} w(x) = v_i, lim_{x → +∞} w(x) = v_j i, j = 0, 1, 2, 3, i ≠ j. The constant c here is the wave speed and its relationship to original model parameters will be established by the analytical treatment. The value of c is a solution to the above eigenvalue problem.

As shown in Fig 2, the model admits the existence of three types of regimes: (i) a small amplitude wave corresponding to a low-level infection spreading at constant c₀ through the tissue, (ii) a large amplitude wave corresponding to a high-level infection spreading at constant c₁ through the tissue, and (iii) a large amplitude two-wavefront solution corresponding to a low-level infection followed by a hight level one evolving at different wave speeds c₁ and c₀, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Schematic representation of the spatial patterns of virus infection dynamics as travelling waves.

Typical wavefront solutions taking a steady state value v_i at the left end and another steady state v_j, i ≠ j at the right end. The travelling waves evolve with the speed c. A qualitative relationship between the initial viral load and the emerging pattern on the spatiotemporal pattern of virus spread is sketched.

https://doi.org/10.1371/journal.pone.0168576.g002

Model problems with time delay

A qualitative analysis of the behaviour of solutions of the delay reaction-diffusion model for a general bell-shaped immune response function f(⋅) presents a formidable challenge. We proceed using simplifying approximations to it.

In this section we consider Eq (3) with a piece-wise constant function f(v): (23) where a, b ≥ 0 and w₀ > 0 are some constants. If a < b, then it is an approximation of the growing branch of the function f(v) in Fig 1, if a > b, of its decreasing branch. We look for conditions of the existence of a travelling wave solution v(x, t) = w(x − ct), which satisfies the equation (24) with the following boundary conditions (25) We further assume that w₀ < 1 − b. We look for a solution of problem Eqs (24) and (25). Since it is invariant with respect to translation in space, we suppose that the discontinuity of the function f(w(x + cτ)) occurs for x = 0, that is w(cτ) = w₀. Hence instead of Eq (23) we can write (26) These equations should be completed by the relations (27) and by condition (25). Since eq (26) do not contain delay, we will use the phase space analysis to find the trajectories which satisfy condition (27).

Bistable case

We begin with the case where 0 < b < 1 < a in which the piece-wise constant function (23) approximates decreasing branch of the function f(w) in Fig 1 where immune response decreases with the increase of infection level. Under these assumptions the stationary points w = 0 and w = 1 − b of the corresponding equation dw/dt = F(w) are stable. Instead of the second-order equations in Eq (26) we will consider two systems of the first-order equations: (28) Analysis of the trajectories of these systems allows one to formulate the following statement.

Proposition 1. Let a > 1, 0 < w₀ < 1 − b and (29) Then for any τ > 0 there is a unique positive value of c for which there exists a solution of problem Eqs (25)–(27).

The proof of this statement is given in (S1 Appendix). Condition (29) ensures the positiveness of the wave speed c. It can be shown that the wave speed decreases if τ increases. To this end we consider the following approximate model obtained via linearization of system (26) about the boundary values at infinity: (30) This approximating system admits an analytical solution and allows one to find the speed of infection wave propagation as a function of the model parameters. We have from Eq (30): (31) where (32) From condition (27) we get the following equations: (33) From these equations we obtain the following relationship between the wave speed and time delay in immune response. We can write it as an explicit analytical dependence of τ on c: (34) Fig 3 shows an example of wave propagation and comparison of the analytical and numerical results for the wave speed (see S2 Appendix for details). The wave speed decreases with time delay. From the biological point of view, this means that infection wave propagates slower if the delay in the clonal expansion of the antiviral immune response gets longer.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Wave propagation in the bistable case (left) for Eq (1) with function f(v) given by formula (23) (a = 1.1, b = 0.1, w₀ = 0.1, D = 10⁻⁴).

The curves show the function v(x, t) at successive moments of time. The wave speed dependence on time delay (right). The lower curve shows the results of numerical simulations, the upper curve is the analytical approximation by formula (34).

https://doi.org/10.1371/journal.pone.0168576.g003

Monostable case

As above, we consider problem Eqs (24) and (25) with function (23) and reduce it to problem Eqs (25)–(27). We will assume in this section that 0 < a < b < 1 and 0 < w₀ < 1 − b. Then the first system in Eq (28) has two stationary points for 0 ≤ w ≤ 1, (0, 0) is a a node or a focus, depending on the value of c, and (1 − b, 0) is a saddle point. The second system in Eq (28) has only one stationary point (0, 0), which is a node or focus. We suppose that . Then the point (0, 0) is a stable node for both systems.

Proposition 2. Suppose that 0 < a < b < 1 and 0 < w₀ < 1 − b. Then for any and τ > 0 there exists a solution of problem Eqs (24) and (25) with function (23).

The proof of this proposition is given in (S1 Appendix). It can be verified that the solution becomes non-monotone for τ sufficiently large.

We now provide an analytical approximation to the solution of the system. Solution of eq (30) changes if a < 1 since there are two bounded exponentials for x > 0: (35) where (36) Representation of solution Eq (35) holds if μ₁ ≠ μ₂. From condition (27) we get the following equations: (37) Expressing k₁ from the first equation, we obtain the system of two equations with respect to k₂ and k₃: (38) Since μ₁ > μ₂ > 0, λ > 0, τ > 0, then, assuming that c > 0, we have (39) Therefore the determinant of system (38) is positive, and it has unique values of k₂, k₃ for any c and τ.

We obtain below an analytical representation of the approximate solution for the minimal speed c₀. Instead of Eq (35) we now have: (40) where μ = c/2, and from Eq (27) it follows that (41) Hence we get (42) which suggest that the solution is increasing for x < 0 if k₁ < 0. Therefore the travelling wave is non-monotone for τ sufficiently large (Fig 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Propagating wave in the monostable case is monotone for small time delay (τ = 1, left) and non-monotone for large time delay (τ = 8, right).

The values of other parameters are a = 0.1, b = 0.3, w₀ = 0.1, D = 10⁻⁴.

https://doi.org/10.1371/journal.pone.0168576.g004

Solution w(x) is positive and monotonically decreasing for x > 0 if μ₁ and μ₂ are real and positive. This condition is satisfied for . Hence we consider the values of speed greater than or equal to the minimal speed . Numerical simulations with the reaction-diffusion model presented in Fig 4 show that, similar to the case without delay, solution v(x, t) with the initial condition v(x, 0) such that v(x, 0) ≥ 0 for all x, v(x, 0)¬ ≡ 0 and v(x, 0)≡0 for x sufficiently large, approaches the travelling wave with the minimal speed. The details of the numerical solution of the delay reaction-diffusion model are given in S2 Appendix.

Time oscillations

In this section we analyse the emergence of travelling wave solutions with a periodic structure established after the wave propagation. We do not assume here that the function f(v) is given by equality Eq (23). Let the following equality hold f(v₀) = 1 − v₀ for some v₀, and f′(v₀) > 0. Then v(x, t) = v₀, x ∈ (−∞, +∞), t ∈ (−∞, +∞) is a stationary solution of Eq (1). In order to study its local asymptotic stability with respect to small perturbations, we look for the solution of this equation in the form (43) where a, ϵ are real numbers, ϵ is a small parameter, and λ is an eigenvalue. Substituting the above function in Eq (1) and equating the terms with the first power of ϵ, we get the following characteristic equation: (44) (we do not assume here that D = 1). The stability boundary of the steady state solution can be computed by considering the characteristic roots in the form of purely imaginary eigenvalues λ = iϕ. Then we have (45) Therefore the following equalities must hold for the real and imaginary parts: (46) (47) Set z = ϕτ. Then from Eqs (46) and (47) we obtain (48) The first equation has a solution if the right-hand side is greater than −1. In particular, if D = 0, then the condition reduces to f′(v₀) ≥ 1. Using z determined from the first equation, we find τ from the second equation.

Proposition 3. If f′(v₀) > 1 + Da²/v₀, then for all τ > z/(v₀ f′(v₀)sin z) the solution v₀ of Eq (1) is unstable. Here z is determined from the first equation in Eq (48).

Our analysis suggests that if the initial condition v(x, 0) does not depend on the space variable, then the solution v(x, t) of Eq (1) is also homogeneous in space. In this case, depending on the values of parameters, the solution of the model either convergence to the stationary solution v₀ or to stable periodic time oscillations both being spatially homogeneous.

However, this behavior can be different in the case of travelling wave propagation. If we fulfil the linear stability analysis of the homogeneous in space stationary solution v₀ in the moving coordinate frame attached to the wave, we obtain the same stability conditions as before. It follows from the first equation in Eq (48) that the onset of stability depends on the wave number a of the spatial perturbation and on the diffusion coefficient. The steady-state solution v₀ appears to be more stable with respect to spatially non-uniform perturbations (a ≠ 0) than with respect to perturbations which are homogeneous in space (a = 0). The frequency of the spatial perturbations depends on the ratio of wave speed and the frequency of the time oscillations, a = ϕ/c.

Hence the model predicts the existence of three spatiotemporal regimes of travelling wave propagation summarized in Fig 5 in the case of the linear function f(v) = rv. In the first one, both types of perturbations, i.e., the spatial perturbations and time perturbations homogeneous in space, decay with time. They manifest themselves as decaying spatial oscillations behind the wave front followed by a spatially constant solution (Fig 5, left). Another regime of wave propagation is characterized by the decaying spatially heterogeneous perturbation and the persisting homogeneous in space perturbations (Fig 5, middle).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Numerical simulations of eq (1) with the function f(v) = rv (r = 2, D = 10⁻⁴).

Wave propagation for three different values of time delay, τ = 1.4, 2, 4, respectively. For small time delay (left) space and time oscillations decay, for intermediate time delay (middle) space oscillations decay while time oscillations persists, for large time delay (right) both of them persist.

https://doi.org/10.1371/journal.pone.0168576.g005

Finally, if the spatial perturbations persist, then the travelling wave propagation takes place with a (moving) periodic structure established behind the wave front (Fig 5, right). Note that the wave speed in all these cases equals c₀.

Full-scale viral regulation of immune response

We proceed with the analysis of the spatiotemporal pattern formation for model 1. In the previous section we studied analytically and numerically three distinct reduced complexity models that correspond to various specific modes of immune response regulation, i.e. activation, suppression. The first one models the virus-driven activation of the immune response with a growing function f(v), the second variant corresponds to the down-regulation of immune response with a decreasing function f(v). In both of the above cases we considered step functions in order to get some analytical results about the generic behavior of solutions. The third model problem considers a growing function f(v) similar to the first version but with a continuous function for which we analyzed emergence of time oscillations. We consider now a full-scale viral regulation of the immune response described by Eq (1), i.e. the whole function f(v), v ≥ 0 as shown in Fig 1. The corresponding behavior of solutions will be assessed using numerical numerical simulations and insight gained via the analyses of the simplified model problems.

Let us recall that equation f(v) = 1 − v has three roots (solutions), v₁, v₂ and v₃. In the model without delay (Section “Virus spread without time delay”), Eq (3) has wave solutions with the limits w(−∞) = v₁, w(∞) = 0 for all values of the speed greater than or equal to the minimal speed . Similar result holds for the model problem in the monostable case (Section “Model problems with time delay”) when the delay τ > 0. Suppose that the stationary point v₁ is stable (cf. Section “Time oscillations”), and the initial condition has limits 0 and v₁ at ±∞. Numerical simulations show that solution of Eq (1) converges to the wave with the minimal speed. If the time delay is sufficiently large, then the wave is not monotone, similar to the model problems considered in Sections “Monostable case” and “Time oscillations”.

The travelling wave with the limits w(−∞) = v₃, w(∞) = v₁ corresponds to the bistable case. It exists for a unique value c₁ of wavefront speed in the case without delay τ = 0 and for the model problem with non-zero delay τ > 0.

Let us consider the case where c₁ < c₀. Then in the case without delay (τ = 0) the wave with the limits w(−∞) = v₃, w(∞) = 0 does not exist. The behavior of solutions of Eq (1) is described by a system of two waves propagating one after another with the speeds c₀ and c₁. A similar behavior is observed when the delay is present in the regulation of the immune response but is sufficiently small (Fig 6). The solution can be monotone behind the second wave or non-monotone at all depending on the value of τ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Numerical simulations of different regimes of infection spreading depending on time delay, τ = 0.4, 0.95, 1.5, 10;D = 0.0001.

For small time delay (two left figures: τ = 0.4, 0.95), there are two consecutive waves of infection propagating with different speeds. The first wave can be non-monotone. For large time delay (two right figures: τ = 1.5, 10), the second wave propagates faster and they finally merge forming a single wave which can be either monotone or non-monotone.

https://doi.org/10.1371/journal.pone.0168576.g006

Contrary to the case considered in Section “Bistable case”, the speed of the bistable traveling wave increases with time delay. This increase occurs because the bistable wave is preceded by the a monostable wave characterized by a lower level of infection. Therefore the values of the function f(v) at the bistable wave also decrease due to time delay. Hence the speed c₁ increases as a function of τ, while the speed c₀ of the monostable wave does not depend on τ. For τ sufficiently large the two waves merge and form a single wave with the limits w(−∞) = v₃, w(∞) = 0 (Fig 6). This effect does not exist in the model problems considered in Section “Model problems with time delay”. It is specific for the system of waves where the bistable wave follows the monostable one. If τ is large enough, this resulting wave becomes non-monotone.

If f′(v₁) is sufficiently large, i.e., the sensitivity of the T cell activation to the antigen is high then the stationary point v₁ becomes unstable (Section “Time oscillations”). The spatiotemporal regimes of infection spreading show a broad spectrum of dynamic patterns. The monostable wave becomes non-monotone with decaying or persisting oscillations behind it (Fig 7). One type of patterns is characterized by a transition zone between decaying space oscillations and the bistable wave with perturbed time oscillations of the homogeneous solution v₁. Space oscillations become more complex for larger values of τ which represent a second type of the spatial dynamics. If time delay is sufficiently large, then the two travelling waves merge, as before, forming a single stable non-monotone wave.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Spatiotemporal regimes of infection spreading.

The monostable wave becomes non-monotone with decaying or persisting oscillations behind it. The type of patterns on the left is characterized by a transition zone between decaying space oscillations and the bistable wave with perturbed time oscillations of the homogeneous solution v₁. Space oscillations become more complex for larger values of τ which represent a second type of the spatial dynamics. If time delay is sufficiently large, then the two travelling waves merge, as before, forming a single stable non-monotone wave. The values of time delay are, respectively, τ = 0.7, 1, 1.5, 2;D = 0.0001.

https://doi.org/10.1371/journal.pone.0168576.g007

Fig 8 shows the last series of simulations in which the slope (sensitivity) of f′(v₁) is larger than in the above set of simulations and the instability of the steady state solution v₁ is more pronounced. In this case, we observe the existence of a monostable wave with spatial oscillations behind it. This wave can be separated from the bistable wave by a zone of irregular oscillations. Further increase of the delay results in a qualitative change of the spatial patterns of the infection spread. The two travelling waves do not merge any longer and the monostable wave is not followed by steady space oscillations. Aperiodic oscillations are observed behind the wave front which propagates, as before, at a constant speed c₀.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Existence of a monostable wave with spatial oscillations behind it.

This wave is separated from the bistable wave by a zone of irregular oscillations. Increase of the delay value results in a qualitative change of the spatial patterns of the infection spread. The two travelling waves do not merge and the monostable wave is not followed by steady space oscillations. Aperiodic oscillations are observed behind the wave front which propagates at a speed c₀. The values of time delay are, respectively, τ = 1, 2, 3, 4;D = 0.0001.

https://doi.org/10.1371/journal.pone.0168576.g008

We summarize the analysis of the existence of various dynamical regimes of infection wave propagation on the parameter plane representing the strength of the immune response, characterized by parameter f_m and the time delay, representing the duration of the clonal expansion of the virus-specific T cells in Fig 9. For small delays and weak immune responses the infection spreads as a stationary wave (region 1, left low domain). An increase in the strength of the immune response leads to the emergence of two consecutive waves propagating with different speeds (region 2). Strong immune responses cause the infection spread in the form of two waves with spatiotemporal patterns between them for delays which are not too large (region 3). However, if the time needed for the immune response to develop exceeds a certain threshold (region 1) then the infection propagates, as before, as a stationary wave.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Different regimes of infection wave propagation on the parameter plane (f_m is the maximum of the function f(v), and τ is the time needed for the development of immune response), see Supplementary materials: 1—stationary wave propagation (cf. Fig 6, two right images, Fig 7, right image), 2—two consecutive waves with different speeds (cf. Fig 6, two left images), 3—two waves with spatiotemporal patterns between them (cf. Fig 7, three left images and Fig 8).

https://doi.org/10.1371/journal.pone.0168576.g009

Discussion

Viral propagation in host organisms during infection includes two processes, the local spreading inside permissive tissues and organs, and the global propagation between different organs. The target organs for HIV and SIV infections are lymphoid tissues such as GALT, lymph nodes, spleen, etc., whereas most other viruses replicate outside organized lymphoid organs, e.g., the hepatitis viruses (B and C) primarily infect and spread in the liver cells. In this study we suggest a novel mathematical model for the spatiotemporal dynamics of virus infection in the face of an antiviral immune response. The model is formulated with a 1D reaction-diffusion equation for the spatially distributed viral population. The reaction term considers the immune response regulated by a bell-shaped dose-response relationship with viral load. It takes into account the time-lag between the activation of the virus-specific T cells and the elimination of the virus by clonally expanded CTLs. The model was used to gain insights into plausible regimes of the viral spread in tissues. To this end a combination of analytical studies of the simplified versions of the model and numerical simulations were conducted to elucidate the role of key infection parameters such as dose, diffusion coefficient characterizing the efficacy of the cell-to-cell or free virus-to-cell spread, and delay in the establishment of the antiviral immune response.

Our mathematical model is based upon the following general phenomenology of virus infections. Viruses reproduce inside host cells and transport to new cells either by direct cell-cell contact or through the extracellular matrix. Infected cells are killed by CTLs. Thus, infection spreading in the tissue is determined by virus multiplication, virus transport and immune responses. The intensity of the immune response depends on the infection level. It has a specific bell shape as shown in Fig 1. The form of the function f(v) shows that the immune response intensifies with the infection level up to some critical value where this function attains its maximum. However, for greater virus levels, it decreases because viruses affect the cells that participate in immune response. Moreover, the immune response is not instantaneous, but requires some time for cell proliferation and maturation. Hence, in order to describe virus multiplication, transport and death we consider reaction-diffusion equations with a time delay. Further extensions of the model include the shift towards larger spatial dimensions, i.e, 2D- or 3D spatial models of viral spread and explicit consideration of the spatial-temporal dynamics of the immune responses in tissues [71, 72]. As the number of studies in the field is still rather limited as outlined in S1 Table, there is a number of challenges in terms of the modelling approaches and their computational treatment.

Although the proposed mathematical model is a one-dimensional non-linear PDE model with a time lag, the behavior of the model solutions is already quite complex. It depends on the intensity of immune response (function f(v)), the initial virus load (initial condition) and the growth rate of the virus in accordance with the existing view of the regulation of virus infections [1, 2, 4, 11, 29, 33, 34, 38]. We assume that the maximal intensity of the immune response is sufficiently high to overcome virus growth and spread in a low-dose infection. Under this assumption, the reaction-term equation f(v) = 1 − v has three positive solutions v₁ < v₂ < v₃. The steady state solution v₁ can be either asymptotically stable or unstable, i.e. exhibit an oscillatory regime, the steady state v₂ is always unstable, v₃ is always asymptotically stable.

The model predicts the existence of a diverse spectrum of spatiotemporal patterns of virus dynamics. The initial infection can convert the virus free steady-state v = 0 into one of the following infection states: a low level infection v₁, a high level infection v₃ or in oscillatory regime.

1. If the initial virus load is sufficiently low, then the model predicts the existence of a mono-stable wave propagating with the speed c₀ determined by the diffusion coefficient D and by the immune response f′(0). The wave propagation depends on the value of the time delay τ. The wave front can be followed by either a constant low level infection v₁ or a spatially periodic viral distribution moving through the tissue with the same speed c₀.
2. If the initial virus load is sufficiently high, then the regimes of travelling wave propagation depend on the speed of the bistable wave c₁. If it is negative, then the initial infection evolves to a low level spatially homogeneous infection. If this speed is positive but less than c₀, then the infection spreads as a combination of two travelling waves. The first one, a low level travelling wave is followed by a transitory region of spatially heterogeneous or irregular infection finally replaced by a high level spatially homogeneous infection.
3. Finally, if c₁ is sufficiently high, then there is only one travelling wave establishing a high infection level v₃. This transition can be non-monotone as far as the viral distribution is concerned. Note that the speed of the bistable wave c₁ increases with the increase of the time delay of the immune response.

The predicted patterns of virus infection spread admit the following biological interpretation. The traveling wave solution corresponds to an infection of the target organ spreading at some tempo resulting in a homogeneous or spatially periodic distribution of the virus in the organ. Depending on the parameters of the immune response and the virus diffusion, oscillatory dynamics can emerge indicating that the virus distribution in the organ can be spatially non-homogeneous and time-dependent. This could manifest as oscillations of the virus density with the amplitude different in different parts of the target organ. It should be noted, that oscillatory dynamics of the virus infection is considered to provide an evolutionary advantage to the virus by allowing it to evade the control by the immune response and to establish virus persistence [73]. In addition, the alternating phases of high- and low-virus replication levels allow the target cells to recover consumed resources and to circumvent cell death. Overall, the oscillatory dynamics seems to increase the robustness of virus persistence in face of the immune response. Virus distribution in lymph nodes can be highly inhomogeneous even if the tissue itself is homogeneous. A clear example is provided by HIV infection. Clinical and experimental studies of HIV and SIV infections indicate that the virus growth occurs in multiple local bursts reflecting local non-equilibrium interactions between HIV and immune activated cells [31]. It was further discussed in [32] that in chronic HIV infection efficient transmission of the virus is limited to microscopic clusters of T cell in lymphoid tissues suggesting that the continuity of virus production is a result of spatially separated bursts. Mathematically, this would correspond to an oscillatory or irregular spatially inhomogeneous distribution of virus revealed in our study.

Spatial aspects of virus replication and the impact of the immune response on virus dynamics were thoroughly studied for SIV infection of monkeys [10]. Viral replication in skin patches of monkeys was clearly shown to be spatially heterogeneous (Fig 2c and 2d) and dependent on the time needed for CTLs to reach the site of local virus replication (i.e., the delay between activation of antigen-specific CTLs and in situ infiltration). It was proposed that the spatial spreading of the virus underlies the dynamical escape from the immune response. A series of latest studies on HIV visualized in cerebellum and lymph node tissues [74, 75] showed different patterns of virus distribution. HIV infected cells were clustered in cerebellum but were dispersed in lymph node tissues. The studies conclude that a high variability in the amount of HIV in tissues within an individual deserves further experimental analysis with a special attention to precise location of tissue sampling as factor affecting the amount of HIV recovered (non-homogeneity). In vivo studies of virus distribution and spreading in target tissues are complemented by the analysis of the impact of virus infection parameters on the infection spread between cells in in vitro systems [76].

The above findings concerning the spatiotemporal dynamics of the virus spread have direct implications for experimental and clinical studies of the biology of infections. The invention of imaging and visualization technologies provided a novel type of spatial information on the infection-host interaction process [7, 11, 77–79]. To interpret the data and gain a quantitative understanding of the determinants of the infection process from the image-derived data, the application of mathematical models is needed [80]. Indeed, the spatial patterns of pathogen dynamics can be complex and far from being uniform [7], and require the invention of appropriate concepts, patterns and parameters for their analysis. The major work in this direction remains to be done. Our analysis suggests that viral spreading in tissues can proceed as a travelling wave (characterized by either a low- or a high viral load), spatially periodic or irregular oscillations, and a combination of those. The evidence for the existence of spatially heterogeneous patterns, called multifocal infection, of virus spread in HIV infection in lymphoid tissue has been published recently [7]. The images of the tissues with fluorescently labelled HIV RNA [7] allow a range of interpretations depending on the choice of the cross-section: a bistable travelling wave (see Fig 1A), monostable low-level infection Fig 1B, or oscillatory mode (Fig 1D).

The possibility that the infection spreads as waves differing in their amplitude suggests the existence of spatiotemporal mechanisms of disease pathogenesis. Indeed, due to a delay τ in the development of the immune response, the magnitude of the immune reaction at time t will be consistent with the level of stimulation at time t − τ. If the level of infection at time t is considerably different from that at the earlier time of the immune stimulation, then the magnitude of the response might appear to be out of consistency with the current level of infection in situ. This in turn would result either in a higher destruction of the infected tissues (which would be otherwise protected by exhaustion mechanisms) or in an insufficient response favoring virus persistence, depending on the kinetics. Therefore, the immune response and the infection kinetics should be highly coordinated to ensure a proper clearance of the infection.

If the virus is not eliminated from the infected host then a chronic or persistent virus infection is established. To diagnose chronic infections with respect to the presence and extent of a disease, and to assess the results of therapeutic interventions, the quantification of the infection process is necessary. A biopsy-based medical testing may be performed involving extraction of tissues for examining the distribution and intensity of the infectious process. A recent example of the study of viral reservoirs in HIV infection based on examination of lymph node and GALT biopsies from subjects under antiretroviral therapy is provided in [7]. It is of note that the sections represent 3 × 10⁻⁴% of the lymph node tissue and the pinch biopsy area is 3 × 10⁻⁶ m² compared to the area of the intestine in an adult human being ∼300 m². This enormous difference in scale rises the issue of how to ensure that the biopsy is informative if the viral infection spread can be spatially essentially heterogeneous as discussed above. Obviously this calls for the need to coordinate interdisciplinary studies of virologists and clinicians with mathematicians in order to clearly define in quantitative terms the spatial patterns of infection. Today infections in tissues are solely classified as diffuse, zonal, patchy, spotty, confluent, etc. Further insights into the virus infection spreading will require extension of the one-dimensional model to the two-dimensional consideration.

As “space” becomes a frontier in immunology [72, 77], further multidisciplinary research is needed towards a mechanistic explanation of these types of pathology signs on a firm basis of the virus ontogeny, the function of specific cells, mediators, and tissue responses analyzed from small blood and tissue samples using computational models with a final aim of identifying the right targets and design more efficient therapies [81].

Supporting Information

Acknowledgments

The research was funded by the Russian Science Foundation (Grant no. 15-11-00029) to G.B., A.M., V.V. A.M. was also partially supported by a grant from the Spanish Ministry of Economy and Competitiveness and FEDER (Grant no. SAF2013-46077-R). S.T. and V.V. were also partially supported by FONDECYT (Chile) project 1150480.