Effect of time delay on pattern dynamics in a spatial epidemic model.

Time delay, accounting for constant incubation period or sojourn times in an infective state, widely exists in most biological systems like epidemiological models. However, the effect of time delay on spatial epidemic models is not well understood. In this paper, spatial pattern of an epidemic model with both nonlinear incidence rate and time delay is investigated. In particular, we mainly focus on the effect of time delay on the formation of spatial pattern. Through mathematical analysis, we gain the conditions for Hopf bifurcation and Turing bifurcation, and find exact Turing space in parameter space. Furthermore, numerical results show that time delay has a significant effect on pattern formation. The simulation results may enrich the finding of patterns and may well capture some key features in the epidemic models.

Introduction

Major public health threat brought about by infectious diseases has drawn increasingly attention in recent years, such as the severe acute respiratory syndromes (SARS) outbreak in 2003  [1], [2], the avian influenza A (H7N9) outbreak in China in 2013  [3], [4]. A great deal of effort has been made toward exploring more realistic mathematical models for the transmission dynamics of infectious diseases. The aim of these models is to both explain the observed epidemiological patterns and predict the consequences of the introduction of public health interventions to control the spread of diseases (see e.g., Diekmann and Heesterbeek  [5] for an overview).

In the past several decades, various epidemic models have been formulated and studied to reveal mechanisms of disease transmission (see, for example,  [6], [7] and the references therein). In modeling of communicable diseases, the incidence rate, i.e. the rate of new infections, often plays a key role in guaranteeing that the model does indeed give a reasonable qualitative description of the disease dynamics  [8], [9]. Great interest has been devoted to investigating the nonlinear dynamics, including Hopf bifurcation, saddle–node bifurcation, Bogdanov–Takens bifurcation, existence of periodic and homoclinic orbits, and so on. For example, Yorke and London  [10] showed that the incidence rate with positive and time-dependent was in excellent agreement with the simulation results for measles outbreaks. Capasso and Serio  [11] adopted a saturated incidence rate of the form , , to represent “crowded effect” or “protection measure” in modeling the cholera epidemics in Bari in 1973. To incorporate the effect of behavioral changes for certain communicable diseases, Liu and coworkers  [12], [13] extended that idea and proposed a nonlinear saturated incidence function , where represented the infection force of the disease, and measured the inhibition effect from the behavioral change of the susceptible individuals when the number of infectious individuals increased, with , , , . There were numerous results about nonlinear incidence rates in the literature, and we refer the reader to Ref.  [8] for a general review.

Apart from the incidence rate, many other factors, such as stochasticity  [14], seasonality  [15], [16], the distribution of latent and infectious periods  [17], [18] and spatial structure  [19], can influence epidemic dynamics. In particular, the study of spatial epidemiology becomes an exciting and important area of research because space may strongly influence lots of important epidemiological phenomena due to the localized nature of transmission or other forms of interaction  [20], [21], [22], [23]. As a result, mathematical models with both time and space are more realistic in considering the process of epidemic spreading. To our knowledge, both discrete patchy models and continuous reaction–diffusion systems are used to study the spatial heterogeneity  [24]. Patchy models are often used to describe directed movement among patches, while reaction–diffusion systems are suitable for random spatial dispersal. For example, we may find epidemic wavefronts in reaction–diffusion models, which corresponded to the real world observations as in the spread of the Black Death in Europe from 1347 to 1350  [25], [19].

Moreover, time delay is ubiquitous in most biological systems like predator–prey models and epidemiological models. We note that in disease transmission models, time delay is an important quantity accounting for many epidemiological mechanisms. In particular, time delays can be introduced to model constant incubation period or sojourn times in an infective state. For example, using an model with a maturation delay and vertical disease transmission Busenberg et al.  [26] obtained some periodic solutions. An epidemic model with exponential demographic structure, disease related deaths and two delays corresponding to the latent and immune constant periods respectively was first analyzed by Cooke and van den Driessche  [27]. They identified a delay-dependent threshold parameter determining the local asymptotical stability of the equilibrium states. For a brief review of delay differential equations arising from disease modeling, we refer the reader to van den Driessche  [28].

In our previous papers  [29], [30], [31], [32], [33], [34], [35], we considered the mechanisms of pattern formation in epidemic models induced by spatial (cross-) diffusion or noise. To the best of our knowledge, little attention  [36] has been paid to the joint effect of delay and diffusion on pattern dynamics, and the mechanism of delay-induced Turing instability is not well understood. Our main focus, however, is on how the incorporation of a constant time lag (an incubation period) into the epidemic model alters the aforementioned qualitative results.

