Coexistence in two-species competition with delayed maturation

Inter- and intraspecific competition is most important during the immature life stage for many species of interest, such as multiple coexisting mosquito species that act as vectors of diseases. Mortality caused by competition that occurs during maturation is explicitly modelled in some alternative formulations of the Lotka–Volterra competition model. We generalise this approach by using a distributed delay for maturation time. The kernel of the distributed delay is represented by a truncated Erlang distribution. The shape and rate of the distribution, as well as the position of the truncation, are found to determine the solution at equilibrium. The resulting system of delay differential equations is transformed into a system of ordinary differential equations using the linear chain approximation. Numerical solutions are provided to demonstrate cases where competitive exclusion and coexistence occur. Stability conditions are determined using the nullclines method and local stability analysis. The introduction of a distributed delay promotes coexistence and survival of the species compared to the limiting case of a discrete delay, potentially affecting management of relevant pests and threatened species.


Introduction
Population control technologies aim to replace or to suppress vectors of diseases, such as arthropods that carry and transmit pathogens to animals or humans.Risk assessment of the large-scale implementation of population control technologies relies in part on mathematical modelling of the competitive dynamics between the species to be replaced or eliminated and the species introduced, such as transgenic constructs (Beeton et al. 2020(Beeton et al. , 2022;;Pagendam et al. 2020;Taghikhani et al. 2020).We are interested in competition between disease vectors such as mosquitoes, where maturation of the offspring is long compared to the lifespan of the adult.In this case, immature stages, hatching eggs and larvae foraging in a water reservoir are subject to far more intense competition for resources than adult mosquitoes, who have an abundant source of food (Noden et al. 2016).
Multi-species models of competition studied in the literature are often extensions of the continuous Lotka-Volterra competition model (Gause and Witt 1935;Hofbauer et al. 1987;Turchin 2013) or its discrete equivalent, the Leslie-Gower model (Leslie and Gower 1958).The Lotka-Volterra competition model incorporates density dependent growth that takes into account both intraspecific competition between individuals of the same species and interspecific competition between individuals of different species.Generalisations of the logistic Lotka-Volterra competition model, since the first successful applications to yeast and protozoans (Gause and Witt 1935), include non-linear competitive interactions (Gilpin and Ayala 1973) to capture growth behaviour for a wider diversity of species.Mechanistic approaches incorporate components of the processes that result in competition (Schoener 1976).For example, Beeton et al. (2020) included the process of choosing a mate when hybridising subspecies are competing.
Sexually mature individuals can be separated from immature individuals by using maturation delay, either in a discrete form assuming a single known delay period, or a distributed form containing a range of different delays in a known distribution.The end of the maturation period signifies that the individual is now able to breed.Modelling a maturation delay provides a simpler alternative to including age-structured subpopulations that require an extensive number of parameters to characterise each subpopulation.The kernel of a gamma distribution for a distributed delay stands as a more realistic option than a discrete delay (Blythe et al. 1984), allowing for more heterogeneity in representing the maturation by accommodating individual variation.Variability of the maturation duration could be a successful strategy for a species in response to a changing environment.
The insertion of a delay into the Lotka-Volterra competition model has been achieved by extending the delayed logistic model, called the Hutchinson-Wright model, to multiple species (Gopalsamy 1980;Gopalsamy and Aggarwala 1980;Smith 1995;Chen et al. 2010).The delay in the Hutchinson-Wright model causes growing, steady or declining oscillations around the equilibrium point (Murray 2002;Berezansky and Braverman 2003).Using the Hutchinson-Wright model as a template for adding a delay to the Lotka-Volterra competition leads to the existence of oscillations under some conditions.Biological populations are not easily described by oscillatory solutions, especially given sparse biological data.The oscillations observed in nature may not be directly related to maturation delay but may be due to external variations such as seasonal changes.
Here, we build on an alternative approach, proposed by Arino et al. (2006), that elaborates the model with a delay from a mechanistic point of view.This alternative to the delayed logistic model (Arino et al. 2006;Lin et al. 2018;Baker and Röst 2020) is obtained by replacing the birth term in the logistic growth by a function that depends on maturation delay.Arino et al. assume that the recruitment rate of immatures to the adult population is based on new births from the adult population N alive at past time t − τ , i.e.N (t − τ ), and on survival until the current time t.Death and intraspecific competition rates are instantaneous, and depend only on the current population.The number of remaining individuals that survive attrition during maturation is obtained by solving the equation without a birth term, dN (s)/ds = −m N (s) − a N 2 (s) at time t, using the known initial condition N (t − τ ).The full model in Arino et al. (2006) takes the following form, including the resulting birth term, The positive equilibrium, named the "delayed carrying capacity", exists only if the maturation delay τ of the species is below a critical threshold, implying that the population goes extinct if the delay is too long.This approach has also been applied to the Lotka-Volterra competition model, with a discrete delay (Lin et al. 2022), to show that evolution favours a short maturation delay in order to increase the carrying capacity of the population.
In this work, we build on the delayed Lotka-Volterra competition model for two species by representing the maturation delay with a truncated Erlang distribution.Including a truncation point allows for the case where a juvenile must either mature within a defined period or die.We consider a population of adults colonising a new environment, like the introduction of flying insects into a cage, or the beginning of the rainy season for mosquitoes in a seasonal environment.In our approach, mature and immature individuals together form a population at equilibrium in an environment where a perturbation is introduced by inserting or removing adults of one or two species.In the next section, we derive the equations of the model with a distributed delay.In Sect.3, we discuss the numerical results obtained, and we study the stability of the system.We complete the presentation of the results by discussing a biological example.
where Ni (t) represents the abundance of species i as a function of time t.All parameters are positive: the mortality rate mi , the reproduction rate ri , the carrying capacity Ki of the breeding environment, the density dependence effect âi of species i on its conspecifics, and the effect b j of species j on species i.It is assumed that ri > mi .The non-trivial equilibrium for the single species model when b j = 0 is Ki (r i − mi )/(r i âi ).
The model in ( 2) is often studied under two types of competition, or "struggle for existence" (Gause and Witt 1935).The case of weak interspecific competition is characterised by the condition where the two species compete for common food or resources while also belonging to different ecological niches, such that one species has particular resources that are not consumed by the other species.Weak interspecific competition can in some cases give rise to coexistence where both species survive.Strong interspecific competition is defined by the complement of (3) such that and represents a situation where two species compete strongly for common resources, leading to the less numerous species losing the competition and being excluded.We introduce the normalised variables where t = t, r i = ri and m i = mi .We now introduce maturation delays s i for i = {1, 2}.We assume that the population size of immature offspring decreases due to inter-and intraspecific competition.We also assume that adults experience density independent mortality and negligible competition.We reformulate the Lotka-Volterra competition model for two species in (5) by generalising the expression of the delayed growth for one species calculated in (1), such that where the expression r i R i (t) is the number of immatures born (e.g.eggs) at time t − s i that survive to maturation and emerge as adults at time t.The birth rate of new immatures at time t − s i is assumed to be proportional to the size of the adult population at time t − s i , where r i is the number of hatching offspring per adult, such that R i (t Density dependent mortality then applies to these offspring from time t − s i until time t.We find R i (t) by evaluating where 1/R i dR i is the per capita rate of change of the immature population born at time t − s i and expected to emerge at t, where μ i > 0 is the rate of mortality of the immature offspring, and a i N i (s)+b j N j (s) is the density dependent mortality specified in Eq. ( 5).The left hand side of Eq. (7 )) and we then obtain We also define the function f i (t) as We can then substitute Eq. ( 8) into Eq.( 6), and use Eq. ( 9) to define a coupled ODE model for two species with a discrete delay, such that where (i, j) = {(1, 2), (2, 1)}.The function f i (t) represents the mean population density during maturation and d f i /dt when negative corresponds to the mean population loss during maturation due to intra-and inter-specific competition across the maturation period.The initial conditions N i (t) = φ i (t) > 0 for t ≤ 0, are bounded continuous functions mapping from (−∞, 0] to R + .The initial conditions, called history, describe a scenario where the population of adults at time t = 0 is positive and where the configuration of the population at t < 0 is at equilibrium, including the trivial equilibrium where N i (t < 0) = 0.A perturbation occurs at t = 0 when the population of adults increases or decreases.We aim to convert the discrete delay described in Eq. ( 11) into a distributed delay, where maturation time varies randomly between individuals according to a defined probability distribution.The kernel used is the probability density function (PDF) of the gamma distribution where p ≥ 0 is the shape parameter of the distribution, β > 0 is the rate parameter of the distribution and ( p) is the gamma function.The area under the curve in ( 12) is ∞ 0 G (s | β, p) ds = 1 and the average delay is ∞ 0 sG (s | β, p) ds = p/β.The limit of the probability density function in (12) when p → ∞ is a shifted Dirac distribution, δ(s − p/β), corresponding to the kernel of the discrete delay (Smith 2011).
It is more relevant computationally to use a bounded support, consistent with a finite interval of integration.A limited support of the distributed kernel signifies biologically that there is an upper limit to the amount of time necessary for an individual to mature and that no maturation is possible beyond this limit.The threshold for the maturation period is an intrinsic characteristic of the species.We use the reflectedshifted-truncated gamma distribution, as studied in Waymyers et al. (2018), to bound the support of the delay kernel.The reflected-shifted-truncated gamma distribution is defined by the probability density function where > 0 is the shift and u = − s, ( p, β ) is the incomplete upper gamma function, such that (x, a) = ∞ a t x−1 exp (−t) dt for any real x > 0. The PDF in (13) is normalised so that similarly the area under the curve The truncated distribution used for a kernel with bounded support appears also in Vittadello et al. (2021) to describe the distributed delay in a heterogeneous biological population of cells.
Figure 1 shows the skewness of the probability density function in (12) for different values of p, represented by a pink shape for a mean delay p/β = 2 and by a green shape for a mean delay p/β = 3.The peak or mode of the distribution is at the location s = 0 when p = 1 (Fig. 1a).The mode moves away from s = 0 towards the position p/β when p = 2, 10 and 20 (Fig. 1b-d).The probability density functions in Fig. 1e-h are obtained from the PDFs in (a)-(d) after a reflection over the y-axis, followed by a shift of = 4 to the right and a left truncation at s = 0. Classical models using delays for biological populations favour shifted distributions (Blythe et al. 1984) as they assume a minimum mandatory time for maturation.The reflectedshifted-truncated gamma distribution can approximate the minimum mandatory time for maturation by using p 1 as in Fig. 1d, h where the probability that the delay occurs between u = and u = − 1 is nearly zero, corresponding to a minimum mandatory time of maturation equal to one.
The discrete delay described in Eq. ( 11) in the form (Smith 2011) taken over all possible delays s, where δ is the Dirac function.We can use this formulation to convert our discrete delay into a distributed delay using the reflected-shifted-truncated gamma distribution in Eq. ( 13).The delay N (t −s i ) can then be replaced by where the use of s i in the denominator of Eq. ( 11) is now replaced with a scaling factor i .Note that the distributed delay has been applied independently to N i (t) and f i (t), so that density dependent mortality of immatures is independent of the maturation period u.Distributed delays are well studied in the literature; for example, the alternative logistic model ( 1) for one species is generalised with a distributed delay in Lin et al. (2018).However, our work is different from previous mechanistic models that consider a maturation delay.To our knowledge, this is the first analysis of a distributed delay model of two species where there is no competition among adults.
The linear chain approximation (MacDonald 1978;Smith 2011;Wolkowicz et al. 1997) allows us to transform a system of delay differential equations (DDEs) with a distributed delay represented by a gamma distribution into a system of ordinary differential equations.We define the term containing the distributed delay in the right hand side of dN i /dt as where 16) is the normalising constant that originates from the PDF of the truncated gamma distribution.The PDF G (u | β i , p i , μ i ) associated with the delay u must be multiplied by the normalising constant as the distribution is truncated.We will then express the derivative of the PDF as a linear combination of PDFs.The kernel function in Eq. ( 17) satisfies the initial value problem where k = 2 . . .p i .We aim to work with a finite number of ordinary differential equations (ODEs), such that p i ∈ N + and p i ≥ 1.The results in this work are obtained by considering G (u | β i , p i , μ i ) as the Erlang distribution, a special case of the gamma distribution, where the shape of the distribution p i is discretised such that The Erlang distribution with a shape p i and a rate β i gives the elapsed time until a chain of p i independent events of maturation have occurred at a rate β i .
The integral in Eq. ( 16) is defined as the function x i, p i (t) in the system of equations where k = 1 . . .p i and x i,k (t) is a continuous bounded function mapping from (−∞, 0] to R. The same process is applied to the term containing the distributed delay in the right hand side of d f i /dt that we define as the function y i, p i (t) from the system of equations where k = 1 . . .p i and y i,k (t) is a continuous bounded function mapping from (−∞, 0] to R. The system of ODEs made up of the derivatives of x i,k (t) with respect to t can be obtained using the Leibniz integral rule, by applying the properties of the derivative of the convolution and by using Eqs.( 18)-( 19).The same approach yields the system of ODEs constituting the derivatives of y i,k (t) with respect to t.We obtain the following system of 2( p 1 + p 2 + 2) ordinary differential equations Fig. 2 Numerical solutions illustrating survival, coexistence and competitive exclusion.The numerical solutions of the system of ODEs ( 22)-( 27) are shown starting with initial conditions for (i, j) = {(1, 2), (2, 1)}, where k = 2 . . .p i and e p i −1 (

Results and discussion
In this section, we present the numerical solutions of the system of ODEs ( 22)-( 27) and the conditions of existence of the equilibrium points.We combine the method of nullclines and some results from local stability analysis to understand the features of the system at equilibrium.

Numerical solutions of the system of ODEs
The numerical solutions of the system of ODEs ( 22)-( 27) are obtained using the Julia programming language; specifically with the solver BS3 in the package Differ-entialEquations.jl (Rackauckas and Nie 2017).This solver is an implementation of the Bogacki-Shampine method (Bogacki and Shampine 1989), an explicit three-stage, third order Runge-Kutta method with an adaptative step size for the time domain.The tolerance is set to 1 × 10 −6 to obtain the numerical results presented in this work.The code in Julia that solves the system of ODEs ( 22)-( 27) is available on GitHub.
Figure 2 shows numerical solutions N 1 (t) and N 2 (t) from ( 22)-( 27) with parameters m 1 = m 2 = μ 1 = μ 2 = 0.3, r 1 = r 2 = 1 and a 1 = a 2 = 1 constant accross all cases.The parameters b 1 , b 2 , β 1 , β 2 , p 1 , p 2 , 1 and 2 are chosen so that the resulting solutions illustrate extinction of one species in Fig. 2a, stable coexistence in Fig. 2b, and competitive exclusion in Fig. 2c.Note that these parameters chosen for Fig. 2 do not represent known biological species, though an example for two known species using biologically plausible parameters is discussed in Sect.3.4.The truncation position is chosen to be at the 99th percentile of each maturation delay distribution (where the cumulative distribution function is equal to 0.99).The survival of species 1 in Fig. 2a is achieved by choosing β 1 = 1 when p 1 = 1 (blue line) and β 1 = 20 when p 1 = 20 (green line), and the extinction of species 2 is obtained by choosing β 2 = 1/8 when p 2 = 1 (red line) and β 2 = 20/8 when p 2 = 20 (yellow line).Comparison of Fig. 2a (where species 2 goes extinct) and Fig. 2b (where species 2 survives) shows that the rate β 2 must be increased from β 2 = 1/8 in Fig. 2a to β 2 = 1/2 in Fig. 2b when p 2 = 1, and from β 2 = 20/8 to β 2 = 10 when p 2 = 20, to ensure survival.The criterion for weak interspecific competition as defined in Eq. ( 3) is used in 2b, where the maturation period is short enough that the population gain due to birth is able to compensate for the loss due to mortality, resulting in stable coexistence of both species.Conversely, competitive exclusion in Fig. 2c is obtained by using the criterion for strong interspecific competition in Eq. ( 4).Species 1 wins the competition in Fig. 2c where the rate β 1 = 1 is greater than β 2 = 1/2 when p 1 = p 2 = 1 and β 1 = 20 is superior to β 2 = 10 when p 1 = p 2 = 20.We do not show the total extinction of both species that correspond to a maturation delay too long for each species to survive.The condition that relates the survival of each species to the distributed maturation delay is given in Sect.3.2.1.
The Chi-squared distribution is a special case of the gamma distribution ( 12) with a rate of 1/2 and a shape of k/2, where k are the degrees of freedom, as illustrated in pink in Fig. 1a for a shape of one.Figure 3 presents the numerical solutions N 1 (t) and N 2 (t) from ( 22)-( 27) where β 1 = β 2 = 1/2 as in the Chi-squared distribution and p i = k/2.We show how the solutions are different when the shifts of the truncation 1 and 2 correspond to the first quartile, the second quartile or the 99th percentile of the Chi-squared distribution, equivalent to the location where the cumulative distribution is 0.25, 0.50 and 0.99.We obtain coexistence for any quantile shown in (a)-(c) when p 1 = 2 and p 2 = 3.We obtain the survival of only species 1 when p 1 = 2 and p 2 = 5 for any quantile shown in (d)-(f).Figure 3 illustrates how the type of equilibrium does not change with the quartile used for the choice of 1 or 2 .The magnitudes of N 1 and N 2 vary significantly with 1 and 2 when the truncation chosen corresponds to a percentile lower than the 99th percentile.The difference indicates that the truncation must be chosen carefully so it corresponds to a high percentile of the gamma distribution.

Stability analysis
It can be easily shown that the solutions N 1 (t) and N 2 (t) are limited to the positive quadrant and bounded.By examining Eqs. ( 22)-( 27), we can also readily show that any value of f i (t) = fi at equilibrium is proportional to N i (t) = Ni at equilibrium such that fi = p i 1 − exp (−β i i ) e p i −1 (β i i ) /(β i i ) Ni , and that xi = ȳi = Ni .We can then directly examine the equilibrium points of Eqs. ( 22)-( 27) in the plane (N 1 , N 2 ).We can further show that the subsystem of Eqs. ( 22)-( 24) being Fig. 3 Numerical solutions at first and second quartiles, and at 99th percentile of the Chi-squared distribution.The solutions of system of ODEs ( 22)-( 27) are shown starting with initial conditions at equilibrium implies that the subsystem of Eqs. ( 25)-( 27) is also at equilibrium, and vice-versa.As a result, only Eqs. ( 22)-( 24) need be studied to determine the equilibrium states, and we thus consider that f i (t) and y i,k (t) are not relevant to describe the history in regards to equilibria.For numerical simulations, we choose initial conditions where either f i (0) = fi and y i,k (0) = ȳi,k or f i (0) = 0 and y i,k = 0, with N i (0) > 0 and x i,k (0) ≥ 0.

Equilibrium points and nullclines
We present the existence conditions of the equilibrium points of Eqs. ( 22)-( 27).Two outcomes can be distinguished visually in Figs. 2 and 3: the boundary equilibrium where one species is extinct while the other survives, and the coexistence equilibrium where both species coexist.The boundary equilibrium where Ni > 0 and N j = 0, for (i, j) = {(1, 2), (2, 1)}, satisfies the condition that corresponds to equilibrium points We can use l'Hôpital's rule to find that the limit of − ln The result is useful to understand what happens to the equilibrium point when the kernel of the distributed delay approaches the kernel of a discrete delay δ(s − p i /β i ).The abundance at the semitrivial equilibrium ( 28) is such that Ni → ∞ when p i /β i 1, exceeding (r i − m i )/(r i a i ), the expected carrying capacity for the single species model without a delay.The reason why the abundance in the delayed model exceeds that of the single species model with no delay is that density dependence and interspecific competition take effect in the space where the offspring are maturing, and are considered negligible in the space where the adults live.The offspring maturing in a short delay would very quickly vacate the space occupied, and would liberate resources for new offspring.A maturation delay approaching zero would thus mean a carrying capacity approaching infinity.
The condition of survival for species i in (29) corresponds to the critical condition ln Extinction occurring when condition (30) is false is not caused by competition.The species is simply extinct because the maturation period is too long for the species to survive at the current rates of reproduction and mortality.We examine the following conditions for the existence of the equilibrium where both N1 and N2 are positive, with where α i is the inverse of the normalising constant of the PDF of the truncated Erlang distribution such that α i → 1 when i → ∞.The equilibrium point N1 , N2 following the conditions (31)-(32) correspond to with Fig. 4 Phase plane with stable coexistence equilibrium.The nullclines obtained from ( 31)-( 32) are shown in blue and in green for N 1 and in red and in yellow for N 2 , when p 1 = p 2 = 1 and when p 1 = p 2 = 20 respectively, using where α i is defined in (33).The existence of (34) assumes that the boundary or semitrivial equilibrium (29) exists or that condition (30) is true for (i, j) = {(1, 2), (2, 1)}.
Again, as we did in (29), we can use l'Hôpital's rule and the assumption that p i /β i is constant, to find that A → (μ i p i /β i ) ln(m i /r i ) in the limit where p i → ∞.
The trivial equilibrium and the semi-trivial equilibria are global asymptotic equilibria if the threshold condition for survival in Eq. ( 30) is not fulfilled for one or two species (see Lyapunov function in Appendix A).
We study the equilibrium point (34) in the cases of weak and strong interspecific competition.We are studying a projection of the solution on the plane (N 1 , N 2 ) from a phase space that represents the dynamical system from Eqs. ( 22)-( 27).The plane (N 1 , N 2 ) is appropriate to study the solution around the equilibrium points as the variables x i,k and y i,k follow directly from the variables N i , and as f i is simply the mean of N i during maturation.The variables x i,k , y i,k and f i represent what is happening during maturation, while N i is the population at the current time.The equilibrium point in (34) exists for weak interspecific conditions if all the following conditions are true: The stable coexistence equilibrium when the interspecific competition is weak is shown in the phase planes of Fig. 4. The nullclines illustrated in blue for N 1 and in red for N 2 , when p 1 = p 2 = 1, are obtained from Eqs. ( 31) and ( 32) respectively.The intersection point of the nullclines in the positive quadrant corresponds to the stable coexistence equilibrium.Any solution starting at N 1 > 0 and N 2 > 0 will move towards the coexistence equilibrium when the conditions for its stability and existence are fulfilled.An example is shown when p 1 /β 1 = p 2 /β 2 = 2 for p 1 = 1 in Fig. 4a where N1 = N2 , and when p 1 /β 1 < p 2 /β 2 for p 1 = 1 in (b) and (c).N1 is being advantaged over N2 in (b) and (c).Figure 4 shows the nullclines in green for N 1 and in yellow for N 2 when p 1 = p 2 = 20.A closer look at the green and yellow lines in (c) shows that the intersection of the two lines is outside the positive quadrant.We conclude that there exists p 1 = p 2 1 such that the species coexist for p 1 = p 2 = 1 5 Phase plane with unstable coexistence equilibrium.The nullclines obtained from ( 31)-( 32) are shown in blue and in green for N 1 and in red and in yellow for N 2 , when p 1 = p 2 = 1 and when p 1 = p 2 = 20 respectively, using r 1 = r 2 = 1 with the ratios p 1 /β 1 = 2 and p 2 /β 2 = 2, 3 and 3.25 in (a)-(c) respectively.The vector field (grey arrows), the trajectories solutions (grey curves) and the equilibrium points (black discs) are shown for p 1 = p 2 = 1 (colour figure online) and such that species 1 excludes species 2 for p 1 = p 2 1, where the ratios p 1 /β 1 and p 2 /β 2 are maintained constant for p 1 = p 2 = 1 and for p 1 = p 2 = 20.
The equilibrium point in (34) exists for strong interspecific competition if all the following conditions are true: The conditions in (37) are obtained by reversing each inequality in (36).Figure 5 shows an example of bistability where the interspecific competition is strong and the coexistence equilibrium point is unstable.The asymptotic equilibrium depends on the starting point of the trajectory.The size of N2 in the equilibrium point ( N1 = 0, N2 > 0) is equal to that of N1 in the equilibrium point ( N1 > 0, N2 = 0), in Fig. 5a, as the unstable coexistence equilibrium is such that N1 = N2 .The size of N2 in the equilibrium point ( N1 = 0, N2 > 0) is smaller than N1 in the equilibrium point ( N1 > 0, N2 = 0), in Fig. 5b and c.The green and yellow lines in Fig. 5 representing the nullclines for N 1 and N 2 respectively show that the unstable coexistence equilibrium also exists in Fig. 5a-b when p 1 = p 2 = 20.Figure 5c shows that there exists p 1 = p 2 1, with the ratios p 1 /β 1 and p 2 /β 2 maintained, such that the unstable coexistence equilibrium exists for p 1 = p 2 = 1, where the winning species depends on the initial conditions, and such that the unstable coexistence equilibrium disappears for p 1 = p 2 1, where species 1 excludes species 2.
The strong interspecific competition conditions also include a 1 a 2 = b 1 b 2 , where the coexistence equilibrium would exist only if the nullclines for N 1 and N 2 are superimposed, meaning that there is an infinity of equilibrium points.

Local stability
The system of Eqs. ( 14)-( 15) is reduced to Eq. ( 14) for i = 1, 2 and linearised around equilibrium points ( N1 , N2 ) from Sect.3.2.1 where we substitute the exponential solution [u(t) v(t)] T = exp (λt) [b c] T .We use the properties of the Erlang function to determine Smith (2011).The characteristic equation can be obtained with the determinant of the following 2-by-2 matrix where , and where f1 and f2 can be expressed in terms of N1 and N2 respectively, such that f1 = τ 1 N1 1 and f2 = τ 2 N2 2 , and assuming that The complete expression of the characteristic equation is given in Appendix B, Eq. ( B16).Here, we give the specific characteristic equations around the trivial and semi-trivial equilibrium points and the corresponding stability conditions.
The eigenvalues when ( N1 , N2 ) = (0, 0) are obtained by the characteristic equation The equilibrium point ( N1 , N2 ) = (0, 0) is stable if ln(m 1 ) − ln(r 1 ) + p 1 [ln(β 1 + μ 1 ) − ln β 1 ] > 0 and ln(m 2 ) − ln(r 2 ) + p 2 [ln(β 2 + μ 2 ) − ln β 2 ] > 0, corresponding respectively to A 1 > 0 and A 2 > 0, and the equilibrium point is unstable otherwise, meaning that if both species are not fit for survival, both species are extinct, and if both species are fit for survival, the solution will move away from the equilibrium ( N1 , N2 ) = (0, 0) to competitive exclusion or coexistence.The characteristic equation becomes when N 1 = 0 and N 2 > 0. Equations ( 40) and ( 41) can be expressed as M i,1 (λ)M i,2 (λ) = 0 where M i,1 (λ) correspond to the first line of the right hand side of the equation and M i,2 (λ) is the second line of the right hand side of the equation.M i,1 (λ) = 0 have roots with negative real parts if the critical condition of survival (30) of the species i is fulfilled.M i,2 (λ) = 0 in Eq. ( 40) have roots with negative real parts if τ 1 a 1 A 2 ≥ τ 2 b 1 A 1 , and M i,2 (λ) = 0 in Eq. ( 41) have roots with negative real parts if τ 2 a 2 A 1 ≥ τ 1 b 2 A 2 .We conclude that the conditions of stability of existing equilibrium points ( N1 , 0) or (0, N2 ) correspond partially to the conditions of existence (37) of the unstable coexistence equilibrium.At the opposite, the conditions of stability that make equilibrium points ( N1 , 0) or (0, N2 ) unstable correspond partially to the conditions of existence (36) of the stable coexistence equilibrium.
We give an example of how to study local stability when The eigenvalues are all real and negative if m 1 (β 1 + μ 1 ) > r 1 β 1 and m 2 (β 2 + μ 2 ) > r 2 β 2 equivalent to A 1 > 0 and A 2 > 0, where A i is defined in (35).The characteristic Eq. ( 40) when p 1 = p 2 = 1 is in the form M 1,1 (λ)M 1,2 (λ), where M 1,1 (λ) is a polynomial of third order and M 1,2 (λ) is a polynomial of second order, that can be studied using the Routh-Hurwitz criterion for stability.The roots of polynomial M 1,1 (λ) all have negative real parts if the conditions in (28) are fulfilled, and the roots of M 1,2 (λ) all have negative real parts if 37) when α 1 and α 2 approach one.The polynomial of order six corresponding to the characteristic equation for the coexistence equilibrium is given in Eq. (B17).The Routh-Hurwitz criterion can be verified with the help of symbolic software.
We now give the characteristic equations when p 1 → ∞ and p 2 → ∞ and when p 1 /β 1 and p 2 /β 2 are constant.We use l'Hôpital's rule to determine that g We use the theorem that states that z has a negative real part if c + d < 0 and d ≥ c, as proved in Hayes (1950) and in Smith (2011).We find that the eigenvalues have a negative real part if ln(r 1 /m 1 ) + p 1 /β 1 μ 1 > 0 and ln(r 2 /m 2 ) + p 2 /β 2 μ 2 > 0, meaning that ( N1 , N2 ) = (0, 0) is stable if the critical conditions for survival for species 1 and 2 are not fulfilled and unstable otherwise.The characteristic equation around ( Ni > 0, N j = 0) can be separated into two transcendental equations.The first transcendental equation is in the form We verified that the three following conditions necessary for absolute stability (Brauer 1987;Smith 2011) are satisfied for every value of the delay β i / p i when N i > 0 and N j = 0.The first condition is that the real part of the roots of P(λ) must be positive or zero.The second condition is The characteristic equation around the coexistence equilibrium point takes a more complicated form (see Appendix B).Such forms of transcendental equations with two discrete delays are approached in An et al. (2019).We remind the reader that the competition between two species with a discrete delay is also studied in Lin et al. (2022).

