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

Martin Strugarek, Laetitia Dufour, Nicolas Vauchelet, Luís Almeida, Benoît Perthame, Daniel Antunes Maciel Villela

► To cite this version:

Martin Strugarek, Laetitia Dufour, Nicolas Vauchelet, Luís Almeida, Benoît Perthame, et al.. Oscillatory regimes in a mosquito population model with larval feedback on egg hatching. Journal of Biological Dynamics, 2019, 13 (1), pp.269-300. 10.1080/17513758.2019.1593524. hal-01674280

HAL Id: hal-01674280

https://hal.science/hal-01674280

Submitted on 10 Jan 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Oscillatory regimes in a mosquito population model with larval
feedback on egg hatching

Martin Strugarek¹,², Laetitia Dufour², Nicolas Vauchelet³, Luis Almeida², Benoît Perthame²,
and Daniel A.M. Villela⁴

¹AgroParisTech, 16 rue Claude Bernard, F-75231 Paris Cedex 05, France

²Université Pierre et Marie Curie, Paris 6, Laboratoire Jacques-Louis Lions, UMR 7598 CNRS, 4 Place Jussieu, 75252 Paris, France

³LAGA - UMR 7539 Institut Galilée Université Paris 13 99, avenue Jean-Baptiste Clément 93430 Villetaneuse - France

⁴PROCC, Fundação Oswaldo Cruz, Av. Brasil, 4365. Manguinhos CEP 21040-360, Rio de Janeiro, Brazil

January 10, 2018

Abstract

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.

**Key-words:** Mathematical biology; Hopf bifurcation; Slow-fast dynamics; Egg hatching; Mosquitoes life cycle;

MSC numbers: 34C23; 34E15; 92D25;

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 aegypti* 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.

# 1 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

$$ egin{cases} 
\frac{d}{dt}E = \beta_E A - E(H(E, L) + \delta_E), \\
\frac{d}{dt}L = EH(E, L) - L(\phi(L) + \delta_L + \tau_L), \\
\frac{d}{dt}P = \tau_L L - \delta_P P - \tau_P P, \\
\frac{d}{dt}A = \tau_P A - \delta_A A.
\end{cases} 
\quad (S_4) $$

We interpret the parameters as follows: $\beta_E > 0$ is the intrinsic oviposition rate; $\delta_E, \delta_L, \delta_P, \delta_A > 0$ are the death rates for eggs, larvae, pupae and adults, respectively; $\tau_L, \tau_P > 0$ are the transition rates from larvae to pupae and pupae to adults, respectively; $\phi$ 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 (S4) 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 = \frac{\tau_L}{\delta_P + \tau_P}L$.

To justify this approximation more rigorously, we assume $\tau_P, \delta_P = O(1/\varepsilon)$ (quantifying the “fast dynamics” for pupae) and define $\bar{\tau}_P = \varepsilon\tau_P, \bar{\delta}_P = \varepsilon\delta_P$. We introduce $P = \varepsilon M$ and then we find the following equations on $M$ and $A$ (those on $E$ and $L$ are untouched)

$$ egin{cases} 
\varepsilon \frac{dM}{dt} = \tau_L L - \varepsilon M(\tau_P + \delta_P), \\
\frac{dA}{dt} = \varepsilon \tau_P M - \delta_A A.
\end{cases} $$

This method follows the classical justification of Michaelis-Menten laws (see [20], [21]). We end up with

$$ egin{cases} 
\varepsilon \frac{dM}{dt} = \tau_L L - M(\bar{\tau}_P + \bar{\delta}_P), \\
\frac{dA}{dt} = \bar{\tau}_P M - \delta_A A,
\end{cases} $$

