Delay-induced instabilities of stationary solutions in a single species nonlocal hyperbolic-parabolic population model

The present work investigates the effects of maturation and dispersal delays on dynamics of single species populations. Both delays have been incorporated in a single species nonlocal hyperbolic-parabolic population model, which admits traveling and stationary wave solutions. We reduce the model into various forms and obtain the corresponding analytical solutions. Analysis of the reduced models indicates that the dispersal delay can result in loss of monotonicity, where the solutions oscillate as they converge to a positive equilibrium. The stability analysis of the general model reveals that the maturation time delay admits a Hopf bifurcation threshold, which is expressed as a function of the dispersal delay. The numerical simulations of the general model suggest that the global stability of the stationary wave solutions is lost when the dispersal delay is increased from zero. In conclusion, population models with maturation and dispersal delays can give new insights into the complex dynamics of single species.


1.
Introduction. Mathematical modeling of population dynamics dates back to the late eighteenth century where exponential growth (also known as Malthusian) models were experimented for different single species populations. Later such models were replaced with logistic growth models which include a decay mechanism due to increased intra-specific competition and depletion of resources. An extended version of the logistic growth model [30] is given by where all constants a, b, r and k are positive, r is the growth rate, k is the carrying capacity, a is the growth regulator and b is the depletion factor. The original logistic growth model is obtained when a = 1 and b = 1. The modified versions of the logistic growth model include the following two: where L > 0 represents the harvested rate, L < 0 is the restocked rate, and αe −βu (t) relates to the migration of individuals to each colony. The other modified version of the logistic growth model is the Gompertz model [10,14,47] du(t) dt = u(t)(α − β ln u(t)) where α and β are per capita birth and death rates, respectively. The abovementioned models are categorized as density dependent birth-death models with the general form where u(t) is the population size at time t, and b(u) and d(u) are the birth and death functions, respectively. The general model (4) is missing several key factors such as spatial movement of individuals, age structure, time delay, and nonlocality.
In particular, the changes in the growth rate du(t)/dt are not necessarily a function of population size at the present time, but at previous times. To this end, several spatially-homogeneous delay models were constructed to study the effects of delay on population dynamics. The books by K. Gopalsamy [15], Y. Kuang [22] and J. Hale and S. Verduyn Lunel [16] are excellent sources of spatially homogeneous delay modeling and analysis, where the last book provides a solid theoretical background for studies of delay differential equations. Inclusion of a diffusion term to discrete delay models was the first approach to study delay effects in processes of population dispersal. Then the birth-death model can be extended to ∂u(x, t) ∂t = D ∂ 2 u(x, t) ∂x 2 + f (u(x, t), u(x, t − τ )), where u ∈ R n and f is a vector field that includes both birth and death functions. The delayed Reaction-Diffusion (RD) equation (5) has been studied in the book by J. Wu [52], which covers existence and compactness, linear stability, invariant manifolds, Hopf bifurcation, and traveling wave solutions. An example of delayed RD equation (5) is given by which represents an extended version of the logistic growth model (also known as the Hutchinson model) with delay and diffusion. Such a model has been well studied [28,37,55], but also has been criticized for the lack of nonlocality [9] and age structure [46,45]. A model that takes into account the maturation time delay, age structure and nonlocality of individuals was first introduced in the seminal paper by So et al. [39]. Their model was later extended to two-dimensional spatial domains of rectangle [24], disk [8], unbounded strip [51], and the entire R 2 [2]. Despite the major breakthrough in modeling single species population, the model developed by So et al. [39] overlooked the travel time. It has been well documented in ecology literature [1,13] that dispersal of many species occurs with time delays. The fact that it takes time for individuals to move from one location to another has not been considered in the model by So et al. This issue was first addressed in the work by Ou and Wu [29]. However, details of their model derivation were not revealed. The present work aims to fill these gaps by providing the details of model derivation and investigating the effects of dispersal and maturation delays on dynamics of mature populations. The dispersal delay reflects the travel time for population movements between the patches, whereas the maturation time delay accounts for the time from birth until the individuals are sexually mature. The concept of dispersal delays has also been incorporated into Lotka-Volterra type models [11,20], which reflects time delays for interpatch migration. The work by Neubert et al. [27] showed that dispersal delay can stabilize the positive equilibrium of a prey-predator model. This is an interesting result, since the delay terms often have destabilizing effects [15,22]. Other models with dispersal delay include the impulsive diffusive single species model proposed by Wan et al. [49,50] and the model by Takeuchi et al. [44]. They showed that in an environment with different food-rich and food-poor patches the permanence of single species is ensured if the dispersion between the patches is sufficiently small. It has been shown that the Fisher RD equation [2] can be extended to a time-delayed hyperbolic-parabolic reaction diffusion (HPRD) model, which provides more accurate results related to the neolithic transition [12,32]. Time-delayed diffusion has also been applied in other studies such as optics in turbid media [18] and x-ray bursters [34].
The model investigated in the present work includes two different types of delay: (i) A maturation time delay which conveys the effects of birth in the past to the growth rate at the present, and (ii) a dispersal delay which takes into account the effects of inter-patch movement on population dynamics.
The rest of this paper is organized as follows. Section 2 provides the details of model derivation and model justification. Section 3 is a collection of reduced models where the analytical and numerical solutions of some reduced models are presented. It is shown that the dispersal delay will give rise to oscillatory convergence to constant equilibrium solutions. Section 4 covers the stability analysis of the general model, where the expression for a Hopf bifurcation threshold is derived. Also the effects of dispersal delay on stability of stationary solutions are numerically investigated. Section 5 provides a discussion of the main implications of the results. 2. The model.

