Next Article in Journal
On the Maximum Likelihood Estimators’ Uniqueness and Existence for Two Unitary Distributions: Analytically and Graphically, with Application
Next Article in Special Issue
Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction
Previous Article in Journal
Analyzing Dynamics: Lie Symmetry Approach to Bifurcation, Chaos, Multistability, and Solitons in Extended (3 + 1)-Dimensional Wave Equation
Previous Article in Special Issue
Molecular-Memory-Induced Counter-Intuitive Noise Attenuator in Protein Polymerization
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation

by
Ikram Medjahdi
1,
Fatima Zohra Lachachi
1,
María Ángeles Castro
1 and
Francisco Rodríguez
1,2,*
1
Department Applied Mathematics, University of Alicante, E-03080 Alicante, Spain
2
Multidisciplinary Institute for Environmental Studies (IMEM), University of Alicante, E-03080 Alicante, Spain
*
Author to whom correspondence should be addressed.
Symmetry 2024, 16(5), 609; https://doi.org/10.3390/sym16050609
Submission received: 29 March 2024 / Revised: 10 May 2024 / Accepted: 11 May 2024 / Published: 14 May 2024
(This article belongs to the Special Issue Mathematical Modeling in Biology and Life Sciences)

Abstract

:
The Klausmeier–Gray–Scott model of vegetation dynamics consists of a system of two partial differential equations relating plant growth and soil water. It is capable of reproducing the characteristic spatial patterns of vegetation found in plant ecosystems under water limitations. Recently, a discrete delay was incorporated into this model to account for the lag between water infiltration into the soil and the following water uptake by plants. In this work, we consider a more ecologically realistic distributed delay to relate plant growth and soil water availability and analyse the effects of different delay types on the dynamics of both mean-field and spatial Klausmeier–Gray–Scott models. We consider distributed delays based on Gamma kernels and use the so-called linear chain trick to analyse the stability of the uniformly vegetated equilibrium. It is shown that the presence of delays can lead to the loss of stability in the constant equilibrium and to a reduction of the parameter region where steady-state vegetation patterns can arise through symmetry-breaking by diffusion-driven instability. However, these effects depend on the type of delay, and they are absent for distributed delays with weak kernels when vegetation mortality is low.

1. Introduction

Vegetation in dryland ecosystems may exhibit complex dynamics, with alternative stable states prone to critical transitions between them and spatial discontinuities, resulting in different types of characteristic vegetation patterns [1,2,3,4,5]. Abrupt transitions, the so-called catastrophic shifts, from stable vegetated states to alternative bare soil or desert states, may occur when a parameter in the system crosses a tipping point due to the presence of a positive amplifying feedback [1,6].
Different properties of vegetation patterns and spatial metrics have been suggested as indicators of approaching tipping points [6,7,8,9,10,11]. It has also been discussed how self-organising vegetation patterns might help evade tipping points and system collapse, thereby maintaining ecosystem functioning and enhancing drylands’ resilience to worsening environmental conditions [12,13].
Many diverse mathematical models have been proposed to simulate and help analyse pattern formation and vegetation dynamics in drylands, including discrete-space cellular automata models (e.g., [14,15,16]) and, most often, continuous-space models based on systems of reaction-diffusion partial differential equations (e.g., [3,17,18,19,20,21]). In these models, patterns may arise through symmetry-breaking by Turing instability, i.e, the amplification in time of small perturbations from the homogeneous vegetated stable state by diffusion-driven instability [22,23,24,25,26,27], leading to the emergence of non-homogeneous steady-state vegetation distributions that may exhibit different types of spatial symmetry.
The Klausmeier model [3] is a classical reference in dryland vegetation modelling. It can reproduce spatial patterns similar to those found in different semiarid ecosystems using a relatively simple system of two partial differential equations that relate plant biomass ( N ) and soil water ( W ) on a bidimensional slope. Water is added uniformly at a rate A and is lost due to evaporation at a rate L W . Plants take up water at a rate R G ( W ) F ( N ) N , where G ( W ) and F ( N ) correspond, respectively, to the functional response of plants to water and how the presence of vegetation increments water infiltration. Linear functions are assumed in the model for simplicity: F ( N ) = N and G ( W ) = W . Water flows downhill in the X-direction at a speed V, and plant spreading is represented by a diffusion term with a diffusion coefficient D,
W T = A L W R W N 2 + V W X , N T = J R W N 2 M N + D 2 X 2 + 2 Y 2 N ,
where M represents density-independent mortality and J is the yield of plant biomass from uptaken water.
This classical model was later analysed, modified, and extended in different aspects (e.g., [28,29,30,31,32,33,34,35,36,37]). In the case of flat terrains, the gradient term in (1) can be replaced with a second diffusion term with a different diffusion constant. The resulting model is known as the Klausmeier–Gray–Scott model since it is equivalent to a chemical reaction model previously proposed by Gray and Scott [38,39,40].
In nondimensionalized form, the model can be written in terms of only four parameters: a and m, representing water input and plant mortality, and two diffusion coefficients, d 1 and d 2 , for water and plants. With new variables w and n, and transformed spatial and temporal variables, it reads
w t = a w w n 2 + d 1 Δ ( w ) , n t = w n 2 m n + d 2 Δ ( n ) ,
with adequate initial and boundary conditions on a given bounded region.
In a recent work [41], a discrete delay was introduced into this model to represent a lag in water infiltration,
w t = a w w n 2 + d 1 Δ ( w ) , n t = w ( x , t τ ) n ( x , t τ ) n m n + d 2 Δ ( n ) ,
which can be interpreted as plant growth being dependent on the water taken up by the plant at a fixed previous τ -lagged time. However, from a functional perspective, it seems more realistic to consider that plant growth would not be exclusively determined by water infiltrating instantly, either at the moment or in a certain previous time, but rather by water availability during a certain period. This could be represented by a weight function indicating how the importance for plant growth of previous water uptake changes with time, vanishing in the long term. Hence, in this paper, we consider a Klausmeier–Gray–Scott model with a distributed delay for soil water availability.
In the next section, we present the different models under study and analyse the stability of the vegetated equilibrium in both the full spatial models and the corresponding mean-field, non-spatial systems of ordinary differential equations. By considering distributed delays based on Gamma kernels and using the so-called linear chain trick [42,43], we will focus on the effects of the different delay types on the Turing space, i.e., the region of parameters where the onset of non-uniform steady-state patterned solutions is possible, with a more detailed analysis of the possibly more realistic case of a weak kernel of exponentially decaying lagging effects. The main conclusions are summarised and discussed in the last section.
We would like to emphasise the main contributions of this work: the introduction of a distributed delay in the dependence of plant growth on water uptake, which represents a more realistic situation than either instantaneous responses or fixed lag effects, and the results obtained in the next section showing how the type of delay may affect the stability of the homogeneous vegetated equilibrium in different ways. Specifically, there are no effects for an exponentially distributed delay at low mortality values, and the only effect on the onset of Turing instability is the shrinking of the Turing space.

2. Models and Results

2.1. Models without Delay

We first consider the stability properties of the original models without delay. While most results have already been discussed in the literature, we mainly summarise them without derivation. However, in different works, the analysis has usually been limited to certain ranges of the mortality parameter m or with slightly different nondimensionalized forms of the original model. For the sake of readability, we present the results for general values of mortality, along with schematic proofs to fix notations and to introduce expressions that are referred to below when analysing the models with delay.

2.1.1. Non-Spatial Model

The mean-field, non-spatial model corresponding to system (2) is given by the ordinary differential equation system
w t = a w w n 2 , n t = w n 2 m n .
This system has one, two, or three equilibria, depending on whether a < 2 m , a = 2 m , or a > 2 m , which are also the spatially uniform equilibria of (2). When a < 2 m , the only equilibrium is P 0 = a , 0 , the desert state, which is always stable. When a > 2 m , there are two more equilibria,
P 1 = a + a 2 4 m 2 2 , a a 2 4 m 2 2 m
and
P 2 = a a 2 4 m 2 2 , a + a 2 4 m 2 2 m ,
where P 1 is always unstable. These two vegetated equilibria coincide when a = 2 m at the point P 12 = m , 1 , which is unstable and where the system undergoes a fold bifurcation. Detailed stability analyses of these equilibria can be found in previous works (e.g., [36,37,41]).
We focus on the conditions for stability of P 2 , the uniformly vegetated state, summarised in the next lemma.
Lemma 1.
Consider system (4) with a > 2 m . The vegetated equilibrium P 2 is stable if m 2 or if m > 2 and a > a c , where
a c = m 2 m 1 ,
and it is unstable otherwise.
Proof. 
The linearised system at the equilibrium ( w * , n * ) is given by
w t = 1 + n * 2 w t 2 w * n * n t , n t = n * 2 w t + 2 n * w * m n t .
At the equilibrium P 2 , one has
n * = a + a 2 4 m 2 2 m
and w * n * = m . Thus, the characteristic equation is
p λ = λ 2 + 1 + n * 2 m λ + m n * 2 1 = 0 .
Since n * > 1 , the only condition for all roots of (7) to have negative real parts is
1 + n * 2 m > 0 ,
which is immediately satisfied when m 2 , whereas for m > 2 , from (6), it is equivalent to a > a c , as given in (5). □
The results of Lemma 1 are presented in Figure 1, showing the regions of stability and instability of P 2 in the ( m , a ) -parameter space.
As recalled in the next lemma, there is a Hopf bifurcation at the critical value a c (see (Theorem 2.1 in [37])).
Lemma 2.
Consider system (4) with a > 2 m . For m > 2 , the system undergoes a Hopf bifurcation at a = a c .
Figure 2 shows a phase portrait and a bifurcation diagram illustrating the results in Lemmas 1 and 2. For m > 2 , local perturbations of P 2 lead the system to either P 2 or to the desert state P 0 , depending on whether a is greater or lower than the critical value a c , respectively.