and in the limit $\varepsilon \to 0$, we recover our claim under the form $M = \frac{\tau_L}{\bar{\tau}_P + \bar{\delta}_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 $\phi(L) = cL$ in order to obtain the simplified system

$$ egin{cases} 
\frac{d}{dt}E = \beta_E A - \delta_E E - h(L)E, \\
\frac{d}{dt}L = h(L)E - \delta_L L - cL^2 - \tau_L L, \\
\frac{d}{dt}A = \frac{\tau_P \tau_L}{\delta_P + \tau_P} L - \delta_A A.
\end{cases} 
\quad (S_3) $$

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 $\delta_A$ and $\beta_E$ are large (and A itself is small), it makes sense to set in this system, at first order, $A = \frac{\tau_P \tau_L}{(\delta_P + \tau_P) \delta_A} L$.

Finally, system (S4) reduces to the following system in dimension 2:

$$ \left\{ \begin{aligned} \frac{d}{dt} E &= b_E L - d_E E - h(L) E, \\ \frac{d}{dt} L &= h(L) E - d_L L - c L^2, \end{aligned} \right. \qquad (1) $$

where $b_E = \beta_E \frac{\tau_P \tau_L}{(\delta_P + \tau_P) \delta_A} > 0$, $d_E = \delta_E > 0$ and $d_L = \delta_L + \tau_L > 0$.

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 environment-related 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:

$$ h \in C^1([0, \infty)), \quad h > 0, \quad \max_L h(L) := h_0 < +\infty. \qquad (2) $$

To ensure the instability of the trivial equilibrium (0,0) and rule out population extinction, we assume

$$ d_E d_L < h(0) (b_E - d_L). \qquad (3) $$

We also assume, for the matter of simplification of later computations

$$ b_E > d_L + d_E. \qquad (4) $$

For several mosquito species, it is actually possible to identify the biological parameters $\tau_L$, $\delta_A$, $\delta_E$, $\delta_L$ and the adult density at equilibrium on the field (i.e A such that $\frac{d}{dt} A = 0$). From the formula $L = A \frac{\delta_A}{\tau_L}$, we deduce larvae density at equilibrium on the field (this density is called $\bar{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

## 2.1 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 ($\bar{E}$, $\bar{L}$) are determined by a non-linear
relation on $\overline{L}$

$$ \left\{ \begin{aligned} \overline{E} &= \frac{b_E \overline{L}}{d_E + h(\overline{L})}, \\ c\overline{L} &= b_E - d_L - \frac{d_E b_E}{d_E + h(\overline{L})}. \end{aligned} \right. \qquad (5) $$

We observe that solutions of (5) are positive if and only if $h(\overline{L}) > \frac{d_E d_L}{b_E - d_L}$. 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) \le K$. Moreover, there exists at least one positive steady state of (1) if and only if

$$ \min_{x \ge 0} \left( cx + \frac{d_E b_E}{d_E + h(x)} \right) \le b_E - d_L. \qquad (6) $$

Furthermore, all steady states $(\overline{E}, \overline{L}) \ne (0,0)$ satisfy $0 < c\overline{L} < b_E - d_L - \frac{d_E b_E}{d_E + h_0}$.

For the first point, we do not use any property of $h$, but merely the fact that $cL^2/L \to +\infty$ as $L \to +\infty$. 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

$$ \frac{d}{dt}(E+L) = b_E L - d_E E - d_L L - cL^2 \le -d_E(E+L) + U_M, $$

where $U_M := \frac{(b_E+d_E-d_L)^2}{4c}$ is the maximum of $L \mapsto (b_E + d_E - d_L)L - cL^2$. Consequently the claim holds with $K = U_M/d_E$.

Let

$$ f(x) = cx + \frac{d_E b_E}{d_E + h(x)} - (b_E - d_L). $$

Then $\overline{L}$ defines a steady state of (1) if and only if $f(\overline{L}) = 0$, by (5).

Continuity of $f$ yields the conclusion since $h_0 = \max h$.

<figure><img src="image_2.png" /></figure>

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.

**Lemma 2.2.** The steady state (0,0) is unstable (locally linearly) if and only if (3) holds.

A non-trivial steady state $(\overline{E}, \overline{L})$ of (1) is unstable (locally linearly) if and only if either

$$ h'(\overline{L})\overline{E} - d_L - 2c\overline{L} - d_E - h(\overline{L}) > 0, \qquad (7) $$

or

$$ cd_E \overline{L} - d_E h'(\overline{L}) \overline{E} + c \overline{L} h(\overline{L}) < 0 \quad \text{and} \quad h'(\overline{L}) \overline{E} - d_L - 2c \overline{L} - d_E - h(\overline{L}) \le 0. \qquad (8) $$

*Proof.* We divide the proof into three steps.

Firstly we linearize system (1) around a steady state $(\overline{E}, \overline{L})$. Setting $E = \overline{E}+e+...$ and $L = \overline{L}+\ell+...$, we find

$$ \left\{ \begin{aligned} \frac{d}{dt}e &= b_E \ell - d_E e - h(\overline{L})e - h'(\overline{L})\overline{E}\ell, \\ \frac{d}{dt}\ell &= h(\overline{L})e + h'(\overline{L})\overline{E}\ell - d_L\ell - 2c\overline{L}\ell. \end{aligned} \right. $$
The eigenvalues $\lambda$ of the above linear system are given by the determinant

$$ \begin{vmatrix} -d_E - h(\bar{L}) - \lambda & b_E - h'(\bar{L})\bar{E} \\ h(\bar{L}) & h'(\bar{L})\bar{E} - d_L - 2c\bar{L} - \lambda \end{vmatrix} = 0. $$

After straightforward computations, we obtain:

$$ \lambda^2 - \lambda(h'(\bar{L})\bar{E} - d_L - 2c\bar{L} - d_E - h(\bar{L})) + d_E(d_L + 2c\bar{L} - h'(\bar{L})\bar{E}) + h(\bar{L})(d_L + 2c\bar{L} - b_E) = 0. \quad (9) $$

Secondly we look at the trivial steady-state. Taking $\bar{E} = \bar{L} = 0$ in equation (9), we obtain:

$$ P(\lambda) := \lambda^2 + \lambda(d_L + d_E + h(0)) + d_E d_L + h(0)(d_L - b_E) = 0. \quad (10) $$

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

$$ \lambda^2 - \operatorname{tr}(A)\lambda + \det(A) = 0, $$

where $A$ is the Jacobian matrix of the linearized system (1). Using (5) we find

$$ d_E(d_L + 2c\bar{L} - h'(\bar{L})\bar{E}) + h(\bar{L})(d_L + 2c\bar{L} - b_E) = cd_E\bar{L} - d_Eh'(\bar{L})\bar{E} + c\bar{L}h(\bar{L}), $$

and thus

$$ \left\{ \begin{array}{l} \operatorname{tr}(A) = h'(\bar{L})\bar{E} - d_L - 2c\bar{L} - d_E - h(\bar{L}), \\ \det(A) = cd_E\bar{L} - d_Eh'(\bar{L})\bar{E} + c\bar{L}h(\bar{L}). \end{array} \right. \quad (11) $$

The discriminant $\Delta$ of this polynomial is $\Delta = (\operatorname{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 $\Delta < 0$ then the real part of the roots is $\frac{\operatorname{tr}(A)}{2}$. The steady state is unstable if and only if $\operatorname{tr}(A) > 0$.

If $\Delta \ge 0$ then the bigger root is $\frac{\operatorname{tr}(A)+\sqrt{\Delta}}{2}$. Hence the steady state is unstable if and only if $\operatorname{tr}(A) > -\sqrt{\Delta}$. This is true if and only if either $\operatorname{tr}(A) > 0$ or if $\det(A) < 0$ and $\operatorname{tr}(A) \le 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 $Q_0 = \frac{b_E h(0)}{d_L(d_E + h(0))}$.

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:

$$ \left\{ \begin{aligned} T(k) &= \frac{1}{\bar{L}}(2k + \frac{k+d_E}{b_E}(k+d_E-d_L)), \\ D(k) &= \frac{1}{\bar{L}}\frac{k+d_E}{b_E d_E}(k(b_E-d_L)-d_E d_L). \end{aligned} \right. \qquad (12) $$

We can rephrase Lemma 2.2 into: Let $(k, k') = (h(\bar{L}), h'(\bar{L}))$ at some equilibrium $\bar{L}$. The state $(\bar{E}, \bar{L})$ is unstable if and only if either $k' > T(k)$ or $T(k) \ge k' > D(k)$. Thanks to (4) we can define

$$ k_{\pm} := \frac{d_E(b_E + 2d_E + d_L) \pm \sqrt{4d_E^3(b_E - d_E - d_L) + d_E^2(b_E + 2d_E + d_L)^2}}{2(b_E - d_E - d_L)}. \qquad (13) $$

**Lemma 2.5.** Assume (6) holds. If $k > k_+$ then $T(k) < D(k)$, and if $k \in (0, k_+)$ then $T(k) > D(k)$.

*Proof.* We are looking for the $k > 0$ such that $T(k) > D(k)$, that is also written from (12)

$$ k^2(b_E - d_E - d_L) - kd_E(b_E + 2d_E + d_L) - d_E^3 < 0. $$

Recalling that $b_E > d_E + d_L$ by (6), the discriminant is:

$$ \Delta = d_E^2(b_E + 2d_E + d_L)^2 + 4d_E^3(b_E - d_E - d_L) > 0. $$

The roots are exactly $k_{\pm}$, so the polynomial is negative when $k \in (k_{-}, k_{+})$.

We note that $k_{-} < 0$, so $T(k) < D(k)$ if and only if $k > k_{+}$, and $T(k) > D(k)$ if and only if $k \in (k_{-}, k_{+})$. Since $k > 0$, this is equivalent to $k \in (0, k_{+})$. □

Collecting our results on the equilibria we can state

**Proposition 2.6.** Assume (6) holds, and let $(\bar{E}, \bar{L})$ be a positive steady state of (1). Then $k_{+} > \frac{d_E d_L}{b_E - d_L}$ and necessarily $h(\bar{L}) > \frac{d_E d_L}{b_E - d_L}$.

If $h(\bar{L}) > k_{+}$, then $(\bar{E}, \bar{L})$ is unstable if and only if $h'(\bar{L}) > T(h(\bar{L}))$. If $\frac{d_E d_L}{b_E - d_L} < h(\bar{L}) < k_{+}$, then it is unstable if and only if $h'(\bar{L}) > D(h(\bar{L}))$.

Finally, the eigenvalues of the linearized of (1) at $(\bar{E}, \bar{L})$ are complex conjugate and pure imaginary if and only if $h(\bar{L}) > k_{+}$ and $h'(\bar{L}) = T(h(\bar{L}))$.

*Proof.* This is a direct consequence of the previous calculations, except for

$$ k_{+} = \frac{d_E(b_E + 2d_E + d_L) + \sqrt{4d_E^3(b_E - d_E - d_L) + d_E^2(b_E + 2d_E + d_L)^2}}{2(b_E - d_E - d_L)} > \frac{d_E d_L}{b_E - d_L}. \qquad (14) $$

Inequality (14) is equivalent to

$$ \frac{b_E - d_L}{b_E - d_L - d_E} (b_E + d_L + 2d_E + \sqrt{(b_E + 2d_E + d_L)^2 + 4d_E(b_E - d_E - d_L)}) > 2d_L. $$
This inequality holds because $b_E > d_L$ (thanks to (3)). Indeed,

$$ \frac{b_E - d_L}{b_E - d_L - d_E} (b_E + d_L + 2d_E + \sqrt{(b_E + 2d_E + d_L)^2 + 4d_E(b_E - d_E - d_L)}) \\ > (b_E + d_L + 2d_E + \sqrt{(b_E + 2d_E + d_L)^2 + 4d_E(b_E - d_E - d_L)}) \\ > b_E + d_L \\ > 2d_L. $$

Then, setting $k = h(\bar{L})$, $k' = h'(\bar{L})$ and using the notations (11), the eigenvalues of the linearized operator are roots of the polynomial

$$ P(\lambda) = \lambda^2 - \lambda \operatorname{tr}(A) + \det(A). $$

Hence the roots are pure imaginary if and only if $\operatorname{tr}(A) = 0$ and $\det(A) > 0$. From the definition of $T, D$ in (12), $\operatorname{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_+$.

□

## 2.2 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 $\phi(L)$, as in Section 1 such that $\phi(L) \to +\infty$ as $L \to +\infty$ 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 $\bar{L}$ (which corresponds to choosing a type of breeding site). Then, to each value $k = h(\bar{L})$ corresponds a non-necessarily unique $c(k)$ that makes $\bar{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₄) the Jacobian matrix at any point X = (E, L, P, A) reads

$$J(X) = \begin{pmatrix} -\delta_E - h(L) & -Eh'(L) & 0 & \beta_E \\ h(L) & h'(L)E - \delta_L - \tau_L - 2cL & 0 & 0 \\ 0 & \tau_L & -\delta_P - \tau_P & 0 \\ 0 & 0 & \tau_P & -\delta_A \end{pmatrix},$$

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

$$P(\lambda) = (\lambda + A_1)(\lambda + A_2)(\lambda^2 + A_3\lambda + A_4) - C,$$

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, +\infty)$). This condition reads

$$\delta_A(\delta_P + \tau_P)((\delta_E + h(L))(-h'(L)E + \tau_L + \delta_L + 2cL) + Eh(L)h'(L)) > \beta_E \tau_L \tau_P h(L).$$

At equilibrium,

$$\delta_L + \tau_L + cL = \frac{h(L)}{h(L) + \delta_E} \frac{\beta_E \tau_L \tau_P}{\delta_A(\delta_P + \tau_P)},$$

therefore $P(0) > 0$ and thus any equilibrium where $h' < 0$ must be (locally) stable, in system (S₄) 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 $\bar{L}$ is actually unique? The second equation in (5) is also written

$$h(\bar{L}) = d_E \frac{d_L + c\bar{L}}{b_E - d_L - c\bar{L}}. \qquad (15)$$

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 $\kappa : L \mapsto d_E \frac{d_L + cL}{b_E - d_L - cL}$ is convex on $(0, (b_E - d_L)/c)$ and goes to $+\infty$ at $(b_E - d_L)/c$. So for instance uniqueness is guaranteed if (3) holds and either, for all $L \in (0, (b_E - d_L)/c)$, $h''(L) < 0$ or

$$h'(L) < \kappa'(L) = \frac{d_E c b_E}{(b_E - d_L - cL)^2}.$$
## 3 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.

### 3.1 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 $\varepsilon > 0$, $\eta: \mathbb{R}_+ \to \mathbb{R}_+$, and assume at first that all parameters (except for $h$) may depend on $\varepsilon$. We transform the variables $(E, L)$ from (1) into $v_\varepsilon := \varepsilon E$ and $u_\varepsilon := \frac{1}{\eta(\varepsilon)}L$. These new variables satisfy

$$ egin{cases} 
\dot{v}_{\varepsilon} = \varepsilon\eta(\varepsilon)b_{E}u_{\varepsilon} - (d_{E} + h(\eta(\varepsilon)u_{\varepsilon}))v_{\varepsilon} =: f_{\varepsilon}(u_{\varepsilon}, v_{\varepsilon}), \\ 
\varepsilon\dot{u}_{\varepsilon} = \frac{1}{\eta(\varepsilon)}h(\eta(\varepsilon)u_{\varepsilon})v_{\varepsilon} - d_{L}\varepsilon u_{\varepsilon} - c\eta(\varepsilon)\varepsilon u_{\varepsilon}^{2} =: g_{\varepsilon}(u_{\varepsilon}, v_{\varepsilon}). 
\end{cases} 
\quad (16) $$

We assume that parameters scale in such a way that the following limits exist, as $\varepsilon \to 0$:

$$ egin{cases} 
f_{\varepsilon} \xrightarrow{L^{\infty}} f, \quad g_{\varepsilon} \xrightarrow{L^{\infty}} g, \\ 
u_{\varepsilon}(t=0) = u_{\varepsilon}^{0} \to u_{0}, \quad v_{\varepsilon}(t=0) = v_{\varepsilon}^{0} \to v_{0}. 
\end{cases} 
\quad (17) $$

In addition, we assume that the zero set of $g$ is “non-degenerate” in the sense:

$$ \forall v \ge 0, \quad \{\sigma \ge 0, g(\sigma, v) = 0\} \text{ does not contain any open interval.} \quad (18) $$

We give below a simple proof of the following fact, in the spirit of Tikhonov's theorem on dynamical systems [10].

**Theorem 3.1.** Consider system (16) with $d_E, d_L, \bar{L}$ and $h$ fixed, $b_E(\varepsilon) = \frac{h(\bar{L})+d_E}{\varepsilon\bar{L}}$, $\eta(\varepsilon) = \frac{\bar{L}^2}{h(\bar{L})-\varepsilon d_L \bar{L}}$, for $\varepsilon$ small enough, and $c_\varepsilon = \frac{1}{\varepsilon\eta(\varepsilon)}$. Let $\bar{E}(\varepsilon) := 1/\varepsilon$. Then $(\varepsilon\bar{E}(\varepsilon), \frac{1}{\eta(\varepsilon)}\bar{L}) = (1, \frac{h(\bar{L})-\varepsilon d_L \bar{L}}{\bar{L}})$ is a steady state of (16) for all $\varepsilon > 0$ and (17) holds.

In addition, solutions of system (16) along with any bounded initial data admits a limit as $\varepsilon \to 0$: there exists $u, v \in L^1 \cap L^\infty(0, T)$ for all $T > 0$ such that $v_\varepsilon \to v$ uniformly and $u_\varepsilon \to u$ in $L^p(0, T)$ for all $p < \infty$.

Moreover, if initial data $u_\varepsilon^0, v_\varepsilon^0$ are such that $(\operatorname{sgn}(g_\varepsilon), \operatorname{sgn}(f_\varepsilon))(u_\varepsilon^0, v_\varepsilon^0)$ is constant for $\varepsilon$ small enough, then $(u, v)$ is periodic, $g(u(t), v(t)) = 0$ for almost every $t > 0$ and the trajectory is uniquely defined from $f$ and $g$ with $\frac{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 $(\bar{E}, \bar{L})$ of (1) are given by (5):

$$ d_L + c\bar{L} = \frac{h(\bar{L})b_E}{h(\bar{L}) + d_E}, \quad \bar{E} = \frac{b_E\bar{L}}{h(\bar{L}) + d_E}. $$
<figure><img src="image_3.png" /><figcaption>Figure 1: u is in x-axis, v in y-axis. Red dashed curves correspond to nullclines $g_{\epsilon} = 0$ ($\dot{u}_{\epsilon} = 0$) and blue dotted curves to nullclines $f_{\epsilon} = 0$ ($\dot{v}_{\epsilon} = 0$). The four figures correspond to decreasing values of $\epsilon$ from top-left to bottom-right (0.5, 0.1, 0.01 and 0.001). In yellow and purple, two trajectories $t \mapsto (u_{\epsilon}(t), v_{\epsilon}(t))$ are shown, for two different initial conditions (respectively (0.35, 1.2) and (0.2, 0.5)).</figcaption></figure>

Thus in all generality (allowing all parameters to depend on $\epsilon$), the scalings fit for our purpose (i.e. with $\bar{E}(\epsilon) = 1/\epsilon$) are exactly those for which $\epsilon = \frac{h(\bar{L}(\epsilon)) + d_E(\epsilon)}{b_E(\epsilon)\bar{L}(\epsilon)}$ and there exists $\eta(\epsilon) = O(1)$ such that

$$d_L(\epsilon)(h(\bar{L}(\epsilon)) + d_E(\epsilon)) + \frac{b_E(\epsilon)\bar{L}^2(\epsilon)}{\eta(\epsilon)} = h(\bar{L}(\epsilon))b_E(\epsilon).$$

It turns out that $\bar{L}(\epsilon)(\epsilon(b_E - d_L) - \frac{1}{\eta}) = d_E$. Hence to guarantee $\eta(\epsilon) = O(1)$ it is required that

$$b_E - d_L = O(1/\epsilon).$$

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 $\epsilon c_{\epsilon}$ and $\epsilon b_E(\epsilon)$ have limits $1/\eta_0, \xi > 0$ respectively as $\epsilon \to 0$ (this is the case with the scaling
used in Theorem 3.1), we choose $\eta(\varepsilon) = O(1)$ such that $c_\varepsilon\eta(\varepsilon)\varepsilon = 1$ and end up with

$$ 
\left\{ 
\begin{array}{l}
\dot{v} = \varepsilon\eta b_E u - (d_E + h(\eta u))v =: f_\varepsilon(u, v), \\
\varepsilon\dot{u} = \frac{1}{\eta}h(\eta u)v - d_L\varepsilon u - u^2 =: g_\varepsilon(u, v).
\end{array}
\right. 
$$

The limits $f$ and $g$ are given by

$$ f(u, v) = \eta_0 \xi u - (d_E + h(\eta_0 u))v, \quad g(u, v) = \frac{1}{\eta_0} h(\eta_0 u) - u^2. \qquad (19) $$

## 3.2 Proof of the main result

We proceed to the proof of Theorem 3.1 in three steps. First, scaled quantities $u_\varepsilon$ and $v_\varepsilon$ remain uniformly bounded independently of $\varepsilon$, 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 $\varepsilon > 0$ and $t > 0$,

$$ |u_{\varepsilon}(t)|, |v_{\varepsilon}(t)|, |f_{\varepsilon}(u_{\varepsilon}(t), v_{\varepsilon}(t))|, |g_{\varepsilon}(u_{\varepsilon}(t), v_{\varepsilon}(t))| \leq C. $$

Hence, up to extraction, $v_\varepsilon$ converges to $v$ uniformly on compact sets $[0, T]$ by the Ascoli theorem. Then, the convergence of an auxiliary quantity gives convergence of $u_\varepsilon$:

**Lemma 3.3.** For all $T > 0$.

$$ \|g_{\varepsilon}(u_{\varepsilon}, v_{\varepsilon})\|_{L^2(0,T)} = O(\sqrt{\varepsilon}). \qquad (20) $$

Moreover, there exists $u, v \in L^1 \cap L^\infty$ such that after extraction of a subsequence $u_\varepsilon \to u$ in $L^p(0, T)$ for all $1 \le p < \infty$, as $v_\varepsilon \to v$ uniformly.

*Proof.* Let $B(t, u) := \int_0^u g^2(\sigma, v(t))d\sigma$, where $v$ is the limit of $v_\varepsilon$ (obtained by the Ascoli theorem) and $g$ is the limit of $g_\varepsilon$ (from (17)). From (18) we deduce that for all $t$, $u \mapsto 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) \in L^p(0, T)$ for all $p < \infty$ and $T > 0$ such that

$$ \int_0^{u_\varepsilon(t)} g_\varepsilon^2(\sigma, v_\varepsilon(t)) d\sigma \xrightarrow{\varepsilon \to 0} w(t), \qquad (21) $$

then defining $u(t) := A(t, w(t))$ we can conclude that $u_\varepsilon = A(\cdot, \int_0^{u_\varepsilon} g^2(\sigma, v) d\sigma) \xrightarrow{\varepsilon \to 0} u = A(\cdot, w)$.

Indeed, we notice that

$$ \int_0^{u_\varepsilon(t)} g_\varepsilon^2(\sigma, v_\varepsilon(t)) d\sigma - \int_0^{u_\varepsilon(t)} g^2(\sigma, v_\varepsilon(t)) d\sigma \to 0, $$

and

$$ \int_0^{u_\varepsilon(t)} g^2(\sigma, v_\varepsilon(t)) d\sigma - \int_0^{u_\varepsilon(t)} g^2(\sigma, v(t)) d\sigma \to 0. $$
Since $u_\epsilon$ is uniformly bounded,

$$ \left| \int_{0}^{u_{\epsilon}(t)} g_{\epsilon}^{2}(\sigma, v_{\epsilon}(t)) d\sigma - \int_{0}^{u_{\epsilon}(t)} g^{2}(\sigma, v(t)) d\sigma \right| \leq u_{\epsilon}(t) (\|g_{\epsilon}^{2} - g^{2}\|_{\infty} + C \|v_{\epsilon} - v\|_{\infty}), $$

for some $C > 0$ which depends only on $\partial_v g$. Hence (21) implies

$$ \int_{0}^{u_{\epsilon}(t)} g^{2}(\sigma, v(t)) d\sigma \xrightarrow{\epsilon \to 0} w(t). $$

Therefore we only need to prove (21) to complete the proof. To do so we first obtain (20) by computing

$$ \begin{aligned} \frac{g_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t))^{2}}{\epsilon} &= g_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t)) \dot{u}_{\epsilon} \\ &= \frac{d}{dt} \int_{0}^{u_{\epsilon}(t)} g_{\epsilon}(\sigma, v_{\epsilon}(t)) d\sigma - f_{\epsilon}(u_{\epsilon}, v_{\epsilon}) \int_{0}^{u_{\epsilon}(t)} \partial_{v} g_{\epsilon}(\sigma, v_{\epsilon}(t)) d\sigma. \end{aligned} $$

Hence

$$ \frac{1}{\epsilon} \int_{0}^{T} (g_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t)))^2 dt = \int_{u_{\epsilon}(0)}^{u_{\epsilon}(T)} g_{\epsilon}(\sigma, v_{\epsilon}(t)) d\sigma - \int_{0}^{T} f_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t)) \int_{0}^{u_{\epsilon}(t)} \partial_{v} g_{\epsilon}(\sigma, v_{\epsilon}(t)) d\sigma dt. $$

Since $f_\epsilon, g_\epsilon$ and $\partial_v g_\epsilon = \frac{1}{\eta(\epsilon)}h(\eta(\epsilon)u_\epsilon)$ are uniformly bounded, we deduce that

$$ \int_{0}^{T} g_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t))^{2} dt = O(\epsilon). $$

This gives (20). Then we introduce

$$ w_{\epsilon}(t) := \int_{0}^{u_{\epsilon}(t)} g_{\epsilon}^{2}(\sigma, v_{\epsilon}(t)) d\sigma. $$

We compute

$$ \dot{w}_{\epsilon}(t) = \frac{1}{\epsilon} g_{\epsilon}^{2}(u_{\epsilon}(t), v_{\epsilon}(t)) \epsilon \dot{u}_{\epsilon} + f_{\epsilon}(u_{\epsilon}(t), v_{\epsilon}(t)) \int_{0}^{u_{\epsilon}(t)} 2g_{\epsilon}(\sigma, v_{\epsilon}(t)) \partial_{v} g_{\epsilon}(\sigma, v_{\epsilon}(t)) d\sigma. $$

By the previous point, $t \mapsto \frac{1}{\epsilon} g_{\epsilon}^{2}(u_{\epsilon}(t), v_{\epsilon}(t))$ is uniformly (in $\epsilon$) bounded in $L^1$. In addition, $t \mapsto \epsilon \dot{u}_{\epsilon}(t)$ is uniformly (in $\epsilon$) bounded in $L^{\infty}$, by the Lemma 3.2. The second term $f_{\epsilon} \int g_{\epsilon} \partial_{v} g_{\epsilon}$ is uniformly bounded as well.

As a consequence, $w_\epsilon$ is uniformly (in $\epsilon$) bounded in $BV_{\text{loc}}$. This implies that up to extraction, $w_\epsilon \to w$ in $L^1$. Because $w_\epsilon$ is also bounded in $L^\infty$, convergence actually takes place in all $L^p$ spaces. $\square$

<figure><img src="image_4.png" /></figure>

Finally, the shapes of $(f, g)$ allow us to describe simply the limit trajectories. We use the following assumptions: for all $\epsilon > 0$ small enough, we assume that the right-hand sides of system (16) satisfy

(R.1) the set $\mathbb{R}^2 \setminus \{f_\epsilon = 0, g_\epsilon = 0\}$ has exactly 4 connected components, whose measures do not vanish as $\epsilon \to 0$,
(R.2) $f(u_0, v_0) \neq 0$, $g(u_0, v_0) \neq 0$ and the couple $(\operatorname{sgn}(f_{\varepsilon}(u_0^{\varepsilon}, v_0^{\varepsilon})), \operatorname{sgn}(g_{\varepsilon}(u_0^{\varepsilon}, v_0^{\varepsilon})))$ is constant and equal to $(\operatorname{sgn}(f(u_0, v_0)), \operatorname{sgn}(g(u_0, v_0)))$.

We also assume that the uniform limits $f, g$ of $f_{\varepsilon}, g_{\varepsilon}$ satisfy

(L.1) the curve $\Upsilon := \{g = 0\}$ is the graph of a function $\phi \in C^1(\mathbb{R}_+, \mathbb{R}_+)$ with $\phi(\infty) = \infty$ and $\phi(0) = 0$,

(L.2) the function $g$ is positive on the epigraph of $\phi$,

(L.3) the function $\phi$ has exactly two local extrema,

(L.4) on the graph of $\phi$, $\operatorname{sgn}(f) = -1$ except for a bounded set.

**Lemma 3.4.** With these assumptions we have:

There exists a unique $\tau > 0$ and a (unique up to translations) $\tau$-periodic function $(u_{\tau}, v_{\tau}) : \mathbb{R}_+ \to \Upsilon$ such that $v_{\tau}$ is Lipschitz-continuous, $u_{\tau}$ is piecewise continuous, for all $t \ge 0$, $v_{\tau} = \phi(u_{\tau})$ everywhere, $\dot{v}_{\tau} = f(u_{\tau}, v_{\tau})$ almost everywhere and the discontinuities of $u_{\tau}$ are located at times $t$ such that $\phi$ has a local extremum at $u_{\tau}(t^-)$.

There exists $\tau_1 \ge 0$ and $\tau_2 \in [0, \tau)$ such that for all $t > \tau_1$, $(u, v)(t) = (u_{\tau}, v_{\tau})(t + \tau_2)$. Moreover, by construction $\tau_1$ and $\tau_2$ are uniquely defined from $u_0$ and $v_0$, so the limit $(u, v)$ is in fact unique and the whole family $(u_{\varepsilon}, v_{\varepsilon})_{\varepsilon}$ converges as $\varepsilon$ goes to 0.

Clearly from (19), Lemma 3.4 applies with the hypotheses of Theorem 3.1 and

$$ \phi(u) = \frac{\eta_0 u^2}{h(\eta_0 u)}, \quad \eta_0 = \frac{\bar{L}^2}{h(\bar{L})}, \qquad (22) $$

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_{\tau}, v_{\tau})$ is classical and can be done by pasting together solutions of Cauchy problems given (locally) by $\dot{v}_{\tau} = f(\phi^{-1}(v_{\tau}), v_{\tau})$, on intervals where $\phi$ is invertible. Uniqueness comes from the crucial fact that discontinuities of $u_{\tau}$ are assumed to be located at local extrema of $\phi$.