Model derivation.
We are concerned with particles (e.g., insects, reptiles, viruses and bacteria) that have two age classes: the immature (i.e., unable to reproduce) and the mature particles. Let u(t, a, x) denote the population density of particles at time t, age a and spatial location x ∈ R. Let τ ≥ 0 be the average time required for maturation of each individual. Then is the total adult population size at time t and spatial position x. Below we shall use the concept of isotropic random walks to derive a deterministic model for studying the dynamics of mature individuals. Using this concept, any particle has the same probability to jump to the left or to the right, and also all jumps have the same length of δ. We also assume that the particles spend the same time interval T to perform any of the jumps. Then we have where d(a) is the death rate for the individuals at age a during the jump. Since T and δ are very small, we can use Taylor's formula to expand the above equation up to the second-order terms and obtain where r = T /2 is the dispersal delay term and D(a) = (1 − d(a)T )δ 2 /2T is the agedependent diffusion rate. The dispersal delay r can be thought of the inter-patch travel time in a patchy environment or simply the time required for the particle to move a length of δ to the left or right.

Remark 1.
Inclusion of the second-order Taylor's terms is the main difference in our modeling approach compared to that of So et al. [39]. The second-order terms allow us to take into account the dispersal delay r ≥ 0. When r = 0, our approach coincides with that of So et al. [39].
We aim to derive a more realistic model that reflects the dynamics of the mature population with respect to the maturation time delay τ and the dispersal delay r. For simplicity we assume that the diffusion and the death rates of the mature population are age independent. In particular, we assume that D(a) = D m and d(a) = d m for all a ∈ [τ, ∞]. Then by taking the integral from both sides of equation (9) with respect to a from τ to ∞, using equation (7), adding a constant source k ≥ 0, and noting that u(t, ∞, x) = 0, we have To obtain the model of the mature population, we need to calculate u(t, τ, x) and (2∂/∂t + ∂/∂a)u(t, τ, x). To this end, fix s ≥ 0 and let for s ≤ t ≤ s + τ. Using the chain rule we have Then from equation (9) we get with Equation (13) is a linear reaction-diffusion equation. Applying the method of separation of variables to (13) and (14), we obtain the general solution where k(s, ξ) is determined as follows: In other words, the birth function b(w(s, x)) is the Fourier transform of k(s, ξ). Therefore k(s, ξ) is the inverse Fourier transform of b(w(s, x)) and hence For a ∈ [0, τ ], denote the diffusion and death rates of the immature population with D I (a) and d I (a), respectively. By setting s = t − τ and we have Let and Then ε ∈ [0, 1], represents the portion of the immature population that survives and becomes mature. The end values ε = 0 and ε = 1 correspond to extreme situations at which none or all of the immature individuals survive, respectively. Substituting (18) and (19) into (17) we get By letting s = t − τ and substituting (20) into equation (11) we get Substituting (21) into equation (10) leads to Multiplying both sides of equation (21) by r and taking the partial derivative with respect to t, we get (23) Assuming that r is small, the term with r 2 can be dropped and substituting (23) into equation (22) we obtain the following model of the mature population which is a second-order Hyperbolic-Parabolic Reaction-Diffusion (HPRD) equation.
We will justify this model in the following subsection.