2.1.2. Spatial Model

To analyse the spatial model (2), the problem has to be defined in a given domain with appropriate boundary conditions. For simplicity, we consider the one-dimensional case, with the domain being the interval [ 0 , π ] and Neumann no-flux boundary conditions,
w x ( 0 , t ) = n x ( 0 , t ) = w x ( π , t ) = n x ( π , t ) = 0 , t 0 .
However, the problem in a general domain Ω , usually in R or R 2 , can be dealt with similarly by considering the corresponding spatial eigenvalue problem [25]. A stability analysis of (2) in a general domain for m < 2 was presented in [36]. In a general interval [ 0 , l π ] for m > 2 , with a slightly different nondimensionalized form of (1), a stability analysis was presented in [37].
Linearising system (2) at the equilibrium ( w * , n * ) and considering perturbations of the form e λ t cos ( k π ) , according to the Neumann boundary conditions, where k is the wavenumber, and writing μ k = k 2 one obtains the set of characteristics equations
p k λ = λ 2 + S k 1 λ + S k 2 = 0 ,
with S k 1 = d 1 + d 2 μ k + 1 + n * 2 m , and
S k 2 = d 1 d 2 μ k 2 + d 2 1 + n * 2 d 1 m μ k + m n * 2 1 .
We focus on the conditions for Turing instability to occur, that is, the conditions for the uniformly vegetated equilibrium P 2 to be stable in the non-spatial model and unstable in the spatial one [25]. These conditions are given in the next lemma.
Lemma 3.
Consider system (2) with a > 2 m , and assume that the equilibrium P 2 is stable in the non-spatial system (4). Then, writing r = d 1 / d 2 , the spatially uniform equilibrium P 2 is unstable if and only if
r > r c = 3 n * 2 1 + 2 n * 2 ( n * 2 1 ) m
and there exists k 1 such that
l k 2 l + ,
where
l = ( r m 1 n * 2 ) ( r m 1 n * 2 ) 2 4 r m ( n * 2 1 ) 2 d 1 .
Proof. 
For k = 0 , the characteristic Equation (10) reduces to that of the non-spatial problem (7). Thus, since we are assuming that P 2 is stable in (4), one has that S 0 1 and S 0 2 are both positive. Since S 0 1 > 0 implies that S k 1 > 0 for all k, the system is unstable at P 2 if and only if S k 2 is negative, and for this to happen, one must have
d 2 1 + n * 2 d 1 m < 0 ,
that is,
r = d 1 d 2 > 1 + n * 2 m ,
and
d 2 1 + n * 2 d 1 m 2 4 d 1 d 2 m n * 2 1 > 0 ,
or, equivalently,
m 2 r 2 m ( 3 n * 2 1 ) r ( 1 + n * 2 ) 2 > 0 ,
which holds when r satisfies (12) or when
0 r < 3 n * 2 1 2 n * 2 ( n * 2 1 ) m ,
but in the latter case, condition (15) is not satisfied. Since n * > 1 for P 2 , one has
3 n * 2 1 2 n * 2 ( n * 2 1 ) m < 1 + n * 2 m < 3 n * 2 1 + 2 n * 2 ( n * 2 1 ) m ,
and therefore both (15) and (16) are fulfilled if and only if condition (12) holds.
Besides condition (12), to ensure that S k 2 is negative, since { μ k } = { k 2 } is a discrete sequence, it must also hold that some μ k lies between the two positive roots of the second-order equation
d 1 d 2 μ 2 + d 2 1 + n * 2 d 1 m μ + m n * 2 1 = 0 ,
which are given by l in (14). □
Remark 1.
If r m 1 n * 2 is sufficiently small, there is no wavenumber satisfying (13), and the spatially homogeneous equilibrium is stable. Thus, as is well known [25], pattern formation by diffusion-driven instability requires the ratio r of diffusion coefficients to be sufficiently high. This is to be expected in real situations since water diffusion, with coefficient d 1 , is expected to be higher than the spread of plants represented by their diffusion term, with coefficient d 2 . When Turing instability occurs, the eigenfunctions corresponding to the set of wavenumbers satisfying (13) determine the spatially heterogeneous, patterned, steady-state solution.

2.2. Models with Distributed Delay

Now, in (2), we introduce a distributed delay in the product w n , corresponding to water uptake by unit plant biomass, to account for the dependence of plant growth on previously available water,
w t = a w ( x , t ) w ( x , t ) n ( x , t ) 2 + d 1 Δ ( w ) , n t = n ( x , t ) t g α p t s w x , s n x , s d s m n ( x , t ) + d 2 Δ ( n ) ,
where g α p is a Gamma-distributed kernel,
g α p t = α Γ p α t p 1 e α t .
This kernel depends on two parameters, p Z + and α > 0 , where the average delay is p / α with variance p / α 2 . Thus, to represent a mean delay τ , we can take α = p / τ , with variance τ / p . In this way, increasing p concentrates the distribution more around τ , approaching the model with discrete delay (3) as p tends to infinity (Figure 3).
We analyse in more detail the case of a weak kernel, p = 1 , where the effect on plant growth of previous water uptake decays exponentially, as it is likely the most realistic case. However, we also discuss the differences with more concentrated weight distributions, strong kernels with p > 1 , with a specific analysis of the case p = 2 .

2.2.1. Non-Spatial Model with a Weak Kernel

