Oscillatory regimes in a mosquito population model with larval feedback on egg hatching

Understanding mosquitoes life cycle is of great interest presently because of the increasing impact of vector borne diseases in several countries. There is evidence of oscillations in mosquito populations independent of seasonality, still unexplained, based on observations both in laboratories and in nature. We propose a simple mathematical model of egg hatching enhancement by larvae which produces such oscillations that conveys a possible explanation. We propose both a theoretical analysis, based on slow–fast dynamics and Hopf bifurcation, and numerical investigations in order to shed some light on the mechanisms at work in this model.


Introduction
Today numerous areas of the world are severely affected by mosquito-borne viral diseases, with notable examples including dengue, chikungunya and Zika (see [5]).Scientists are hard at work to find new and efficient ways to mitigate the impact of, or even eradicate these arboviral infections, and especially target vector control.
A beneficial implementation of any of the vector control methods requires a good understanding of the local vector population's bio-ecology, and a reliable monitoring of its dynamics.To achieve better knowledge, this monitoring needs not only be demographic (using trap counts), but can also use genetic data -for the example of Wolbachia see [12] and [26].However, studies in Rio de Janeiro throughout the past decade have shown that monitoring urban populations of Aedes aegypti is a difficult task (see [13], [7], [23]), largely because of environmental variations (spatial heterogeneity, seasonality, etc.).
A first -to the best of our knowledge -systematic comparison of two complex models of Aedes aegypti population dynamics, relevant for a control program, was done recently in [16].Another application of proper modeling of mosquito's life-cycle is the risk estimation for disease emergence (see [11]).
We believe that the intrinsic life cycle of Aedes aegpyti may still be improperly modeled, and effort should be put in the direction of integrating several key features in the models.Among these features, we have in mind the transitions between the stages (egg, larva, pupa, adult) or even within these stages (larval instars, etc.) because in theory, any of these transitions (ovipositing behavior, hatching, pupation, mating, etc.) can give rise to nonlinearity.Nonlinearities ought to be taken into account when using collected data, so that they do not blur the picture we get of the actual population's dynamics.In addition, synchronizing or de-synchronizing effects, either in time or space, are possible outputs of these nonlinearities, and can result in variations in crucial traits of the mosquito populations, such as vector capacity (see [14], [4]) We focus exclusively in this work on one single aspect of the evolution of the mosquito population, setting the hypothesis that the larval density in breeding sites directly impacts the hatching rate.Previous works on hatching and larvae dynamics include [3], [2], where stochastic models with food dynamics were used.However, to the best of our knowledge, no mathematical work has been published on the very topic of hatching enhancement through larval density since the experimental findings of [17].Observations on this phenomenon are uneasy to obtain in the field but can be assessed in the lab (see [8]).Further research in this field could benefit from mathematical modeling tools able to take it into account and this may help monitoring the dynamics of mosquito populations.
We develop a mathematical model of the dynamics of mosquito population, with the requirements that this model be sufficiently generic to match experimental observations across various conditions and sufficiently simple so that it is possible to handle it theoretically and interpret it.Therefore we choose to develop a deterministic model based on a system of ordinary differential equations, as was done, for example, in [15].Our simplistic model involves the positive influence of larvae on the system, acting on hatching rate.We show that this feature can explain oscillations.
We draw a general picture of the system's properties in Section 2, and justify rigorously the use of a two-population model as a further simplification for the identification of the qualitative properties induced by hatching feedback.Then we focus on two parameter regimes of particular interest.Firstly (Section 3) when the quantity of eggs is large compared to the quantity of larvae, oscillations can appear and we are faced to a slow-fast oscillatory regime giving rise to oscillation profiles comparable to those of the FitzHugh-Nagumo system (Theorem 3.1).We can compute the amplitude of the oscillations in this case, where they are typically large, and also their period.Secondly (Section 4), we show that our model presents a Hopf bifurcation at any positive equilibrium of the system, assuming the quantity of larvae promotes hatching.The bifurcation occurs as the feedback becomes stronger (Theorem 4.2).In this case we can compute the period of the oscillations at the bifurcation point.We provide numerical results for the system parametrized (roughly) for a tropical area such as Rio de Janeiro, showing that the range of possible oscillations is wide.