The paper is organized as follows. In Section  2, we introduce a spatial epidemic model with Neumann boundary conditions and nonzero initial conditions. In Section  3, we summarize the dynamics of spatial epidemic model without time delay. Then, we derive the conditions for Turing instability occurring in Section  4 through the method of linear stability analysis. In Section  5, we illustrate the effect of time delay on the emergence of Turing patterns by performing extensive numerical simulations. Finally, conclusions and discussion are presented in Section  6.

A spatial epidemic model

In this paper, we use a simple - epidemic model with nonlinear incidence rate to investigate the effect of time delay on the spatial epidemic model. Let and be the number of susceptible, infected and infectious individuals at time . To incorporate saturation or multiple exposures before infection, Liu et al.  [12], [13] proposed a nonlinear incidence rate with , . This form of nonlinear incidence rate, without a periodic forcing, could produce much wider range of interesting dynamical phenomena in comparison to bilinear incidence rate . Owing to their simple form, they cannot involve many of the complex biological factors. However, they often shed insightful lights to help us understand some complex processes  [30], [34], [35]. In the present paper, let and . Then we have the following spatial epidemic model with nonlinear incidence rate where is the recruitment rate of the susceptible, is the natural death rate of the population, and is the disease-related death rate from the infected. Here, represents the space, and is the usual Laplacian operator in two-dimensional space. The susceptible and infected individuals diffusion coefficients are denoted by and , respectively. From the biological point of view, we suppose that all the parameters are positive throughout the paper. More details about the model can be found in our previous paper  [35].

In general, we are interested in the self-organization of patterns and choose the following nonzero initial conditions and Neumann (zero-flux) boundary conditions where denotes the size of the system in the directions of and , is the outward unit normal vector of the boundary . Neumann boundary conditions imply that the boundary of the model domain is simply reflective, and that the domain is isolated or insulated from the external environment  [37].

If the incubation period is assumed to be a constant , we obtain a delayed spatial epidemic model with nonlinear incidence rate as follows, with the initial conditions

Some results about the model without delay

To give insights into the effect of time delay or diffusion on system (3), it is of significance to investigate the local dynamics of its corresponding non-delay and non-diffusion model. Now, we consider the case of spatially homogeneous states, i.e. non-diffusion model. From a biological point of view, we are interested in the non-negative steady states , . We summarize the local dynamics of the model around the equilibrium states, see our previous paper  [35] for more details. From Ref.  [35], we know that the system has three equilibrium states

, which corresponds to extinction of the disease;

Coexistence of the and population which is unstable and a saddle by direct calculations;

Coexistence of the and population which is a stable node. Hereafter, we denote by this local stable equilibrium. To ensure the positivity of and , one has the following condition

In the sequential, we only focus on the locally asymptotically stable equilibrium . For the sake of establishing our main results in Section  4, we give a brief review on diffusion-induced Turing instability based on  [35] for system (1).

Linearizing system (1) around the spatially homogeneous equilibrium point (, ) for small space- and time-dependent fluctuations and expanding them in Fourier space we have the following characteristic equation where and the Jacobian matrix at is given by where and Here, or corresponds to the one- or two-dimension space, and are spatiotemporal perturbations, and is the wave number.

Solving Eq. (5) yields the characteristic polynomial of the original problem (1) where and Therefore, the roots of (6) yield the dispersion relation

At the bifurcation point, two equilibrium states of the model intersect and exchange their stability. Biologically speaking, this bifurcation corresponds to a smooth transition between equilibrium states. The Hopf bifurcation is space independent and breaks the temporal symmetry of a system, which gives rise to space-uniform and periodic-time oscillations, while the Turing bifurcation breaks spatial symmetry, leading to the formation of stationary-time and space-oscillation patterns[29], [32].

Hopf instability occurs if a pair of imaginary eigenvalues cross the real axis from the negative to the positive axis and the diffusion is absent  [29], [32]. From a mathematical point of view, the Hopf bifurcation appears when