It is clear that (17) and (2) have the same constant equilibria, as well as (4) and the non-spatial model corresponding to (17),
w t = a w ( t ) w ( t ) n ( t ) 2 , n t = n ( t ) t g α p t s w s n s d s m n ( t ) .
To analyse the stability of P 2 , we use the so-called linear chain trick [42,43]. For p = 1 , we introduce one new variable,
z ( t ) = t g α 1 ( t s ) w ( s ) n ( s ) d s ,
and consider the extended system
w t = a w ( t ) w ( t ) n ( t ) 2 , n t = n ( t ) z ( t ) m n ( t ) , z t = α w ( t ) n ( t ) z ( t ) .
If P = ( w * , n * ) is an equilibrium of (18), then P = ( w * , n * , w * n * ) is the corresponding equilibrium in (19), and we drop the third component in what follows. The next theorem states the conditions for the uniform equilibrium P 2 to be stable in (19).
Theorem 1.
Consider system (19) with a > 2 m . The vegetated equilibrium P 2 is stable independently of α if m 2 or if m > 2 and a > a c α = m 2 m , and it is unstable if m > 2 and a < a c , as defined in (5).
When m > 2 and a c < a < a c α , the stability depends on the value of α, being stable for α > α c and unstable for α < α c , where
α c = 2 m n * 2 n * 2 + 1 2 1 + n * 2 m .
Proof. 
The linearised system of (19) at the equilibrium P = ( w * , n * , w * n * ) is given by
w t = ( 1 + n * 2 ) w ( t ) 2 w * n * n ( t ) , n t = ( w * n * m ) n ( t ) + n * z ( t ) , z t = α n * w ( t ) + α w * n ( t ) α z ( t ) .
Thus, at P 2 , where w * n * = m , one obtains the characteristic equation
p α λ = λ 3 + 1 + n * 2 + α λ 2 + α 1 + n * 2 m λ + α m n * 2 1 = 0 .
Note that writing q α λ = 1 α p α λ , one has lim α + q α λ = p λ , as given in (7) for the non-spatial model without delay.
Since n * > 1 at P 2 , the conditions for all the roots of (22) to have negative real parts are 1 + n * 2 m > 0 , i.e., condition (8), which holds when m 2 or when m > 2 and a > a c , and
1 + n * 2 + α α 1 + n * 2 m α m n * 2 1 > 0 .
When m 2 , condition (23) is satisfied independently of α since
1 + n * 2 + α 1 + n * 2 m m n * 2 1 1 + n * 2 + α m n * 2 1 > 0 .
For m > 2 , condition (23) might be dependent on α . Since α > 0 , it is equivalent to
α > m n * 2 1 1 + n * 2 m 1 + n * 2 = 2 m n * 2 1 + n * 2 2 1 + n * 2 m ,
that is, α > α c , as defined in (20), which holds for any value of α if α c 0 .
Taking into account the expression of n * in terms of a and m in P 2 , (6), with some algebraic manipulation, it can be checked that α c > 0 a < a c α = m 2 m . □
Figure 4 shows the region in the ( m , a ) -parameter space where the stability of P 2 depends on the value of α , C ( α ) = { ( m , a ) | m > 2 , a c < a < a c α } , as given in Theorem 1.
As shown in the next theorem, when the stability of P 2 depends on α , there is a Hopf bifurcation at α = α c .
Theorem 2.
Consider system (19) with m > 2 and a c < a < a c α . The system undergoes a Hopf bifurcation at α = α c .
Proof. 
When α = α c the characteristic Equation (22) has one pair of conjugate pure imaginary roots, λ 1 , 2 = ± i α c 1 + n * 2 m , and one real negative root, λ 3 = 1 + n * 2 + α c . To check the transversality condition for a Hopf bifurcation to exist, write (22) in the form
λ 3 + 1 + n * 2 λ 2 + α q λ = 0 ,
where
q λ = λ 2 + n * 2 m + 1 λ + m n * 2 1 .
Differentiating with respect to α in (25), one has
3 λ 2 + 1 + n * 2 2 λ + α q λ d λ d α = q λ ,
so
d λ d α 1 = 3 λ 2 + 1 + α + n * 2 2 λ + α n * 2 m + 1 q ( λ ) .
For λ = λ 1 = i α c n * 2 m + 1 , one has
q λ 1 = n * 2 m + 1 n 2 + 1 + i 2 m n * 2 n * 2 + 1 2 ,
and
q λ 1 n 2 + 1 i 2 m n * 2 n * 2 + 1 2 = n * 2 m + 1 2 m n * 2 .
Writing s λ for the numerator of (26), for α = α c , one has
s λ 1 = 2 n * 2 + 1 2 2 m n * 2 + i 2 m n * 2 1 n * 2 m + 1 2 m n * 2 n * 2 + 1 2 ,
so
s λ 1 n 2 + 1 i 2 m n * 2 n * 2 + 1 2 = 2 n * 2 + 1 2 2 m n * 2 2 n * 2 m + 1 .
Hence, from (26)–(28), one obtains
d λ d α α = α c 1 = n * 2 + 1 2 2 m n * 2 2 m n * 2 n * 2 m + 1 2 < 0 .
Figure 5 illustrates the results of Theorems 1 and 2. For the example in this figure, with m = 6 , one has a c 16.10 and a c α 20.78 . Thus, a c < a = 18 < a c α , and the stability of P 2 depends on α being greater or lower than the critical value defined in Theorem 1, α c 11.09 . The equilibrium values at P 2 are n * 2.62 and w * 2.29 , with the system converging to the equilibrium under small perturbations when α > α c (Figure 5a). For α < α c , the equilibrium P 2 is unstable, and the system oscillates with increasing amplitudes for decreasing α values (Figure 5b,c), converging to the desert state P 0 for lower α values (Figure 5d), corresponding to larger mean delays.
Remark 2.
Considering the linearised system (21) at the corresponding equilibria, it is not difficult to see that P 0 is always stable and P 1 is always unstable, with no effect of the distributed delay on their stability.

2.2.2. Spatial Model with a Weak Kernel

Now, we introduce into (2) the new variable
z ( x , t ) = t g α 1 ( t s ) w ( x , s ) n ( x , s ) d s ,
and the linearised extended system at the equilibrium P = ( w * , n * , w * n * ) is given by
w t = ( 1 + n * 2 ) w ( x , t ) 2 w * n * n ( x , t ) + d 1 Δ ( w ) , n t = ( w * n * m ) n ( x , t ) + n * z ( x , t ) + d 2 Δ ( n ) , z t = α n * w ( x , t ) + α w * n ( x , t ) α z ( x , t ) .
Considering perturbations as in the spatial model without delay, and with μ k = k 2 , one obtains the set of characteristics equations
p k α λ = λ 3 + l k 1 λ 2 + l k 2 λ + l k 3 ,
with
l k 1 = d 1 + d 2 μ k + 1 + n * 2 + α ,
l k 2 = d 1 d 2 μ k 2 + d 2 1 + n * 2 μ k + α S k 1 ,
and l k 3 = α S k 2 , where S k 1 and S k 2 are the coefficients of the spatial model without delay (18). We note that by letting q k α λ = 1 α p k α λ , one has lim α + q k α λ = p k λ , and the characteristic equation of the spatial model without delay is recovered.
As shown in the next theorem, when the non-spatial model with delay is stable, the conditions for Turing instability to occur are the same as in the spatial model without delay.
Theorem 3.
Consider system (17) with a > 2 m , and assume that the equilibrium P 2 is stable in the non-spatial system (18). Under these conditions, P 2 is unstable if and only if it is unstable in the spatial model without delay (2).
Proof. 
First, we show that if P 2 is unstable in the model without delay (2), it is also unstable in (17). Since we are assuming that it is stable in the non-spatial model with delay (18), it holds that S 0 1 > 0 . Hence, since S k 1 is increasing with k, it also holds that S k 1 > 0 for all k. Thus, for P 2 to be unstable in (2), S k 2 has to be negative, which implies that l k 3 = α S k 2 is also negative. Therefore, P 2 is unstable in (17).
Assume now that P 2 is stable in both (18) and (2). Then, S k 1 and S k 2 are both positive, so all the coefficients in (30) are also positive. The only condition for stability that has to be checked is the positivity of l k 1 l k 2 l k 3 , which can be written in the form
b 0 ( μ k ) α 2 + b 1 ( μ k ) α + b 2 ( μ k ) ,
where b 0 ( μ k ) = S k 1 ,
b 1 ( μ k ) = d 1 + d 2 μ k + n * 2 + 1 S k 1 + d 1 d 2 μ k 2 + d 2 μ k n * 2 + 1 S k 2 = d 1 + d 2 2 μ k 2 + d 1 m + ( 2 + 2 n * 2 m ) d 1 + d 2 μ k + n * 4 + 2 1 m n * 2 + 1 ,
and
b 2 ( μ k ) = d 2 μ k d 1 + d 2 μ k + n * 2 + 1 d 1 μ k + n * 2 + 1 .
Thus, b 0 ( μ k ) and b 2 ( μ k ) are always positive, and b 1 ( μ k ) is increasing with k, so it will be positive for any k if
b 1 ( 0 ) = n * 4 + 2 1 m n * 2 + 1 > 0 .
It is easy to see that (35) is positive when m 2 , and it can be shown, with some algebraic manipulation, that it is also positive when m > 2 and a > a c α as defined in Theorem 1. Hence, since we are assuming that P 2 is stable in (18), the only remaining case where the positivity of (33) has to be checked is when a c < a < a c α and α > α c , as given in Theorem 1.
If a < a c α , then b 1 ( μ k ) < 0 , and (33) can only be negative if
b 1 μ k 2 4 b 0 μ k b 2 μ k > 0 ,
so there are two real roots, α 1 < α 2 and α 1 < α < α 2 . However,
α 2 = b 1 μ k + b 1 μ k 2 4 b 0 μ k b 2 μ k 2 b 0 μ k ,
and
b 1 μ k 2 4 b 0 μ k b 2 μ k < b 1 μ k = b 1 μ k ,
so since 0 < b 1 ( μ k ) < b 1 ( 0 ) and 0 < b 0 ( 0 ) < b 0 ( μ k ) ,
α 2 < b 1 μ k b 1 μ k 2 b 0 μ k = b 1 μ k b 0 μ k < b 1 0 b 0 0 = α c .
Hence, when α > α c , (33) is also positive. Therefore, P 2 is stable in (17) if it is stable in (18) and (2). □
Remark 3.
Similarly to Theorem 2 for the non-spatial model with delay (18), it can be shown that for m > 2 and a c < a < a c α , there is a Hopf bifurcation at α = α c in the spatial system (17), with oscillatory unstable and stable solutions for α < α c and α > α c , respectively.
Figure 6 shows an example where the spatially homogeneous solution is stable in the model without delay since r = d 1 / d 2 = 1 < r c 3.95 , so the instability condition (12) in Lemma 3 is not satisfied. For m = 2.95 , one obtains a c 6.23 and a c α 7.16 , so a c < a = 6.5 < a c α , and α c 5.15 . The equilibrium P 2 is unstable for α < α c (Figure 6a) with mean delay 1 / α = 0.4 and stable for α > α c (Figure 6b) with mean delay 1 / α = 0.125 .
Remark 4.
As shown in Theorem 3, when the homogeneous vegetated equilibrium is stable in the non-spatial system (18), spatially non-homogeneous steady-state solutions may arise from Turing bifurcations under the same conditions as in the system without delay. The presence of a delay does not affect these steady-state solutions, in the same way that the presence of a delay does not affect the homogeneous equilibria.
Figure 7 shows two examples of non-homogeneous solutions for m 2 (left panels), where the onset of Turing instability is not affected by the delay in system (17), and for m > 2 and a c < a < a c α , where the Turing space is restricted by the condition α > α c . In the first example, with m = 1 , a = 2.2 , d 1 = 50 , and d 2 = 0.1 , the limits given in (14) are l 0.03 and l + 9.90 , so the wavenumbers with positive eigenvalues are μ k = k 2 for k = 1 3 . In the second example, with m = 2.5 , a = 5.25 , d 1 = 50 , and d 2 = 0.1 , the limits are l 0.02 and l + 24.92 , so the wavenumbers with unstable modes are μ k = k 2 for k = 1 4 .

