DELAY REACTION-DIFFUSION EQUATION FOR INFECTION DYNAMICS

. Nonlinear dynamics of a reaction-diﬀusion equation with delay is studied with numerical simulations in 1D and 2D cases. Homogeneous in space solutions can manifest time oscillations with period doubling bifurcations and transition to chaos. Transition between two regions with homogeneous os- cillations is provided by quasi-waves, propagating solutions without regular structure and often with complex aperiodic oscillations. Dynamics of space dependent solutions is described by a combination of various waves, e.g., bistable, monostable, periodic and quasi-waves.


1.
Introduction. The dynamics of virus infection in organs and tissues depends on the balance of the reaction-type processes, i.e. the virus replication and the development of antiviral immune response, and the spatial propagation of viruses and cells, e.g. via diffusion. The delay nonlinear reaction-diffusion equations provide a fundamental tool for modelling the spatiotemporal infection dynamics. We consider the equation where u τ (x, t) = u(x, t − τ ), f (u) is a sufficiently smooth function which will be specified below. This equation is introduced in [3] as a model of viral infection spreading in tissues. Here u is the local virus concentration, the diffusion term describes virus random motion in the extracellular matrix or between the host cells. Virus reproduction is described by the logistic term ku (1−u), and its elimination by the immune cells is given by the term kuf (u τ ), where f (u τ ) specifies the efficiency of immune response-mediated virus elimination. The value of function f (·), i.e. the strength of the antiviral immune response, depends on virus concentration at time t − τ . The time delay τ is determined by the duration of proliferation and maturation of immune cells. It can be of the order of several days. The form of the function f (u) is determined by the regulation of immune response as described in detail in [3]. For small u, it is a growing function since the virus infection stimulates immune response. For large u, it decreases since high virus concentration can induce functional or physical exhaustion of virus-specific immune cells. The properties of this function are essential for the analysis of behavior of solutions of equation (1). Let us also note that a more complete model with an additional equation for the concentration c of immune cells is introduced in [2]. It can be reduced to the simplified model with the single equation under some assumptions.
Set (4) We will denote the first one [0, v 0 ]-wave and the second one [v 0 , v 2 ]-wave (the values in the brackets indicate the limits at infinity). The first wave corresponds to the so-called monostable case (since v = 0 is unstable). Such waves exist for all values of the speed c greater than or equal to the minimal speed c 0 . The bistable [v 0 , v 2 ]-wave exists for a unique value of speed c = c 1 . If c 1 > c 0 then there exists a [0, v 2 ]-wave, otherwise it does not exist and there are two waves propagating one after another with different speeds. Existence of such waves and systems of waves is well known as well as convergence of solutions of the Cauchy problem to them (see [13,14] and the references therein).
Thus, there are three regimes of infection spreading in the tissue in the case without delay. It can be transition to a low infection level described by the [0, v 0 ]wave, to high infection level by the [0, v 2 ]-wave, and a two-stage transition, first, by the [0, v 0 ]-wave and then by the [v 0 , v 2 ]-wave.
The situation becomes more complex in the case with delay. In [3] we studied some model examples of the delay reaction-diffusion equations. Here we will present a more systematic investigation. Let us note, first of all, that the stationary solution v 0 of the delay equation can become unstable resulting in the appearance of time oscillations. In Section 2 we study the global stability of the homogeneous in space solution and its instability. The instabilities in DDEs are well known and they are studied in numerous works for various models. If the function f (u) is monotone (Figure 1, left), then only periodic oscillations are observed. In the case of a non-monotone function (Figure 1, right), we also observe period doubling bifurcations and transition to chaos (Section 2). The homogeneous in space solutions can lose their stability with respect to spatial perturbations. In order to study them we consider equation (1) on the interval 0 < x < L with periodic boundary conditions. Many different spatiotemporal patterns can occur depending on the values of parameters. They include standing and travelling waves, and also other periodic and irregular spatiotemporal patterns.
If we consider equation (1) with a piecewise constant initial condition, then we observe two regions with homogeneous in space oscillations and a transition zone between them. This transition zone spreads in space as a wave front separating two regions where behaviour of solutions is different. In particular, these can be homogeneous oscillations from one side of the wave front and various spatiotemporal patterns from the other side. We call such propagating solutions quasi-waves since their structure can evolve with time. Emergence of quasi-waves is related to the instability of homogeneous solutions and propagation of spatial perturbations. We discuss them in Section 3.
Monostable and bistable waves can also exist for the delay equation. The proof of their existence represents a complex mathematical question studied in some particular cases for this and for some other equations [11]. If time delay is sufficiently large, then instead of stationary waves with constant limits at infinity we observe stationary monostable waves periodic at infinity and periodic in time bistable waves. Existence of stationary monostable waves periodic at infinity is studied in [4,5,6,7,8]. The existence of generalized travelling waves is shown for nonlocal equations and for all speeds greater or equal than the minimal speed [1]. Similar results hold for the monostable delay equation. Numerical simulations of monostable waves in the delay equations are presented in [9].
Dynamics of solutions of the delay equation is described by different combinations of monostable and bistable waves with homogeneous oscillations, periodic waves and quasi-waves. We study this dynamics in Section 4. Section 5 is devoted to 2D simulations. We conclude that 1D waves are stable with respect to 2D perturbations while quasi-waves can become unstable. The summary of the observed regimes is presented in Section 6.
2. Homogeneous in space oscillations and their stability. We consider functions f (u) shown in Figure 1. Let us note that the function in the left image is a particular case of the function in the right image where f 2 = f 3 . We consider this particular form of the function f (u) in order to study the properties of solutions starting with this simpler case.
is a stationary solution of equation (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 v(x, t) = v 0 + ǫe λt+iax , where a, ǫ are real numbers, ǫ is a small parameter, and λ is an eigenvalue. Substituting the above function in (1) and equating the terms with the first power of ǫ, we get the following characteristic equation: (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 Therefore the following equalities hold for the real and imaginary parts: Set z = φτ . Then from (7), (8) we obtain 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 0 ) ≥ 1. For the value of z determined from the first equation, we find τ from the second equation. (1) is unstable. Here z is determined from the first equation in (9).
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 equation (1) is also homogeneous in space. In this case, depending on the values of parameters, the solution of the model either converges to the stationary solution v 0 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 0 in the moving coordinate frame attached to the wave, we obtain the same stability conditions as before. It follows from the first equation in (9) 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 0 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.
We use numerical simulations to study large amplitude oscillations of solutions. In the case of monotone function f (u) only periodic oscillations are observed. More complex behavior of solutions occur for the non-monotone f (u) (Figure 2). We vary the value f 3 for all other parameters fixed. We observe period doubling bifurcations with periods 1,2,4,8 for f 3 between 2.5 and 1.898, and then with periods 7,5,3 for f 3 between 1.987 and 1.88. This sequence of periods possibly indicates the presence of chaotic oscillations.

2.2.
Their stability with respect to spatial perturbations. Time oscillating solutions homogeneous in space can lose their stability with respect to spatial perturbations. Near the stability boundary, behavior of solution is determined by the corresponding eigenfunction. Further in the instability region, spatiotemporal patterns can be more complex including irregular spatiotemporal oscillations. Some examples are shown in Figure 3 (monotone f (u)) and Figure 4 (non-monotone f (u)) in the case of periodic boundary conditions. The choice of initial conditions is not essential here.

2.3.
Global stability of the homogeneous in space solution. In this subsection, we study spatially homogeneous solutions of equation (1). Simplifying, we assume here that the coefficients D, k are equal to 1. Then each such solution should satisfy the functional differential equation If we assume that f : [0, +∞) → (0, +∞) is a continuous function and φ is a continuous non-negative function such that φ(0) > 0, then it can be easily verified that this problem has a unique positive solution existing for all t ≥ 0. Furthermore, a straightforward analysis shows that the graph of each positive solution of equation (10) eventually enters and then remains in the horizontal strip R × (0, 1). We suppose that there exists a positive constant u 0 such that Consider first the case where the function f (s) is increasing. We will say that Theorem 2.1. Suppose that f is a strictly increasing function over [0, 1] and either The proof of this theorem is given in the Appendix.
Remark 1. Taking S(x, y) in the forms xy, x + y or x + y + xy, we obtain the following sufficient conditions for the global attractivity of the positive equilibrium u 0 : We will suppose now that f (u) has a single maximum, f (M ) = max 0≤u≤1 f (u) and f is increasing over (0, M ) and decreasing over (M, 1).

Lemma 2.2.
For each positive solution u of (10) there exists T ≥ 0 such that u(t) ≤ M for all t ≥ T if one of the following conditions is satisfied The proofs of these statements are given in the Appendix.
3. Transition between homogeneous oscillations: Quasi-waves. If the stationary point v 0 is unstable and the initial condition u(x, 0) does not depend on x, then the solution of equation (1) is periodic in time and constant in space. Consider now a piecewise constant initial condition, u(x, t) = u 0 for x ≤ L/2, u(x, t) = u 1 for x > L/2 for some u 0 and u 1 , −τ ≤ t ≤ 0. In this case, we observe periodic time oscillations at the left part of the interval and at the right part of the interval with a phase shift. These two intervals where solutions are space independent are separated by a transition zone. In this section we study how the width of the transition zone grows in time and the behavior of solutions inside it.
If the values of parameters are close to the stability boundary of the homogeneous in space solution, and small amplitude oscillations occur, then a maximum of the solution emerges at the center of the transition zone periodically in time. This maximum splits and propagates towards the borders of the transition zone where it disappears. Such solution is shown in Figure 5 (left): periodic oscillations of the solution at the center of the interval (upper row), a snapshot of solution u(x, t) as a function of x (middle row), the positions of the maxima of the solution on the (x, t)-plane (lower row). The number of maxima of the solution remains the same from one period to another one, but they become more flat with time because of the slow growth of the width of the transition zone.
Let us note that maxima of solutions are determined numerically by the following condition: u(x i ) is considered as a maximum if u(x i ) > u(x i±1 ) + δ for some small positive δ (δ = 10 −6 in all simulations). Here x i−1 , x i , x i+1 are three neighboring discretization points. Therefore, if the maximum is sufficiently flat, such that the difference with the neighbouring points is less than δ, then the point is not identified as a maximum. Such blanch zones can be seen in the figure for large simulation time.
The width of the transition zone growth in time as √ t due to diffusion, see Figure  5 (lower, left). Such growth rate is similar to that for the linear diffusion equation with a piecewise constant initial condition.
If we increase τ and move further in the instability region ( Figure 5, middle and right columns), then the width of the transition zone grows linearly in time. New maxima emerge in the transition zone during its growth and move inside it. Modulated periodic oscillations occur at the center of the interval, while the oscillations remain periodic at the boundaries.
Oscillations inside the transition zone can have some specific time and space scales and the speed of spatial displacement. However they are not periodic. Moreover, growth rate of the transition zone and the oscillation inside it can depend on the initial condition. We call such solutions quasi-waves in order to emphasize that they are irregular and possibly transient.  Figure 6 shows examples of quasi-waves for the non-monotone function f (u). Depending on the values of parameters, the speed of growth can be decaying (left) with a single peak of solution moving periodically from the left to the right of the transition zone, or it can be linear in average (middle and right) with multiple peaks appearing, moving and merging inside it. Patterns of maxima location have some characteristic time and space scales. The former is related to the period of homogeneous oscillations, the latter can be related to the wave number of the corresponding maximum eigenvalue.
Let us note that homogeneous oscillations (at the boundary) are periodic in all three cases shown in Figure 6 with periods 1, 2 and 4, respectively. The oscillations at the center of the interval preserve their periodicity in the first two cases and they become irregular in the last case.
In the case of large time delay τ , the amplitude of homogeneous oscillations becomes also large, and the solution approaches the stable stationary point. Though the space homogeneous solution does not converge to the stable stationary point, the space dependent solution converges to it in some space intervals. Such intervals can either gradually grow (Figure 7, left) or they can merge with the other intervals (Figure 7, right). Thus, we observe propagation of an irregular quasi-wave with transition to the stable stationary solution.

Waves without time delay.
Monostable waves. We begin the analysis of one-dimensional (1D) waves in the model without delay:  (12) has a travelling wave solution u(x, t) = w(x − ct) with the limits at infinity. It exists for a unique value of the wave speed c whose sign is determined by the sign of the integral Let us note that the speed c enters the nonlinearity due to time delay. Depending on stability of the stationary point v 0 , there are three types of waves. The first one is observed in the case where the point v 0 is stable (Figure 8, left). Then the wave has the limits at infinity: Stability conditions (9) depend on the wave number a of the perturbation. The most unstable is a = 0. Space periodic modes stabilize the solution. Therefore, solution v 0 can be unstable with respect to spatially homogeneous solutions and stable with respect to oscillatory modes with certain wave numbers.
The wave number of the oscillations behind the wave is determined by the period of space homogeneous oscillations and by the wave speed. The minimal speed of propagation in the case without delay is given by the formula c 0 = 2 DF ′ (0). We expect that it remains the same in the case with delay. We get the analytical value c 0 = 0.006, and numerically c 0 = 0.0059.  Similar waves are observed in the case with delay if the point v 0 remains stable. Existence of such waves is studied in [11,12]. Behavior of solutions becomes different if it is unstable. Then homogeneous oscillations occur in front of the bistable wave. These two regimes are separated in space by various periodic or transient patterns. These intermediate patterns can be either periodic waves, quasi-waves or their combination. Figure 9 (left) shows the case where the speed of the periodic wave is greater than the speed of the bistable wave, and Figure 9 (middle) the case where the relation between the wave speeds is opposite. Combination of a periodic wave and quasi-wave is shown in Figure 9 (right). Different quasi-waves separating the bistable wave and the homogeneous oscillations are shown in Figure 10.
In the monostable-bistable case the dynamics of solution is determined by the relation between the speed of the bistable wave c 1 and the minimal speed of the monostable wave c 0 . If the former is greater than the latter, then one [0, v 2 ]-wave is formed. Otherwise, there are two waves propagating one after another, first the monostable wave and after it the bistable one with a lesser speed. These two waves can be separated by a constant in space solution v 0 , by a periodic wave or by a combination of a periodic wave with a quasi-wave ( Figure 11).
5. 2D problem. In this section we consider the same problem in two space dimensions in a rectangular domain with the no-flux boundary conditions. The initial condition is considered as a step-wise function, positive in some small region and zero elsewhere. Behavior of solutions in this case is similar in many aspects to the solutions in 1D. The homogeneous in space oscillations can lose their stability with respect to spatial oscillations. Resulting patterns depend on the initial perturbation. If the initial perturbation is 1D, that is the initial condition depends, for example, on x and does not depend on y, then this remains true for the solution. 2D perturbations do not develop, at least during the time interval determined by possible duration of numerical simulations. If we consider a radially symmetric initial condition, then the solution remains also radially symmetric in the beginning of the simulation (Figure 12). However, after some time, when the perturbation reaches the boundaries of the square domain, the solution loses its symmetry and, depending on the values of parameters, periodic or irregular oscillations develop. This symmetry loss is related to the geometry of the domain. We can expect that the solution preserves its radial symmetry in a circular domain.
Quasi-waves described above in the 1D case are also observed in the 2D case. In the square domain with radially symmetric initial conditions, the 2D oscillations without radial symmetry develop with time. Figure 13 shows two examples of quasiwaves. In the first case, the region of 2D oscillations is separated from the region of homogeneous in space oscillations by a radial wave propagating in space (cf. Figures  5, 6). In the second case, the 2D quasi-wave separates the stable stationary value v 2 and the region with the homogeneous in space periodic in time oscillations (cf. Figure 7). Thus, radially symmetric quasi-waves can lose their stability and 2D oscillations emerge. It can happen because square grid provides perturbations of radial solutions. We finish this section with the 2D simulations of wave propagation in rectangular domains. If the initial condition is 1D, then the solution is also one-dimensional, and 2D perturbations do not develop during the time interval limited by the duration of numerical simulations. If we introduce a small perturbation in the initial condition ( Figure 14), then we observe propagation of the 1D wave. However, the quasi-wave, which propagates slower, develops 2D oscillations. On the other hand, the quasiwave remains 1D without initial 2D perturbations. Thus, 1D waves are stable with respect to 2D perturbations and 1D quasi-waves are unstable. 6. Conclusions: Summary of the regimes. The following regimes are observed.
• Homogeneous in space: periodic, period doubling, transition to chaos.
• If homogeneous in space oscillations lose their stability with respect to spatial perturbations, then various periodic or aperiodic patterns emerge.  • Quasi-waves propagating with a constant (in average) speed or decaying speed and separating the regions with homogeneous oscillations or travelling waves.
They can be periodic in time and close to periodic in space, or they manifest complex aperiodic spatiotemporal dynamics without clearly established structure. • Monostable waves. They have limit v = 0 at +∞ and v 0 or a periodic in space solution at −∞. These are stationary waves, that is solutions of equation (15). The monostable wave can also have limit v 2 at −∞ when it results from merging of the monostable and bistable waves. • Bistable waves. They have limit v 2 at −∞ and v 0 at +∞. If the latter is unstable, then a time periodic wave is observed. • Combinations of different waves: -monostable -bistable: two waves propagate with different speeds one after another. Periodic or quasi-waves can separate them, -monostable -homogeneous oscillations: the unstable equilibrium v 0 is established after the monostable wave. The region with homogeneous oscillations is separated by an interval growing in time, -periodic bistable -homogeneous oscillations, -periodic bistable -periodic wave -homogeneous oscillations. Quasi-waves can separate the first or the second pair.
• 2D waves. One-dimensional waves remain stable with respect to 2D perturbations. Quasi-waves are unstable and they can develop 2D irregular spatiotemporal patterns. In this work we have identified various new regimes in the dynamics of delay reaction-diffusion equation. Biologically, our study suggests a richness of spatiotemporal patterns of infection dynamics in organs and tissues. The novel insights should stimulate further interdisciplinary research in this important field. Some questions remain beyond the scope of this study. We do not construct detailed bifurcation diagrams and do not study multiplicity of solutions (for the same values of parameters) specific for such kind of nonlinear dynamics. These are the topics for future investigations. 1150480. The work was supported by the "RUDN University Program 5-100" (to V.V.) and the Russian Science Foundation grant number 18-11-00171 (to N.B. and G.B.).

Appendix.
Proof of Theorem 2.1. For a positive solution u(t) of equation (10), consider the following values We are going to prove that actually u ∞ > 0.
Next we set v(t) = u(t) − α and we claim that v(t) > 0 for all t ≥ τ. By contradiction we assume that there exists t 0 > τ such that v(t 0 ) = 0 and v(t) > 0 for all t ∈ [0, t 0 ], so we have v ′ (t 0 ) ≤ 0. By substituting v(t 0 ) in (10) we get the following contradiction proving the claim.
Assume now that f (0) > 1 − θ where f (θ) = 1. Take ε > 0 such that We claim that there exists T ≥ τ such that u(t) < θ − ε for all t ≥ T. We first suppose that u(t) ≥ θ − ε for all t ≥ τ. Then for t ≥ 2τ Consequently there exists t 0 ≥ τ such that u(t 0 ) < θ −ε. Next, if there exists t 1 ≥ t 0 such that u(t 1 ) = θ − ε and u ′ (t 1 ) ≥ 0, then (17) yields the following contradiction proving the claim. Finally, since f (θ − ε) < 1, and u(t) < θ − ε, for all t ≥ T, then we can use the same idea as in the first part of this proof in order to show the persistence of the solution of (10). This shows that u ∞ > 0. Next, it is easy to see that there exist sequences {t n } and {s n } converging to +∞ and such that u ′ (t n ) → 0, and u(t n ) → u ∞ ≥ lim u(t n − τ ) as n → ∞ and u ′ (s n ) → 0, and u(s n ) → u ∞ ≤ lim u(s n − τ ) as n → ∞.
Hence, g * (u ∞ ) = u ∞ > 0 so that u ∞ = u ∞ = u 0 . This means that each solution u(t) of equation (10) converges to the positive equilibrium. The theorem is proved.
Proof of Lemma 2.2. We suppose that (K1) is satisfied. We first claim that there exists t 0 ≥ τ such that u(t 0 ) ≤ M. In fact, suppose on the contrary that u(t) > M for all t ≥ τ. Then in view of (10) and (K1) we have, hence lim t→∞ u(t) = l ≥ M. Furthermore, there exists a sequence (t n ) such that u(t n ) → l and u ′ (t n ) → 0 and t n → ∞. Substituting u(t n ) in (10) and passing to the limit we get either l = u 0 < M or l = 0 which is a contradiction. Consequently there exists t 0 ≥ τ such that u(t 0 ) ≤ M. Suppose now that there exists t 1 ≥ t 0 such that u(t 1 ) = M and u ′ (t 1 ) ≥ 0. Thus, again by the equation of (10) we have we reach a contradiction. Then there exists T ≥ 0 such that u(t) ≤ M for all t ≥ T. Now suppose that (K2) is satisfied. We begin by proving that there exists T ≥ 0 such that u(t) < A for all t ≥ T. To do that, we first assume that u(t) ≥ A for all t ≥ τ. Using the equation of (10) and (K2) we obtain for t ≥ 2τ, Combining this observation with (20) we reach a contradiction. Therefore there exists T ≥ 0 such that u(t) < A for all t ≥ T.
Further observe that f (s) > 1 − M for all s ∈ [0, A). We claim that there exists T 1 ≥ T + τ such that u(t) ≤ M for all t ≥ T 1 . Indeed, we first suppose that M ≤ u(t) < A for all t ≥ T + τ then since f is decreasing over [M, A], we have Proof of Theorems 2.3 and 2.4. The proof of these theorems is based on Theorem 2.1 and Lemma 2.2. Indeed, from Lemma 2.2, u(t) ≤ M for all t ≥ T. Further, according to Theorem 2, the condition f (M ) < 1 or f (0) > 1−θ with f (θ) = 1 gives the persistence of solutions of problem (10). Finally, since f is an increasing function over (0, M ), we can repeat the arguments developed in the proof of Theorem 2.1 in order to establish the convergence of u(t).