A positive equilibrium state is said to be Turing instability if it is stable for the corresponding non-spatial model of (1), but becomes unstable with respect to homogeneous perturbation due to diffusion. A general linear analysis  [38], [39] shows that the necessary conditions for yielding Turing patterns for model (1) are given by The first two conditions guarantee that the equilibrium (, ) is stable for the non-diffusion model of (1), and becomes unstable for model (1) if transits from a negative value to a positive one (corresponding to last two conditions).

Mathematically speaking, the Turing bifurcation occurs when and the wavenumber satisfies

Analysis of time-delayed model (3)

Linear stability analysis

Now in this subsection, we will concentrate on the stability of system (1). Obviously, system (3) has the same equilibria as that of corresponding non-diffusion model of system (1). Following the approach in Ref.  [40], [41] and assuming to be small we replace in Eq. (3) to write as

Expanding in Taylor series and neglecting the higher-order nonlinearities, then Eq. (10) becomes where , and . Noting that , we obtain the following equations

The homogeneous steady state (i.e. fixed point) of the dynamical system satisfies and . Considering small spatiotemporal perturbations and on a homogeneous steady state (, ), then we have

By expanding the reaction terms around the steady state in Taylor series up to first order and rearranging the terms, we finally obtain

Assume that spatiotemporal perturbations and take the following form where is the growth rate of the perturbation in time , and and are the wavenumbers of the solutions. Upon inserting (15) into Eq. (14), we obtain the following matrix equation for eigenvalues where .

From Eq. (16) we get the following characteristic equation for the eigenvalues of the associated stability matrix where

Our aim here is to establish the threshold or critical value of the delay time for which the system with time delay, which is otherwise stable with respect to homogeneous perturbation, becomes unstable. Now we are in a position to investigate the effects of time delay and diffusion on the dynamical system (3), and under what conditions for time delay to destabilize the steady state of the system and bring about spatiotemporal instability. We know that the onset of instability requires that at least one of and is violated. Note that and , we always have . Hence, we consider the potential growth instability in the following two cases: (i) is violated; (ii) is violated.

The conditions for stability of the homogeneous steady state for system (3) with delay are and . Since is always positive, the condition can be obtained if . This allows the range of values for determined by the following condition

if , i.e.  , for all small time delays , then is always satisfied, which implies that homogeneous steady state is locally asymptotically stable independent of delay;

if , we obtain the following range of values for ensuring the stability

Diffusion-induced instability

In this subsection, we consider the second case, i.e.  is violated. In other words, the condition of instability is determined by the sign of . Since the numerator of determines the Turing critical line which is found to be independent of , the emergence of the instability in this situation is induced only by the diffusion. Direct calculation shows that the critical value of bifurcation parameter for Turing critical line equals where

The effect of delay on system (3)

In this subsection, we consider the first case, in which is violated. Delay-diffusion-induced instability means that a positive equilibrium is uniformly asymptotically stable in the reaction–diffusion model without a delay effect (e.g., system (1)), while it becomes unstable with respect to homogeneous steady state for the reaction–diffusion model with a delay (e.g., system (3)).

In this case, the condition of instability is . Moreover, since the denominator of is always positive, the condition of instability reduces to the following: which implies that the lower bound of must satisfy for some .

There may be a large variety of distinct patterns that can be observed by varying the parameters slightly. To well see the bifurcation condition, we let , , and . In Fig. 1 , we show the Turing space in plane for model (3). is positive equilibrium existence line (the black line), is Hopf bifurcation line corresponding to (the green one), is Turing bifurcation line (the red one). The Turing space marked is bounded by and . For parameters in this domain, above , the positive equilibrium of corresponding non-delay non-diffusive model is stable; below , the corresponding solution of system (1) is unstable. In other words, Turing instability occurs, therefore Turing patterns emerge. Introduction of delay may give us another useful handle for further manipulation of the instability region between Hopf and Turing curves. The Hopf bifurcation curve (the blue line) corresponding to is also shown in Fig. 1. Therefore, for a fixed parameter set and delay , the transmission rate beyond a critical threshold asserts a condition of instability of the homogeneous steady state of the system even if is below the Hopf bifurcation line corresponding to (the green one). We proceed to explore this in the next section.

Fig. 1

Bifurcation diagram for parameter region for , , , and . The Turing space of model (1) is marked by . Region is located above the positive equilibrium existence line (the black line ), and bounded by the Hopf bifurcation line (the green line ) and the Turing bifurcation line (the red line ). Also, the Hopf bifurcation line corresponds to , while or 0.05, the Hopf bifurcation curve for system (3) (the blue line) is below .