2.2.3. Models with Strong Kernels

Next, we discuss how the type of delay may differentially affect the stability of P 2 . We provide a detailed analysis only for the model with p = 2 since models with strong kernels are likely to be less realistic than an exponentially decaying kernel for modelling the dependence of plant growth on previous water uptake. Notwithstanding, a Gamma kernel with p = 2 and a small mean delay might also represent a plausible modelling option, combining an essentially exponential-like decay with avoiding completely instantaneous responses (compare, for instance, the strong kernel with p = 2 and α = 2 with the weak kernel with p = 1 and α = 1 in Figure 3a).
Now, consider model (18) with p = 2 . Introducing the two new variables
z i ( t ) = t g α i ( t s ) w ( s ) n ( s ) d s , i = 1 , 2 ,
we can write the extended system
d w ( t ) d t = a w ( t ) w ( t ) n ( t ) 2 , d n ( t ) d t = n ( t ) z 2 ( t ) m n ( t ) , d z 1 ( t ) d t = α w ( t ) n ( t ) z 1 ( t ) , d z 2 ( t ) d t = α z 1 ( t ) z 2 ( t ) .
In this model, z 2 represents the effect of previous water uptake on plant growth, similar to z in model (29), but with a different weighted function, g α 2 instead of g α 1 . Meanwhile, z 1 is an auxiliary variable.
The next theorem states the conditions for the uniform equilibrium P 2 to be stable in (37).
Theorem 4.
Consider system (37) with a > 2 m . The vegetated equilibrium P 2 is unstable if m > 2 and a < a c , as defined in (5). When m 2 or m > 2 and a > a c , the stability depends on the value of α, being stable for α > α c , 2 and unstable for α < α c , 2 , where α c , 2 is the only positive root of
α ( 1 + n * 2 m ) 2 ( 1 + n * 2 + α ) 2 + α m m ( n * 2 1 ) ( 1 + n * 2 + 2 α ) 2 = 0 .
Proof. 
From the linearised system at P 2 , one obtains the characteristic equation p α , 2 = 0 , where
p α , 2 λ = ( λ + α ) 2 λ λ + 1 + n * 2 α 2 m λ + 1 n * 2 ,
that is,
p α , 2 λ = λ 4 + ( 1 + n * 2 + 2 α ) λ 3 + 2 α ( 1 + n * 2 ) + α 2 λ 2 + α 2 1 + n * 2 m λ + m α 2 ( n * 2 1 ) .
Since n * 2 > 1 at P 2 , all coefficients of p α , 2 are positive if 1 + n * 2 m > 0 . This condition, as given in (8), determines the stability of P 2 in the corresponding model without delay. This condition fails, and the system is unstable for m > 2 and a < a c .
There is one more condition for all the roots of (39) to have negative real parts, which can be obtained using the Hurwitz criterion, Q ( α ) > 0 , where
Q ( α ) = α ( 1 + n * 2 m ) 2 ( 1 + n * 2 + α ) 2 + α m m ( n * 2 1 ) ( 1 + n * 2 + 2 α ) 2 ,
or, equivalently,
Q ( α ) = A 0 α 3 + A 1 α 2 + A 2 α + A 3 ,
where A 0 = 1 + n * 2 m , A 1 = ( 1 + n * 2 m ) m + 4 ( 1 + n * 2 ) 4 m ( n * 2 1 ) ,
A 2 = ( 1 + n * 2 ) 2 ( 1 + n * 2 m ) ( 1 + n * 2 ) 4 m ( n * 2 1 ) ,
and A 3 = m ( 1 + n * 2 ) 2 ( n * 2 1 ) .
When m 2 or m > 2 and a > a c , condition (8) holds. Since n * > 1 at P 2 , one has A 0 > 0 and A 3 < 0 , which implies that Q ( α ) has at least one positive root. If A 2 < 0 , there is only one sign change in the coefficients of Q ( α ) , and, according to Descartes’ rule of signs, there is only one positive root. If A 2 > 0 , then 4 m ( n * 2 1 ) < 2 ( 1 + n * 2 m ) ( 1 + n * 2 ) , which implies
A 1 > ( 1 + n * 2 m ) m + 4 ( 1 + n * 2 ) 2 ( 1 + n * 2 ) > 0 .
Thus, in this case, there is also only one sign change, and therefore Q ( α ) has only one positive root, which defines the critical value α c , 2 . □
We note that since α is a factor in the first term of (41) and the second term is negative, condition Q ( α ) > 0 may fail for any values of m and a when α is sufficiently small. Thus, unlike in the model with a weak kernel, the delay can make the vegetated equilibrium P 2 unstable in the whole range of the ( m , a ) -parameter space and, in particular, when m 2 .
The next theorem shows the result corresponding to Theorem 2 for system (37).
Theorem 5.
Consider system (37) with m 2 or m > 2 and a > a c . The system undergoes a Hopf bifurcation at α = α c , 2 .
Proof. 
Let λ = i w , w > 0 , be a pure imaginary root of p α , 2 . Then, separating the real and imaginary parts, one obtains
w 4 α 2 n * 2 + 2 + α w 2 + α 2 m n * 2 1 = 0
and
w 2 α + 1 + n * 2 w 2 + α 2 1 + n * 2 m = 0 .
From (43), one has
w 2 = α 2 1 + n * 2 m 2 α + 1 + n * 2
and substituting this expression into (42) and simplifying, one obtains
2 1 + n * 2 m α 3 + 4 n * 4 + 8 7 m n * 2 m m 1 + 4 α 2        + 2 n * 2 + 1 n * 4 + 2 3 m n * 2 + m + 1 α m n * 2 1 n * 2 + 1 2 = 0 ,
which can be shown to be equivalent to Q ( α ) = 0 , that is, α = α c , 2 .
The transversality condition for a Hopf bifurcation to exist at α = α c , 2 can be checked similarly to the case p = 1 , although it involves more cumbersome computations.
Writing p α , 2 ( λ ) in the form
p α , 2 λ = λ 4 + 1 + n * 2 + 2 α λ 3 + 2 α 1 + n * 2 λ 2 + α 2 q λ ,
with
q λ = λ 2 + 1 + n * 2 m λ + m n * 2 1 ,
one obtains
λ α 1 = 4 λ 3 + 3 1 + n * 2 + 2 α λ 2 + 4 α 1 + n * 2 λ + α 2 q λ 2 λ 3 + 2 1 + n * 2 λ 2 + 2 α q λ .
Substituting q λ = 2 λ + 1 + n * 2 m into this expression, and setting λ = i w while evaluating at α = α c , 2 , one obtains, after some rearrangements and simplifications and using (43) and the expression defining α c , 2 , that is, Q ( α c , 2 ) = 0 , the following:
λ α α = α c , 2 1 = 1 2 R 1 + i I 1 R 2 + i I 2 = 1 2 R 1 R 2 + I 1 I 2 + i I 1 R 2 R 1 I 2 R 2 2 + I 2 2 ,
where
R 1 = 3 1 + n * 2 + 2 α c , 2 w 2 + α c , 2 2 1 + n * 2 m = 2 α c , 2 2 1 + n * 2 m , I 1 = 2 α c , 2 2 1 + n * 2 + α c , 2 4 w 2 w = 2 α c , 2 2 n * 4 + 4 + 3 α c , 2 n * 2 + 2 α c , 2 2 + 3 + 2 m α c , 2 + 2 2 α c , 2 + 1 + n * 2 w , R 2 = α c , 2 m n * 2 1 1 + n * 2 + α w 2 = α c , 2 2 1 + n * 2 m 2 1 + n * 2 + α c , 2 2 + α c , 2 m 1 + n * 2 + 2 α c , 2 1 + n * 2 + α 2 α c , 2 + 1 + n * 2 2 = α c , 2 2 1 + n * 2 m 1 + n * 2 + m α c , 2 + n * 2 + 1 2 2 α c , 2 + 1 + n * 2 2 , I 2 = α c , 2 1 + n * 2 m w 2 w = α c , 2 1 + n * 2 m α c , 2 + 1 + n * 2 2 α c , 2 + 1 + n * 2 w .
Since
R 1 R 2 + I 1 I 2 = 2 α c , 2 4 1 + n * 2 m 2 1 + n * 2 + m α c , 2 + n * 2 + 1 2 2 α c , 2 + 1 + n * 2 2 + 2 α c , 2 4 2 n * 4 + 4 + 3 α c , 2 n * 2 + 2 α c , 2 2 + 3 + 2 m α c , 2 + 2 1 + n * 2 m 2 α c , 2 + 1 + n * 2 2 α c , 2 + 1 + n * 2 3 = 2 α c , 2 4 1 + n * 2 m 2 2 α c , 2 + 1 + n * 2 3 2 n * 4 + 4 + 3 α c , 2 n * 2 + 2 α c , 2 2 + 3 + 2 m α c , 2 + 2 α c , 2 + 1 + n * 2 2 α c , 2 + 1 + n * 2 1 + n * 2 + m α c , 2 + n * 2 + 1 2 = 2 α c , 2 4 1 + n * 2 m 2 2 α c , 2 + 1 + n * 2 3 n * 6 + 3 + 2 α c , 2 n * 4 + 3 α c , 2 2 + 4 + m α c , 2 + 3 n * 2 + 2 α c , 2 3 + 3 α c , 2 2 + 2 + m α c , 2 + 1 > 0 ,
it holds that
λ α α = α c , 2 1 < 0 .
Figure 8 illustrates the results of Theorems 4 and 5, showing changes in the stability of P 2 through Hopf bifurcations at α = α c , 2 , in parameter regions where P 2 is stable independently of α , in the corresponding model with p = 1 . The top panels correspond to the region m 2 ( m = 1.5 , a = 4 ), with α c , 2 0.83 . The lower panels correspond to the region m > 2 and a > a c α 7.35 ( m = 3 , a = 8 ), with α c , 2 3.62 . In both cases, the left panels correspond to α = 1.1 α c , 2 , with the system converging to P 2 equilibrium values, while the right panels correspond to α = 0.9 α c , 2 , with the system converging to periodic solutions. Lower α values, that is, larger mean delays, would drive the system to the desert state P 0 , similar to the behaviour shown in Figure 5.
The next theorem is analogous to Theorem 3 for the strong kernel case with p = 2 .
Theorem 6.
For p = 2 , consider system (17) with a > 2 m , and assume that the equilibrium P 2 is stable in the corresponding non-spatial system (18). Under these conditions, P 2 is unstable if and only if it is unstable in the spatial model without delay (2).
Proof. 
Similar to the non-spatial case, introducing the new variables
z i ( x , t ) = t g α i ( t s ) w ( x , s ) n ( x , s ) d s , i = 1 , 2 ,
we transform (17) into the extended system
d w ( t ) d t = a w ( x , t ) w ( x , t ) n ( x , t ) 2 + d 1 Δ ( w ) , d n ( t ) d t = n ( x , t ) z 2 ( x , t ) m n ( x , t ) + d 2 Δ ( n ) , d z 1 ( t ) d t = α w ( x , t ) n ( x , t ) z 1 ( x , t ) , d z 2 ( t ) d t = α z 1 ( x , t ) z 2 ( x , t ) .
From the linearised system at equilibrium, and considering perturbations as in the spatial models without delay and with distributed delay with a weak kernel, writing μ k = k 2 , one obtains the set of characteristics equations
p k α , 2 λ = λ 4 + Φ k 1 λ 3 + Φ k 2 λ 2 + Φ k 3 λ + Φ k 4 ,
where
Φ k 1 = ( d 1 + d 2 ) μ k + 1 + n * 2 + 2 α , Φ k 2 = d 1 d 2 μ k 2 + d 2 n * 2 + 1 μ k + 2 α S k 1 + α 2 m + α , Φ k 3 = 2 α d 2 μ k n * 2 + d 1 μ k + 1 + α 2 S k 1 ,
and Φ k 4 = α 2 S k 2 , with S k 1 and S k 2 being the coefficients of the characteristic equation for the spatial model without delay (10).
The arguments are now similar to those used in the proof of Theorem 3. If P 2 is stable in the non-spatial system (18) and unstable in the spatial model without delay (2), it must hold that S k 2 is negative and, consequently, Φ k 4 < 0 , so P 2 is also unstable in (47).
Assume now that P 2 is stable in both (18) and (2). Then, S k 1 and S k 2 are both positive, so all the coefficients in (48) are also positive. The only condition for stability that has to be checked is the positivity of Φ k 1 Φ k 2 Φ k 3 Φ k 3 2 Φ k 1 2 Φ k 4 , which can be written in the form
Φ k 1 Φ k 2 Φ k 3 Φ k 3 2 Φ k 1 2 Φ k 4 = α Q ( α ) + C 0 α 4 + C 1 α 3 + C 2 α 2 + C 3 α + C 4 ,
with Q ( α ) given in (41), and where C 0 = 2 d 1 + d 2 μ k ,
C 1 = μ k 4 d 1 + d 2 2 μ k + d 1 8 n * 2 + 8 + m + d 2 8 + 8 n * 2 3 m = μ k 4 d 1 + d 2 2 μ k + d 1 8 n * 2 + 1 + m + d 2 5 1 + n * 2 + 3 1 + n * 2 m , C 2 = μ k 2 d 1 + d 2 d 1 2 + 4 d 1 d 1 + d 2 2 μ k 2 + 2 d 1 2 3 1 + n * 2 + m + 2 d 1 d 2 10 1 + n * 2 + m + 8 d 2 2 1 + n * 2 + 2 d 2 2 1 + n * 2 m μ k + 2 d 1 3 + 4 n * 2 + n * 4 + 2 m + 2 d 2 5 + 7 n * 2 + 2 n * 4 + m + 2 n * 2 3 d 2 + 2 d 1 1 + n * 2 m , C 3 = μ k 4 d 1 d 2 d 1 + d 2 2 μ k 3 + d 1 + d 2 d 1 2 m + 4 d 2 3 d 1 + d 2 1 + n * 2 μ k 2 + 2 d 2 2 4 + 7 n * 2 + 3 n * 4 + 2 d 1 d 2 6 + m + 5 n * 4 + 11 n * 2 + d 1 2 m n * 2 + 3 + 2 d 2 n * 2 d 2 + d 1 1 + n * 2 m μ k + n * 2 + 1 d 2 n * 4 + 5 n * 2 + 4 + m + d 1 3 n * 2 m + 3 n * 2 d 2 1 + n * 2 m ,
and C 4 = 2 μ k 2 d 2 d 1 + d 2 μ k + 1 + n * 2 d 1 μ k + n * 2 + 1 2 .
Since P 2 is stable in the non-spatial model, it holds that Q ( α ) and ( 1 + n * 2 m ) are both positive, so all the terms on the right-hand side of (49), except perhaps C 3 , are clearly positive. But for C 3 , the only possibly non-positive terms is
d 2 n * 4 + 5 n * 2 + 4 + m + d 1 3 n * 2 m + 3 n * 2 d 2 1 + n * 2 m ,
which is greater than
n * 2 d 2 1 + n * 2 + d 1 m 3 n * 2 .
Since we are assuming that P 2 is stable in the spatial system without delay, from (15), it must hold that d 2 ( 1 + n * 2 ) d 1 m , so
n * 2 d 2 1 + n * 2 + d 1 m 3 n * 2 d 1 m n * 2 + 3 n * 2 > 0 .
Remark 5.
We do not pursue the analysis here, but it can be shown that for any p Z + , model (18) can be transformed into an extended model with p + 2 equations by introducing p new variables, similar to what was done in (36) for p = 2 . Then, by linearising the extended model at P 2 , we can obtain the characteristic equation p α , p λ = 0 , where
p α , p λ = ( λ + α ) p λ λ + 1 + n * 2 α p m λ + 1 n * 2 ,
which coincides with p α λ in (22) for p = 1 , and with p α , 2 in (39) for p = 2 . Also, by letting
q α , p λ = 1 α p p α λ
one has lim α + q α , p λ = p λ , and the characteristic equation (7) for the non-spatial model without delay is recovered.
For increasing p, there are new extra conditions for P 2 to be stable. Thus, by keeping the mean delay τ = p / α and increasing p, the effect of the delay in plant growth becomes more concentrated around the mean delay, and the stability of P 2 may only be compromised. Under these conditions, since α is not a root of (50), p α , p λ = 0 is equivalent to
λ λ + 1 + n * 2 m λ + 1 n * 2 1 + λ τ p p = 0 ,
and in the limit, as p tends to infinity, one obtains the characteristic equation of the model with discrete delay (3) [41],
λ λ + 1 + n * 2 m λ + 1 n * 2 e λ τ = 0 .