Distributed delay promotes coexistence and survival
It is useful to understand how the equilibrium for a species moves from extinction to persistence depending on the parameters of the model.We use the following method to determine the values of N1 and N2 depending on the parameters of the model.We determine if the species i can survive by verifying if condition (30) is true.The species that does not pass the test in (30), independently of competition, is extinct and the species that does pass the test survives if the other is extinct.If condition (30) is true for both species, we establish the existence of the coexistence equilibrium point (34), then we determine if the coexistence equilibrium is stable (36) or unstable (37).It is possible that the coexistence equilibrium does not exist in the positive quadrant.We evaluate the intersection of the nullclines (31) and (32) to determine the outcome of competitive exclusion.An intersection in the quadrant where N1 is positive and N2 is negative means that species 1 excludes species 2 if a 1 a 2 < b 1 b 2 , and that species 2 excludes species 1 otherwise.An intersection in the quadrant where N1 is negative and N2 is positive means that species 2 excludes species 1 if a 1 a 2 < b 1 b 2 , and that species 1 excludes species 2 otherwise.
Figure 6 gives examples of how coexistence, competitive exclusion and extinction depend on the shape p i .Figure 6 is subdivided in 9 different cases where p 1 = 1, 2 and 20, and p 2 = 1, 2 and 20.Each case shows the region of coexistence in brown, the region where species 1 wins in blue, the region where species 2 wins in red and the region where both species are extinct in black.The regions are illustrated on a map where p 1 /β 1 is on the horizontal coordinate and where p 2 /β 2 is on the vertical coordinate.Using the functional form p i /β i for our unfolding parameters allows us to compare the results of the model with a distributed delay to the model with a discrete delay, as p i /β i represents the maturation duration in the discrete delay kernel when p i → ∞.
As p 1 /β 1 and p 2 /β 2 both vary from 1/10 to 8, as a result β 1 and β 2 vary from 10 to 1/8 in Fig. 6a, and both species are extinct in the region where β 1 < 1/7.8 and β 2 < 1/7.8, or where p 1 /β 1 > 7.8 and p 2 /β 2 > 7.8.The region of coexistence shown in brown in Fig. 6a is approximately four times larger than the region of coexistence delimited by the white curve representing the region of coexistence for the discrete delay.In Fig. 6i, where β 1 and β 2 vary from 200 to 20/8, both species are extinct in the region where p 1 /β 1 > 4 and p 2 /β 2 > 4, or where β 1 < 5 and β 2 < 5.The region of coexistence shown in brown in Fig. 6i approximates the corresponding region delimited by the white curve for a discrete delay.The species with the shortest maturation is advantaged when a species excludes another species in Fig. 6.The parameter p i plays a role in the persistence of the species, as it corresponds to the number of events expected to occur at fixed rate of maturation.Species 1, for example, is advantaged at the expense of species 2 in Fig. 6b where p 1 = 1 and p 2 = 2.
We did not present the equivalent of Fig. 6 for strong interspecific competition.The regions of competitive exclusion (red or blue) and total extinction (black) shown in Fig. 6 for a 1 = a 2 = 1 and b 1 = b 2 = 0.5 would have the same shapes and areas if b 1 = b 2 = 2.The region of coexistence in brown would become a region of bistability.The effect of species 2 on species 1, represented by the parameter b 1 , and the effect of species 1 on species 2, represented by the parameter b 2 , play a role in persistence and extinction, as shown for a 1 a 2 < b 1 b 2 in Fig. 7, where p 1 = 1.The region of stable coexistence, coloured in brown, is larger when b 1 = b 2 = 0.1 in Fig. 7a whereas the region is much smaller when b 1 = b 2 = 0.8 in Fig. 7c.Species 2, shown in red, is advantaged at the expense of species 1 where b 1 = 0.1 and b 2 = 0.8 in Fig. 7c.Species 1 is advantaged over species 2 where p 2 = 20 in (g)-(i) when compared to the results in (a)-(c) where p 2 = 1, as the rates β 2 used in (g)-(i) are much smaller than the rates used in (a)-(c).
We can also interpret the results in Fig. 7 by comparing the regions of coexistence (brown) and survival or competitive exclusion (red and blue) to the regions delimited by the white curves representing the discrete delay.All the cases show that the distributed delay promotes coexistence of both species and the survival of a species.The region of total extinction for the discrete delay, delimited by the white lines, would be larger than the narrow black rectangle that represents the region of extinction for the distributed delay.The regions in Fig. 7 f) and (i).The region in brown would represent the bistable equilibrium of strong interspecific competition.Figures 6 and 7 highlight the flexibility offered by the parameters of the distributed delay allowing for more coexistence and survival compared to the discrete delay.The results are shown for a restricted set of parameters in Figs. 6 and 7. Our analysis and an extensive exploration of the results for a larger domain of parameters show that the following observation remains true: the threshold for survival depending partially on the expression [(β i + μ i )/β i ] p i is always smaller for p i = 1 than for p i → ∞ if p i /β i remains constant.