Models and their reduction
The life cycle of a mosquito (male and female) consists of two main stages: the aquatic stage (egg, larva, pupa), and the adult stage.We adopt a population biology point of view, which means that we describe the mosquitoes life-cycle thanks to a system of ordinary differential equations.For the purpose of studying the impact of larval density on hatching, we introduce the number densities of each population: A(t) (adults), E(t) (eggs), L(t) (larvae) and P (t) (pupae).
In a compartmental model, one can suppose the following type of dynamics We interpret the parameters as follows: β E > 0 is the intrinsic oviposition rate; δ E , δ L , δ P , δ A > 0 are the death rates for eggs, larvae, pupae and adults, respectively; τ L , τ P > 0 are the transition rates from larvae to pupae and pupae to adults, respectively ; φ tunes an extra-death term due to intra-specific competition (this term is non-linear and we assume that it depends only on the larval density); finally, H(E, L) is the hatching rate, which may in general depend on larval density L and on egg density E, neglecting a possible effect of pupae.
In order to reduce (S 4 ) to a simpler model we suppose pupa population at equilibrium.This boils down to assuming that the time dynamics for pupae is fast compared to the other compartments and thus P = τ L δ P + τ P L.
To justify this approximation more rigorously, we assume τ P , δ P = O(1/ε) (quantifying the "fast dynamics" for pupae) and define τ P = ετ P , δ P = εδ P .We introduce P = εM and then we find the following equations on M and A (those on E and L are untouched) This method follows the classical justification of Michaelis-Menten laws (see [20], [21]).We end up with and in the limit ε → 0, we recover our claim under the form M = τ L τ P +δ P L. This simplification enables us to reduce the model to dimension 3. From now on we also assume H(E, L) = h(L) and φ(L) = cL in order to obtain the simplified system We can proceed to a further reduction by supposing adult population at equilibrium.This boils down to assuming that the time dynamics for adult mosquitoes is fast compared to the other compartments.
Exactly as above with the pupae, in the approximation when δ A and β E are large (and A itself is small), it makes sense to set in this system, at first order, A = τ P τ L (δ P +τ P )δ A L. Finally, system (S 4 ) reduces to the following system in dimension 2: where We perform this model reduction because it is sufficient to take into account the larval effect.Indeed, we show and quantify how the larval density-dependent hatching rate effectively generates oscillations, without any other source of instability (like time-delay, temperature variations or other environmentrelated effects).However, for future practical applications, further studies including the use of a more biologically realistic model will be mandatory.
According to experimental data (results from [8]) and mainly guided by a biological intuition we assume the hatching undergoes saturation for large values of L: To ensure the instability of the trivial equilibrium (0, 0) and rule out population extinction, we assume We also assume, for the matter of simplification of later computations For several mosquito species, it is actually possible to identify the biological parameters τ L , δ A , δ E , δ L and the adult density at equilibrium on the field (i.eA such that d dt A = 0).From the formula L = A δ A τ L , we deduce larvae density at equilibrium on the field (this density is called L throughout this paper).
We warn the reader about what we call "equilibrium density on the field" and about parameter values.We do not claim they can precisely reproduce population variations as observed in field experiments.We simply use rough estimation of their orders of magnitude so as to prove the concept of population oscillations due to density-dependent hatching rate.See paragraph 2.2 for additional comments.
This warning made, from now on we consider that parameters b E , d E , d L , adults and larvae density at equilibrium on the field are known; the competition parameter c and the hatching function h are unknown.The known parameters are set at a given place and temperature (see [23], [25]) and we work with a fixed temperature, so the previous biological parameters are fixed and time-independent.
Our general goal is thus to assert the possible range of remaining parameters c and h(L) depending on the qualitative properties of solutions.
2 Study of the reduced model

Basic properties, equilibria and their stability
With the assumption (2) we know that solutions remain non-negative.Furthermore, the trivial equilibrium (0, 0) is a steady state of (1) and all the other steady states (E, L) are determined by a non-linear We observe that solutions of (5) are positive if and only if In addition: Lemma 2.1.Assume (2) holds.Then there is a constant K > 0 such that for all non-negative t, L(t) + E(t) ≤ K.Moreover, there exists at least one positive steady state of (1) if and only if Furthermore, all steady states (E, L) = (0, 0) For the first point, we do not use any property of h, but merely the fact that cL 2 /L → +∞ as L → +∞.Note that with estimates on h, more restrictive properties can be obtained, in the sense that one could construct strictly smaller positively stable and attractive sets.
Proof.We notice that where Then L defines a steady state of (1) if and only if f (L) = 0, by (5).Continuity of f yields the conclusion since h 0 = max h.
From now on we always assume that (6) holds, so that there exists at least one positive steady state of (1).Then we analyze the stability of those steady states.
A non-trivial steady state (E, L) of (1) is unstable (locally linearly) if and only if either Proof.We divide the proof into three steps.Firstly we linearize system (1) around a steady state (E, L).Setting E = E+e+. . .and L = L+ +. .., we find The eigenvalues λ of the above linear system are given by the determinant After straightforward computations, we obtain: Secondly we look at the trivial steady-state.Taking E = L = 0 in equation ( 9), we obtain: We are looking for the condition such that (0, 0) is linearly unstable (we are interested in the conditions when the mosquito population does not tend to zero in nature).In other words, we expect that the polynomial P has a root with positive real part.Since the first order coefficient is positive we end up with condition (3) and the first point of the lemma is proved.
Finally we consider non-trivial steady states.We rewrite (9) as where A is the Jacobian matrix of the linearized system (1).Using (5) we find and thus The discriminant ∆ of this polynomial is ∆ = tr(A) 2 − 4 det(A) and the steady state is unstable if and only if there exists a root with positive real part.
There are two cases: If ∆ < 0 then the real part of the roots is tr(A) 2 .The steady state is unstable if and only if tr(A) > 0.
. Hence the steady state is unstable if and only if tr(A) > − √ ∆.This is true if and only if either tr(A) > 0 or if det(A) < 0 and tr(A) ≤ 0.