Main results

To compare the analytical predictions from the aforesaid analysis, we have to perform extensive numerical simulations of the system under the influence of time delay, using Eq. (3), by the explicit Euler method with parameters defined above. The continuous problem defined by the reaction–diffusion system in two-dimensional space is solved in a discrete domain with lattice sites. All our numerical simulations employ the Neumann boundary conditions with a time step size of and space step size . In the present paper, we set , and approximate the Laplacian describing diffusion by differences over   [42]. The simulations are initiated with spatially random perturbations of ∼0.01% around the steady state . In this paper, we want to know the effect of delay on the distribution of infected individuals. As a result, we only show results of pattern formation about one distribution of . We conduct the simulations until they reach a stationary state or suggest a behavior that does not seem to change its characteristics anymore.

In the following part, we choose suitable values of parameters for simulation: For the different values of located in the Turing space (the domain in Fig. 1), different categories of Turing patterns for the distribution of can emerge. In each pattern, the blue (red) represents the low (high) density of the infected . However, we concern the effect of time delay on the pattern formation of .

For the non-delay spatial model (1), i.e.  in model (3), when for the above parameter set, Turing instability occurs. Figs. 2 –3 present the evolution of the spatial pattern of infected population at 0, 1000, 50,000, 100,000 iterations, for and 42 in the Turing space, respectively. In Fig. 2, the regular stripe patterns prevail over the whole domain at last, and the dynamics of the system does not undergo any further changes; while in Fig. 3, the spotted spatial patterns prevail the whole domain at last. We use these two case as baseline, and investigate the effect of delay on the pattern formation.

Fig. 2

Snapshots of contour pictures of the time evolution of the infected at different instants in heterogeneous environment. Numerically simulated (in two-dimensional space; grid size 200×200 with and ) delayed system (3) for the parameters , , , , , , and , which are in the Turing region. (a) 0 iterations; (b) 1000 iterations; (c) 50,000 iterations; (d) 100,000 iterations.

Fig. 3

Snapshots of contour pictures of the time evolution of the infected at different instants in heterogeneous environment. Numerically simulated (in two-dimensional space; grid size 200×200 with and ) delayed system (3) for the parameters , , , , , , and , which are in the Turing region. (a) 0 iterations; (b) 1000 iterations; (c) 50,000 iterations; (d) 100,000 iterations.

The parameter values of Fig. 4 are the same as those in Fig. 2, but with time delay, . All of the figures show the evolution of the spatial patterns at 0, 1000, 50,000, 100,000 iterations, with small random perturbation of the stationary state and of the spatially homogeneous systems. From these figures, we can see that the regular stripe patterns also prevail over the whole domain at last, and the dynamics of the system does not undergo any further changes. However, it should be noted that there is an essential difference between Fig. 2 (non-delayed spatial system) and Fig. 4 (delayed spatial system), both starting from the same initial conditions. In particular, the large regular stripe patterns in Fig. 2 break down to small pieces in Fig. 4, and the direction of the stripe of Fig. 4 is also changed.

Fig. 4

Snapshots of contour pictures of the time evolution of the infected at different instants in heterogeneous environment. Here, the parameters are the same as in Fig. 2 except . (a) 0 iterations; (b) 1000 iterations; (c) 50,000 iterations; (d) 100,000 iterations.

The parameter values of Fig. 5 are the same as those in Fig. 3, but with delay . Different from Fig. 3, in which the spotted patterns prevail the whole domain finally, we find that stationary stripe and spot patterns coexist in the distribution of the infected population density at last in Fig. 5 and the dynamics of the system does not undergo any further changes. Moreover, we can see that the time delay largely postpone the pattern formation of the infected population , for example, comparing Fig. 3(c) with Fig. 5(c).

Fig. 5

Snapshots of contour pictures of the time evolution of the infected at different instants in heterogeneous environment. Here, the parameters are the same as in Fig. 3 except . (a) 0 iterations; (b) 1000 iterations; (c) 50,000 iterations; (d) 100,000 iterations; (e) 150,000 iterations; (f) 200,000 iterations.

Fig. 6 shows an interesting phenomena. For the delayed spatial system (3), if the parameter values are set as those in Fig. 4, Turing stability occurs for . In other words, delay can widen the Turing space. Letting , we find that stationary stripe and spot patterns emerge mixed in the distribution of the infected population density, and the dynamics of the system does not undergo any further changes finally. It should be noticed that the patterns are induced only by delayed-spatial system since there are no patterns for non-delayed spatial system, i.e.  .