Biological example
We give an applied example of the model to show possible outcomes with real biological parameters.Two species of mosquitoes larvae, Anopheles gambiae sensu stricto and Anopheles arabiensis, vectors of malaria in Africa, can be found in the same reservoirs of water where they mature.We use assumptions about the biology of both species from the literature (Kirby and Lindsay 2009;Paaijmans et al. 2009).Species 1, An. gambiae s.s., has a shorter maturation than species 2, An.arabiensis.Species 1 has also the biggest carrying capacity, as the female of species 1 reaches a smaller size than the female of species 2 at full growth.The maturation delay could be shortened for both species, in warmer water ponds for example (Kirby and Lindsay 2009;Agyekum et al. 2022).
Figure 8 shows the populations of An. gambiae s.s.(N 1 (t)) and A. arabiensis (N 2 (t)), where the non dimensional N i (t) = 1 corresponds to Ni (t) = Ki .The maturation in Fig. 8 is delayed and the distribution of individual delays within the population of one species is heterogeneous, at opposed to the classic Lotka-Volterra model where the maturation is instantaneous and homogeneous.We selected maturation delays and mortality rates from Kirby and Lindsay (2009) that are physically plausible for both species.The parameters β 1 and β 2 in Fig. 8 are chosen so that τ 1 = 13 and τ 2 = 14 in (a), τ 1 = 11 and τ 2 = 12 in (b).The shifts corresponding to the truncation are chosen so that α 1 = α 2 = 0.99: this gives shifts 1 = 59.9 and 2 = 64.5 in (a) for p 1 = p 2 = 1, and 1 = 20.7 and 2 = 22.3 in (a) for p 1 = p 2 = 20.
The potential for coexistence that arises from weak interspecific competition is illustrated in Fig. 8a, b when p 1 = p 2 = 1: for longer delays in (a), and for shorter delays but higher mortality rates in (b).Coexistence fails in both cases when p 1 = p 2 = 20: in (a), A. arabiensis becomes extinct, and in (b) both species reach extinction.Coexistence of An. gambiae (blue) and A. arabiensis (red) is possible in Fig. 8b when The solutions in Fig. 8c and d are obtained with the same parameters as in (a) and (b) respectively, except for a 1 , a 2 , b 1 and b 2 that we modified to satisfy the condition for strong interspecific competition.An. gambiae s.s.excludes A. arabiensis in Fig. 8c. A. arabiensis wins the competition with appropriate initial conditions in Fig. 8d only when p 1 = p 2 = 1.