Remark 2.3.
There is a link with the basic offspring number Q 0 (defined in [6]).This dimensionless number is the average number of offspring generated by a single fertilized mosquito: from the method in [22], we can compute We remark that the first statement in Lemma 2.2 boils down to the classical property: trivial equilibrium point is unstable if and only if Q 0 > 1.
Remark 2.4.As in nature we can observe oscillations of eggs and larvae density [13], we pay attention in this work to oscillations around the positive steady states described in Lemma 2.2.We show in Section 4 that these solutions exhibit oscillations, by applying the Hopf bifurcation theorem.This behavior occurs only if the non-trivial steady state is unstable.
For the sake of conciseness we define the following functions: We can rephrase Lemma 2.
Thanks to (4) we can define Proof.We are looking for the k > 0 such that T (k) > D(k), that is also written from ( 12) Recalling that b E > d E + d L by ( 6), the discriminant is: The roots are exactly k ± , so the polynomial is negative when k Collecting our results on the equilibria we can state Proposition 2.6.Assume (6) holds, and let (E, L) be a positive steady state of (1).
Finally, the eigenvalues of the linearized of (1) at (E, L) are complex conjugate and pure imaginary if and only if h(L) > k + and h (L) = T h(L) .
Proof.This is a direct consequence of the previous calculations, except for Inequality ( 14) is equivalent to

This inequality holds because b
Then, setting k = h(L), k = h (L) and using the notations (11), the eigenvalues of the linearized operator are roots of the polynomial Hence the roots are pure imaginary if and only if tr(A) = 0 and det(A) > 0. From the definition of T, D in (12), tr(A) = 0 if and only if k = T (k).As det(A) > 0 if and only if k < D(k), by Lemma 2.5 this holds whenever k > k + .

