Minimum Vector Control Intensity to Get a Stable Fixed Point in a Mosquito Dynamic Model

ABSTRACT

Vector-borne diseases are a cause of concern all around the world, especially in Brazil. In the past few years, the Brazilian health system faced recurrent epidemics such as Dengue and Malaria as well as new cases of Chikungunya, Zika and Yellow Fever. Vector control continues to be one of the most important counter measures against these types of diseases. Mathematical models are important tools for planning vector control strategies. In this work we present an approach in order to calculate what is the minimum vector control intensity to obtain stability in a simple population’s dynamics model of mosquitoes. We combined numerical simulations with analytic results. Transcritical bifurcations appear in our analysis considering different control’s parameters values for the eggs, larvae, pupae and adults mosquitoes populations. A discussion about combined strategies of vector control was also showed.

Keywords:
modelling; population dynamic; stability; simulation

1 INTRODUCTION

The mosquito Aedes Aegypti is the vector for the transmission of diseases such as Dengue Fever, Chikungunya, Zika and Yellow Fever. In recent years, there has been a great increase in the number of cases of those diseases ^{4) (9) (11) (14}, considered some of the major challenges to health systems in several countries, especially Brazil. This situation occurs in large part due to disordered population growth, since the mosquito develops in places where there is a constant accumulation of water, often due to large concentrations of people ³. The vector control is one of the most important strategies for helping control the mosquito population.

There are a lot of papers in the literature bringing mathematical model for those diseases ^{1) (2) (6) (8) (10) (12) (15}. In this work, we constructed simple models of mosquito population dynamics, based on those previous works, reaching to one version with different vector control mechanisms to be applied isolated and in a combined strategy. Based on qualitative studies of differential equations system (stability analysis), and numerical simulations, our results classify the linear stability of the system’s equilibrium points and reveal the conditions for which the population of mosquitoes can grow indefinitely, die out or converge to a certain number. In a simple but general/reproducible approach, we illustrate the process of quantification of the minimum vector control intensity to keep the mosquito population stable.

2 MODELS AND METHODS

The Aedes aegypti mosquito life cycle has four stages: egg, larva, pulp and adult (mosquito itself). In the adult phase the mosquito is able to transmit a virus to the human during a bite, in order to get blood to develop its eggs. In this section we show the model construction and the computational and mathematical tools for solving and analyzing the model and the solutions.

2.1 Models

Let E(t), L(t), P(t), M(t) be the functions that represent the number of Aedes aegypti’s eggs, larvae, pupae and mosquitoes’ population respectively in a instant t of time.

Consider the following models:

