Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessReview articles

Travelling fronts in time-delayed reaction–diffusion systems

Jan Rombouts

Jan Rombouts

Laboratory of Dynamics in Biological Systems, Department of Cellular and Molecular Medicine, University of Leuven (KU Leuven), 3000 Leuven, Belgium

Google Scholar

Find this author on PubMed

Lendert Gelens

Lendert Gelens

Laboratory of Dynamics in Biological Systems, Department of Cellular and Molecular Medicine, University of Leuven (KU Leuven), 3000 Leuven, Belgium

Google Scholar

Find this author on PubMed

Thomas Erneux

Thomas Erneux

Optique Nonlinéaire Théorique, Université Libre de Bruxelles, Campus Plaine C.P. 231, 1050 Bruxelles, Belgium

[email protected]

Google Scholar

Find this author on PubMed


    We review a series of key travelling front problems in reaction–diffusion systems with a time-delayed feedback, appearing in ecology, nonlinear optics and neurobiology. For each problem, we determine asymptotic approximations for the wave shape and its speed. Particular attention is devoted to their validity and all analytical solutions are compared to solutions obtained numerically. We also extend the work by Erneux et al. (Erneux et al. 2010 Phil. Trans. R. Soc. A 368, 483–493 (doi:10.1098/rsta.2009.0228)) by considering the case of a slowly propagating front subject to a weak delayed feedback. The delay may either speed up the front in the same direction or reverse its direction.

    This article is part of the theme issue ‘Nonlinear dynamics of delay systems’.

    1. Introduction

    Time delays naturally appear in our everyday life. The reaction time of a driver is about 2 s, although it depends on various factors such as muscle memory, anticipation and age. The stability of quietly standing up depends on a neural controller that creates corrective actions after a time lag of 70–90 ms. Pharmacodynamic effects of drugs are typically delayed in relation to plasma drug concentrations. If the delay is short (minutes) then the mechanism is probably a distribution process, whereas if the delay is long (hours or longer) then the mechanism is more likely to be physiological [19]. In ecology, delays arise from maturation times, and in epidemiology, there is a time lag between host infection and disease. When time-delayed systems are coupled in space, diffusion comes into play, allowing a wide range of spatio-temporal phenomena, such as travelling structures, waves or localized states.

    Time-delayed reaction–diffusion problems typically appear in population dynamics where the delay τ is the time from birth to maturity [10,11]. They are also used to model infectious diseases where the delay τ represents the time needed for a cell that is infected by a virus to release new viruses [1216]. Often, a solution consists of a travelling front. These fronts are invasion waves where a spatially uniform steady state is propagating with a constant speed into another state.

    Two distinct approaches to modelling the delayed reaction–diffusion problem have been proposed. One direction consists of supplementing the classical reaction–diffusion equation by terms proportional to the delay [1315]. The migration and growth of the density of population p is described by

    where F(p) represents the instantaneous reaction or growth term. This equation contains Cattaneo modifications of the Fick diffusion equation as proposed in 1948 to overcome the infinitely fast propagation [17]. Different derivations of this equation are possible. One of them consists of expanding Fick's Law with a time delay: the flux of particles depends on the gradient at a previous time. By then expanding the future flux J(t + τ) around τ = 0, we obtain equation (1.1). Another way of deriving it is by starting from microscopic considerations. A modern derivation of the Cattaneo equation can be found in [18]. By implementing the time delay in this way, the time-delayed problem is reduced to the determination of a front solution of a conventional reaction–diffusion problem. This greatly facilitates comparisons with experimental observations [15].

    A different and more radical approach to modelling the effects of a time delay is to deal with the functional equation itself. Fickian diffusion is combined with a nonlinear reaction term that depends on a time delay [10,16]. The delayed reaction–diffusion equation then has the form

    In ecology, the function F(p, p(t − τ)) is typically inspired by reliable growth models such as the delayed logistic equation [19]. In some specific cases, equation (1.2) may be rigorously derived from an age-structured model [10]. Needless to say, the analysis of equation (1.2) is mathematically difficult.

    In this work, we limit our review to scalar reaction–diffusion equations and concentrate on monostable and bistable front solutions. More complex spatio-temporal structures have been studied for higher-order reaction–diffusion systems (see, for example, [2024]). Our main objective is to explore several analytical techniques appropriate for constructing front solutions of equation (1.2).

    These methods take advantage of either a particular form of F(p, p(t − τ)) or limiting values of some parameters. The goal is to derive valuable approximations for both the front shape and its propagation speed. Although these methods have been successfully used for classical reaction–diffusion problems [25], results with time delay remain rare. Here, we consider a series of key problems in ecology, nonlinear optics and neurobiology. For most of these problems, we determine asymptotic approximations and verify their validity by numerical simulations of the original equations.

    The outline of the paper is as follows. The case of monostable fronts is the most documented problem in the population biology literature. In §2, we review that the speed of propagation of these fronts is a decreasing function of the delay in a typical case, where the time delay is due to maturation. When the population dynamics includes an Allee effect, we obtain bistable fronts linking two stable steady states. The interest in bistable fronts is more recent, but it is also very relevant for biochemical reaction–diffusion systems. In §3, we determine the propagation speed for two different problems and contrast their behaviours in the limit of large delays. Finally, we discuss some future directions and possible applications in the field of cell biology.

    2. Monostable travelling fronts

    The ecology of invasions by animals and plants has received more and more attention in recent years, mainly because nearly every ecosystem has been invaded by exotic organisms with potentially drastic consequences for the native fauna or flora [26]. Mathematical reaction–diffusion models are proposed to describe or predict the fate of some particular invasions. Early attempts can be traced back to 1937, when R.A. Fisher was interested in the spatial spread of new genes that appeared in a population. He proposed the now famous FKPP (Fisher–Kolmogoroff–Petrovsky–Piscounoff) equation [25, p. 439]

    where u is the local population density. In this model, population expansion arises from, on the one hand, a balance between the local growth with maximum population density α/β, where α is the linear growth rate, and, on the other hand, the diffusion of individuals with diffusion coefficient D. A linear stability analysis of the homogeneous (non-spatial) system reveals that the zero steady state is unstable, and the state with u = α/β is stable (figure 1a). In the spatial system, a front connects the unstable uniform steady state u = 0 and the uniform non-zero steady state u = α/β (figure 1a,b). This front propagates into the unstable state with a speed c.
    Figure 1.

    Figure 1. (a) In a non-spatial system with logistic dynamics, there is one stable steady state. The state with zero population is unstable. (b) A front connects the two states and moves such that the populated state overtakes the unpopulated state. (c) The final speed depends on the steepness of the initial profile. Dashed lines: initial conditions, with different k. Solid lines: profile after t = 1.5 time units. Other parameters are α = 10, β = D = 1. A solution with steeper initial profile is slower. Any starting profile that is steeper than the blue curve will converge to a front with minimal speed cm=2αD. (d) Speed as a function of k, the exponential decay of the initial profile. The blue line denotes the analytical solution; the dots are obtained by explicit simulation using an Euler forward difference scheme. The numbers 1 and 2 correspond to the values of k used as initial profiles in (c). (Online version in colour.)

    A property of this system is that an infinite number of such fronts exist. In the 1970s, a number of authors studied types of initial data that can evolve to fronts with different speeds c. If the initial conditions satisfy [27]

    then, if k α/D, the system evolves to a front with a minimal speed
    However, if k<α/D, then the evolution is to a front having a higher speed c = Dk + α/k (figure 1c,d).

    Al-Omari & Gourley [10] considered the evolution of mature (adult) populations where the immature (juvenile) members of the population do not disperse. For many insect species, the juvenile phase is a larval phase in which the individuals do not move. Using an age structure model, they derive the following delay reaction–diffusion equation:

    where u represents the number of mature members of the population and the delay τ is the time taken from birth to maturity. Compared to equation (2.1), equation (2.3) provides a more realistic model of a single population by recognizing that individuals cannot reproduce right away from birth, but only after they have matured, which takes a time τ. It is assumed that the birth rate at any time is proportional to the number of adults at that time, so that the birth rate at time t is αu(t). The second term on the right-hand side of equation (2.3) represents the rate of recruitment into the adult population. This rate is essentially the birth rate, τ time units ago, decreased by the factor exp(γτ) accounting for mortality during the juvenile phase. The parameter γ thus measures the juvenile death rate while β measures the adult death rate. As for equation (2.1), monotone travelling front solutions connecting the two uniform equilibria exist for any speed c exceeding some critical (τ-dependent) minimum value c = cm.

    The transcendental equation for cm is documented by eqn (3.9) in [10] or eqn (2.14) in [28] and is given by

    The above equation has been solved numerically. The minimal speed cm = cm(τ) is shown in figure 2 for two different values of γ together with their local approximation
    Figure 2.

    Figure 2. Minimal speed above which a front connecting the two uniform steady states is possible. It has been obtained by solving numerically equation (2.4). Fixed parameters are α = β = D = 1 and γ = 0 or 0.5 (indicated in the figure). The straight lines are the small-τ approximations given by equation (2.5). The dots denote results obtained by numerical simulation of the equations using an Euler forward difference scheme, with as initial condition a step function. (Online version in colour.)

    The expression (2.5) clearly illustrates the fact that the minimal speed is decreasing with the delay. Equation (2.4) is obtained as follows: first, we switch to travelling wave coordinates by substituting z = x + ct in equation (2.3), then we linearize the resulting equation around u = 0 and examine the decay to zero of u(z) as z. A front that decays to zero in an oscillatory fashion is of no interest since u is a population and thus needs to remain positive. The interesting result from an ecological point of view is that the minimum speed decreases compared to equation (2.2) as the delay increases.

    3. Bistable travelling fronts

    (a) Population dynamics

    The situation is different if the population growth is damped by a strong Allee effect (after the American zoologist W.C. Allee). The Allee effect has to do with the fact that the fitness of small populations is sometimes negative, i.e. if the population density is too small, the species or group of individuals will not survive. In other words, the zero steady state is stable (figure 3a). To model this effect, an extra factor is added to the reaction term of the FKPP equation. Let us consider a population whose growth is described by a cubic polynomial

    where the ranges 0 < a < 1 and −1 < a < 0 correspond to the strong and weak Allee effects, respectively. For the strong Allee effect, there is a front connecting the two stable uniform steady states u = 0 and u = 1 (figure 3a,b). In contrast to monostable fronts, where an infinite number of speeds is possible for a front depending on the initial condition, this bistable front has a unique speed given by [29]
    Figure 3.

    Figure 3. (a) When a cubic function is used in the reaction part of the equation (equation (3.1)), three steady states can appear, of which two are stable. This picture corresponds to a strong Allee effect. (b) A front exists which connects the two stable steady states. The direction of the front depends on the parameter a. (Online version in colour.)

    This expression is obtained by inserting u=(1+exp(Az))1 into equation (3.1) reformulated in terms of the travelling wave coordinate z = x + ct. By solving this equation using as boundary conditions u( − ∞) = 0 and u(∞) = 1, we obtain A and c (equation (3.2)). By reversing the boundary conditions, A and c change sign.

    If a < 1/2, then c > 0, which implies that the invasion is successful: the populated steady state invades the unpopulated state. If instead a > 1/2, then c < 0. This means that the front travels in the other direction, and the low population density state wins out in the end. Interpreted in ecological terms, this means that an invasion fails.

    When there is a weak Allee effect, which can be modelled by taking −1 < a < 0, the zero steady state is unstable, but has a growth rate which is less than that at intermediate population densities. This case is a crossover between the strong Allee effect and logistic growth. There now is a front connecting the unstable zero steady state to the stable populated state, but determining the final propagation speed is more subtle. See [29, Sec. 3.3] for a discussion.

    (b) Optics

    The obvious question now is how these results for bistable fronts change when a time delay is introduced. We will now broaden our selection of topics, since travelling wave solutions of simple bistable reaction–diffusion equations subject to a time-delayed feedback are attracting the interest of many scientists with different objectives and experiences. In optics, the control of localized structures appearing in broad-area lasers is realized by optical feedback from a distant mirror. Small-area bistable lasing spots (10 μm) were first observed within the aperture of a broad-area vertical-cavity surface emitting laser (VCSEL) (80 μm) with frequency-selective feedback [30,31] and have immediately prompted theoretical studies [32]. The impact of time-delayed feedback on the dynamics of localized structures or cavity solitons has been further investigated theoretically for the cases of a driven nonlinear optical resonator [33,34] and broad-area VCSELs [3537]. Such patterns may be stationary or oscillatory, static or moving (see, for example, figure 4b). In [33], the following delayed Swift–Hohenberg (SH) equation

    was investigated numerically and analytically in two space dimensions. The quantities a, b, c are constant parameters, and η and τ are defined as the gain and the delay of the feedback, respectively. The time-delayed feedback does not modify the stationary time-independent state, but instead changes its stability properties [38,39]. An important result is that localized stationary solutions (spots) can be set into motion by the delayed feedback (figure 4b). It is the first delay-induced instability for the SH equation but other bifurcation scenarios are possible for higher-order reaction–diffusion systems [23]. The authors emphasized a necessary condition for moving spots given by the inequality
    This condition, which will later reappear, comes from the linear stability analysis around a stationary state that leads to the following characteristic equation for the growth rate λ:
    Figure 4.

    Figure 4. (a) The condition for a localized solution to start moving is given by τ > − 1/η and η < 0 (filled region). (b) Blue: localized stationary state. When τ is big enough, the bump starts moving (orange curve). The shape also becomes asymmetric depending on the direction of movement. Curves were obtained by simulating equation (3.3) using an Euler forward difference scheme, with a = 0, b = 1, c = 1.6 and η = − 1. (Online version in colour.)

    A necessary condition for stability is obtained by first determining the real roots. Keeping τ fixed, we analyse λ = λ(η) from equation (3.5) in the implicit form ητ=λτ(exp(λτ)1)1. We find that λ is unique and negative if

    Equation (3.5) admits no purely imaginary roots, and it can be shown that all the complex roots have a negative real part (the roots of λ=abexp(λτ) are studied in [40, p. 60]). The critical point
    is a bifurcation point that corresponds to a double-zero eigenvalue of the characteristic equation (3.5).

    (c) Neurobiology

    Bistable fronts appear in the neurosciences when one needs to describe the spreading of neuronal activities. One model that has been studied in this context is described by the following equation:

    where the last term is a feedback term, which depends on the delayed variable uτ. In [41], the feedback is a delayed spatial average defined by
    The length of the one-dimensional system is L = 1 and the constants a, η and τ are parameters. For ε = 0, this system admits a front solution
    which is stationary (figure 5).
    Figure 5.

    Figure 5. Hyperbolic tangent solution to equation (3.8), with ε = 0. The fact that a simple explicit formula exists for the leading-order problem facilitates the analysis. See also equation (3.12). (Online version in colour.)

    The influence of the time-delayed feedback on the shape and speed of the front given by equation (3.10) can be studied by taking advantage of the small parameter ε. In this way, Boubendir et al. found that the steady front may undergo a Hopf bifurcation to oscillatory travelling fronts [41].

    The exact form of the feedback term g(uτ) greatly influences the dynamics. A different global delayed feedback for a Ginzburg–Landau equation exhibiting a subcritical bifurcation is considered in [42]. They noted that, for a large delay, localized solutions can start oscillating. Equation (3.8) has been used as an early model for brain depression waves appearing for both migraine and stroke [43,44]. In those models, the delayed feedback was of the form

    and its effects for small feedback strength have been investigated in [45] for the symmetric case a = 0. Note that, in our present study, the feedback strength is given by εη, whereas in equation (3.3) it is η. In the presence of a time-delayed feedback (ε≠0), the steady front (3.10) starts to slowly move provided that εη < − 1/τ, the same conditions as for the Swift–Hohenberg case (equation (3.7), figure 4a). In this case, a strong enough feedback will set the front into motion. Since the system is symmetric, the front can travel in either direction with the same speed. We now generalize the results of [45] and explore the asymmetric case a≠0. Specifically, we seek a solution u = u(z) where z = x − ct that satisfies
    Prime means differentiation with respect to z. By taking advantage of the small parameter ε, we find that the leading solution is again given by (3.10). As in [45], we formulate the O(ε) problem and apply a solvability condition. The latter provides an equation for the speed c = c(τ) given by
    In this equation αexp(cτ). Figure 6a,b shows the solutions of equation (3.13) with η = − 1. The pitchfork bifurcation in orange corresponds to the case a = 0. The imperfect bifurcation in blue is for a = 0.01. In the language of the theory of singularities, the effect of εa≠0 is to unfold the pitchfork bifurcation. A local analysis near the bifurcation point at τ = ε−1 similar to the one in [45] indicates that the upper and lower branches are stable while the middle branch is unstable for τ > 1/ε. From (3.13), we note that the speed c approaches constant values as τ given by
    There are two important results that are worth stressing. First, since we take η = − 1 here, the two upper branches (c > 0) are coming closer together as a1 (figure 6c). Furthermore, their left limit point moves to large τ in the same limit a1. Consequently, if a > 1, the lower branch (c < 0) is the only branch of solutions. We conclude that a front that is too asymmetric may only be observed to be propagating in one direction (i.e. one specific sign of c). The delayed feedback then speeds up the front in that same direction. Second, the fact that the speed approaches constant values as τ implies that we may replace u(z + ) in equation (3.12) by its asymptotic limit for z large, i.e. u = 1. The reduced problem then becomes an ordinary differential equation which is easier to solve. Using (3.12) with u(z + ) = 1 we found analytically the correct speed limits.
    Figure 6.

    Figure 6. (a) Time evolution from a step-like initial profile for equation (3.12). Parameters: a = 0.01, η = − 1, τ = 25, ε = 0.1. For the moving front, the profile shape becomes asymmetric. (b) Bifurcation diagram obtained by solving equation (3.13). The orange curve is the symmetric case a = 0; the blue curve is obtained for a = 0.01. The parameter ε = 0.1. The dots denote solutions obtained by direct numerical simulation (forward Euler, a = 0.01). (c) Bifurcation curve for different values of a with ε = 0.1. The higher the asymmetry, the larger the delay needs to be in order to obtain a solution moving in the other direction. (Online version in colour.)

    The final example we discuss is again motivated by simple bistable models in neurobiology. The following equation exhibits a delayed threshold nonlinearity:

    and has been investigated in [46]. Here, H(u) denotes the Heaviside function, and the function H(u − a) − u can be seen as a piecewise linear version of a cubic polynomial with zeros at 0 and 1 (figure 7a). The advantage of using piecewise linear functions is that this approach allows one to find an exact solution for a stable front connecting u = 0 and u = 1. The speed c = c(τ) is again unique and satisfies the implicit equation
    Figure 7.

    Figure 7. (a) The cubic function used before for bistable fronts is replaced by a piecewise linear function. There are still two stable steady states. This form allows for an analytical solution. (b) Speed of the propagating front of equation (3.15). The broken line is the large-τ limit given in (3.17). The dots are the results obtained by direct numerical simulation (forward Euler). (Online version in colour.)

    where λ±=(c±c2+4)/2. In the limit of very small delays (τ0) and very large delays (τ), the front speeds go to cc0 and cc, respectively. These speeds are given by the following expressions:

    Figure 7b shows c as a function of τ. We note that the speed is decaying to zero as τ, which contrasts with our previous result using (3.11) as the feedback. There, the speed approaches constant values in the limit of large delays.

    4. Discussion

    Since the 1970s, there have been systematic studies of the impact of spatial diffusion on population dynamics when time lags are taken into account. More recently, new reaction–diffusion problems have appeared in neurobiology and nonlinear optics with the aim to understand how delays may control the progression of travelling waves. In this review, we emphasized the need for simple one-variable models in order to have an alternative to extensive numerical simulations. Key problems have been selected and various analytical techniques have been used to obtain results of physical significance.

    For the population model (2.3) and the delayed threshold model (3.15), the delay always slows down moving fronts. If the delayed feedback η(u(t − τ) − u) is added to the reaction–diffusion equation as for the Swift–Hohenberg equation (3.3) and the bistable equation (3.8) with (3.11) and a = 0, different responses are possible. For negative feedback strength (η < 0) and sufficiently large time delays τ, stationary solutions are set into motion. If the basic reaction–diffusion equation admits a slowly propagating front as equation (3.8) with (3.11) and a≠0, the delay may either speed up the front in the same direction or, interestingly, reverse its direction.

    While looking for the linear stability properties of the unstable uniform solution provides the propagation speed in the case of the monostable front, we need to solve the nonlinear problem for the bistable case. Preference has been given to asymptotic techniques. One novel observation is the unfolding of the pitchfork bifurcation illustrated in figure 6 for the asymmetric bistable equation (3.8) with (3.11) and a≠0. A second novel contribution of our work are the comparisons between analytical and numerical solutions which allow us to evaluate the quantitative validity of our approximations.

    Fronts are the building blocks of more complex structures such as pulses and wave trains, and asymptotic studies on the effects of a time delay on propagation speeds are currently being investigated. As explained by Gourley et al. [10], diffusion and time delays are not independent of each other since individuals have not been at the same point in space at previous times. One way to correct this is to introduce a non-local term in the reaction–diffusion equation, i.e. a weighted average in space. The effects of spatial non-locality without delay have recently been studied in different settings [4751]. Non-locality allows one to enhance or decrease front interaction and to create new structures.

    Besides population dynamics, epidemiology and neurobiology, there are systems in cell biology where it is relevant to study time-delayed reaction–diffusion equations. In these systems, a time delay can be used to model a chain of intermediate, and usually unknown, chemical reactions. One example of such a system is the early embryonic cell cycle of the frog Xenopus laevis. In these cells, travelling waves have been observed experimentally [52]. Similar to models in population biology, these waves could be interpreted as invasion waves, where one state (the mitotic state of the cell in which the cell divides) invades another state (interphase). This early cell cycle is controlled by a time-delayed negative feedback loop [53,54]. In future work, we plan to use and extend the techniques described here to study how such waves organize cell division in space and time.

    Data accessibility

    This article has no additional data.

    Authors' contributions

    J.R. performed the numerical simulations, helped with calculations and participated in writing. L.G. participated in interpreting results, structuring the content and writing the manuscript. T.E. initiated the research, did the literature review, made calculations and drafted the manuscript.

    Competing interests

    We have no competing interests.


    The authors acknowledge support by the Research-Foundation Flanders (FWO-Vlaanderen) for individual support (J.R.) and project support (L.G. grant G0A5217N).


    One contribution of 14 to a theme issue ‘Nonlinear dynamics of delay systems’.

    Published by the Royal Society. All rights reserved.