Fig. 6

Snapshots of contour pictures of the time evolution of the infected at different instants in heterogeneous environment. Here, the parameters are the same as in Fig. 4 except . (a) 0 iterations; (b) 1000 iterations; (c) 50,000 iterations; (d) 100,000 iterations.

Conclusions and discussion

To conclude, a spatial epidemic model with both nonlinear incidence rate and time delay is investigated. The numerical results correspond perfectly to our theoretical findings. Specifically, there is a range of parameters in plane where the different spatial patterns can be obtained. On the one hand, for the non-delayed spatial system, different spatial patterns, such as regular stripe and spot patterns, can emerge for a range of parameter values  [35]. On the other hand, time delay has a significant impact on the pattern formation. More specifically, there are three aspects. Firstly, time delay can change the pattern essence, such as the arm length and the direction of the stripe, see Fig. 2, Fig. 4. Secondly, time delay can largely postpone the pattern formation of the infected population density , comparing Fig. 3, Fig. 5. Thirdly, time delay can widen the Turing space. In other words, for the non-delayed spatial epidemic model, there is no pattern occurring; while for the delayed one, stationary stripe and spot patterns may coexist in the distribution of the infected population density, and the dynamics of the system does not undergo any further changes finally (see Fig. 6 and main text for more details).

By the above analysis, we can find that the qualitative dynamics of the delayed spatial epidemic model (3) is fundamentally different from the non-delayed one (1) when time delay is slightly changed. In the past few years, a great deal of attention has been paid to transitions between different dynamical regimes as a result of perturbation of the system’s parameters  [43]. For simplicity, we let other parameter values remain fixed and vary only one parameter, such as or in the present paper. In Ref.  [29], Sun et al. studied a spatial model with logistic growth and nonlinear incidence rates with , and obtained not only a stripe-like pattern but also a spot pattern, or coexistence of the two. Then, in Ref.  [35] they assumed the recruitment of population with constant rate and set and . However, little attention has been paid to the time delay accounting for constant incubation period or sojourn times in an infective state, which may have significant impact on pattern formations of infected population density. Hence, in this paper, we unfold it by extensive numerical simulations.

The methods and results in the present paper may enrich the research of pattern formation in the spatial epidemic models and may well explain the field observations in some areas. For example, Jewell et al.  [44] investigated the spatial and temporal dynamics of foot-and-mouth disease outbreak during the 2007, and suggested undetected (occult) infections in the UK. Su et al.  [45] used a reaction–diffusion system to explore the relationship between malaria fever and parasite replication cycles. Here, we would like to remark that conclusive evidence of spatial patterns that have the peculiarities in interacting epidemiological systems is still to be found and there are an increasing number of indications of patterns in realistic ecosystems, such as vegetation distribution  [46], planktonic interaction  [47] and a few work on prey–predator type interaction  [48]. Our spatial epidemic model with time delay may be more realistic to capture some key features of the complex variation and to explain the observation in spatial structure to most species  [49], [50].

However, it should be noted that the method in this paper is particularly suitable for short time delay based on the Taylor series expansions, whereas the time delay is much larger, one should investigate the stability matrix  [51], [52], [53] of a reaction–diffusion system and give necessary or sufficient conditions which guarantee that its uniform steady state undergoes a Turing bifurcation. Usually, one needs taking care of the joint effect of diffusion and time delay  [54] to obtain stability conditions which depends on transcendental equation associating characteristic eigenvalue with a function . In addition, the presence of noise  [55], [56], [57], [58] or other terms may give rise to a rich variety of dynamical effects  [59], including noise-enhanced stability  [60], noise-delayed extinction  [61], [62] and noise-induced transitions  [63]. For example, pattern formation induced by the noise in two competing species was analyzed by Valenti et al.  [55]. Nonmonotonic behavior of spatiotemporal pattern formation in noisy population dynamics has been found by Fiasconaro et al.  [56]. A model for epidemic dynamics was analyzed by Chichigina et al.  [57], using a pulse noise model with memory. The effect of multiplicative noise, always present in population dynamics, in the form of a pulse train with regulated periodicity, on a parametric instability was investigated by Chichigina et al.  [58]. These and some other related issues are left for further investigation and discussion.