$$\begin{gathered} \left\{\begin{array}{l}\frac{dE}{dt}={\varphi }_m.M\left(t\right)-{\sigma }_e.E\left(t\right) \\ \frac{dL}{dt}={\sigma }_e.E\left(t\right)-{\sigma }_l.L\left(t\right) \\ \frac{dP}{dt}={\sigma }_l.L\left(t\right)-{\sigma }_p.P\left(t\right) \\ \frac{dM}{dt}={\sigma }_p.P\left(t\right)-{\mu }_m.M\left(t\right)\end{array}\right. \end{gathered}$$ (2.1)

$$\begin{gathered} \left\{\begin{array}{l}\frac{dE}{dt}={\varphi }_m.M\left(t\right)-{\sigma }_e.E\left(t\right) \\ \frac{dL}{dt}={\sigma }_e.E\left(t\right)-{\sigma }_l.L\left(t\right) \\ \frac{dP}{dt}={\sigma }_l.L\left(t\right)-{\sigma }_p.P\left(t\right) \\ \frac{dM}{dt}={\sigma }_p.P\left(t\right)\left(1-\frac{M\left(t\right)}{k_m}\right)-{\mu }_m.M\left(t\right)\end{array}\right. \end{gathered}$$ (2.2)

$$\begin{gathered} \left\{\begin{array}{l}\frac{dE}{dt}={\varphi }_m.M\left(t\right)-({\sigma }_e+{\mu }_e).E\left(t\right) \\ \frac{dL}{dt}={\sigma }_e.E\left(t\right)-({\sigma }_l+{\mu }_l).L\left(t\right) \\ \frac{dP}{dt}={\sigma }_l.L\left(t\right)-({\sigma }_p+{\mu }_p).P\left(t\right) \\ \frac{dM}{dt}={\sigma }_p.P\left(t\right)\left(1-\frac{M\left(t\right)}{k_m}\right)-{\mu }_m.M\left(t\right)\end{array}\right. \end{gathered}$$ (2.3)

Note that (2.1) is a linear ODE system therefore, it’s expected that it will have only one equilibrium point its analysis and results are expected to be quite simple and not realistically possible. The system in itself is fairly coupled: the entries of the compartments, besides the first, are the transitions from the previous one. The σ parameter represents this transition. The entry of the first compartment is the eggs that each mosquito puts (represented by the parameter φ _m ) and the sink of the last one is the death of its individuals (µ _m ) thus it’s reasonable for one to think that these two parameters will be pivotal to the system’s stability. Meanwhile, (2.2) is a nonlinear ODE system due to the logistical term on the last equation which will cause the appearance of a second equilibrium point, and while the system’s overall stability is expected to depend on the same parameters as the first model because of their resemblance, this second equilibrium point will cause a different behaviour as it will be seen seen in (2.1). Last but not the least, we have the (2.3) system which will be the main focus of this work due to the vector control presented in each one oh the its compartments (represented by the µ parameters). Its system is also nonlinear due to the logistical term but it does not resemble the previous ones since it’s not coupled. So , its stability analysis is expected to depend of much more parameters when compared to (2.1) and (2.2).

In Table 1, the ”-” in the Value column represents that the parameter will not have a fixed numerical value, i.e we will vary it in different simulations.

Table 1
Summary of the parameters.

Symbol	Parameter	Unit	Domain	Value
σ e	Per-capita transition rate, from egg to larva	day −1	[0,1]	$\frac{1}{5}$
σ l	Per-capita transition rate, from larva to pupa	day −1	[0,1]	$\frac{1}{7}$
σ p	Per-capita transition rate, from pupa to mosquito	day −1	[0,1]	$\frac{1}{3}$
µ e	Per-capita mortality rate of eggs	day −1	[0,1]	-
µ l	Per-capita mortality rate of larvae	day −1	[0,1]	-
µ p	Per-capita mortality rate of pupae	day −1	[0,1]	-
µ m	Per-capita mortality rate of mosquitoes	day −1	[0,1]	$\frac{1}{7}$
φ m	Oviposition rate per mosquito	day −1	(0,∞)	-
k m	Carrying capacity of the environment	-	(0,∞)	100

2.2 Methods

To study the linear stability of the equilibrium points (where all the equations are equal to zero) of a nonlinear ODE system, we can restrict ourselves to a local analysis around the fixed point by linearization. Let X’ = F(X ) be a autonomous system where F(X ) = (F ₁(X ), ..., F _n (X )) and X = (x ₁ , . . . , x _n ), P ∈ ${\mathrm{ℝ}}^n$ an equilibrium point. We call P a stable equilibrium point if for a given ε > 0, there exists δ > 0 such as the solution X = f (t), satisfying ∥ f (0) − P∥ < δ implies that ∥ f (0) − P∥ < ε. If P is not stable then we call it a unstable equilibrium point.

The system has a linear approximation X’ = J(P)· X , where J(P) is the system’s jacobian matrix calculated at P:

$J\left(P\right)=\frac{\partial (F_1,\dots ,\partial F_n)}{\partial (x_1,\dots ,x_n)}=\left[\begin{array}{ccc}\frac{\partial F_1}{\partial x_1}\left(P\right)& \dots & \frac{\partial F_1}{\partial x_n}\left(P\right)\\ ⋮& & ⋮\\ \frac{\partial F_n}{\partial x_1}\left(P\right)& \dots & \frac{\partial F_n}{\partial x_n}\left(P\right)\end{array}\right]$ (2.4)

Thus, the eigenvalues of (2.4) are calculated and we conclude that if none of them has a positive real part, then the equilibrium point P is locally stable. If at least one has a positive real part, then the equilibrium point is unstable. All of the affirmations above are formally demonstrated at ^{7), (5}.

As for the system’s bifurcations, they are determined by calculating the eigenvalues of the jacobian matrix for a given set of parameters while varying one specific parameter in its domain. For our models, there are two cases for a jacobian matrix calculated at a fixed point: all of its eigen- values have a negative real part or there is one with a positive real part. Besides that we found out that the equilibrium points, when there are more than one, have opposite stability: when the first is stable the other is unstable and vice-versa. So, the bifurcation happens at the value in which a root of a fixed point’s jacobian matrix changes its signal, resulting in a exchange of stability between the equilibrium points.

A fourth order Runge Kutta algorithm was applied for numerical solving and the MATLAB software for spectral analysis.

3 RESULTS AND DISCUSSION

In this section it is shown the stability analysis for each model, reaching to a version with different vector control mechanisms, to be applied isolated and in a combined strategy.

3.1 Model 1

The linear system (2.1) possess a single equilibrium point and it’s the trivial one: 0 = (E, L, P, M) = (0, 0, 0, 0). For any given point, its Jacobian matrix is:

$J=\left[\begin{array}{cccc}-{\sigma }_e& 0& 0& {\varphi }_m\\ {\sigma }_e& -{\sigma }_l& 0& 0\\ 0& {\sigma }_l& -{\sigma }_p& 0\\ 0& 0& {\sigma }_p& -{\mu }_m\end{array}\right]$ (3.1)

whose characteristic polynomial is

$p\left(\lambda \right)=(\lambda +{\mu }_m)(\lambda +{\sigma }_e)(\lambda +{\sigma }_l)(\lambda +{\sigma }_p)-{\varphi }_m{\sigma }_e{\sigma }_l{\sigma }_p$ (3.2)

The signal of its roots depends of two parameters: φ _m and µ _m . If the second is greater than the first (µ _m > φ _m ), all the eigenvalues of (3.1) have a negative real part (two complex and two real), therefore the equilibrium point is stable. Otherwise, one real eigenvalue has a positive real part and the remaining ones have a negative real part (a complex and its complex conjugate, and a real number), therefore in this case the equilibrium point is unstable. This specific behaviour is due to the fact that there exists a bifurcation in φ _m = µm as illustrated by the following diagram:

Figure 1
Bifurcation diagram with the parameter values shown in the beginning of the section.

The vertical axis represents the fourth coordinate of the equilibrium point and the horizontal axis, the values of the φ _m parameter. The highlighted curve is the equilibrium point’s fourth coordinate. A arrow pointing upwards (away from the curve) indicates that the equilibrium point is unstable and the one pointing downwards (to the curve) indicates it’s stable. It’s notable that the bifurcation occurs when φ _m =µ _m .

Since there’s only one equilibrium point, if it’s unstable, then the population of the mosquitoes (and subsequently the eggs, larvae and pupae’s population) grows indefinitely. Otherwise all the populations converge to zero for a sufficiently great value of t.

3.2 Model 2

The nonlinear system (2.2) possess two equilibrium points: the trivial one, which we will refer as 0 from now on and $P^{*}=(E,L,P,M)=k_m({\varphi }_m-{\mu }_m)\left(\frac{1}{{\sigma }_e},\frac{1}{{\sigma }_l},\frac{1}{{\sigma }_p},\frac{1}{{\varphi }_m}\right)$. Its jacobian matrix is

$J=\left[\begin{array}{cccc}-{\sigma }_e& 0& 0& {\varphi }_m\\ {\sigma }_e& -{\sigma }_l& 0& 0\\ 0& {\sigma }_l& -{\sigma }_p& 0\\ 0& 0& {\sigma }_p-\frac{{\sigma }_{p}M\left(t\right)}{k_m}& -{\mu }_m-\frac{{\sigma }_{p}P\left(t\right)}{k_m}\end{array}\right]$ (3.3)

Calculating the matrix (3.3) at 0 gives us the exact same matrix as the one of the first model, therefore it has the same characteristic polynomial as (3.2). Thus, the conclusion about this equilibrium point remains the same: if φ _m > µ _m , it’s unstable, otherwise it’s stable. However the logistical term that differs this model from the first one and makes this system nonlinear prevents the populations from growing indefinitely because of the nontrivial equilibrium point’s appearance.

Calculating the matrix (3.3) at P ^∗ gives us:

$J\left(P^{*}\right)=\left[\begin{array}{cccc}-{\sigma }_e& 0& 0& {\varphi }_m\\ {\sigma }_e& -{\sigma }_l& 0& 0\\ 0& {\sigma }_l& -{\sigma }_p& 0\\ 0& 0& {\sigma }_p\frac{{\mu }_m}{{\varphi }_m}& -{\varphi }_m\end{array}\right]$

whose characteristic polynomial is

$p\left(\lambda \right)=(\lambda +{\varphi }_m)(\lambda +{\sigma }_e)(\lambda +{\sigma }_l)(\lambda +{\sigma }_p)-{\mu }_m{\sigma }_e{\sigma }_l{\sigma }_p$

The signal of its roots depends on φ _m and µ _m but their dependence is the exact opposite of the showed in (3.2): if φ _m > µ _m then every eigenvalues has a negative real part, therefore the nontrivial equilibrium point is stable. If φ _m < µ _m , one eigenvalue has a positive real part, therefore the point is unstable.

Figure 2
Bifurcation diagram with the parameter values shown in the beginning of the section.

Note that, a bifurcation occurs when φ _m = µ _m , since if φ _m < µ _m the trivial equilibrium point is stable and the nontrivial one is unstable and if φ _m < µ _m both points change their respective stability. By ¹³, we say that an exchange of stability between both equilibrium points takes place at the bifurcation point which constitutes a transcritical bifurcation.

In this model in particular, the population of mosquitoes will not grow indefinitely, since it will converge to one of the equilibrium points: if φ _m < µ _m it will eventually die out and if φ _m > µ _m , it will stabilize at the nontrivial fixed point.

There is a very important observation here: the typically values of φ _m are really greater than µ _m . This mathematical relation shows us the importance of including vector control in this model.

All the simulations will use (E, L, P, M) = (1, 0, 0, 0) as the initial point from now on. The simulations run show exactly what happens if we set a smaller or a bigger value of φ _m than the bifurcation point. Using the parameters showed in Table 1 and setting φ _m as 10 in Figure 3a and as $\frac{1}{25}$ in Figure 3b note that the bifurcation point is ${\mu }_m=\frac{1}{7}$ and we can see that in the first case, the population of mosquitoes does not converge to zero and in case of Figure 3b, it dies out.

Figure 3
Temporal series of the mosquito population

3.3 Model 3

The nonlinear system (2.3) possess two equilibrium points: the trivial one and

$P`*=(E,L,P,M)=k_m(\alpha -\beta )\left(\frac{{\varphi }_m}{\alpha ({\sigma }_e+{\mu }_e)},\frac{{\varphi }_m{\sigma }_e}{\alpha ({\sigma }_e+{\mu }_e)({\sigma }_l+{\mu }_l)},\frac{{\mu }_m}{\beta {\sigma }_p},\frac{1}{\alpha }\right)$ (3.4)

where α = σ _e σ _l σ _p φ _m and β = µ _m (µ _e + σ _e )(µ _l + σ _l )(µ _p + σ _p ). Its jacobian matrix is:

$J=\left[\begin{array}{cccc}-({\sigma }_e+{\mu }_e)& 0& 0& {\varphi }_m\\ {\sigma }_e& -({\sigma }_l+{\mu }_l)& 0& 0\\ 0& {\sigma }_l& -({\sigma }_p+{\mu }_p)& 0\\ 0& 0& {\sigma }_p-\frac{{\sigma }_{p}M\left(t\right)}{k_m}& -{\mu }_m-\frac{{\sigma }_{p}P\left(t\right)}{k_m}\end{array}\right]$ (3.5)

Applying J at 0 gives us:

$J\left(\mathit{0}\right)=\left[\begin{array}{cccc}-({\sigma }_e+{\mu }_e)& 0& 0& {\varphi }_m\\ {\sigma }_e& -({\sigma }_l+{\mu }_l)& 0& 0\\ 0& {\sigma }_l& -({\sigma }_p+{\mu }_p)& 0\\ 0& 0& {\sigma }_p& -{\mu }_m\end{array}\right]$ (3.6)

whose characteristic polynomial is

$p_1\left(\lambda \right)=(\lambda +({\sigma }_e+\mu e\left)\right)(\lambda +({\sigma }_l+\mu l\left)\right)(\lambda +({\sigma }_p+\mu p\left)\right)(\lambda +{\mu }_m)-\alpha$

Its roots depend not only on the oviposition and mortality parameters but also on the other parameters, more specifically, α and β . If α < β each eigenvalue has a negative real part, while if α > β there exists one with a positive real part. And so, we conclude that in the first case, the trivial equilibrium point is stable and in the second one, it’s unstable.

Applying the matrix (3.5) at P ^∗ gives us:

$J\left(P^{*}\right)=\left[\begin{array}{cccc}-({\sigma }_e+{\mu }_e)& 0& 0& {\varphi }_m\\ {\sigma }_e& -({\sigma }_l+{\mu }_l)& 0& 0\\ 0& {\sigma }_l& -({\sigma }_p+{\mu }_p)& 0\\ 0& 0& {\sigma }_p\frac{\beta }{\alpha }& -{\mu }_m\frac{\alpha }{\beta }\end{array}\right]$ (3.7)

whose characteristic polynomial is:

$p_2\left(\lambda \right)=(\lambda +({\sigma }_e+{\mu }_e\left)\right)(\lambda +({\sigma }_l+{\mu }_l\left)\right)(\lambda +({\sigma }_p+{\mu }_p\left)\right)\left(\lambda +\frac{{\mu }_m\alpha }{\beta }\right)-\beta$

The stability of this equilibrium point also depends of α and β and it is the opposite of the trivial point: if α > β it will be stable and if α < β it will be unstable.

Now, analyzing each control parameter (µ _e , µ _l , µ _p ) separated, i.e setting different values for µ _e and the other two being zero, and repeating this process for each µ _i , i ∈ e, l, p we get the bifurcation diagrams seen in Figure 4.

Figure 4
Bifurcation diagrams for each control parameter.

We can calculate the exact point of bifurcation by using the relation α = β . Therefore, using the parameters values of Table 1 and setting ${\varphi }_m=\frac{1}{2}$, the bifurcation point seen in Figure 4a is ${\mu }_e=\frac{1}{2}$, in Figure 4b is ${\mu }_l=\frac{5}{14}$ and in Figure 4c is ${\mu }_p=\frac{5}{6}$.

Simulations were run to illustrate the conclusions. In Figure 5a it was used µ _e = 0.30 and in Figure 5b, µ _e = 0.55. In Figure 6a, we set µ _l = 0.30 and in Figure 6b, µ _l = 0.45. In Figure 7a, we used µ _p = 0.75 and in Figure 7bµ _p = 0.95. Finally, using the values set in the simulations that the population does not converge to 0 at the same time, i.e µ _e = 0.30, µ _l = 0.30, µ _p = 0.75, Figure 8 shows that those strategies combined cause the population of mosquitoes to die out.

Figure 5
Temporal series for the mosquito population considering only the egg’s control parameter. (a) µ _e = 0.30. (b) µ _e = 0.55.

Figure 6
Temporal series for the mosquito population considering only the larvae’s control parameter. (a) µ _l = 0.30. (b) µ _l = 0.45.

Figure 7
Temporal series for the mosquito population considering only the pupae’s control parameter. (a) µ _p = 0.75. (b) µ _p = 0.95.

Figure 8
Temporal series for the mosquito population considering control parameters in all populations.

4 CONCLUSIONS AND PERSPECTIVES

Our results, based on analogical stability analysis and numerical simulations, highlighted the importance of including a carrying capacity quantity (transition between model one to model two) and a vector control intensity (transition between model two to model three) in this kind of modeling. For the non-trivial models (second and third ones), typically values of φ _m are really greater than µ _m . The vector control allows us to get stable non-trivial solutions even for more realistic values of oviposition φ _m .

Singly, the vector control acting on the different stages of mosquito population demands high intensity, so much effort. However, when we combine all vector control strategies, considering the population of eggs (mechanic or chemistry), larvae (mechanic or chemistry) and pupae (mechanic), it is possible to get stable non-trivial solutions even for some high values of oviposition φ _m even with none control on the mosquito population.

For the following works we are going to perform our approach to a set of more refined models in the literature, comparing among different models the minimum intensity of vector control to get stable non-trivial solutions.

Acknowledgments

This work is partially supported by the project FAPESP 2015/50122-0 and DFG-IRTG-1740/2, and CNPq 139180/2018-4 and 420338/2018-7.

REFERENCES

• ¹
  G. Adamu, M. Bawa, M. Jiya & U. Chado. A Mathematical Model for the Dynamics of Zika Virus via Homotopy Perturbation Method. Journal of Applied Sciences and Environmental Management, 21 (2017), 615.
• ²
  F.B. Agusto, S. Bewick & W.F. Fagan. Mathematical model of Zika virus with vertical transmission. Infectious Disease Modelling, 2 (2017), 244-267.
• ³
  L. Barsante, R. Cardoso, J. Luiz Acebal, C. Silva & A. Eiras. Modelo Entomológico Determinístico sob Efeito da Pluviosidade para o Aedes aegypti e o Aedes albopictus. TEMA: Tendências em Matemática Aplicada e Computacional, 19 (2018), 289-303.
• ⁴
  C.T. Codeço, D.A. Villela & F.C. Coelho. Estimating the effective reproduction number of dengue considering temperature-dependent generation intervals. Epidemics, 25 (2018), 101-111.
• ⁵
  R.L. Devaney. “An Introduction to Chaotic Dynamical System”. Westview Press (2003).
• ⁶
  L.Z. Goldani. Yellow fever outbreak in Brazil, 2017. Brazilian Society of Infectious Diseases, 21(2) (2017), 123-124.
• ⁷
  D.L. Kreider, R.G. Kuller & D.R. Ostberg. “Elementary Differential Equations”. Edgard Blucher (1972).
• ⁸
  R.M. Lana, M.M. Morais, T.F.M.d. Lima, T.G.d.S. Carneiro, L.M. Stolerman, J.P.C. dos Santos, J.J.C. Cortés, A.E. Eiras & C.T. Codeço. Assessment of a trap based Aedes aegypti surveillance program using mathematical modeling. PLOS One, 13 (2018), 1-16.
• ⁹
  R. Lowe, C. Barcellos, P. Brasil, O.G. Cruz, N.A. Honório, H. Kuper & M.S. Carvalho. The Zika Virus Epidemic in Brazil: From Discovery to Future Implications. Int J Environ Res Public Health, 15 (2018), 1-18.
• ¹⁰
  G. Phaijoo & D. Gurung. Mathematical Model of Dengue Disease Transmission Dynamics with Control Measures. Journal of Advances in Mathematics and Computer Science, 23 (2017), 1-12.
• ¹¹
  S.M. Raimundo, M. Amaku & E. Massad. Equilibrium Analysis of a Yellow Fever Dynamical Model with Vaccination. Computational and Mathematical Methods in Medicine, 2015 (2015).
• ¹²
  L.B.L. Santos, S.T.R. Pinho, R.F.S. Andrade & et al. Periodic forcing in a three-level cellular automata model for a vector-transmitted disease. Physical Review E, 80 (2009).
• ¹³
  S.H. Strogatz. “Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry, and engineering”. Addison-Wesley Publishing Company (1994).
• ¹⁴
  B. Wahid, A. Ali, S. Rafique & M. Idrees. Global expansion of chikungunya virus: mapping the 64-year history. International Journal of Infectious Diseases, 58 (2017), 69-76.
• ¹⁵
  H.M. Yang, J.L. Boldrini, A.C. Fassoni & et al. Fitting the Incidence Data from the City of Campinas, Brazil, Based on Dengue Transmission Modellings Considering Time-Dependent Entomological Parameters. PLOS One, (2016).

Publication Dates

• Publication in this collection
  28 July 2023
• Date of issue
  Jul-Sep 2023

History

• Received
  10 Dec 2018
• Accepted
  27 Jan 2023