Discussion on the nonlinearities and the equilibrium values
We discuss in this paragraph the nonlinearities of system (1), and the role they play.
First we justify the use of a competition term.Solutions of (1) are bounded (Lemma 2.1), but this holds only thanks to the nonlinear competition term −cL 2 in the equation describing the larvae dynamics.More generally, any competition term φ(L), as in Section 1 such that φ(L) → +∞ as L → +∞ yields the same result.However, in the absence of such a competition, a priori bound on the solutions cannot be obtained, and no phenomenon keeps the population finite.For Aedes mosquitoes, the amount of available food in the breeding sites is an actual resource limitation that can trigger massive death of larvae if the amount of food per larva drops down too low (see [2]).Therefore, we choose the simplest (i.e.quadratic) competition term to represent this competition for resources, and this ensures mathematically that solutions remain bounded.
Still, the competition parameter c is extremely hard to assess from experimental data, and the values we use in this work should be handled with care.Usually, we fix a value for a positive equilibrium L (which corresponds to choosing a type of breeding site).Then, to each value k = h(L) corresponds a non-necessarily unique c(k) that makes L an equilibrium of (1).We treat k as a free parameter in this study.It has been observed that the hatching rate indeed is extremely dispersed (see for instance the experimental results of [17]), depending not only on the mosquito population and the environmental conditions but also on the egg batches themselves.In future works expanding on the simplest oscillatory behavior we describe here, this variability in the actual value of k should be taken into account if the model outputs are to be linked with experimental data.
Second, we discuss the hatching rate function h, which is crucial to our study.From now on, we require h to be increasing.Indeed, Proposition 2.6 shows that a steady state is always stable if h is decreasing.Hence only an increasing h can produce stable oscillations.This mathematical assumption is supported by a simple biological hypothesis: larvae promote hatching.
An interesting feature of this intuition is that it can be subsequently extended to higher-dimensional systems such as (S 4 ).In other words, it is not an artifact produced by considering only a 2-dimensional system but a robust qualitative property for these systems.Indeed, for (S 4 ) the Jacobian matrix at any point X = (E, L, P, A) reads hence if h (L) < 0 then J(X) is a Metzler matrix (it has positive extra-diagonal coefficients): the system is cooperative in this case.Its characteristic polynomial may be written where A i , C > 0. Being a Metzler matrix, J has a real dominant eigenvalue.This matrix is stable if and only if this eigenvalue is negative; in other words, if and only if P (0) > 0 (since P is increasing on (0, +∞)).This condition reads At equilibrium, , therefore P (0) > 0 and thus any equilibrium where h < 0 must be (locally) stable, in system (S 4 ) as well as in system (1).Adding "neutral" compartments keeps this property true and we can be confident in concluding that only a positive effect of larvae on hatching rate can destabilize the equilibrium and lead to (local) oscillations.Some preliminary experiments ran by one of the authors seem to indicate that the larval impact on hatching may depend on larval development stage.Taking this into account would require model complexification.For instance, to model hatching impact discrepancies between first instar (positive) and last instar larvae (negative) we could add at least one compartment in (1).However, we focus here on the simplest oscillations-producing mechanism.The hatching function being increasing and bounded, it is reasonable to assume that h is S-shaped and smooth, which is what we use in the rest of the paper.
Third, having discussed the two nonlinearities in (1), we are left with an important question about steady states: how to ensure that L is actually unique?The second equation in ( 5) is also written The number of positive steady states depends strongly on function h.Being a S-shaped function does not guarantee uniqueness.Therefore, it should be checked case by case except for some simple function families.We illustrate this fact in Appendix A with Hill functions.Still, we notice that uniqueness is guaranteed if (3) holds and either, for all L ∈ (0, (b

The slow-fast oscillatory regime
In order to understand periodic solutions to (1), we examine a possible regime with a small parameter and then prove the oscillation result (Theorem 3.1).We have in mind here the analysis of the FitzHugh-Nagumo system.Numerical illustration, amplitude and period computation in some particular cases can be found in Appendix B.

Parameter regime and main result
Here, we assume that the egg stock is large, and its dynamics slow compared with the larvae stock.This identifies a small parameter leading to a slow-fast system.
More precisely, let ε > 0, η : R + → R + , and assume at first that all parameters (except for h) may depend on ε.We transform the variables (E, L) from (1) into v ε := εE and We assume that parameters scale in such a way that the following limits exist, as ε → 0: In addition, we assume that the zero set of g is "non-degenerate" in the sense: We give below a simple proof of the following fact, in the spirit of Tikhonov's theorem on dynamical systems [10].
In addition, solutions of system (16) along with any bounded initial data admits a limit as ε → 0: for almost every t > 0 and the trajectory is uniquely defined from f and g with dv dt = f (u, v).
Figure 1 illustrates the slow-fast dynamics.Before proving Theorem 3.1, we justify the particular scaling choices in its statement.Non-trivial equilibrium (E, L) of (1) are given by (5):  Thus in all generality (allowing all parameters to depend on ε), the scalings fit for our purpose (i.e. with E(ε) = 1/ε) are exactly those for which and there exists η(ε) = O(1) such that Therefore the scaling choice made in Theorem 3.1 is in some sense "generic".Note that for every possible parameter scaling we get a (possibly different) limit in (17).For instance, assuming εc ε and εb E (ε) have limits 1/η 0 , ξ > 0 respectively as ε → 0 (this is the case with the scaling used in Theorem 3.1), we choose η(ε) = O(1) such that c ε η(ε)ε = 1 and end up with The limits f and g are given by

Proof of the main result
We proceed to the proof of Theorem 3.1 in three steps.First, scaled quantities u ε and v ε remain uniformly bounded independently of ε, as can be proved from direct computation using the bound K from Lemma 2.1.
Lemma 3.2.There exists C > 0 such that for all ε > 0 and t > 0, Hence, up to extraction, v ε converges to v uniformly on compact sets [0, T ] by the Ascoli theorem.Then, the convergence of an auxiliary quantity gives convergence of u ε : Moreover, there exists u, v ∈ L 1 ∩ L ∞ such that after extraction of a subsequence u ε → u in L p (0, T ) for all 1 ≤ p < ∞, as v ε → v uniformly.
Proof.Let B(t, u) := u 0 g 2 (σ, v(t))dσ, where v is the limit of v ε (obtained by the Ascoli theorem) and g is the limit of g ε (from ( 17)).From (18) we deduce that for all t, u → B(t, u) is increasing.Hence there exists a smooth function A(t, u) such that for all t, u, A(t, B(t, u)) = u.
If there exists w(t) ∈ L p (0, T ) for all p < ∞ and T > 0 such that then defining u(t) := A(t, w(t)) we can conclude that Indeed, we notice that Since u ε is uniformly bounded, for some C > 0 which depends only on ∂ v g.Hence (21) implies Therefore we only need to prove (21) to complete the proof.To do so we first obtain (20) by computing Hence This gives (20).Then we introduce We compute By the previous point, t As a consequence, w ε is uniformly (in ε) bounded in BV loc .This implies that up to extraction, w ε → w in L 1 .Because w ε is also bounded in L ∞ , convergence actually takes place in all L p spaces.

Lemma 3.4. With these assumptions we have:
There exists a unique τ > 0 and a (unique up to translations) τ -periodic function (u τ , v τ ) : R + → Υ such that v τ is Lipschitz-continuous, u τ is piecewise continuous, for all t ≥ 0, v τ = φ(u τ ) everywhere, vτ = f (u τ , v τ ) almost everywhere and the discontinuities of u τ are located at times t such that φ has a local extremum at u τ (t − ).
Clearly from (19), Lemma 3.4 applies with the hypotheses of Theorem 3.1 and thus proving the remaining part of the theorem.
Proof of Lemma 3.4.Thanks to assumptions (R.1), (L.1), (L.3) and (L.4), the construction of (u τ , v τ ) is classical and can be done by pasting together solutions of Cauchy problems given (locally) by vτ = f (φ −1 (v τ ), v τ ), on intervals where φ is invertible.Uniqueness comes from the crucial fact that discontinuities of u τ are assumed to be located at local extrema of φ.
From the previous lemmas we know that (u, v) ∈ Υ almost everywhere.In addition, uniform boundedness of f ε (u ε , v ε ) ensures that v is Lipschitz continuous.
Then, we claim that if t > 0 is such that φ has no local extremum at u(t + ), then there exists τ 0 > 0 such that (u, v) is continuous on (t, t + τ 0 ).This point is the key of the proof.To prove it, let u i be such that φ (u i ) < 0. We solve only the simpler problem , where the inverse of φ is taken locally (this is possible for ε small enough since φ (u i ) < 0 and v ε is uniformly Lipschitz-continuous), we obtain for some r ε (t) between u ε (t) and φ −1 ( v ε (t)).We have ∂ 1 g ≤ −α < 0 on a neighborhood of (u i , φ(u i )), so on this neighborhood w ε remains small (it is a o(ε)), which in turn proves that ( r ε , v ε ) remains in this neighborhood.In particular, u ε converges to some function u which is continuous at t = 0 (since it is equal to φ −1 ( v(t)) on a positive neighborhood of 0).We do not write the full proof because the derivation we use here extends readily at the price of tedious notations.A full proof should use f ε , g ε rather than f, g, and rise some analogue φ ε of φ at level ε > 0, for ε small enough, which is locally invertible on a neighborhood of the initial data.It does not require more assumptions than the ones we stated.This is enough to get all the results of Lemma 3.4, except for the initial layer which we treat now.To fix the notations, we assume that φ has a local minimum equal to φ m at u m and a local maximum equal to φ M > φ m at u M < u m .Moreover, let u 0 m < u M such that φ(u 0 m ) = φ(u m ).For α, β ∈ {1, −1}, we also introduce Z β α := {sgn(f ) = α, sgn(g) = β}.We define a mapping π : R 2 → Υ by π = Id on Υ and if (u, v) ∈ Z β α then π(u, v) = (u 1 , v) such that φ(u 1 ) = v and sgn(u 1 − u) = β.The projection π is well-defined thanks to the assumptions on φ and g, except on (u 0 m , +∞) × {φ m }, on which we let π ≡ (u 0 m , φ m ).Then (u, v)(0 + ) = π(u 0 , v 0 ).To prove this, one simply has to check the behavior of u ε (since v ε and v are Lipschitz continuous).As above, we claim that the first-order behavior is simply given by the "layer equation" which makes u ε converge exponentially fast to π(u 0 , v 0 ) 1 , thanks to assumptions (R.2) and (L.2).Up to tedious notations and thanks to ( 17) and (R.2), this result extends to u ε 0 , v ε and g ε .
) After the initial layer, the trajectory of (u, v) remains on Υ s .This follows from the sign of f on Υ u : because of the continuity property, the trajctory cannot exit Υ s but at (u m , φ m ) (or (u M , φ M ), respectively).At these points however, Υ u is repulsive since v must be continuous, v < 0 ( v > 0, respectively) and Υ u lies locally in {v > φ m } (respectively in {v < φ M }).
Still, the initial data does not need to be projected directly by π on It remains to check that τ 1 < +∞.For all T > 0, as long as φ has no local extremum at u(t) for t ∈ (0, T ), u is continuous.Thanks to our assumption (R.1), there are two connected components in Υ s , on each one of whom sgn(f ) is constant.Because of assumption (L.4), f must be negative on the unbounded connected component.Therefore (u, v) remains on (0, T ) in a part of Υ where |f | is positively bounded from below (one of the two connected components of Υ s ) and has the appropriate sign.This yields the existence of τ 1 < +∞.
Then for all t ≥ τ 1 we have v(t) ∈ [φ m , φ M ], and the trajectory is uniquely defined onwards.
Remark 3.5.We did not treat the case when the limit of (u ε 0 , v ε 0 ) belongs to Υ (relaxing assumption (R.2)).In this case indeed, no general result can be obtained, unless the various convergence speeds (of f ε , g ε , u ε 0 and v ε 0 ) are quantified.
Remark 3.6.The last point of Theorem 3.1 implies that the amplitude of the oscillations (in u, v) at the limit ε → 0 can be computed if one knows these parameter scale in ε thanks to only f and g.Their period τ can also be computed directly from f and φ.As in the proof of Lemma 3.4 we denote the intervals of values taken by u(t) where it is continuous (and thus C ∞ ) as [u 0 m , u M ] and [u m , u 0 M ] respectively, and let du .
Then we have

Hopf bifurcation
Numerical observations (see Section A and Appendix A) show that the system (1) has a stable periodic solution oscillating around the non-zero steady state, even far from the slow-fast asymptotics we explored in the previous section.We now prove the local existence of this periodic solution using the Hopf bifurcation theorem (Theorem 8.8 from [19], with a classical proof in [18]; see also [10]) for 2 × 2 systems of differential equations.

The function class H L
To find out a possible bifurcation parameter, we choose the hatching function h within a special class, for which we fix the value of one specific steady state L. With this setting, we can state a bifurcation theorem using the simple bifurcation parameter h (L), which represents the sensitivity of hatching rate to larval density at equilibrium.However, it is worth noting that our argument does not rely on the structure of this class of functions, and may be adapted, for instance, to the Hill functions considered in Appendix A.
For a fixed L the class of functions under consideration that fits our purposes is Graphs of these functions are shown in Figure 2. We use the immediate properties that these functions  Also we can solve the equation in L, aπ . Hence L is positive under the stated condition.
Remark 4.1.From Lemma 2.2, for h of class H L , the state (0, 0) is unstable if and only if .
Remark 4.3.k + is given by (13), and ã is such that the normal form coefficient α N (see [19]) of our system is negative if a > ã.We simply give a numerical justification of the existence of ã as the computations appear to be very long (see the proof below).
Remark 4.4.The value of a must be greater than a crit to ensure that the linearized operator has complex eigenvalues.
The bifurcation diagram for S γ (a 0 , b 0 ) in Figure 3 is obtained by XPPAUT software [9].Proof of Theorem 4.2.We set λ 1,2 (γ) = α(γ) ± iβ(γ) (with γ a real parameter), the two eigenvalues of J P (γ) the Jacobian matrix associated to our system and computed in (0,0).We call γ c a bifurcation value, and α N (γ) the normal form coefficient of the system (see [19]).Firstly, we only need to study complex conjugate and pure imaginary eigenvalues of J P (γ) to find the bifurcation value γ c , which means also to look for γ c such that α(γ c ) = 0 and β(γ c ) = 0. Thanks to Proposition 2.6 we know that this is the case when k > k + i.e. aπ 2 > k + or equivalently a > a crit (by definition, a crit = 2k + π ).Moreover since h (L) = ab (direct computation from ( 25)), we know that the bifurcation value is located at the level of the graph G of function b (defined in (26)) And we can set γ c = 0. Secondly, we have to see if dα dγ (γ c ) > 0, this means to check that tr(J P (γ) ) changes sign at the bifurcation value γ c .Let γ −→ z(γ) = α(a 0 , b(a 0 ) + γ).We recall that α(γ) is a function of a and b. Since and we obtain that z (γ c ) = z (0) = aE 2 and it is always positive.
Thirdly, we have to study the normal form coefficient of the system computed in γ c = 0 and find when α N (γ c ) = 0. To get the normal form coefficient, we have to transform the system (S γ (a 0 , b 0 )) and we use the steps from [19].In a first step we reduce the initial system (S γ (a 0 , b 0 )) to a system where the equilibrium (E, L) becomes the origin.By the change of variables x = E − E and y = L − L, (S γ (a 0 , b 0 )) becomes: Then as (E, L) is an equilibrium, we can simplify (28) into which we write as where The system (30) can also be written under the matrix form We call M the first (2 × 2) matrix in the right-hand-side.Now, to obtain the normal form coefficient, one way is to perform a linear change of variables so as to get In our case, we can have an idea of the normal coefficient only in a neighborhood of γ = 0.Because we want to make a simple linear change of variables, we are looking for a matrix P such that P M P −1 = N and that at the bifurcation value γ = 0, tr(M ) = 0 = tr(N ) and det Next we obtain the matrix system (31) where In a final step we compute the normal form coefficient using the previous formulas and the expression that exists in two dimensions given in [19] which is: The coefficient is easy but very tedious to compute, and we used the computer algebra system Maple [1] to get its expression.
In our case the coefficient is equal to zero for some value ã > 0, and is always negative for a > ã (as it appears that ã < a crit , this is sufficient by definition of ( 27)).Then α N (γ c ) = 0 for a = ã.
Finally, we want to have for all real γ in a neighborhood of 0, α N (γ)α(γ) < 0. Thanks to Maple we have α N (0) < 0, in a neighborhood of γ = 0, for a > ã with ã small.So we can apply the Hopf bifurcation theorem that ensures there exists a limit cycle (periodic solution) when α(γ) > 0 (i.e tr(J P (γ)) > 0), and moreover this cycle is stable as α(γ) > 0: we are faced to a supercritical bifurcation.

Discussion on the period of the oscillations
The period of the oscillating solutions are relevant to the biological problem in consideration, because they can be compared with observations in nature.Proposition 4.5.As γ → 0 + , the periodic solution of the system (S γ (a, b(a))) has a frequency ω and a period T 0 = 2π/ω given by the expression Proof.As γ → 0 + , the oscillations frequency is given by the imaginary part of the root of the polynomial equation (9) in the case of non-trivial steady state.The frequency is ω γ = det(J P (γ) ), where the expression of det(J P (γ) ) is Then the expression of ω = ω 0 follows.
Remark 4.6.At the bifurcation value, the parameter k can be linked with the period T 0 .Let T 0 a given period observed experimentally, then we find a corresponding value for k as the positive root of the following characteristic polynomial: Away from the bifurcation value, the real part of the eigenvalues is greater than zero and the period of the oscillations can only be obtained numerically.Unfortunately, this case is more relevant as the Hopf bifurcation theorem asserts that the amplitude is increasing with the parameter γ.In other words, for fixed a the amplitude of the oscillations is an increasing function of b.

Conclusion
We show that introducing internal regulation in the form of a larval-density-mediated hatching rate in a compartmental model for mosquito population dynamics induces stable oscillations.These oscillations can be rather simply understood from the mathematical point of view either as cycles produced by a Hopf bifurcation (Theorem 4.2), in a first parameter regime, or as the typical slow-fast behavior (close to FitzHugh-Nagumo model, Theorem 3.1) in a second parameter regime.
Our study supports the idea that understanding internal life-cycle regulation can effectively help modeling and simulating population dynamics properly.Ongoing experiments of some of the authors try to reproduce the larval density impact on hatching which was observed in [8] and may shed some light on this misunderstood phenomenon.In particular, restricting the parameters and possible oscillations range could only be reached by assessing as precisely as possible the actual hatching feedback.
In this paper we neglect environmental variations.Therefore it leaves open for future studies the deep question of linking internal life-cycle regulation and external variations (induced, for instance, by rainfall and temperature) in order to get a better description of the mosquito populations dynamics.However, it was observed that population oscillations may happen on periods much shorter than seasonal variations, and this justifies the study of internal regulations as possible triggers.
Another possible extension of our works is the adaptive dynamics of hatching regulation trait.Indeed, synchronizing the egg hatching may be beneficial for a population in a given environment, but also be detrimental if rare and extreme events can annihilate larval population, for instance.The egg stage can be seen indeed as a quiescent, refuge state for the species (this approach was studied in [24]).Here we prove that positive feedback of larvae on egg hatching tends to make the population size oscillate, creating distinct generations (synchronizing effect) while negative feedback tends to stabilize the population size, which may be detrimental on the long run if, for example, the favorable period for larvae and adult development is typically short.The above computations extend to the 3-dimensional system (S 3 ), and numerical observations are similar.Indeed, the condition (35) guaranteeing positivity of the trace of the Jacobian at the unique positive equilibrium, rewrites for system (S 3 ) as (p − 2) In this case we define and the above condition is equivalent to X(k) < 1.
Exactly as in the two-dimensional case, we explore the full range of (34) by choosing the parameters (ι, ζ) ∈ (0, 1 − X(k)) × (0, 1) and defining α ι (k) and a ι,ζ (k) by the same formulas as before.For all the numerical values we took for ι and ζ, we always found oscillating solutions.Examples (dynamics of larvae and of (E, L, A) in the three dimensional space) are shown in Figure 5. B Amplitude and period computation in the slow-fast regime In the slow-fast approach, system (1) exhibits oscillations with known amplitude and period at the limit ε → 0. We show here how to compute this amplitude analytically.To do so, we simply compute the local extrema of u → ηu 2 h(ηu) .The first-order necessary condition is xh (x) = 2h(x), where x = ηu.This provides with a general method to determine the limit trajectories.With the previous example from (32), h(x) = h m + a x p αL p + x p , this boils down to Letting y = x p , we end up with a second-order polynomial, for which the analytical computation can be pushed a few steps further.In particular, its discriminant is Hence there are exactly two positive local extrema if and only if (2−p) 2 a > 8ph m and (2−p)a+4h m < 0.
The first condition implies the second one if p > 2, and the second one is impossible if p ≤ 2. Therefore the only case when there are two local extrema is when Under assumption (36) we find that the extrema (y M < y m ) are located at With the notations of Lemma 3.4, Then we can compute u 0 r for r ∈ {m, M } by solving η•(u 0 r ) 2 h(ηu 0 r ) = φ r .Unfortunately this cannot be done analytically.However, the amplitude of the oscillations in terms of v is equal to With E = 1/ε, we expect that the oscillations of E have amplitude 1+α p , by (22).Hence the amplitude of egg oscillations is equal to We can simplify this expression one step further by letting ρ := h m /a.Then we notice that q In particular we notice that the amplitude depends only on the function h through ρ, α (hence L) and p, and not on any other biological parameter, under the constraints (36).
An interesting case is when p → +∞, where h approaches a step function from h m to h m + a, with its jump located at αL.In this limit we can compute the amplitudes in u and v: If we assume d E = 0 (for simplicity), using formula (23), we can also obtain in this case an analytical expression for the period of the oscillations: C Numerical oscillations, period and amplitude close to the bifurcation We illustrate the statements from Section 4 with numerical examples.Biological parameters of (1) are taken at a temperature around 25 • C which leads to A = 3.4 mosquitoes per 100 square meters (taken from a physical situation described in [23]) and b E = 20.94,d L = 0.15.(taken from [25]).To fit the condition d E d L , d E is fixed arbitrarily at 1 180 .We note that condition ( 4  The parameters a and b are chosen so that Theorem 4.2 applies, which proves the existence of periodic solutions close to the non-trivial steady states.We perform numerical test by letting a parameter j vary in a set J of 18 values between 0.05 and 4 in order to obtain 162 couples (a i , b j i ) i=1,... where b i,min is the minimal b that can be chosen for a i to obtain oscillations (if b < b i,min the solutions can not oscillate), i.e. for which the trace of the linearized operator is equal to 0.
The hatching functions are: In We provide numerical results for i ∈ {1, 4, 9} and j ∈ {0.1, 0.25, 0.5} initial data close to the steady state (which is drawn in dashed line).Two sets of initial data are chosen, (E(0), L(0)) = (E, L) (green) and (E(0), L(0)) = (E, L + 0.02) (blue), which gives oscillations that appear to be periodic in time.Table 3: Steady states, period and amplitude of oscillations for a = .5 ) and a time variable evaluated in [0, 150] days.Considering the blue curves, we sum up in the Tables 1, 2 and 3 what we obtain for the period and the oscillations' amplitude taken by the solutions.In the last line of the tables we give a value of b that can be chosen to obtain a period of about 10 days.Relative amplitude of the oscillations is expressed as a percentage of the (constant) value L.
It is possible to achieve the same period T 0 for different couples of parameters (a, b).For a fixed a, when b is increasing, the period T 0 and the amplitude of larvae are increasing too.The amplitude, on the contrary, mainly depends on b.This is illustrated in Figure 9.

2 π
. Indeed, for given values (k, k ) ∈ R 2 + , the choice of a = 2k π and b = k a gives the solution since h(L)

Figure 3 :
Figure 3: Supercritical Hopf bifurcation diagram with a 0 = 0.2.The bifurcation parameter b is in x-axis, the diagram shows extreme values of the periodic solution for L (the L scale is in y-axis).The steady state is stable (red line) until the bifurcation point (point number 2) is reached.A periodic solution appears and is stable (green points) until a bigger value of b, where it becomes unstable (blue circles).The amplitude of the periodic solution grows with the parameter b.

Figure 9 :Figure 8 (
Figure 9: Larvae dynamics period T 0 in days (left) and larvae dynamics amplitude (Amp) in percentage of L (right), for different couples (a, b).

Table 1 :
our tests the steady state changes with i (for example E 1 = 145.92,E 4 = 59.59 and E 9 = 30) but we always have L = A δ A τ L = 1.13.Steady states, period and amplitude of oscillations for a = .1

Table 2 :
Steady states, period and amplitude of oscillations for a = .25