2.2.
Model justification. The second-order delayed HPRD model can be rewritten as ∂ ∂t is the survival-birth convolution term [9]. Table 1 provides a description of the parameters and functions used in model (26). Equation (26) shows that the growth rate ∂w(t, x)/∂t of the mature population is a function of population dispersal, convoluted birth function, the rate of change in the growth rate and death of the mature population. By discretizing equation (26) with respect to x we can justify each component of the model. As shown in Figure 1, the single species population is spread out to adjacent habitats that are on a straight line. Each habitat is centered at a location x i and has a different quality of life. Individuals randomly move between the habitats and the movement is associated with the dispersal delay r > 0. The diffusion term is discretized as which shows how the random walk of the mature population between the habitats relates to the growth rate ∂w(t, x)/∂t. Specifically the diffusion is likely to reduce Figure 1. A schematic representation of the delayed HPRD model (26). Individuals can randomly move between infinitely many habitats on a straight line. Considering habitats centered at x i−1 , x i and x i+1 , the effects of birth function b(w), mortality rate d m , survival rate ε, dispersal delay r, migration between habitats f α * b(w), and migration from other regions k i are indicated in the diagram. Notes: b(w) and ε are dimensionless; see (18) and (19) for ε and f α (x), respectively.
the growth rate, when w(t, x i ) > max{w(t, x i−1 ), w(t, x i+1 )}. The diffusion rate D m > 0 is a measure of how fast the population dispersal between the habitats takes place. As individuals randomly move, birth and mortality occur within each habitat. Mortality of the mature population is density dependent and reduces the growth rate in each habitat with the standard term d m w(t, x i ). Despite mortality, the birth of individuals has a long term effect on the mature population. The key point is that individuals who were born in one habitat may randomly spread into different habitats and influence the growth rate of mature populations in those habitats. The survival-birth convolution term (f α * b(w)) (t−τ,xi) is a spatial weighted sum of all such influences on a habitat centered at x i . In particular, the offspring of individuals, who were born at time t − τ and location y have now reached to location x i and have become mature at the present time t. Therefore, they can influence the growth rate of the mature population in the habitat centered at x i .
is a decreasing function of y. As indicated in equation (19), the function f α (x i − y) is the Gaussian distribution showing that the further the distance x i − y, the less likelihood of reaching from habitat y to habitat x i . The coefficient ε of the survival-birth convolution f α * b(w) is the survival rate of the immature population. Since the dispersal of individuals has a delay, the flow of individuals reaching from y to x i does not occur instantaneously. Also, the term shows that the rate of the change in the spatial sum influences the growth rate of mature population in habitat x i . As shown in Figure 1, we have assumed that all habitats are on a straight line. Nevertheless mature individuals may migrate from other locations to each habitat. For simplicity we assume that the migration rate to each habitat is the same (i.e., k i = k for all i ∈ Z). In addition to the above-mentioned factors, the growth rate of mature population ∂w(t, x)/∂t is regulated by the rate of change in the growth rate ∂ 2 w(t, x)/∂t 2 . In particular, the term −r∂ 2 w(t, x)/∂t 2 inversely relates the growth rate with the concavity of the mature population due to dispersal delay r ≥ 0. When the mature population is growing and the corresponding curve is concave up, the term −r∂ 2 w(t, x)/∂t 2 is negative and reduces the growth rate. Whereas a concave down growth makes the term −r∂ 2 w(t, x)/∂t 2 positive and increases the growth rate. Hence, the term −r∂ 2 w(t, x)/∂t 2 captures the impact of dispersal delay, which makes the growth rate of the mature population bounce back and forth within a certain interval. In the next section we will see that the dispersal delay will cause the mature population w(x, t) to oscillate around the constant ratio k/d m .
3. Model reduction and analysis. The delayed HPRD model (26) can be reduced to several models, some of which have already been studied. In the following we provide a summary of the reduced models and the related analysis. In the most naive case we may assume that the population is spatially homogeneous (i.e., D I = D m = 0), there is no maturation or dispersal delays (i.e., τ = r = 0) and the immature population does not survive (i.e., ε = 0). Then model (26) is reduced to which has a globally stable equilibrium at w * = k/d m with k ≥ 0. In this case, the survival of the mature population totally depends on the migration rate k. When k > 0, the population w(t) monotonically converges to the value k/d m , whereas k = 0 results in population extinction (i.e., lim w(t) = 0 as t → ∞). By taking into account the dispersal delay r > 0, model (28) is elevated to which can result in oscillatory solutions. Namely, the general solution of model (29) is given by The eigenvalues λ 1 and λ 2 are complex when rd m > 1/4. Hence the population w(t) oscillates as it converges to the value k/d m when rd m > 1/4. Otherwise the convergence to k/d m is monotonic. In the presence of diffusion (i.e., when D m > 0), the mature population is spatially inhomogeneous and model (29) is elevated to As shown in the next section, equation (30) plays a major role in the analysis of the general model (26) and investigating the effects of dispersal delay. Assuming that a portion of the immature population survives (i.e., when ε > 0), model (28) is changed to d dt which includes the contribution of the immature population εb(w(t)) by becoming mature and reproducing. Note that for all birth functions b(w) we have b(0) = 0 and b (0) > 0. Hence, depending on the shape of the birth function b(w) (see Figure  1 of [6]), model (31) may have up to two positive equilibria w * 2 and w * 3 , and the trivial equilibrium w * 1 , provided that k = 0. With two positive equilibria w * 2 < w * 3 , model (31) exhibits founder-control (i.e., lim t→∞ w(t) = w * 3 when w(0) > w * 2 or lim t→∞ w(t) = 0 when w(0) < w * 2 ). The founder-control state is also known as Extinction-Survival (ES) outcome, where the initial population size determines the fate of the single species. With only one positive equilibrium w * 2 , the model exhibits When k > 0, model (31) cannot have the trivial equilibrium w * 1 = 0. However, the establishment at the positive equilibrium is still possible. By including the maturation time delay (τ > 0) to model (31), we get to another reduced form of the general model (26) given by d dt The dynamics of model (32) are the same as those of model (26) provided that εb (w * i ) + d m > 0 for i = 1, 2, 3. If the direction of the last inequality is reversed, then the stability of w * i becomes delay dependent, where w * i loses its stability when τ is greater than a critical valueτ (the proof is similar to the proof of Theorem 3 in [6]). When τ >τ , oscillatory and periodic solutions bifurcate from w * i > 0, which refers to a different quality of population establishment (i.e., the population size w(t) never reaches to zero but it oscillates over time). By including the diffusion term to the model we get where D m > 0 is the diffusion rate of the mature population. Also, by including the maturation time delay to model (33) we obtain the following delayed RD model Given that the diffusion has a stabilizing effect, the constant equilibria of models (33) and (34) preserve their stability and the model dynamics remain the same with respect to constant equilibria. The new dynamics of models (33) and (34) compared to those of the previous models (28) -(32) are the existence of traveling and stationary wave solutions. The traveling wave solutions of these models have been studied by several researchers (see [52] for a review). Also, the existence of stationary wave solutions of models (33) and (34) has been proven in [3,6] and their stability has been numerically examined. Formation of a traveling wavefront may resemble the spatial dispersal of invasive species, whereas a stationary wavefront or wave pulse relates to space-dependent extinction-survival (see Table 3 of [6]). Further analysis of model (34) includes the asymptotic behavior of the solutions [26], Hopf bifurcation [38], approximated wave solutions using a differential transform method [4], solutions of the model (34) with Dirichlet boundary conditions [40], effects of maturation time delay in the absence of diffusion [37] and the global stability of the solutions [37].
When the immature population is immobile (i.e. D I = 0), the value of α will be zero and the heat kernel f α (x) changes to Dirac delta function. This will eliminate the nonlocality (i.e., the convolution term) from the general HPRD model (26) and the effect of birth on population growth will be localized. In particular, we have which reduces model (26) to This model represents dynamics of the mature population without including the effects of those who are at location x, but were born at a location y = x. On the other hand, when the dispersal delay is negligible (i.e., when r = 0), the general HPRD model (26) is reduced to which assumes that the random walk occurs instantaneously. Model (36) was first introduced by So et al. [39], where they proved the existence of traveling wave solutions using the method of upper-lower-solution. Later Liang and Wu [23] numerically investigated the behavior of model (36) with an extra advection term. They found that the traveling wavefront exhibits humps when the monotonicity conditions are violated. The behavior of the traveling wavefront was further studied in [7] where it was shown that the monotonic wave solution can become oscillatory when all eigenvalues of the nontrivial equilibrium have nonzero imaginary parts.
The wave solution of model (36) was also approximated using a boundary layer method combined with the asymptotic expansion [2,5]. Moreover, the existence of stationary wave solutions and the convergence of the w(x, t) solution of model (36) to the stationary wave solutions were investigated in [3,6,39]. The numerical solutions of the traveling and stationary solutions strongly suggest that these solutions are locally and globally stable [3,6]. The existence, stability and behavior of the stationary and wave solution of the model (36) have several biological implications which have been discussed in [3]. For instance, it has been shown that the threshold property of single species with Allee effect does not hold when the diffusion is present (i.e., when D m > 0; see figures 3(d) and 4(d) of [3]). The rest of this section is focused on the reduced model (30). Before constructing the analytic τ -periodic solution of model (30) we prove the following theorem.   For model (30), the behavior of the solution depends on the corresponding characteristic equation (113) provided in Appendix A. When r > 0, the eigenvalues λ 1 and λ 2 become complex numbers for some ζ ∈ R with a magnitude large enough. However, in the absence of dispersal delay (i.e. when r = 0), the characteristic equation (113) leads to a single eigenvalue which is real and negative for all ζ ∈ R. We therefore have the following remark.  To visualize the oscillatory versus monotonic convergence, we used the software package Comsol Multiphysics 5.3 to solve the reduced model (30) for two cases of r > 0 and r = 0. When r = 0, Figure 2 (a) shows that the solution w(t, x) monotonically converges to k/d m = 2, as t → ∞. Whereas Figure 2 (b) shows that for r = 2 the convergence becomes non-monotonic and the solution w(t, x) oscillates as it converges to k/d m = 2. The video clips "monotonic.gif" and "oscillatory.gif" available in the supplementary document show the monotonic and oscillatory convergences, respectively. The specific parameter values used for the simulation are d m = 2, D m = 0.05 and k = 4. We assumed that the mature population initially has a density of w(0, x) = 2+2e −0.01(x+30) 2 −2e −0.01(x+10) 2 +e −0.05x 2 −2e −0.03(x−15) 2 + 2e −0.05(x−20) 2 − 1.9e −0.01(x−45) 2 with the velocity ∂w(0, x)/∂t = 0. The simulation results remained the same when different initial and parameter values were utilized.
A τ -periodic solution of model (30) is a solution w(t, x) such that w(t − τ, x) = w(t, x) for some τ > 0. Similar to [33], the following theorem shows that model (30) admits a τ -periodic solution.
Theorem 2. The reduced model (30) has a τ -periodic solution of the form w(t, exp(λ n x)[C n cos(β n t + γ n x) + D n sin(β n t + γ n x)], A n , B n , C n and D n are constants and λ n = βn 2Dmγn . Theorem 2 is proved by the method of separation of variables. See Appendix A for details. (118) in Appendix A, we get that γ n is an increasing function of r > 0, whereas β n is a decreasing of τ > 0. Hence the dispersal delay r > 0 and the maturation time delay τ > 0 both decrease the value of λ n = β n /2D m γ n which make the solution W (t, x) converge to k/d m at a slower rate. This implies that the dispersal delay prolongs the oscillatory behavior of solution W (t, x) of the reduced model (30).

Remark 4. From equation
The following theorem shows that the reduced model (30) may have inverted τ -periodic solutions.
The explicit form of W 2 (x, t) is the same as equations (37) and (38) with A 0 = 0, B 0 = 0 and β n = 2π(n − 1)/τ for n = 1, 2, . . . . The proof is omitted here since it is similar to the proof of Theorem 2.
When the immature population is immobile (i.e., when D I = 0), the value of α changes to zero and the general model (26) reduces to (39) Let w(t − τ, x) = v(t, x). Also we may extend the birth function b(v) to b(v, w). This will allow us to investigate birth functions such as the one used in the Hutchison equation (6). Assume that b(v, w) has the property that for some real-valued function G. Then the reduced model (39) can be rewritten We have the following theorem. (a n θ n (x) + b n γ n (x))e cnt + k/d m , where the constants a n , b n and c n are determined using the Fourier series corresponding to the functions p(x) and q(x). The eigenfunctions θ n (x) and γ n (x) are in the form of sin Substituting (43) into (41) and separating the terms according to variables x and t, the reduced model (41) is transformed to and where β is an arbitrary constant. Equation (44) is a second-order delay differential equation which can be solved using functions g(t) = e σ1t and g(t) = e σ2t for t ∈ [−τ, t]. It can be easily verified that g(t) = Ae crt is a particular solution of (44), where A is an arbitrary constant and c r is determined by solving Equation (45) may have three types of solutions which completes the proof. If |x| < L, then the boundary conditions can determine the specific form of the function h(x). Typically Dirichlet boundary condition leads to the case β < 0 whereas the zero-flux boundary condition gives rise to the case β > 0. When the domain is semi-finite (i.e., x ≥ 0), the right boundary condition is replaced with lim x→∞ w(x, t) = 0. This can correspond to the traveling wave solutions of the reduced model (41), which are solutions of the form The constant c is the wave speed and the solution φ(z) travels in either the left or the right direction without changing its shape. Substituting solution (48) into model (41) we get to the wave equation The next two theorems show that the reduced model (41) has a general solution based on the τ -periodic solution W 1 (t, x) of the reduced model (30).
The reduced model (41) has the general solution where the τ -periodic solution W 1 (t, x) is given by (37) with the only difference of γ n being replaced with γ n . The constant c is obtained by solving Proof. Substituting solution (51) into the reduced model (41) and using equation (30), we get to equation (52). This completes the proof.

Remark 5.
Using the boundedness and the asymptotic property of W 1 (t, x), the solution (51) may converge to k/d m as x → ∞ provided that C n = D n = 0 for all n = 1, 2, .... Theorem 6. The reduced model (41) has the general solution Proof. The constant c is obtained by solving The explicit form of W 2 (x, t) is the same as equations (37) and (38) with A 0 = 0, B 0 = 0 and β n = 2π(n − 1)/τ for n = 1, 2, . . . . Also γ n is replaced with γ n provided in (50). This proof is similar to the proof of the last theorem and therefore it is omitted.
We look for solutions of the reduced model (41) that are in the form of nonlinear superposition of two different wave solutions. In particular we have the following theorem.
Theorem 7. The reduced model (41) admits a wave solution in the form of a multi-scaled traveling wave product given by where z i = k i x + c i t, i = 1, 2 and k 1 and k 2 are positive constants. The function θ(z 2 ) can take the particular form of θ(z 2 ) = A exp(z 2 ).
Proof. Substituting the solution (55) into the reduced model (41) yields the following second order delay differential equation: and The function θ(z 2 ) is the general solution of equation (56) and it has the particular solution θ(z 2 ) = A exp(z 2 ), where A is a constant determined by solving The next section concentrates on the local solutions of the general model (26) and the numerical simulations of the wave solutions. Note that the existence of the traveling wave solutions of the general model (26) has already been proven in the work by Ou and Wu [29]. 4. Analysis of the general model.

4.1.
Linear stability analysis. In this section we first investigate the local behavior of the solution in a neighborhood of the constant equilibrium solution w j of the general model (26). Let the solution w(x, t) of model (26) be approximated by in a sufficiently small neighborhood of w j . Our goal is to find the approximate solution by calculating w(t, x). Substituting (60) into model (26), using the Taylor expansion and dropping the nonlinear terms we get to Assuming that w(x, t) = e λt+ikx the linearized model is rewritten which leads to the characteristic equation and k ∈ R is the wave number.

Remark 6. Assuming the initial history functions
and the general solution of the linearized equation (62) is given by where the A(k) is obtained with the Fourier transforms of g(x), h(x) and the eigenvalue λ depends on k, τ and r.
The following remark indicates the relation between the analysis of the reduced model (30) and the local solutions of the general model (26).

Remark 7.
In equation (30), by letting k = 0 and allowing d m to take positive or negative values, we get that for each λ, Theorems 1-4 are applicable to the linearized model (62).
The eigenvalues λ(k, τ, r) can determine the local behavior of the solutions based on the dispersal delay r and maturation delay τ. If λ(k, τ, r) is real and negative then for each x ∈ R, the approximate solution w(x, t) will converge monotonically to the constant equilibrium w j . When the eigenvalues are all pure imaginary numbers, the approximate solution w(x, t) nearby w j is periodic with respect time t for each x ∈ R. If there is an eigenvalue with positive real part then the solution will diverge from w j . A combination of the above-mentioned cases will result in various spatial patterns (see for example figures 4 and 5 of [6] for the case r = 0). We have the following two theorems, regarding the stability of w j . Theorem 8. Let w j be a constant equilibrium of the general model (26) and let τ = 0. Then w j is locally asymptotically stable if and only if Proof. When τ = 0, the characteristic equation (63) is reduced to Using the Ruth-Hurwitz criterion [21], the necessary and sufficient conditions for and min Given that e −αk 2 is a decaying function of k, the minimum is attained at k = 0 which leads to 1 − εb (w j )r > 0 (71) and d m − εb (w j ) > 0 (72) This completes the proof.
The next theorem shows that stability of a constant equilibrium is lost when the maturation time delay exceeds the Hope bifurcation threshold, which is expressed as a function of the dispersal delay. Theorem 9. Let w i be a constant equilibrium of model (26). If εb (w j ) + d m < 0, then w j loses its stability when τ exceeds the Hopf bifurcation thresholdτ r where Proof. We show that there exists a pair of pure imaginary eigenvalues λ = ±vi when τ exceedsτ .
Considering the Hopf bifurcation of the spatially homogeneous model, we let k = 0. Then substituting λ = vi into the characteristic equation (63) we get Equating the imaginary and real parts of equation (77) to zero, we get v − εb (w j ) sin(vτ ) − εrb (w j ) cos(wτ )v = 0, and respectively. Equations (78) and (79) are rewritten as v = εb (w j )c sin(vτ + θ), and respectively where c = 1 + (rv) 2 and θ = arctan(rv). By squaring both sides of (80) and (81) and adding them we get to By the assumption εb (w j ) + d m < 0, equation (82) is guaranteed to have a real positive root, w s as indicated in equation (74). Substituting w s into equation (81) and solving for τ , we getτ .
Remark 8. Theorem 9 is a generalized version of Theorem 3, part (iii) in [6]. When r = 0, the expression for the Hopf bifurcation threshold is reduced tô (83) Figure 3 shows the graph of the Hopf bifurcationτ (r) as an increasing function of the dispersal delay r. This may suggest that the dispersal delay may stabilize the constant equilibrium w j , which is in agreement with [27].   The remaining part of this subsection is focused on finding conditions for the local stability of w j using the Lambert W function. When r = 0, the characteristic equation (63) is reduced to Let λ r be a root of (63) for r > 0. Then λ r = λ 0 + λ 1 , where λ 0 is the root for the case r = 0, and λ 1 = λ r − λ 0 . Let and g(λ) = λ + D m k 2 + d m . Then the characteristic equation (63) is rewritten as Substituting λ r = λ 0 + λ 1 and noting that Equation (86) is rewritten as or equivalently Dividing both sides by r and rearranging the terms we get or equivalently λ r + 1 2r Let z = τ λ r + τ /2r. Then equation (91) is rewritten as where and Note that A is real-valued, whereas B can be a complex number depending on values of g(λ 0 ). When B = 0, the eigenvalue λ r is given by However, such a solution occurs only for certain values of the wave number k. When B = 0, equation (92) represents a special type of the equation that originally appeared in studies of the double-well Dirac delta potential model [18,25] and later studied in [54]. The solution of (92) is given by which is a quadratic Lambert W q function [19,35]. As discussed in Section 2.2 of [25], solution z can also be described as an ordinary Lambert W function provided that A and B satisfy Using these explanations we have the following theorem about the stability of w j .
Theorem 10. Let w j be a constant equilibrium of the general model (26). If A and B satisfy equation (98) and then w j is locally asymptotically stable.
Proof. Equation (97) is equivalent to The equilibrium w j is stable if (λ r ) < 0. When A < −1, the quadratic equation (98) has complex roots in the form of Hence, The condition A < −1 is equivalent to τ 2 1 r (D m k 2 + d m ) − 1 4r2 > 0, which must be satisfied for all k ∈ R. This is guaranteed by inequality (99), which completes the proof.

4.2.
The numerical algorithm. Here, we outline the numerical algorithm used for solving the initial value problem corresponding to the general model (26). To discretize the time variable, there are various numerical schemes. In this paper, we use the Houbolt method [17,43,54]. This scheme transforms a given timedependent PDE to a series of elliptical differential equations. Consider the following general time-dependent problem: with the boundary condition Bu(x, t) = h(x, t), x ∈ ∂Ω, and the initial condition u(x, 0) = u 0 , x ∈ Ω ∪ ∂Ω. To discretize the time variable, the equation above is reformulated as: The time domain of the PDE is discretized by the Houbolt finite difference method [17,43,54], which requires the Taylor series expansions of u k , u k−1 , and u k−2 at t k+1 as follows: where , and . After solving the system (105) -(107) for ∂u k+1 ∂t and ∂ 2 u k+1 Thus, (104) can be written as In order to fully implement the Houbolt method, we need to approximate the u values at the first three time steps. As a result, the Euler method with a very tiny step was implemented to obtain the initial values of these three time steps. The error of the above discretization is O((∆t) 3 ). If we desire to use an even higher-order time discretization scheme, we can expand a higher-order Taylor series expansion of u k , u k−1 , u k−2 , · · · , u k−l at t k+1 .
In every time-step, a general elliptical equation needs to be solved. The modified local method of approximated particular solutions (LMAPS) can be used to solve equations of this type. In the modified method, the spatial domain is discretized into a set of points in vector format {x i } N i=1 , where N represents the total number of discrete points in the domain and on the boundary. The improvement made is on choosing the polyharmonic splines as the basis. In particular, φ is the particular solution of Laplacian operator ∆ with respect to the polyharmonic splines ψ(r). Therefore, the basis φ is found by integrating polyharmonic spline ψ(r) = r 2m ln(r). With the particular solutions as the basis functions in the collocation technique on local domains, the original elliptic differential equations at each time step can be discretized into an N × N system of algebraic equations, and each equation in the system contains n non-zero entries. The system is a sparse system with N unknowns {û(x j )} N j=1 . Even though an additional polynomial basis needs to be added when we create local small linear systems in LMAPS using polyharmonic splines, the resulting global sparse system remains the same size as in the original LMAPS using other positive definite RBFs. We refer readers to [53] for more details.

Numerical simulations.
While several numerical simulations can be performed, we paid special attention to the stationary wave solutions of the general model (26). Starting with the stationary wave solutions of the reduced model (36), we investigated the effects of the dispersal delay r on the global stability of the stationary waves. Note that a stationary wave solution of the reduced model (36) is also a stationary wave solution of the general model (26). To be consistent with the previous simulations, we used the same birth function b(w) = pw 2 e −aw and the parameter values as those used in [3,6].
As shown in Figure 4, the stationary wave pulse of the HPRD model (26) loses its global stability when the value of dispersal delay r is increased from zero. The specific parameter values used to generate Figure 4   Similar to the previous case, Figure 5 shows that the stationary wavefront loses its global stability when the value of dispersal delay r is increased from zero. The specific parameter values used to generate Figure 5  In both cases of pulse and front, we see that for r > 0 large enough, lim w(t, x) = 0 as t → ∞. The biological implication of these simulations is that dispersal delay may directly influence the fitness of a single species population described by the stationary solutions of the general model (26). An increased dispersal delay in the general model (26) will increase the likelihood of population extinction.
5. Discussion. In this paper, we analyzed the dynamics of the delayed HPRD model (26) and its reduced forms. Table 2 provides a summary of the descriptions  and the main outcomes of the reduced models (see Section 3 and the references provided in the Table). These outcomes include, stability of the constant equilibria, delay-induced bifurcation, existence of traveling and stationary waves and spatiotemporal patterns of the population densities. Both reduced models (29) and (30) predict population survival for k > 0 and extinction for k = 0. Nevertheless, model (30) represents an extinction or survival that is spatially dependent, which is more realistic. As mentioned in Theorem 2, the reduced model (30) has a τ -periodic solution of the form w(t, x) = W 1 (t, x) + k/d m . Also, Theorem 4 provides the exact solution of the reduced model (41) with certain initial history functions. The general model (26) covers the outcomes of the reduced models and provides more insights as outlined below.
(a) Linear stability and Hopf bifurcation. A constant equilibrium of the general model (26) is stable under the conditions of Theorems 8 and 10. The biological implication of Theorem 8 is that in the absence of maturation time delay, the dispersal delay may destabilize the trivial equilibrium. Hence the dispersal delay can be used as a mechanism to prevent the extinction of single species. This coincides with the previous results [1] that the new sites are not populated instantaneously by invasive species. Instead, there is a delay between the time the species is introduced to the new environment until it starts populating the new sites. Under the conditions of Theorem 9, the constant equilibrium may lose its stability if τ >τ (r). The threshold valueτ (r) for maturation time delay is expressed as a function of the dispersal delay r. This is an improvement to the threshold value obtained in Theorem 3 of [6]. However, the functional dependency provided in equation (76) should be empirically validated. Ecological studies often relate the inter-patch movements of male individuals with inbreeding avoidance and those of female individuals with mother-daughter competition avoidance over the food resources (see [31] for a review). Hence, for certain species, equation (76) describing maturation delay as a function of dispersal delay seems to be more meaningful for male individuals.
(b) Global stability the stationary waves. The numerical simulations of the HPRD (26) indicate that the stability of stationary wave pulses and wavefronts are lost when the dispersal delay r is increased. The stationary wave solutions converge to zero as the dispersal delay increases. This result may explain the extinction of single species due to cold weather or other natural factors that increase the dispersal delay of single species.While we only focused on the numerical simulations and the local stability, the general model (26) has rich dynamics. Identifying the sufficient conditions for the global stability of wave solutions is the aim of our future studies. The existence of monotone wavefronts [3,29,48] and the global asymptotic stability of traveling waves [26,36,41,42] have been investigated for a number of delayed population models. Nonetheless, it remains an open problem to prove the global stability of wavefronts corresponding to the general model (26). The stability conditions may be interpreted as the conditions for conserving the endangered single species. In conclusion, the present work shows that population models with maturation and dispersal delays can give new insights into the complex dynamics of single species.
For each ξ, the eigenvalues λ 1 (ξ) and λ 2 (ξ) are the roots of the characteristic equation where D m , d m > 0 and ξ ∈ R. The functions c 1 (ξ) and c 2 (ξ) are determined by solving the system c 1 (ξ) + c 2 (ξ) = F (ξ) (114) where F (ξ) and G(ξ) are the Fourier transforms of f (x) and g(x), respectively. Applying the inverse Fourier transform to solution (112), the general solution of the model (30) with k = 0 is given by When k ≥ 0, the general solution of model (30) is given by w(t, x) = w h (t, x) + k/d m . From the characteristic equation (113) we get that Re(λ i (ξ)) < 0 for i = 1, 2, D m , d m > 0 and ξ ∈ R, which completes the proof.
Proof of Theorem 2. Using the method of separation of variables, we look for τ -periodic solutions of the form v n (t, x) = e −λnx g n (β n t − γ n x).
From (121) we get that Substituting (122) into (121) we get to which leads to solution (37). To make the solution (117) τ -periodic we need to have β n = 2πn/τ . Using the principle of superposition we get to the first series sum in solution (37). In equation (117), by changing −λ n to λ n and −γ n to γ n we may follow the same approach to generate the second series sum in solution (37).