Dynamics of solutions of a reaction-diffusion equation with delayed inhibition

Reaction-diffusion equation with a logistic production term and a delayed inhibition term is studied. Global stability of the homogeneous in space equilibrium is proved under some conditions on the delay term. In the case where these conditions are not satisfied, this solution can become unstable resulting in the emergence of spatiotemporal pattern formation studied in numerical simulations.


1.
Introduction. In this work we study the reaction-diffusion equation with delay suggested in [5] as a model of viral infection spreading in tissues. Here u(x, t) is a normalized virus concentration depending on the space variable x and on time t; the first term in the right-hand side of this equation describes random motion of virus. Its reproduction in the host cells is described by the logistic term u(1 − u), while its elimination by the immune cells by the term −uf (u τ ). The function f (u τ ) describes the concentration of immune cells as a function of the virus concentration at time t− τ , u τ (x, t) = u(x, t − τ ). Time delay takes into account clonal expansion of immune cells stimulated by antigen (virus). This model represents a simplification of a more complete model of immune response described by a system of two reaction-diffusion equations for the concentration of virus and immune cells [6]. This equation can also be considered as a generic model of population growth with delayed inhibition.
In the applications to immune response, the function f (u) grows for small u and decays for large u. There are several stationary points of equation (1), and the behavior of its solutions can be surprisingly complex. Even without diffusion, there are complex oscillations and transition to chaos [2]. In the case with diffusion, there are also waves connecting different equilibria [25]. In this work we prove global stability of the homogeneous in space stationary solutions under some conditions on the function f (u). In the case where these conditions are not satisfied, complex nonlinear dynamics and pattern formation are studied numerically.
In the case of quasi-monotone functional differential equations behavior of solutions can be studied by the comparison principle which excludes complex dynamaics [15], [20], [27]. Investigation of non-monotone functional differential equations is often more involved. In the case without diffusion, there is a number of works devoted to Nicholson blowflies and Mackey-Glass equations [7], [1], [11] - [19], [24], [35]. In particular, equation is studied in [7] for a unimodal function f (u) (it has exactly one maximum). Under some conditions on the initial data, the authors prove the global convergence of solution to the unique positive equilibrium of (2). Functional differential equation with distributed delay, u (t) = −αu(t) + f 1 (u(t)) + τ1 τ0k (s)f 2 (u(s))ds.
is studied in [35]. The authors prove the global convergence of solutions to the unique positive equilibrium of (3) by the method of fluctuations [22]. Global asymptotic and exponential stability of the unique positive equilibrium of the equation is proved by the combination of Lyapunov functionals, monotone semiflow theory and fluctuations method [24] There is a vast literature devoted to delay reaction-diffusion equations (see [8], [23], [28], [29], [30], [32], [33], [36] and the references therein). Most of them consider models in population dynamics such as Nicholson blowflies, Mackey-Glass, and predator-prey models. There are numerous other examples of pattern formation in reaction-diffusion equations (see, e.g., [10], [21]). Let us mention the work [31] where the authors studied the following delay reaction-diffusion equation with the Neumann boundary condition and a unimodal function f . It is proved that the convergence to the unique positive steady state of (5) is strongly linked the convergence of the sequence u n , determined by the difference equation u n+1 = f (u n ). We also mention the work [34] where the authors investigated a nonlocal reactiondiffusion equations in a semi-infinite interval and obtained the global attractivity of the nontrivial equilibrium. Investigation of delay reaction-diffusion equations is often "personalized" in the sense that the methods of their analysis should be adapted to the particular form of the equation. In this work, we develop such methods to prove global convergence of solutions to the positive homogeneous in space equilibrium under certain conditions on the function f (u) (Section 2). In the case where these conditions are not satisfied and this solution becomes unstable, we use numerical simulations to study complex nonlinear dynamics and pattern formation (Section 3).
2. Global convergence to the homogeneous solution. In this work, we study the following class of functional differential equations (6) Throughout this paper, we will make the following assumption: (T1) f is locally Lipschitz continuous.
We define the ordered intervals and for any χ ∈ R, we write χ * for the element of X satisfying χ * (x, θ) = χ for all (x, θ) ∈ [0, 1] × [−τ, 0]. The segment u t ∈ X of a solution is defined by the relation u t (x, θ) = u(x, t + θ) for x ∈ [0, 1] and θ ∈ [−τ, 0]. The family of maps such that (t, φ) → u t (φ) defines a continuous semiflow on X + , [27]. The map U (t, .) is defined from X + to X + which is the semiflow U t , denoted by The set of equilibria of the semiflow generated by (6) is given by Let T (t) (t ≥ 0) be a strongly continuous semigroup of bounded linear operators on C generated by the Laplace operator ∆ under periodic conditions. It is well known that T (t) (t ≥ 0) is an analytic, compact and strongly positive semigroup (see example 1.13, page 26 in [27]). Define F : X → C by We consider the following integral equation with the given initial data For each φ ∈ X, u(., t) with values in C on its maximum interval [0, σ φ ) is called a mild solution of (6) (for existence and uniqueness of this solution, see for instance [9], [16], [17], [27]), and it is called classical if it is in C 2 with respect to x and in C 1 with respect to t. Definition 2.1. Suppose that f is a non-decreasing function. We call a pair of smooth functions (û,ũ) sub-and super-solution of (6) if it satisfies the following properties: ii) The functions (û,ũ) verify the following problems: and Lemma 2.2. Suppose that a smooth function w satisfies the following conditions for some bounded function c(x, t).
The proof of this result is well known. Set L c u = u t − u xx + cu and consider the following problems: and Lemma 2.3. Suppose that f is non-decreasing and a pair (û,ũ) exists. Then we have the following inequalities: Proof. We set w =ũ −ū 1 . Then w satisfies the following problem: From Lemma 2.2, we get w ≥ 0. Similarly we can show thatû ≤ u 1 . Next, we set w =ū 1 − u 1 , and g(u) = u(1 − u). Then we have: where θ is a function betweenû andũ. Since f is a non-decreasing function and u ≥û, then we get Choosing c ≥ f (1) − g (θ(x, t)) and using Lemma 2.2, we obtain w ≥ 0. The assertion of the lemma can now be obtained by induction. Theorem 2.4. Suppose that f is a non-decreasing function, and letũ,û be coupled super and sub-solutions of problem (6). Then the sequences {ū k }, {u k } given by (12)-(13) converge monotonically to a unique solution u of (8) andû ≤ u ≤ũ in Proof. From Lemma 2.3 we conclude that the sequences {ū k }, {u k } converge to some limitsū and u, respectively, andû ≤ u ≤ū ≤ũ. By a classical regularity theorem for parabolic equation, these limits satisfy the equations and the boundary and initial conditions in (6). Hence the limits of {ū k }, {u k } are solutions of (6) ifū = u inD T : Subtracting the equations in (14) and using the mean value theorem, we see that the function w =ū − u satisfies the relation and the same boundary conditions as in (6).
, and the same boundary conditions as in (6), with (after a suitable choice ofc). Next, we multiply the above equation by v and integrate over (0, 1). We obtain: Thus, in view of the equality v(x, 0) = 0 we conclude that , and we show the existence of solution. We proceed similarly as above to prove the uniqueness of the solution betweenû and u.
Lemma 2.5. If φ ∈ X + , then problem (8) admits a unique solution u. In addition, we have the following results: The semiflow U t admits a compact attractor, is a classical solution of (6) for x ∈ [0, 1] and t > τ.
Proof. We first prove that each solution of (8) is bounded. Indeed, the function u solution of (6) is a sub-solution of the following problem Hence, by the classical maximum principle we have u(x, t) ≤ v(t) for all x ∈ [0, 1]. Since v(t) is bounded and converges to 1 as t tends to infinity, then lim sup t→∞ u(x, t) ≤ 1 for all x ∈ [0, 1]. This complete the proof of statement (i). The statements (ii) and (iii) follow from (i) and Theorem 2.2.6 in [27].
We further have the following results.
where Int(X) is the interior of the set X.
). Since φ is positive and not identically null, then in view of Lemma 2.2 we have u(x, t) > 0 for all x ∈ [0, 1] and t > 0.
We consider the following problem Suppose that F satisfies the one-sided Lipschitz condition where c is a nonnegative function and (u,ū) is any pair of ordered lower and upper solution of (15). Proof. The existence of the minimal solution u 1 and the maximal solution u 2 of problem (15) such that u ≤ u 1 ≤ u 2 ≤ū can be deduced by Theorem 3.2.2 in [18]. Now we claim that u 1 = u 2 . Indeed, we have ). Subtracting these two equations and integrating the result over (0, 1), we get By integration by parts and from the periodic conditions, we obtain It follows from u 1 u 2 > 0 and the monotone property of F (u) u that and thus u 1 = u 2 .
Let U (x) be a solution of the stationary problem of (6), namely x ∈ (0, 1), The following result is a direct consequence of Lemma 2.7. Suppose that there exists 0 < u * < 1 such that Next, we study persistence of solutions of problem (6). Lemma 2.9. Suppose that f is an increasing function and f (1) < 1, then the solution of (6) is strongly persistent.
Proof. Since f (0) > 1 − θ, then there exists ε > 0 such that f (0) > 1 − θ + ε. Next, observe that the couple (0, θ−ε) represents sub-and super-solutions of (6), provided In the remaining part of the proof, we will show that each positive solution of problem (6) enters the interval [0, θ − ε] X and stays there. For a nonnegative, non-identically zero function φ , we introduce the following problem: Observe that (0, v(t)) is a pair of sub and super-solutions of problem (6). Therefore, we get 0 ≤ u(x, t) ≤ v(t). Now, we claim that there exists T > 0 such that v(t) ≤ θ − ε for all t ≥ T. Indeed, suppose, first, that v(t) > θ − ε for all t ≥ τ . Then from (19) and the fact that and we obtain a contradiction. Next, suppose that there exists t 1 > t 0 such that v (t 1 ) ≥ 0 and v(t 1 ) = θ − ε. Then This contradiction proves the assertion. Therefore, there exists T > 0 such that u(x, t) ≤ θ − ε for all t ≥ T and x ∈ [0, 1]. The lemma is proved.
Let x n and y n be such that u + (t n ) = u(x n , t n ), and u − (s n ) = u(y n , s n ). Then Hence, we get from (6): Passing to the limit in this equality and using the fact that f is increasing over (M, A), we obtain: Employing the same arguments for s n we get: Next, multiplying the expression (20) by (1 − u ∞ ) and combining the result with (21), we get Therefore, from (H1) we have u ∞ = u ∞ . In the case of (H2), the proof is practically similar to (H1). Indeed, to get the required result, we add (1 − u ∞ ) in (20) and combine it with (21). The same argument is employed if (H3) is satisfied. Thus, in all cases, u(t) converges to a steady state. In view of (18), Lemma 2.9 and Lemma 2.10, each solution of problem (6) converges to the positive steady state. Now, we suppose that f has a single maximum f (M ) = max 0≤u≤1 f (u) and f is increasing over (0, M ) and decreasing over (M, 1). The following lemma will allow us to apply the results of Theorem 2.12 to this case. Proof. First, we suppose that (K1) holds. This implies that there exists ε > 0 such that f (s) ≥ 1 − M + ε for all s ∈ [0, 1]. From this, each solution of (6) satisfies the inequality: We consider the following problem φ(x, s).
By the maximum principle, we have u( . On the other hand, by a straightforward computation we can show that there exists t 0 > 0 such that v(t) ≤ M for all t ≥ t 0 . The assertion (K1) is proved.

Consider the following problem
φ(x, s).
By the maximum principle, we have u( . On the other hand, by a straightforward computation we can show that there exists t 0 > 0 such that v(t) < α for all t ≥ t 0 . The claim is proved.
Further, observe that f (s) > 1 − M for all s ∈ [0, α). Finally, employing the same arguments as in the first part of this proof, we show that there exists T 1 > T + τ such that u(t) ≤ M for all t ≥ T 1 .

3.
Spatiotemporal pattern formation. We study in this section spatiotemporal pattern formation described by the equation in the interval 0 ≤ x ≤ L with periodic boundary conditions and some initial conditions u(x, t) = u 0 (x), −τ ≤ x ≤ 0 defined in the appendix. We will characterize solutions by the positions of their maxima.

Maximum maps.
For each t fixed, denote by x i (t) the coordinates of the strict maxima of the solution u(x, t) of equation (22), that is, of the points for which u(x i (t), t) > u(x, t) for all x in some neighborhood of x i (t). The number of such points is finite or countable, and it can depend on t.
We will call the set of points (x i (t), t), i = 1, ..., N (t) in the domain 0 ≤ x ≤ L, t ≥ 0 the maximum map of the solution u(x, t). We characterize the solution u(x, t) by its maximum map M (u). Some examples are presented in the appendix.
We will see below that this representation of solution can be convenient to characterize complex dynamics. In particular, if the solution can be represented in the form u(x, t) = Φ(t) + v(x, t), where Φ(t) is a periodic function with large amplitude and high frequency of oscillations, and the function v(x, t) has small amplitude and low frequency. Example of such solutions is shown in Figure 1. The maximum map gives a clearer representation of solution than the two-dimensional plot of the function u(x, t) for which high frequency time oscillations hide low frequency amplitude modulations.

Bifurcations and branches of solutions.
In this section we study nonlinear dynamics of solutions of equation (22) in the case where the homogeneous in space periodic solution u 0 (x, t) = Φ(t) loses its stability with respect to spatial perturbations. We fix time delay τ = 1.1, the length of the interval L = 0.5, and vary the value of the diffusion coefficient D.  We will consider the homogeneous in space solution u 0 (t) in the form u 0 (t) = u 0 + φ 1 (t), whereū 0 is its average value and φ(t) a periodic function of variable sign. The first bifurcation occurs for D = 0.001 when the homogeneous solution loses its stability. The bifurcating solution can be approximated as follows: where b = 2π/L, and φ 2 (t) is a periodic function. The value of spatial shift h can be arbitrary since the solution is invariant with respect to translation in space. The spatial perturbation has a single maximum and a single minimum that alternate their positions in time due to the time dependent factor φ 2 (t). The maximum map represents two vertical lines (Figure 2, left; cf. Figure 8, left). By analogy with combustion theory, we call such solutions symmetric modes [26]. In the case of a single maximum, these are one-maximum (1M) symmetric modes. When the principal spatial perturbation u 1 (x, t) = φ 2 (t) cos(bx − h)) vanishes at zeros of the function φ 2 (t), the next perturbations can become visible at the maximum maps (Figure 2, middle; cf. Figure 8, right). These secondary perturbations form dashed lines since they appear at the maximum maps only when they are not hidden by the main perturbations.
Further decrease of the parameter D decreases spatial perturbation of the solution. We expect that this branch of spatiotemporal patterns merges with the branch u 0 (t) of periodic space oscillation for D = 0.00036. For the close value D = 0.0004 (Figure 3, left) the main spatiotemporal perturbation u 1 (x, t) remains as before. However, the maximum map of the next perturbations is now different. It forms a cloud of points with variable densities.
The next branch of solutions is formed by the perturbations u 2 (x, t) = cos(x−ct) specific for travelling waves (Figure 8, middle). Here the constant c is the wave speed. It is observed for the values of D from 0.00035 to 0.000125 ( Figure 4). As before, by analogy with combustion theory [26], such solutions are called onemaximum (1M) spinning modes.
It is interesting to note that the solution with a constant speed shown in Figure 4 (left) loses its stability resulting in appearance of solutions with an oscillation wave speed (Figure 4, middle).
The third branch of solutions appears for D = 0.000195, and it is observed till D = 5 · 10 −4 (smaller values are not considered). It is a 2M symmetric mode, which is similar to 1M symmetric mode but it has two maxima ( Figure 5, left). This solution loses its stability resulting in the emergence of oscillating solutions where the coordinates of the maxima represent periodic functions ( Figure 5, right). Further decrease of D leads to the appearance of more complex aperiodic regimes. The 1M spinning mode and 2M symmetric mode coexist in some interval of D. The response function f (u) is positive for u ≥ 0, and it can be increasing for all positive u, decreasing, or non-monotone. In the latter case, it is increasing for small u and decreasing for large u describing stimulation of immune response by a small viral load and downregulation of immune response by a high virus density due to immunodeficiency and exhaustion of immune cells.
The zeros of the function F (u) = u(1 − u − f (u)) provide homogeneous in space stationary solutions of equation (1). Besides the trivial solution u = 0 it can have up to three positive solutions. Transition between these different stationary points are provided by travelling waves. Their existence and stability is studied in detail in the case without time delay [26]. This question becomes much more involved for the delay reaction-diffusion equation [25].
In this work we study global stability of the homogeneous in space stationary solutions of equation (1). As before, in the case without delay this question is quite simple but it becomes much more complex and interesting in the case with time delay. Even without diffusion, the corresponding delay differential equation can manifest complex dynamics with a sequence of period doubling bifurcations and transition to chaos [2]. Local stability of the stationary solutions can be studied by conventional linear stability analysis. Their global stability requires the development of new approaches. In this work we prove global stability of the homogeneous in space stationary solutions for some classes of functions f (u) (Theorem 2.12).
This solution can lose its stability with respect to time-periodic perturbations resulting in the emergence of homogeneous in space and periodic in time solutions. These solutions can lose their stability in their turn with respect to spatial perturbations leading to the secondary bifurcation of spatiotemporal structures.
In order to study spatiotemporal dynamics of solutions, we introduce the maximum maps which show positions of the local maxima of solutions at the (x, t)-plane. This method appears to be convenient to characterize behavior of solutions taking into account that spatiotemporal structures represent a relatively slow dynamics superposed on high frequency homogeneous in space time oscillations.
The biological meaning of maximum maps depends on the considered applications. In the evolution theory, it describes the dynamics of biological species including their emergence, extinction and drift in the morphological space (see [3], [4]). In the model considered here, maximum maps characterize the dynamics of virus quasi-species in the genotype space or, in another interpretation, the infection development in organs and tissues.