From the previous lemmas we know that $(u, v) \in \Upsilon$ almost everywhere. In addition, uniform boundedness of $f_{\varepsilon}(u_{\varepsilon}, v_{\varepsilon})$ ensures that $v$ is Lipschitz continuous.

Then, we claim that if $t > 0$ is such that $\phi$ has no local extremum at $u(t^+)$, then there exists $\tau_0 > 0$ such that $(u, v)$ is continuous on $(t, t + \tau_0)$. This point is the key of the proof. To prove it, let $u_i$ be such that $\phi'(u_i) < 0$. We solve only the simpler problem

$$ \left\{ \begin{array}{ll} \dot{\hat{v}}_{\varepsilon} = f(\hat{u}_{\varepsilon}, \hat{v}_{\varepsilon}), & \hat{v}_{\varepsilon}(0) = \phi(u_i) + O(\varepsilon), \\ \varepsilon \dot{\hat{u}}_{\varepsilon} = g(\hat{u}_{\varepsilon}, \hat{v}_{\varepsilon}), & \hat{u}_{\varepsilon}(0) = u_i + O(\varepsilon). \end{array} \right. $$

Introducing $\hat{w}_{\varepsilon} := \hat{u}_{\varepsilon} - \phi^{-1}(\hat{v}_{\varepsilon})$, where the inverse of $\phi$ is taken locally (this is possible for $\varepsilon$ small enough since $\phi'(u_i) < 0$ and $\hat{v}_{\varepsilon}$ is uniformly Lipschitz-continuous), we obtain

$$ \dot{\hat{w}}_{\varepsilon} = \frac{\hat{w}_{\varepsilon}}{\varepsilon} \partial_1 g(\hat{r}_{\varepsilon}(t), \hat{v}_{\varepsilon}(t)) - \frac{f(\hat{u}_{\varepsilon}, \hat{v}_{\varepsilon})}{\phi'(\phi^{-1}(\hat{v}_{\varepsilon}))}, \quad \hat{w}_{\varepsilon}(0) = O(\varepsilon), $$

for some $\hat{r}_{\varepsilon}(t)$ between $\hat{u}_{\varepsilon}(t)$ and $\phi^{-1}(\hat{v}_{\varepsilon}(t))$. We have $\partial_1 g \le -\alpha < 0$ on a neighborhood of $(u_i, \phi(u_i))$, so on this neighborhood $\hat{w}_{\varepsilon}$ remains small (it is a $o(\varepsilon)$), which in turn proves that $(\hat{r}_{\varepsilon}, \hat{v}_{\varepsilon})$ remains in
this neighborhood. In particular, $\hat{u}_{\varepsilon}$ converges to some function $\hat{u}$ which is continuous at $t = 0$ (since it is equal to $\phi^{-1}(\hat{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_{\varepsilon}, g_{\varepsilon}$ rather than $f, g$, and rise some analogue $\phi_{\varepsilon}$ of $\phi$ at level $\varepsilon > 0$, for $\varepsilon$ 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 $\phi$ has a local minimum equal to $\phi_m$ at $u_m$ and a local maximum equal to $\phi_M > \phi_m$ at $u_M < u_m$. Moreover, let $u_m^0 < u_M$ such that $\phi(u_m^0) = \phi(u_m)$. For $\alpha, \beta \in \{1, -1\}$, we also introduce $Z_{\alpha}^{\beta} := \{\operatorname{sgn}(f) = \alpha, \operatorname{sgn}(g) = \beta\}$.

We define a mapping $\pi : \mathbb{R}^2 \to \Upsilon$ by $\pi = \mathrm{Id}$ on $\Upsilon$ and if $(u, v) \in Z_{\alpha}^{\beta}$ then $\pi(u, v) = (u_1, v)$ such that $\phi(u_1) = v$ and $\operatorname{sgn}(u_1 - u) = \beta$. The projection $\pi$ is well-defined thanks to the assumptions on $\phi$ and $g$, except on $(u_m^0, +\infty) \times \{\phi_m\}$, on which we let $\pi = (u_m^0, \phi_m)$. Then $(u, v)(0^+) = \pi(u_0, v_0)$. To prove this, one simply has to check the behavior of $u_{\varepsilon}$ (since $v_{\varepsilon}$ and $v$ are Lipschitz continuous). As above, we claim that the first-order behavior is simply given by the “layer equation”

$$ \varepsilon \dot{\tilde{u}}_{\varepsilon} = g(\tilde{u}_{\varepsilon}, v_0), \quad \tilde{u}_{\varepsilon}(0) = u_0, $$

which makes $\tilde{u}_{\varepsilon}$ converge exponentially fast to $\pi(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^{\varepsilon}, v_{\varepsilon}$ and $g_{\varepsilon}$.

Let $\Upsilon_u = \Upsilon \cap ([u_M, u_m] \times \mathbb{R}_+)$ and $\Upsilon_s = \Upsilon - \Upsilon_u$. (Note that $\pi(\mathbb{R}_+^2 - \Upsilon_u) = \Upsilon_s$.) After the initial layer, the trajectory of $(u, v)$ remains on $\Upsilon_s$. This follows from the sign of $f$ on $\Upsilon_u$: because of the continuity property, the trajectory cannot exit $\Upsilon_s$ but at $(u_m, \phi_m)$ (or $(u_M, \phi_M)$, respectively). At these points however, $\Upsilon_u$ is repulsive since $v$ must be continuous, $\dot{v} < 0$ ($\dot{v} > 0$, respectively) and $\Upsilon_u$ lies locally in $\{v > \phi_m\}$ (respectively in $\{v < \phi_M\}$).

Still, the initial data does not need to be projected directly by $\pi$ on $\Upsilon_s \cap \mathbb{R}_+ \times [\phi_m, \phi_M]$. Therefore, we introduce $\tau_1 \ge 0$ as

$$ \tau_1 := \max(0, \sup\{t \ge 0, v(t) \notin [\phi_m, \phi_M]\}). $$

It remains to check that $\tau_1 < +\infty$. For all $T > 0$, as long as $\phi$ has no local extremum at $u(t)$ for $t \in (0, T)$, $u$ is continuous. Thanks to our assumption (R.1), there are two connected components in $\Upsilon_s$, on each one of whom $\operatorname{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 $\Upsilon$ where $|f|$ is positively bounded from below (one of the two connected components of $\Upsilon_s$) and has the appropriate sign. This yields the existence of $\tau_1 < +\infty$.

Then for all $t \ge \tau_1$ we have $v(t) \in [\phi_m, \phi_M]$, and the trajectory is uniquely defined onwards. $\square$

**Remark 3.5.** We did not treat the case when the limit of $(u_0^{\varepsilon}, v_0^{\varepsilon})$ belongs to $\Upsilon$ (relaxing assumption (R.2)). In this case indeed, no general result can be obtained, unless the various convergence speeds (of $f_{\varepsilon}, g_{\varepsilon}, u_0^{\varepsilon}$ and $v_0^{\varepsilon}$) 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 $\varepsilon \to 0$ can be computed if one knows these parameter scale in $\varepsilon$ thanks to only $f$ and $g$. Their period $\tau$ can also be computed directly from $f$ and $\phi$. 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^{\infty}$) as $[u_m^0, u_M]$ and $[u_m, u_M^0]$ respectively, and let

$$ \Psi(u, v) = \int_{u}^{v} \frac{\phi'(u')}{f(u', \phi(u'))} du'. $$
Then we have

$$ \tau = \Psi(u_m^0, u_M) + \Psi(u_M^0, u_m). \qquad (23) $$

# 4 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.

## 4.1 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 $\bar{L}$. With this setting, we can state a bifurcation theorem using the simple bifurcation parameter $h'(\bar{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 $\bar{L}$ the class of functions under consideration that fits our purposes is

$$ H_L := \{h(L) = a \left( \arctan(b(L - \bar{L})) + \frac{\pi}{2} \right), a, b \in \mathbb{R}^+ \}. \qquad (24) $$

Graphs of these functions are shown in Figure 2. We use the immediate properties that these functions

<figure><img src="image_5.png" /><figcaption>Figure 2: Function $h \in H_L$ with $\bar{L} = 2$. Left: $a = 1$, $b = \{0.5, 1, 2\}$. Right: $a = \{1, 1.5, 2\}$, $b = 2$. Curve styles with increasing values in $a$ and $b$: dotted blue, solid red, dashed yellow.</figcaption></figure>

are positive and increasing. For any couple $(k, k') \in \mathbb{R}_+^* \times \mathbb{R}_+^*$, there exists a unique function $h$ of class $H_L$ with $h(\bar{L}) = k$ and $h'(\bar{L}) = k'$. Finally, for all $c > 0$, the steady state relation $h(\bar{L}) = \frac{d_E(d_L+c\bar{L})}{b_E-d_L-c\bar{L}}$ has a positive solution in $\bar{L}$ if $a > \frac{d_E d_L}{b_E-d_L} \frac{2}{\pi}$. Indeed, for given values $(k, k') \in \mathbb{R}_+^2$, the choice of $a = \frac{2k}{\pi}$
and $b = \frac{k'}{a}$ gives the solution since

$$h(\bar{L}) = a\frac{\pi}{2} = \frac{2k}{\pi}\frac{\pi}{2} = k \text{ and } h'(\bar{L}) = ab = \frac{2k}{\pi}\frac{k'}{\frac{2k}{\pi}} = k'$$

Also we can solve the equation in $\bar{L}$, $\frac{a\pi}{2} = \frac{d_E(d_L+c\bar{L})}{b_E-d_L-c\bar{L}}$, which yields $\bar{L} = \frac{a\pi(b_E-d_L)-d_Ed_L}{c(a\pi+cd_E)}$. Hence $\bar{L}$ is positive under the stated condition.

**Remark 4.1.** From Lemma 2.2, for $h$ of class $H_{\bar{L}}$, the state $(0,0)$ is unstable if and only if

$$a > \frac{d_E d_L}{(b_E - d_L)(\frac{\pi}{2} + \arctan(-b\bar{L}))}.$$

## 4.2 Transformation into a canonical form

Let $P = (a,b) \in \mathbb{R}_+^2$ and the function $h_P$ of class $H_{\bar{L}}$

$$h_P(L) = a\left(\arctan(b(L - \bar{L})) + \frac{\pi}{2}\right). \qquad (25)$$

We use the notation $k := h_P(\bar{L}) = a\frac{\pi}{2}$. Let $P: \gamma \mapsto P(\gamma) = (a_0, b_0 + \gamma)$ where $(a_0, b_0) \in \mathbb{R}_+^2$. Then we can associate $P(\gamma)$ to a new system $(S_\gamma(a_0, b_0))$ obtained from (1)

$$\left\{ \begin{array}{l} \dot{E} = b_E L - d_E E - h_{P(\gamma)}(L)E, \\ \dot{L} = h_{P(\gamma)}(L)E - d_L L - cL^2. \end{array} \right. \qquad (S_\gamma(a_0, b_0))$$

This system has a positive equilibrium $(\bar{E}, \bar{L})$ and the Jacobian matrix of the system evaluated in $(\bar{E}, \bar{L})$ is:

$$J_{P(\gamma)} = \begin{pmatrix} -d_E - h_{P(\gamma)}(\bar{L}) & b_E - h'_{P(\gamma)}(\bar{L})\bar{E} \\ h_{P(\gamma)}(\bar{L}) & h'_{P(\gamma)}(\bar{L})\bar{E} - d_L - 2c\bar{L} \end{pmatrix},$$

We set $\lambda_{1,2}(\gamma) = \alpha(\gamma) \pm i\beta(\gamma)$ the eigenvalues of $J_{P(\gamma)}$, when the discriminant of the characteristic polynomial of $J_{P(\gamma)}$ is negative.

## 4.3 Main result

Using function $T$ from (12), we define

$$b(a) := \frac{T(a)}{a}, \quad a_{crit} := \frac{2k^+}{\pi} > 0. \qquad (26)$$

**Theorem 4.2.** There exists $\tilde{a} > 0$ such that: If $a > \max(\tilde{a}, a_{crit})$, $(S_\gamma(a, b(a)))$ has a supercritical Hopf Bifurcation in $\gamma = 0$. In particular:

1. there exists $\gamma_1 < 0$ such that for all $\gamma \in (\gamma_1, 0]$, $(\bar{E}, \bar{L})$ is a stable focus,

2. for all $U$ neighborhood of $(\bar{E}, \bar{L})$, there exists $\gamma_2 > 0$ such that for all $\gamma \in [0, \gamma_2)$, $(\bar{E}, \bar{L})$ is an unstable focus surrounded by a stable limit cycle contained in $U$, which has an amplitude that grows when $\gamma$ grows.
**Remark 4.3.** $k_+$ is given by (13), and $\tilde{a}$ is such that the normal form coefficient $\alpha_N$ (see [19]) of our system is negative if $a > \tilde{a}$. We simply give a numerical justification of the existence of $\tilde{a}$ 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_\gamma(a_0, b_0)$ in Figure 3 is obtained by XPPAUT software [9].

<figure><img src="image_6.png" /><figcaption>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$.</figcaption></figure>

*Proof of Theorem 4.2.* We set $\lambda_{1,2}(\gamma) = \alpha(\gamma) \pm i\beta(\gamma)$ (with $\gamma$ a real parameter), the two eigenvalues of $J_P(\gamma)$ the Jacobian matrix associated to our system and computed in (0,0). We call $\gamma_c$ a bifurcation value, and $\alpha_N(\gamma)$ the normal form coefficient of the system (see [19]).

Firstly, we only need to study complex conjugate and pure imaginary eigenvalues of $J_P(\gamma)$ to find the bifurcation value $\gamma_c$, which means also to look for $\gamma_c$ such that $\alpha(\gamma_c) = 0$ and $\beta(\gamma_c) \neq 0$. Thanks to Proposition 2.6 we know that this is the case when $k > k_+$ i.e. $\frac{a\pi}{2} > k_+$ or equivalently $a > a_{crit}$ (by definition, $a_{crit} = \frac{2k_\pm}{\pi}$). Moreover since $h'(\bar{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))

$$G := \{(a, b) \in \mathbb{R}^2, a > a_{crit}, T(a) = ab = h'_{P(\gamma)}(\bar{L})\}. \qquad (27)$$

And we can set $\gamma_c = 0$.

Secondly, we have to see if $\frac{d\alpha}{d\gamma}(\gamma_c) > 0$, this means to check that $\text{tr}(J_P(\gamma))$ changes sign at the bifurcation value $\gamma_c$. Let $\gamma \mapsto z(\gamma) = \alpha(a_0, b(a_0) + \gamma)$. We recall that $\alpha(\gamma)$ is a function of $a$ and $b$. Since

$$\alpha = \frac{\text{tr}(J_P(\gamma))}{2} = \frac{1}{2} \left( -d_E - \frac{a\pi}{2} + ab\bar{E} - d_L - 2c\bar{L} \right),$$

we have $z'(\gamma) = \partial_b\alpha = \frac{a\bar{E}}{2}$ and we obtain that $z'(\gamma_c) = z'(0) = \frac{a\bar{E}}{2}$ and it is always positive.
Thirdly, we have to study the normal form coefficient of the system computed in $\gamma_c = 0$ and find when $\alpha_N(\gamma_c) \neq 0$. To get the normal form coefficient, we have to transform the system ($S_\gamma(a_0, b_0)$) and we use the steps from [19]. In a first step we reduce the initial system ($S_\gamma(a_0, b_0)$) to a system where the equilibrium ($\bar{E}, \bar{L}$) becomes the origin. By the change of variables $x = E - \bar{E}$ and $y = L - \bar{L}$, ($S_\gamma(a_0, b_0)$) becomes:

$$ \begin{cases} \dot{x} = b_E(y + \bar{L}) - d_E(x + \bar{E}) - a(\arctan(by) + \frac{\pi}{2})(x + \bar{E}), \\ \dot{y} = a(\arctan(by) + \frac{\pi}{2})(x + \bar{E}) - d_L(y + \bar{L}) - c(y + \bar{L})^2. \end{cases} \quad (28) $$

Then as ($\bar{E}, \bar{L}$) is an equilibrium, we can simplify (28) into

$$ \begin{cases} \dot{x} = b_E y - d_E x - \frac{a\pi}{2} x - a(\arctan(by))(x + \bar{E}), \\ \dot{y} = a(\arctan(by) + \frac{\pi}{2})(x + \bar{E}) + \frac{a\pi}{2} x - d_L y - c y^2 - 2cy\bar{L}, \end{cases} \quad (29) $$

which we write as

$$ \begin{cases} \dot{x} = b_E y - d_E x - \frac{a\pi}{2} x - aby\bar{E} + f(x, y), \\ \dot{y} = \frac{a\pi}{2} x - d_L y - 2cy\bar{L} + aby\bar{E} + g(x, y), \end{cases} \quad (30) $$

where

$$ f(x, y) = aby\bar{E} - a \arctan(by)(x + \bar{E}), \quad g(x, y) = -aby\bar{E} + a \arctan(by)(x + \bar{E}) - cy^2. $$

The system (30) can also be written under the matrix form

$$ \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} -\frac{a\pi}{2} - d_E & b_E - ab\bar{E} \\ \frac{a\pi}{2} & -d_L - 2c\bar{L} + ab\bar{E} \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} + \begin{pmatrix} f(x, y) \\ g(x, y) \end{pmatrix}. $$

We call $M$ the first $(2 \times 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

$$ \begin{pmatrix} \dot{X} \\ \dot{Y} \end{pmatrix} = N \begin{pmatrix} X \\ Y \end{pmatrix} + \begin{pmatrix} F(X, Y) \\ G(X, Y) \end{pmatrix}, \quad N := \begin{pmatrix} 0 & -\omega \\ \omega & 0 \end{pmatrix}. \quad (31) $$

In our case, we can have an idea of the normal coefficient only in a neighborhood of $\gamma = 0$. Because we want to make a simple linear change of variables, we are looking for a matrix $P$ such that $PMP^{-1} = N$ and that at the bifurcation value $\gamma = 0$, $\text{tr}(M) = 0 = \text{tr}(N)$ and $\det(M) = -A^2 - BC = \omega^2 > 0$.

$$ \text{We set } M = \begin{pmatrix} A & B \\ C & -A \end{pmatrix} \text{ and we can choose } P = \begin{pmatrix} \frac{\omega+A}{2B\omega} & \frac{1}{2\omega} \\ \frac{\omega-A}{2B\omega} & -\frac{1}{2\omega} \end{pmatrix}, P^{-1} = \begin{pmatrix} B & B \\ \omega-A & -A-\omega \end{pmatrix}. $$

Next we obtain the matrix system (31) where

$$ \begin{pmatrix} X \\ Y \end{pmatrix} = P \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} \frac{x(\omega+A)}{2B\omega} + \frac{y}{2\omega} \\ \frac{x(\omega-A)}{2B\omega} - \frac{y}{2\omega} \end{pmatrix}, \quad \begin{pmatrix} x \\ y \end{pmatrix} = P^{-1} \begin{pmatrix} X \\ Y \end{pmatrix} = \begin{pmatrix} (X+Y)B \\ (\omega-A)X + (-\omega-A)Y \end{pmatrix}, $$
$$ \begin{pmatrix} F(X, Y) \\ G(X, Y) \end{pmatrix} = P \begin{pmatrix} f(x, y) \\ g(x, y) \end{pmatrix} = \begin{pmatrix} \frac{\omega+A}{2B\omega} f(x, y) + \frac{1}{2\omega} g(x, y) \\ \frac{\omega-A}{2B\omega} f(x, y) - \frac{1}{2\omega} g(x, y) \end{pmatrix} = \begin{pmatrix} f(x, y) \left( \frac{\omega+A}{2B\omega} - \frac{1}{2\omega} \right) - \frac{1}{2\omega} c y^2 \\ f(x, y) \left( \frac{\omega-A}{2B\omega} + \frac{1}{2\omega} \right) + \frac{1}{2\omega} c y^2 \end{pmatrix}. $$

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:

$$ \alpha_N(\gamma = 0) = \frac{1}{16} \left( F_{XXX} + F_{YYY} + G_{XXY} + G_{YYY} \right) - \frac{1}{16\omega} \left( G_{XY}(G_{XX} + G_{YY}) - F_{XY}(F_{XX} + F_{YY}) + F_{XX}G_{XX} - F_{YY}G_{YY} \right). $$

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 $\tilde{a} > 0$, and is always negative for $a > \tilde{a}$ (as it appears that $\tilde{a} < a_{crit}$, this is sufficient by definition of (27)). Then $\alpha_N(\gamma_c) \neq 0$ for $a \neq \tilde{a}$.

Finally, we want to have for all real $\gamma$ in a neighborhood of 0, $\alpha_N(\gamma)\alpha(\gamma) < 0$. Thanks to Maple we have $\alpha_N(0) < 0$, in a neighborhood of $\gamma = 0$, for $a > \tilde{a}$ with $\tilde{a}$ small.

So we can apply the Hopf bifurcation theorem that ensures there exists a limit cycle (periodic solution) when $\alpha(\gamma) > 0$ (i.e $\det(J_P(\gamma)) > 0$), and moreover this cycle is stable as $\alpha(\gamma) > 0$: we are faced to a supercritical bifurcation. $\square$

## 4.4 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 $\gamma \to 0^+$, the periodic solution of the system $(S_\gamma(a, b(a)))$ has a frequency $\omega$ and a period $T_0 = 2\pi/\omega$ given by the expression

$$ \omega = \frac{1}{\sqrt{d_E + k}} \left[ k^2(b_E - d_L - d_E) + k(-2d_E^2 - b_E d_E - d_E d_L) - d_E^3 \right]^{\frac{1}{2}}. $$

*Proof.* As $\gamma \to 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 $\omega_\gamma = \sqrt{\det(J_{P(\gamma)})}$, where the expression of $\det(J_{P(\gamma)})$ is

$$ \det(J_{P(\gamma)}) = \frac{1}{d_E + k} \left[ k^2(b_E - d_L - d_E) + k(-2d_E^2 - b_E d_E - d_E d_L) - d_E^3 \right]. $$

Then the expression of $\omega = \omega_0$ follows. $\square$

**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:

$$ k^2 \left[ T_0^2 (b_E - d_L - d_E) \right] + k \left[ T_0^2 (-2d_E^2 - b_E d_E - d_E d_L) - 4\pi^2 \right] - T_0^2 d_E^3 - 4\pi^2 d_E. $$

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 $\gamma$. In other words, for fixed a the amplitude of the oscillations is an increasing function of b.
# 5 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.

**Acknowledgements.** BP has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623). MS, NV and DAM acknowledge partial funding from Inria, France and CAPES, Brazil (processo 99999.007551/2015-00), in the framework of the STIC AmSud project MOSTICAW and from CAPES/COFECUB project Ma-833 15 "Modeling innovative control method for Dengue fever". MS and NV acknowledge partial funding from the ANR blanche project Kibord: ANR-13-BS01-0004 funded by the French Ministry of Research.

## A Observations on a class of hatching functions

Among the many possible choices for a S-shaped hatching function $h$, we numerically and theoretically explore the typical family of Hill functions. We assume the following form with parameters $a, \lambda, p > 0$

$$h(L) = h_m + a \frac{L^p}{\lambda^p + L^p}, \quad h_m > \frac{d_E d_L}{b_E - d_L} \quad (32)$$

Steady states $(\bar{E}, \bar{L})$ of (1) are such that $\bar{L}$ is a solution of $Q(L) = 0$, where

$$Q(L) = -cL^{p+1}(h_m + a + d_E) + L^p((h_m + a)(b_E - d_L) - d_E d_L) - c\lambda^p L(h_m + d_E) + \lambda^p(h_m(b_E - d_L) - d_E d_L).$$

The following lemma is a straightforward consequence of this computation
**Lemma A.1.** When $p = 1$ and $h$ is of type (32), there is a unique steady state of (1).

When $p = 1$ and $h$ is of type (32), then $h'' < 0$, a property that is lost when $p > 1$. Therefore, to simplify the choice of the parameters, now we assume

$$d_E = 0. \qquad (33)$$

Then, condition (6) is fulfilled, the steady state of (1) is unique and is given by

$$\bar{L} = \frac{b_E - d_L}{c}, \quad \bar{E} = \frac{b_E \bar{L}}{h(\bar{L})}.$$

**Proposition A.2.** Let $h$ be of type (32) and assume condition (33) holds. Then (1) has a unique positive steady state $(\bar{E}, \bar{L})$ and its linearization has eigenvalues with negative real parts if and only if

$$k > \frac{a}{1 + \alpha^p} > \frac{\alpha^p + 1}{p\alpha^p} k\left(2 + \frac{k - d_L}{b_E}\right), \quad \alpha = \frac{\lambda}{\bar{L}}, \quad k = h(\bar{L}). \qquad (34)$$

*Proof.* Necessarily $k > k_+ = 0$ (where $k_+$ is defined in (13)). Hence the eigenvalues of the linearized system at $(\bar{E}, \bar{L})$ are in $C\backslash\mathbb{R}$ and the condition for instability of the steady state from Proposition 2.6 simply reads $h'(\bar{L}) > T(h(\bar{L}))$, where $T$ is defined in (12). Then we compute

$$h(\bar{L}) = h_m + \frac{a}{1 + \alpha^p}, \quad h'(\bar{L}) = \frac{ap\alpha^p}{\bar{L}(\alpha^p + 1)^2}.$$

The right-hand side inequality in (34) comes from $h'(\bar{L}) > T(k)$ and the left-hand side from $h_m > 0$. □

If all parameters but $\alpha$ and $a$ are fixed, then condition (34) can be fulfilled if and only if

$$k < (p - 2)b_E + d_L. \qquad (35)$$

Indeed, we need to find $\alpha > 0$ such that $2 + \frac{k-d_L}{b_E} < p \frac{\alpha^p}{1+\alpha^p}$. Note that in particular, this is impossible when $p \le 1$ (since $c\bar{L} = b_E - d_L > 0$ by hypothesis).

We provide below numerical results showing consistent oscillations under condition (35), for 2- and 3-dimensional systems (1) and (S$_3$). To explore the possible behaviors depending on the function $h$ of type (32), we fix the biological parameters (including $\bar{L}$), $p > 1$ and $k = h(\bar{L}) > 0$ such that (35) holds. We introduce the notation $X(k) := \frac{1}{p}(2 + \frac{k-d_L}{b_E}) < 1$ and use two parameters: $\iota \in (0, 1 - X(k))$ and $\zeta \in (0, 1)$, in order to represent the full range of (34). More precisely, we will parametrize $a$ and $\alpha$ with $\iota, \zeta$, as functions of $k$, and then we can go back to a function in (32) by letting $\lambda = \alpha\bar{L}$ and $h_m = k - a/(1 + \alpha^p)$.

We choose

$$\alpha_{\iota}(k) = \left( \frac{X(k) + \iota}{1 - (X(k) + \iota)} \right)^{1/p}$$

and

$$a_{\zeta, \iota}(k) = \zeta k(1 + \alpha_{\iota}^p) + (1 - \zeta)k \frac{(1 + \alpha_{\iota}^p)^2}{\alpha_{\iota}^p} X(k) = k \frac{1 - \zeta(1 - (X(k) + \iota))}{(X(k) + \iota)(1 - (X(k) + \iota))}.$$

For any choice of $\iota$ and $\zeta$, we end up with system (1), $h$ given by (32), featuring a unique, (locally linearly) unstable positive steady state. At least numerically, solutions always exhibit periodic oscillations, as can be seen in Figure 4 for egg dynamics.
<figure><img src="image_7.png" /><figcaption>Figure 4: Egg dynamics from (1) for $h$ of Hill function type. All parameters being fixed, including $p = 3$ and $k = 0.5$, $\iota = 0.05$ (top) or $\iota = 0.2$ (bottom) and $\zeta = 0.2$ (left) or $\zeta = 0.8$ (right).</figcaption></figure>

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)\frac{\beta_E \tau_L}{\delta_A} + \delta_L + \tau_L - \delta_A > k$. In this case we define

$$X(k) := \frac{\delta_A}{p\beta_E\tau_L}(2\frac{\beta_E}{\delta_A}\tau_L - \delta_L - \tau_L + \delta_A + k),$$

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 $(\iota, \zeta) \in (0, 1 - X(k)) \times (0, 1)$ and defining $\alpha_\iota(k)$ and $a_{\iota,\zeta}(k)$ by the same formulas as before. For all the numerical values we took for $\iota$ and $\zeta$, we always found oscillating solutions. Examples (dynamics of larvae and of $(E, L, A)$ in the three dimensional space) are shown in Figure 5.
<figure><img src="image_8.png" /><figcaption>Figure 5: Numerical solutions of (S3) for h defined by two different Hill functions. All parameters being fixed, including p = 3, k = 0.5, $\iota = \frac{1-X(k)}{10}$ and choosing $\zeta = 0.1$ (top) or $\zeta = 0.9$ (bottom).</figcaption></figure>

## 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 $\varepsilon \to 0$. We show here how to compute this amplitude analytically. To do so, we simply compute the local extrema of $u \mapsto \frac{\eta u^2}{h(\eta u)}$. The first-order necessary condition is $xh'(x) = 2h(x)$, where $x = \eta u$.

This provides with a general method to determine the limit trajectories. With the previous example from (32), $h(x) = h_m + a \frac{x^p}{(\alpha \bar{L})^p + x^p}$, this boils down to

$$2(h_m + a)x^{2p} + (\alpha \bar{L})^p((2-p)a + 4h_m)x^p + 2h_m(\alpha \bar{L})^{2p} = 0.$$

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

$$ \begin{aligned} \Delta &= (\alpha \bar{L})^{2p} ((2-p)a + 4h_m)^2 - 16h_m(h_m + a) \\ &= (\alpha \bar{L})^{2p} a((2-p)^2a - 8ph_m). \end{aligned} $$

Hence there are exactly two positive local extrema if and only if $(2-p)^2a > 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 \le 2$. Therefore the only case when there are two local extrema is when

$$ p > 2 \text{ and } \frac{h_m}{a} < \frac{(p-2)^2}{8p}. \qquad (36) $$

Under assumption (36) we find that the extrema ($y_M < y_m$) are located at

$$ (\alpha \bar{L})^p \frac{(p-2)a - 4h_m \pm \sqrt{a^2(p-2)^2 - 8aph_m}}{4(h_m + a)} $$

Let $\xi_{\pm} = (p-2)a - 4h_m \pm \sqrt{a^2(p-2)^2 - 8aph_m}$. With the notations of Lemma 3.4,

$$ u_m = \frac{\alpha}{\eta} \bar{L} \left( \frac{\xi_+}{4(h_m + a)} \right)^{1/p}, \quad \phi_m = \frac{\eta u_m^2}{h(\eta u_m)}, $$

$$ u_M = \frac{\alpha}{\eta} \bar{L} \left( \frac{\xi_-}{4(h_m + a)} \right)^{1/p}, \quad \phi_M = \frac{\eta u_M^2}{h(\eta u_M)}. $$

Then we can compute $u_r^0$ for $r \in \{m, M\}$ by solving $\frac{\eta(u_r^0)^2}{h(\eta u_r^0)} = \phi_r$. Unfortunately this cannot be done analytically. However, the amplitude of the oscillations in terms of $v$ is equal to

$$ A_v := \phi_M - \phi_m. $$

With $\bar{E} = 1/\epsilon$, we expect that the oscillations of $E$ have amplitude

$$ \frac{\phi_M - \phi_m}{\epsilon} = \frac{\alpha^2 \bar{L}^2}{\eta \epsilon} \left( \frac{\left(\frac{\xi_-}{4(h_m+a)}\right)^{2/p}}{h_m + a \frac{\xi_-}{h_m + \xi_-}} - \frac{\left(\frac{\xi_+}{4(h_m+a)}\right)^{2/p}}{h_m + a \frac{\xi_+}{h_m + \xi_+}} \right), $$

where $\eta = \frac{\bar{L}^2}{h(\bar{L})} = \frac{\bar{L}^2}{h_m + \frac{a}{1+\alpha^p}}$, by (22). Hence the amplitude of egg oscillations is equal to

$$ \frac{1}{\epsilon} A_v = \frac{1}{\epsilon} \frac{\alpha^2 (h_m + \frac{a}{1+\alpha^p})}{(4(h_m + a))^{2/p}} \left( \frac{\xi_-^{2/p}}{h_m + a \frac{\xi_-}{h_m + \xi_-}} - \frac{\xi_+^{2/p}}{h_m + a \frac{\xi_+}{h_m + \xi_+}} \right). $$

We can simplify this expression one step further by letting $\rho := h_m/a$. Then we notice that $q_{\pm} := \xi_{\pm}/a = p - 2 - 4\rho \pm \sqrt{(p-2)^2 - 8p\rho}$ and deduce

$$ A_v = \frac{\alpha^2}{1 + \alpha^p} \frac{1 + \rho + \alpha^p}{(4(1 + \rho))^{2/p}} \left( \frac{(\rho q_-)^{2/p}}{1 + \frac{\rho^2 q_-}{1 + \rho q_-}} - \frac{(\rho q_+)^{2/p}}{1 + \frac{\rho^2 q_+}{1 + \rho q_+}} \right). \qquad (37) $$

In particular we notice that the amplitude depends only on the function $h$ through $\rho$, $\alpha$ (hence $\bar{L}$) and $p$, and not on any other biological parameter, under the constraints (36).
An interesting case is when $p \to +\infty$, where $h$ approaches a step function from $h_m$ to $h_m + a$, with its jump located at $\alpha \bar{L}$. In this limit we can compute the amplitudes in $u$ and $v$:

$$ egin{cases} A_u = \frac{\alpha}{\bar{L}}(h_m + a\mathbb{1}_{\alpha<1} + \frac{a}{2}\delta_{\alpha=1})(\sqrt{\frac{\rho+1}{\rho}} - \sqrt{\frac{\rho}{\rho+1}}), \\ A_v = \alpha^2(h_m + a\mathbb{1}_{\alpha<1} + \frac{a}{2}\delta_{\alpha=1})\frac{1}{\rho(h_m+a)}. \end{cases} $$

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:

$$ \tau = \frac{2}{h_m} \log \left( \frac{h_m + a/2 - \alpha \bar{L}}{h_m + a/2 - \alpha \bar{L} \sqrt{\frac{\rho}{1+\rho}}} \right) + \frac{2}{h_m + a} \log \left( \frac{h_m + a/2 - \alpha \bar{L} \sqrt{\frac{1+\rho}{\rho}}}{h_m + a/2 - \alpha \bar{L}} \right) $$

Indeed, $h(u) = h_m$ if $u < \alpha \bar{L}$ and $h(u) = h_m + a$ if $u > \alpha \bar{L}$ so that $f(u, \phi(u)) = \eta_0 u(\xi - u)$ and $\phi'(u) = \frac{2\eta_0}{h_m}u$ if $u < \alpha \bar{L}$ and $\phi'(u) = \frac{2\eta_0}{h_m+a}u$ if $u > \alpha \bar{L}$.

## 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 $\bar{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 \ll d_L$, $d_E$ is fixed arbitrarily at $\frac{1}{180}$. We note that condition (4) is satisfied: $b_E = 20.94 > 0.15 + \frac{1}{180} = d_L + d_E$.

<figure><img src="image_9.png" /><figcaption>Figure 6: Time dynamics of eggs (left) and larvae (right) for $a_1 = 0.1$, $b_1^{0.05} = 2.91$.</figcaption></figure>

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_i^j)_{i=1,...,9; j \in J}$ by

$$ a_i = 0.1 + 0.05(i - 1) \text{ and } b_i^j = b_{i,min} + j \times b_{i,min}, $$
<figure><img src="image_10.png" /><figcaption>Figure 7: Time dynamics of eggs (left) and larvae (right) for $a_2 = 0.25$, $b_2^{0.5} = 4.18$.</figcaption></figure>

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:

$$h_i^j(L) = a_i \left( \arctan(b_i^j(L - \bar{L})) + \frac{\pi}{2} \right).$$

In our tests the steady state changes with $i$ (for example $\bar{E}_1 = 145.92$, $\bar{E}_4 = 59.59$ and $\bar{E}_9 = 30$) but we always have $\bar{L} = \bar{A} \frac{\delta_A}{\tau_L} = 1.13$.

<table><caption>Table 1: Steady states, period and amplitude of oscillations for $a = .1$</caption><thead><tr><th>a = 0.1</th><th>Period (days)</th><th>E&#773;</th><th>L&#773;</th><th>Larvae amplitude (%L&#773;)</th></tr></thead><tbody><tr><td>b = 2.91</td><td>5.18</td><td rowspan="4">145.92</td><td rowspan="4">1.13</td><td>23.1</td></tr><tr><td>b = 4.16</td><td>15.06</td><td>51.3</td></tr><tr><td>b = 8.32</td><td>61.54</td><td>110.22</td></tr><tr><td>b = 3.52</td><td>9.6</td><td>42.08</td></tr></tbody></table>

<table><caption>Table 2: Steady states, period and amplitude of oscillations for $a = .25$</caption><thead><tr><th>a = 0.25</th><th>Period (days)</th><th>E&#773;</th><th>L&#773;</th><th>Larvae amplitude (%L&#773;)</th></tr></thead><tbody><tr><td>b = 2.93</td><td>2.68</td><td rowspan="4">59.59</td><td rowspan="4">1.13</td><td>20.35</td></tr><tr><td>b = 4.18</td><td>6.96</td><td>50.67</td></tr><tr><td>b = 8.37</td><td>22.64</td><td>110.2</td></tr><tr><td>b = 4.99</td><td>9.98</td><td>63.03</td></tr></tbody></table>

We provide numerical results for $i \in \{1, 4, 9\}$ and $j \in \{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)) = (\bar{E}, \bar{L})$ (green) and $(E(0), L(0)) = (\bar{E}, \bar{L} + 0.02)$ (blue), which gives oscillations that appear to be periodic in time.
<figure><img src="image_11.png" /><figcaption>Figure 8: Time dynamics of eggs (left) and larvae (right) for $a_3 = 0.5$, $b_3^2 = 8.44$.</figcaption></figure>

<table><caption>Table 3: Steady states, period and amplitude of oscillations for $a = .5$.</caption><thead><tr><th>a = 0.5</th><th>Period (days)</th><th>E&#773;</th><th>L&#773;</th><th>Larvae amplitude (%L&#773;)</th></tr></thead><tbody><tr><td>b = 2.96</td><td>1.72</td><td rowspan="4">30</td><td rowspan="4">1.13</td><td>18.03</td></tr><tr><td>b = 4.22</td><td>3.8</td><td>49.8</td></tr><tr><td>b = 8.44</td><td>11.5</td><td>109.9</td></tr><tr><td>b = 7.6</td><td>10.1</td><td>99.43</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>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)$.</figcaption></figure>

Simulations are made with $a = a_1$ in Figure 6 (with $b = b_1^{0.05}$) and a time variable evaluated in [0, 100] days; $a = a_2$ in Figure 7 (with $b = b_2^{0.5}$) and a time variable evaluated in [0, 100] days; $a = a_3$ in
Figure 8 (with $b = b_3^2$) 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.

## References

[1] Maple 18. Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario.

[2] V. R. Aznar, M. S. D. Majo, S. Fischer, D. Francisco, M. A. Natiello, and H. G. Solari. A model for the development of *Aedes* (Stegomyia) *aegypti* as a function of the available food. *Journal of Theoretical Biology*, 365:311–324, 2015.

[3] V. R. Aznar, M. Otero, M. S. D. Majo, S. Fischer, and H. G. Solari. Modeling the complex hatching and development of *Aedes aegypti* in temperate climates. *Ecological Modelling*, 253:44–55, 2013.

[4] J. Bara, Z. Rapti, C. E. Cáceres, and E. J. Muturi. Effect of larval competition on extrinsic incubation period and vectorial capacity of *Aedes albopictus* for dengue virus. *PLoS ONE*, 10(5):1–18, 2015.

[5] S. Bhatt, P. W. Gething, O. J. Brady, J. P. Messina, A. W. Farlow, C. L. Moyes, J. M. Drake, J. S. Brownstein, A. G. Hoen, O. Sankoh, M. F. Myers, D. B. George, T. Jaenisch, G. R. W. Wint, C. P. Simmons, T. W. Scott, J. J. Farrar, and S. I. Hay. The global distribution and burden of dengue. *Nature*, 496(7446):504–507, 2013.

[6] O. Diekmann, J. Heesterbeek, and J. Metz. On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations. *J. Math. Biol.*, 28:365–382, 1990.

[7] G. L. C. Dutra, L. M. B. dos Santos, E. P. Caragata, J. B. L. Silva, D. A. M. Villela, R. Maciel-de Freitas, and L. Andrade Moreira. From Lab to Field: the influence of urban landscapes on the invasive potential of *Wolbachia* in Brazilian *Aedes aegypti* mosquitoes. *PLoS Neglect Trop D*, 9(4), 2015.

[8] J. Edgerly and M. Marvier. To hatch or not to hatch? Egg hatch response to larval density and to larval contact in a treehole mosquito. *Ecological entomology*, 17:28–32, 1992.

[9] B. Ermentrout. *Simulating, Analyzing, and Animating Dynamical Systems*. Society for Industrial and Applied Mathematics, 2002.

[10] J.-P. Françoise. *Oscillations en biologie, Analyse qualitative et modèle*. Springer, 2005.

[11] G. Guzzetta, F. Montarsi, F. A. Baldacchino, M. Metz, G. Capelli, A. Rizzoli, A. Pugliese, R. Rosà, P. Poletti, and S. Merler. Potential risk of dengue and chikungunya outbreaks in northern Italy
based on a population model of *Aedes albopictus* (Diptera: Culicidae). PLoS Neglect Trop D, 10(6):1-21, 06 2016.

[12] A. A. Hoffmann, I. Iturbe-Ormaetxe, A. G. Callahan, B. L. Phillips, K. Billington, J. K. Axford, B. Montgomery, A. P. Turley, and S. L. O'Neill. Stability of the *w*Mel *Wolbachia* infection following invasion into *Aedes aegypti* populations. PLoS Neglect Trop D, 8(9):1-9, 09 2014.

[13] N. Honorio, C. Codeço, F. Alves, M. Magalhães, and R. Lourenço-de Oliveira. Temporal distribution of *Aedes aegypti* in different districts of Rio De Janeiro, Brazil, measured by two types of traps. J Med Entomo, 46 (5):1001-1014, 2009.

[14] S. Juliano, R. G.S., R. Maciel-de Freitas, M. Castro, C. Codeço, R. Lourenço-de Oliveira, and L. Lounibos. She's a femme fatale: low-density larval development produces good disease vectors. Memórias do Instituto Oswaldo Cruz, 109(8):1070-1077, dec 2014.

[15] J. Koiller, M. A. da Silva, M. O. Souza, C. T. Codeço, A. Iggidr, and G. Sallet. Aedes, Wolbachia and dengue. Project-Team MASAIE, 2014.

[16] M. Legros, M. Otero, V. Romeo Aznar, H. Solari, F. Gould, and A. L. Lloyd. Comparison of two detailed models of *Aedes aegypti* population dynamics. Ecosphere, 7(10), 2016. e01515.

[17] T. P. Livdahl, R. K. Koenekoop, and S. G. Futterweit. The complex hatching response of *Aedes* eggs to larval density. Ecological Entomology, 9(4):437-442, 1984.

[18] J. Marsden and M. McCracken. *The Hopf Bifurcation and its Applications*, volume 19 of *Applied mathematical sciences*. Springer-Verlag, 1976.

[19] J. Meiss. *Differential Dynamical Systems*. SIAM, 2007.

[20] J. D. Murray. *Mathematical biology. I. An introduction*. Interdisciplinary applied mathematics. Springer, New York, 2002.

[21] B. Perthame. *Parabolic equations in biology*. Lecture Notes on Mathematical Modelling in the Life Sciences. Springer International Publishing, 2015.

[22] P. van den Driessche and J. Watmough. A simple SIS epidemic model with a backward bifurcation. J Math Biol, 40:525-540, 2000.

[23] D. Villela, C. Codeço, F. Figueiredo, G. Garcia, R. Maciel-de Freitas, and C. Struchiner. A Bayesian hierarchical model for estimation of abundance and spatial density of *Aedes aegypti*. PLoS ONE, 10(4), 2015. e0123794.

[24] H. Yang. Assessing the influence of quiescence eggs on the dynamics of mosquito *Aedes aegypti*. *Applied Mathematics*, 5:2696-2711, 2014.

[25] H. Yang, M. Macoris, K. Galvani, M. Andrighetti, and D. Wanderley. Assessing the effects of temperature on the population of *Aedes aegypti*, the vector of dengue. *Epidemiol Infect*, 137:1188-1202, 2009.

[26] H. L. Yeap, G. Rasic, N. M. Endersby-Harshman, S. F. Lee, E. Arguni, H. L. Nguyen, and A. A. Hoffmann. Mitochondrial DNA variants help monitor the dynamics of *Wolbachia* invasion into host populations. *Heredity*, 116(3):265-276, 2016.