Conclusion
We generalised the alternative delayed Lotka-Volterra model by reformulating the model with a distributed delay.The kernel of the Erlang distribution was used for a more realistic representation of the heterogeneity found in the length of maturation in a species.The distributed delay separated the mature from the immature individuals, to represent a species where competition is more important for immature individuals and where maturation time is long compared to lifetime.We used the linear chain trick to approximate the DDEs of the model by a system of ODEs.We showed how the survival of a species depends on the rate of maturation being able to compensate for the rate of loss due to mortality of adults and immature individuals.A species fit for survival enters into competition with another species, leading to competitive exclusion, to coexistence of both species, or to a bistable equilibrium.The necessary conditions for stability were shown for the equilibrium points in function of the parameters of the model.We also gave sufficient conditions to determine a global asymptotic equilibrium.We highlighted differences of the model with a distributed delay as compared to the discrete delay.The model has the overall effect of promoting coexistence and survival of the species.The system is able to induce stable coexistence for a skewed distribution of the maturation delay, whereas a competitive exclusion would be expected for a discrete maturation delay, such that a model that inaccurately assumes a discrete delay would give a completely incorrect result.Another difference of the model is that the chances for species with less fitness to "win" in competition are increased in the bistable equilibrium for a skewed distribution, whereas the extinction of the species would be expected for a discrete maturation delay, as showed in the biological example in Sect.3.4.The model could be extended for three or more species, so we could study whether the model property of allowing more coexistence is maintained in higher dimensions.Future work would also require study of the transient oscillations to see how well they describe the transient behaviour of a population at equilibrium after a perturbation is introduced.
where g(x) is a nonlinear function with equilibrium also at x = 0, defined by Then if we define our Lyapunov function V as given a positive-definite matrix B, then Eq. (A3) above is satisfied by definition.As A is a lower bidiagonal matrix, its eigenvalues λ are equal to its diagonal values, i.e. and is negative-definite as required (See the Mathematica file on GitHub).The full equation for V is thus: For a periodic sequence of real numbers {c 1 , . . ., c n , . ..}where c n+1 ≡ c 1 , n i=1 c i c i+1 , (A15) so using Eq.A15 with the sequences {N i − Ni , x i,k − Ni , . . ., x i, p i − Ni , . ..} for i = 1, 2 gives V (x) < 0 as required if

Fig. 6
Fig. 6 Extinction, survival and coexistence as a function of p 1 and p 2 .The region of competitive exclusion is represented in blue if N1 > 0 and in red if N2 > 0. The region of coexistence where N1 > 0 and N2 > 0 is represented in brown.The region where both populations are extinct is represented in black.The conditions for existence and the stability of each region are obtained with μ 1 = μ 2 = m 1 = m 2 = 0.3, a 1 = a 2 = 1, b 1 = b 2 = 0.5, r 1 = r 2 = 1, p 1 = 1 in (a)-(c), p 1 = 2 in (d)-(f), p 1 = 20 in (g)-(i), p 2 = 1 in (a), (d) and (g), p 2 = 2 in (b), (e) and (h), and p 2 = 20 in (c), (f) and (i).The axis representing p 1 /β 1 and p 2 /β 2 are shown from 0.1 to 8. The white curve corresponds to the frontiers of each region in the limit where p 1 → ∞ and p 2 → ∞, equivalent to the kernel δ(s i − p i /β i ) of a discrete delay (colour figure online)