3. Conclusions and Discussion

In this work, we have analysed the stability of the spatially uniform vegetated state in a Klausmeier–Gray–Scott model of dryland vegetation with a distributed delay in the effect of water availability on plant growth. Our focus was on the conditions for the onset of Turing instability, comparing them with the corresponding original model without delay.
For the models without delay, we have given clear conditions for symmetry-breaking pattern formation by diffusion-driven instability to occur (Lemma 3) for the whole range of model parameters. In previous works, where other aspects of the models were the focus of the analysis, similar conditions were not as clearly defined and were usually restricted to a certain range of the mortality parameter (e.g., Theorem 2.5 in [36] for m < 2 , and Theorem 2.2 in [37] for m > 2 ).
Most previous works on Klausmeier-type models restricted the analysis to the case m < 2 , assumed to represent normal values for undisturbed natural conditions, but higher mortality values could be feasible and expected under harsh environmental conditions. Even under the assumed normal conditions, the upper limit seems to be somewhat arbitrary. For instance, the data presented in [19], referred to in [32] for the range [ 0.05 , 2.0 ] for m assumed there, would give an upper limit of 2.5 for that parameter, with values of m greater than 2 considered extreme but not unrealistic. Although m is usually referred to as mortality, it represents general biomass loss, which could also be the result of external disturbances, such as grazing. Grazing can significantly increase biomass loss (see, e.g., [44,45]) and, in general models where grazing is not specifically modelled, the effect of non-specific or local grazing is an increase in the biomass loss or mortality parameter [46]. General plant biomass loss is also directly related to respiration, and it is explicitly included as such in some vegetation models [47]. Plant respiration increases, through a power relation, with increasing temperature, and hence biomass loss due to respiration is expected to be much higher under future climate change scenarios than historically assumed [48,49].
We have introduced in this work a distributed delay in the product w n , which, in the original formulation of the Klausmeier model, corresponds to the product of the functional response of plants to water, G ( W ) , and the increase in infiltration due to the plants, F ( N ) . For simplicity, linear responses are usually assumed: G ( W ) = W and F ( N ) = N . Thus, the product w n in the scaled model represents water uptake by unit biomass, and the presence of a delay in this product is the consequence of the non-instantaneous effect of water uptake on plant growth. A Klausmeier–Gray–Scott model with a discrete delay in the same product w n was considered in [41], based on the necessary lag for precipitation to infiltrate into the soil and be available as soil water for plant growth. There are more complex dryland vegetation models where surface and soil water are differentiated (e.g., [19,20,50,51]), so a discrete delay for water infiltration would seem appropriate. However, in Klausmeier-type models, there is no such distinction, and the loss of water, besides evaporation, is directly the result of water taken up by plants. Hence, a discrete delay representing infiltration lag does not seem realistic in these models.
The effect on growth of water taken up by plants is certainly not instantaneous, nor does it disappear immediately. Plant growth at a certain time depends on water availability at that moment and also at previous ones, with some decaying effect persisting over time. Hence, a distributed delay, most likely in the form of an exponentially decaying kernel, as considered in this work, seems much more realistic than both a discrete delay and no delay at all. Distributed delays with Gamma kernels, as used in this work, are versatile enough to represent different possibilities of plant growth dependence on previous water uptake. They range from exponential decays to more localized delays and, in particular, weak kernels and strong kernels with p = 2 and a small mean delay represent exponential-like decays well, as would be expected in a base scenario where the effect of water taken up by plants on their rate of growth is constant.
We have shown that the presence of a distributed delay with a weak, exponential kernel does not affect vegetation stability at low mortality values (Theorem 1), in contrast with the effect of a discrete delay [41] or a distributed delay with a strong kernel (Theorem 4). Additionally, we have demonstrated that the only impact on diffusion-driven pattern formation is a shrinking of the Turing space, which occurs by reducing the parameter space where the spatially uniform vegetated state is stable in the non-spatial system (Theorems 3 and 6).
Generalised Klausmeier models with inertial effects, both on slopes and flat terrains, have also been proposed in the form of hyperbolic models [52,53,54]. These models consider non-Fickian diffusion to account for the finite speeds of propagation of disturbances. In this framework, similarly to some of our results, the occurrence of Turing instabilities is not affected by the hyperbolicity of the model, and there are no changes in bifurcation thresholds or in the wavenumbers of the stationary patterns. However, there may be changes in the transient dynamics, a question not addressed in our work.
The main contributions of this work are the introduction of an ecologically realistic form of delay in the classical Klausmeier–Gray–Scott model and the results showing how the effects on the stability of the homogeneous vegetated equilibrium and the onset of Turing instability depend on the type of delay, with it being absent in certain cases. In this paper, we focused on the conditions for Turing instability to occur. Once a patterned vegetation arises, it may evolve through different non-uniform steady-state solutions, helping to maintain vegetation beyond the limit point where a uniform vegetated state would experience a critical transition to the bare soil, desert state [12,13,35,36]. Future work will address the effect of different types of realistic delays on the stability and evolution of patterned vegetation in Klausmeier-type models.

Author Contributions

Conceptualisation, M.Á.C. and F.R.; formal analysis, I.M., F.Z.L., M.Á.C. and F.R.; writing—original draft preparation, F.R.; writing—review and editing, I.M., F.Z.L., M.Á.C. and F.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by “MCIN/AEI/10.13039/501100011033 (Ministerio de Ciencia e Innovación/Agencia Estatal de Investigación) grant number PID2021-125517OB-I00” and “Conselleria de Innovación, Universidades, Ciencia y Sociedad Digital, Generalitat Valenciana grant number CIPROM/2021/001”. I. Medjhadi and F.Z. Lachachi were supported by national grants from the Algerian Ministry of Higher Education and Scientific Research.

Data Availability Statement

Data are contained within the article.

Acknowledgments

We thank the anonymous reviewers for their comments and suggestions, which helped to improve the quality of this paper and the presentation of the results.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Rietkerk, M.; van de Koppel, J. Alternate stable states and threshold effects in semi-arid grazing. Oikos 1997, 79, 69–76. [Google Scholar] [CrossRef]
  2. Aguiar, M.R.; Sala, O.E. Patch structure, dynamics and implications for the functioning of arid ecosystems. Trends Ecol. Evol. 1999, 14, 273–277. [Google Scholar] [CrossRef]
  3. Klausmeier, C.A. Regular and irregular patterns in semiarid vegetation. Science 1999, 284, 1826–1828. [Google Scholar] [CrossRef]
  4. Deblauwe, V.; Barbier, N.; Couteron, P.; Lejeune, O.; Bogaert, J. The global biogeography of semi-arid periodic vegetation patterns. Glob. Ecol. Biogeogr. 2008, 17, 715–723. [Google Scholar] [CrossRef]
  5. Bastiaansen, R.; Jaïbia, O.; Deblauwe, V.; Eppinga, M.B.; Siteur, K.; Siero, E.; Mermoz, S.; Bouvet, A.; Doelman, A.; Rietkerk, M. Multistability of model and real dryland ecosystems through spatial self-organization. Proc. Natl. Acad. Sci. USA 2018, 115, 11256–11261. [Google Scholar] [CrossRef]
  6. Rietkerk, M.; Dekker, S.C.; de Ruiter, P.C.; van de Koppel, J. Self-organized patchiness and catastrophic shifts in ecosystems. Science 2004, 305, 1926–1929. [Google Scholar] [CrossRef]
  7. Kéfi, S.; Rietkerk, M.; Alados, C.; Pueyo, Y.; Papanastasis, V.P.; ElAich, A.; de Ruiter, P.C. Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems. Nature 2007, 449, 213–217. [Google Scholar] [CrossRef] [PubMed]
  8. Guttal, V.; Jayaprakash, C. Spatial variance and spatial skewness: Leading indicators of regime shifts in spatial ecological systems. Theor. Ecol. 2009, 2, 3–12. [Google Scholar] [CrossRef]
  9. Dakos, V.; van Nes, E.H.; Donangelo, R.; Fort, H.; Scheffer, M. Spatial correlation as leading indicator of catastrophic shifts. Theor. Ecol. 2010, 3, 163–174. [Google Scholar] [CrossRef]
  10. Dakos, V.; Kéfi, S.; Rietkerk, M.; van Nes, E.H.; Scheffer, M. Slowing down in spatially patterned ecosystems at the brink of collapse. Am. Nat. 2011, 177, E153–E166. [Google Scholar] [CrossRef]
  11. Kéfi, S.; Guttal, V.; Brock, W.A.; Carpenter, S.R.; Ellison, A.M.; Livina, V.N.; Seekell, D.A.; Scheffer, M.; van Nes, E.H.; Dakos, V. Early warning signals of ecological transitions: Methods for spatial patterns. PLoS ONE 2014, 9, e92097. [Google Scholar] [CrossRef] [PubMed]
  12. Rietkerk, M.; Bastiaansen, R.; Banerjee, S.; van de Koppel, J.; Baudena, M.; Doelman, A. Evasion of tipping in complex systems through spatial pattern formation. Science 2021, 374, eabj0359. [Google Scholar] [CrossRef] [PubMed]
  13. Kéfi, S.; Génin, A.; García-Mayor, A.; Guirado, E.; Cabral, J.S.; Berdugo, M.; Guerber, J.; Solé, R.; Maestre, F.T. Self-organization as a mechanism of resilience in dryland ecosystems. Proc. Natl. Acad. Sci. USA 2024, 121, e2305153121. [Google Scholar] [CrossRef] [PubMed]
  14. Kéfi, S.; Rietkerk, M.; van Baalen, M.; Loreau, M. Local facilitation, bistability and transitions in arid ecosystems. Theor. Popul. Biol. 2007, 71, 367–379. [Google Scholar] [CrossRef]
  15. Mayor, Á.G.; Kéfi, S.; Bautista, S.; Rodríguez, F.; Cartení, F.; Rietkerk, M. Feedbacks between vegetation pattern and resource loss dramatically decrease ecosystem resilience and restoration potential in a simple dryland model. Landsc. Ecol. 2013, 28, 931–942. [Google Scholar] [CrossRef]
  16. Mayor, A.G.; Bautista, S.; Rodríguez, F.; Kéfi, S. Connectivity-mediated ecohydrological feedbacks and regime shifts in drylands. Ecosystems 2019, 22, 1497–1511. [Google Scholar] [CrossRef]
  17. Couteron, P.; Lejeune, O. Periodic spotted patterns in semi-arid vegetation explained by a propagation- inhibition model. J. Ecol. 2001, 89, 616–628. [Google Scholar] [CrossRef]
  18. von Hardenberg, J.; Meron, E.; Shachak, M.; Zarmi, Y. Diversity of vegetation patterns and desertification. Phys. Rev. Lett. 2001, 87, 198101. [Google Scholar] [CrossRef]
  19. Rietkerk, M.; Boerlijst, M.C.; van Langevelde, F.; HilleRisLambers, R.; van de Koppel, J.; Kumar, L.; Prins, H.H.T.; de Roos, A.M. Self-organization of vegetation in arid ecosystems. Am. Nat. 2002, 160, 524–530. [Google Scholar] [CrossRef] [PubMed]
  20. Kéfi, S.; Eppinga, M.B.; de Ruiter, P.C.; Rietkerk, M. Bistability and regular spatial patterns in arid ecosystems. Theor. Ecol. 2010, 3, 257–269. [Google Scholar] [CrossRef]
  21. von Hardenberg, J.; Kletter, A.Y.; Yizhaq, H.; Nathan, J.; Meron, E. Periodic versus scale-free patterns in dryland vegetation. Proc. R. Soc. B Biol. Sci. 2010, 277, 1771–1776. [Google Scholar] [CrossRef]
  22. Turing, A.M. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B 1952, 237, 37–72. [Google Scholar]
  23. Prigogine, I.; Lefever, R. Symmetry breaking instabilities in dissipative systems. II. J. Chem. Phys. 1968, 48, 1695–1700. [Google Scholar] [CrossRef]
  24. Murray, J.D. Mathematical Biology. I. An Introduction, 3rd ed.; Springer: New York, NY, USA, 2002. [Google Scholar]
  25. Murray, J.D. Mathematical Biology. II. Spatial Models and Biomedical Applications, 3rd ed.; Springer: New York, NY, USA, 2003. [Google Scholar]
  26. Rietkerk, M.; van de Koppel, J. Regular pattern formation in real ecosystems. Trends Ecol. Evol. 2008, 23, 169–175. [Google Scholar] [CrossRef] [PubMed]
  27. Borgogno, F.; D’Odorico, P.; Laio, F.; Ridolfi, L. Mathematical models of vegetation pattern formation in ecohydrology. Rev. Geophys. 2009, 47, RG1005. [Google Scholar] [CrossRef]
  28. Sherratt, J.A. An analysis of vegetation stripe formation in semi-arid landscapes. J. Math. Biol. 2005, 51, 183–197. [Google Scholar] [CrossRef] [PubMed]
  29. Sherratt, J.A. Pattern solutions of the Klausmeier model for banded vegetation in semi-arid environments I. Nonlinearity 2010, 23, 2657–2675. [Google Scholar] [CrossRef]
  30. Sherratt, J.A. Pattern solutions of the Klausmeier model for banded vegetation in semi-arid environments II: Patterns with the largest possible propagation speeds. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 2011, 467, 3272–3294. [Google Scholar] [CrossRef]
  31. van der Stelt, S.; Doelman, A.; Hek, G.; Rademacher, J.D.M. Rise and fall of periodic patterns for a generalized Klausmeier–Gray–Scott Model. J. Nonlinear Sci. 2013, 23, 39–95. [Google Scholar] [CrossRef]
  32. Sherratt, J.A. Pattern solutions of the Klausmeiermodel for banded vegetation in semi-arid environments III: The transition between homoclinic solutions. Physica D 2013, 242, 30–41. [Google Scholar] [CrossRef]
  33. Sherratt, J.A. Pattern solutions of the Klausmeier model for banded vegetation in semiarid environments IV: Slowly moving patterns and their stability. SIAM J. Appl. Math. 2013, 73, 330–350. [Google Scholar] [CrossRef]
  34. Sherratt, J.A. Pattern solutions of the Klausmeier model for banded vegetation in semiarid environments V: The transition from patterns to desert. SIAM J. Appl. Math. 2013, 73, 1347–1367. [Google Scholar] [CrossRef]
  35. Siteur, K.; Siero, E.; Eppinga, M.B.; Rademacher, J.D.M.; Doelman, A.; Rietkerk, M. Beyond Turing: The response of patterned ecosystems to environmental change. Ecol. Complex. 2014, 20, 81–96. [Google Scholar] [CrossRef]
  36. Wang, X.; Shi, J.; Zhang, G. Bifurcation and pattern formation in diffusive Klausmeier-Gray-Scott model of water-plant interaction. J. Math. Anal. Appl. 2021, 497, 124860. [Google Scholar] [CrossRef]
  37. Sun, G.-Q.; Zhang, H.-T.; Song, Y.-L.; Li, L.; Jin, Z. Dynamic analysis of a plant-water model with spatial diffusion. J. Differ. Equ. 2022, 329, 395–430. [Google Scholar] [CrossRef]
  38. Gray, P.; Scott, S.K. Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms of multistability. Chem. Eng. Sci. 1983, 38, 29–43. [Google Scholar] [CrossRef]
  39. Gray, P.; Scott, S.K. Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Oscillations and instabilities in the system A +2B→3B; B→C. Chem. Eng. Sci. 1984, 39, 1087–1097. [Google Scholar] [CrossRef]
  40. Gray, P.; Scott, S.K. Sustained oscillations and other exotic patterns of behavior in isothermal reactions. J. Phys. Chem. 1985, 89, 22–32. [Google Scholar] [CrossRef]
  41. Li, J.; Sun, G.-Q.; Guo, Z.-G. Bifurcation analysis of an extended Klausmeier–Gray–Scott model with infiltration delay. Stud. Appl. Math. 2022, 148, 1519–1542. [Google Scholar] [CrossRef]
  42. Cushing, J.M. Integrodifferential Equations and Delay Models in Population Dynamics; Springer: Heidelberg, Germany, 1977. [Google Scholar]
  43. Smith, H. An Introduction to Delay Differential Equations with Applications to the Life Sciences; Springer: New York, NY, USA, 2011. [Google Scholar]
  44. Verwijmeren, M.; Smit, C.; Bautista, S.; Wassen, M.J.; Rietkerk, M. Combined grazing and drought stress alter the outcome of nurse: Beneficiary interactions in a semi-arid ecosystem. Ecosystems 2019, 22, 1295–1307. [Google Scholar] [CrossRef]
  45. Gong, X.; Wang, Y.; Zhan, T.; Wang, C.; Li, C.; Liu, Y. Advances in meta-analysis of the effects of grazing on grassland ecosystems in China. Agriculture 2023, 13, 1084. [Google Scholar] [CrossRef]
  46. Siero, E. Nonlocal grazing in patterned ecosystems. J. Theor. Biol. 2018, 436, 64–71. [Google Scholar] [CrossRef] [PubMed]
  47. Kéfi, S.; Rietkerk, M.; Katul, G.G. Vegetation pattern shift as a result of rising atmospheric CO2 in arid ecosystems. Theor. Popul. Biol. 2008, 74, 332–344. [Google Scholar] [CrossRef] [PubMed]
  48. Chen, Z.; Wu, Y.P.; Feng, G.L.; Qian, Z.H.; Sun, G.Q. Effects of global warming on pattern dynamics of vegetation: Wuwei in China as a case. Appl. Math. Comput. 2021, 390, 125666. [Google Scholar] [CrossRef]
  49. Sun, G.-Q.; Li, L.; Li, J.; Liu, C.; Wu, Y.-P.; Gao, S.; Wang, W.; Feng, G.-L. Impacts of climate change on vegetation pattern: Mathematical modeling and data analysis. Phys. Life Rev. 2022, 43, 239–270. [Google Scholar] [CrossRef] [PubMed]
  50. HilleRisLambers, R.; Rietkerk, M.; van den Bosch, F.; Prins, H.H.T.; de Kroon, H. Vegetation pattern formation in semi-arid grazing systems. Ecology 2001, 82, 50–61. [Google Scholar] [CrossRef]
  51. Gilad, E.; Shachak, M.; Meron, E. Dynamics and spatial organization of plant communities in water limited systems. Theor. Popul. Biol. 2007, 72, 214–230. [Google Scholar] [CrossRef] [PubMed]
  52. Consolo, G.; Currò, C.; Valenti, G. Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments. Appl. Math. Model. 2017, 43, 372–392. [Google Scholar] [CrossRef]
  53. Consolo, G.; Currò, C.; Valenti, G. Supercritical and subcritical Turing pattern formation in a hyperbolic vegetation model for flat arid environments. Physica D 2023, 398, 141–163. [Google Scholar] [CrossRef]
  54. Grifó, G.; Consolo, G.; Currò, C.; Valenti, G. Rhombic and hexagonal pattern formation in 2D hyperbolic reaction–transport systems in the context of dryland ecology. Physica D 2023, 449, 133745. [Google Scholar] [CrossRef]
Figure 1. Regions of stability (S) and instability (U) for the vegetated equilibrium P 2 in the non-spatial model (4).
Figure 1. Regions of stability (S) and instability (U) for the vegetated equilibrium P 2 in the non-spatial model (4).
Symmetry 16 00609 g001
Figure 2. Phase portrait (a, left) and bifurcation diagram (b, right) of system (3) with m = 4 , where a c 9.24 . (a) Evolution of the system from w ( 0 ) = 3.5 and n ( 0 ) = 1.2 with a = 9.3 (blue) and a = 9.1 (red). (b) Stable (blue) and unstable (grey) positive vegetation equilibria and bifurcation points (H: Hopf bifurcation at a = a c ; F: Fold bifurcation at a = 2 m ).
Figure 2. Phase portrait (a, left) and bifurcation diagram (b, right) of system (3) with m = 4 , where a c 9.24 . (a) Evolution of the system from w ( 0 ) = 3.5 and n ( 0 ) = 1.2 with a = 9.3 (blue) and a = 9.1 (red). (b) Stable (blue) and unstable (grey) positive vegetation equilibria and bifurcation points (H: Hopf bifurcation at a = a c ; F: Fold bifurcation at a = 2 m ).
Symmetry 16 00609 g002
Figure 3. Gamma kernels with different parameter values. (a) Weak ( p = 1 ) and strong ( p = 2 ) kernels with different variances p / α 2 . (b) Gamma kernels with equal means and decreasing variances for increasing p values.
Figure 3. Gamma kernels with different parameter values. (a) Weak ( p = 1 ) and strong ( p = 2 ) kernels with different variances p / α 2 . (b) Gamma kernels with equal means and decreasing variances for increasing p values.
Symmetry 16 00609 g003
Figure 4. Regions of stability (S) and instability (U) for the vegetated equilibrium P 2 in the non-spatial model with distributed delay (18), with p = 1 . In the region C ( α ) , the stability is conditional on the value of α , as given in Theorem 1.
Figure 4. Regions of stability (S) and instability (U) for the vegetated equilibrium P 2 in the non-spatial model with distributed delay (18), with p = 1 . In the region C ( α ) , the stability is conditional on the value of α , as given in Theorem 1.
Symmetry 16 00609 g004
Figure 5. Evolution of the system with distributed delay (18) with p = 1 , from initial values close to the vegetated equilibrium P 2 , with parameters m = 6 and a = 18 , for decreasing values of α . (a) α = 1.1 α c . (b) α = 0.9 α c . (c) α = 0.6 α c . (d) α = 0.4 α c .
Figure 5. Evolution of the system with distributed delay (18) with p = 1 , from initial values close to the vegetated equilibrium P 2 , with parameters m = 6 and a = 18 , for decreasing values of α . (a) α = 1.1 α c . (b) α = 0.9 α c . (c) α = 0.6 α c . (d) α = 0.4 α c .
Symmetry 16 00609 g005
Figure 6. Evolution of the system with distributed delay (17) with p = 1 , from initial values close to the vegetated equilibrium P 2 , with parameters m = 2.95 , a = 6.5 , and d 1 = d 2 = 0.1 for two values of α . (a) α = 2.5 . (b) α = 8 .
Figure 6. Evolution of the system with distributed delay (17) with p = 1 , from initial values close to the vegetated equilibrium P 2 , with parameters m = 2.95 , a = 6.5 , and d 1 = d 2 = 0.1 for two values of α . (a) α = 2.5 . (b) α = 8 .
Symmetry 16 00609 g006
Figure 7. Spatially non-homogeneous steady-state solutions resulting from Turing bifurcations. (a, left panels) m = 1 , a = 2.2 , d 1 = 50 , d 2 = 0.1 . (b, right panels) m = 2.5 , a = 5.25 , d 1 = 50 , d 2 = 0.1 .
Figure 7. Spatially non-homogeneous steady-state solutions resulting from Turing bifurcations. (a, left panels) m = 1 , a = 2.2 , d 1 = 50 , d 2 = 0.1 . (b, right panels) m = 2.5 , a = 5.25 , d 1 = 50 , d 2 = 0.1 .
Symmetry 16 00609 g007
Figure 8. Evolution of the system with distributed delay (18) with p = 2 , from initial values close to the vegetated equilibrium P 2 , with parameters in the region m < 2 (top) and m > 2 (bottom), for values of α greater (left) and lower (right) that the critical value α c , 2 . (a) m = 1.5 , a = 4 , α = 1.1 α c , 2 . (b) m = 1.5 , a = 4 , α = 0.9 α c , 2 . (c) m = 3 , a = 8 , α = 1.1 α c , 2 . (d) m = 3 , a = 8 , α = 0.9 α c , 2 .
Figure 8. Evolution of the system with distributed delay (18) with p = 2 , from initial values close to the vegetated equilibrium P 2 , with parameters in the region m < 2 (top) and m > 2 (bottom), for values of α greater (left) and lower (right) that the critical value α c , 2 . (a) m = 1.5 , a = 4 , α = 1.1 α c , 2 . (b) m = 1.5 , a = 4 , α = 0.9 α c , 2 . (c) m = 3 , a = 8 , α = 1.1 α c , 2 . (d) m = 3 , a = 8 , α = 0.9 α c , 2 .
Symmetry 16 00609 g008
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Medjahdi, I.; Lachachi, F.Z.; Castro, M.Á.; Rodríguez, F. Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation. Symmetry 2024, 16, 609. https://doi.org/10.3390/sym16050609

AMA Style

Medjahdi I, Lachachi FZ, Castro MÁ, Rodríguez F. Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation. Symmetry. 2024; 16(5):609. https://doi.org/10.3390/sym16050609

Chicago/Turabian Style

Medjahdi, Ikram, Fatima Zohra Lachachi, María Ángeles Castro, and Francisco Rodríguez. 2024. "Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation" Symmetry 16, no. 5: 609. https://doi.org/10.3390/sym16050609

APA Style

Medjahdi, I., Lachachi, F. Z., Castro, M. Á., & Rodríguez, F. (2024). Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation. Symmetry, 16(5), 609. https://doi.org/10.3390/sym16050609

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop