A Model for Aedes aegypti Infestation According to Meteorological Variables: Case of Caratinga (Minas Gerais - Brazil)

ABSTRACT

Aedes aegypti mosquitoes are the vector of diseases such as dengue, zika, chikungunya, and yellow fever among others. All the stages of development, eggs, larvae, pupa, and the adult of the species have its population modulated by meteorological variables, such as precipitation and temperature through affecting the productivity of breeding sites, metabolic processes, and others. Since adult females are responsible for transmitting the virus, the population of females becomes a direct indicator of the risk of infection. For this reason, some ongoing vector surveillance programs are based on adult female capture. In turn, all the stages of development have its population modulated by meteorological variables, such as precipitation and temperature, through productivity of breeding sites, metabolic processes and others. In this work, field data of capture of females was used to evaluate if a population dynamics model of Aedes aegypti under the effect of weather would be able to forecast field population. The nonlinear dynamic system model comprises: (1) four equations for the populations of the stages of development of the mosquito, designed for the ongoing surveillance program; (2) parametric dependencies of the rates of development on mean temperature and weekly accumulated precipitation. The dependencies on temperature and precipitation are modelled with aim of simplicity with the fewer number of parameters as possible. Temperature dependence is modelled based on values of the related literature under the assumption of existence of a optimum temperature for the rates, getting worse for extreme temperatures. The dependence on precipitation which is barely treated in experiments is modelled under the assumption of a monotonic dependence described by a power law with values estimated in orders of magnitude from data in the literature. By comparison with field data of an entomological indicator based on the number of Ae. aegypti females captured by a public health program in the city of Caratinga (Minas Gerais, Brazil), the model showed a significant correlation (R = 0.75). The result shows that the approach, if refined, can provide forecasting for of the population size.

Keywords: Aedes aegypti; meteorological variables; dynamical systems

RESUMO

Os mosquitos Aedes aegypti são vetores de doenças como dengue, zika, chikungunya, febre amarela, entre outras. Como as fêmeas adultas são responsáveis pela transmissão do vírus, sua população torna-se um indicador direto do risco de infecção. Por esse motivo, alguns programas de vigilância de vetores baseiam-se na captura de fêmeas adultas. Por sua vez, todas as fases de desenvolvimento (ovos, larvas, pupas e o vetor adulto) têm sua população modulada por variáveis meteorológicas, como precipitação e temperatura, através da produtividade dos criadouros, processos metabólicos e outros. Neste trabalho, dados de captura de fêmeas em campo foram utilizados para avaliar se um modelo de dinâmica populacional de Aedes aegypti , sob o efeito de variáveis meteorológicas, seria capaz de prever a população de campo. O modelo de sistema dinâmico não linear compreende: quatro equações para as populações dos estágios de desenvolvimento do mosquito, projetadas para um programa de vigilância em andamento; (2) dependências paramétricas das taxas de desenvolvimento em relação a` temperatura média e precipitação acumulada semanalmente. As dependências de temperatura e precipitação são modeladas, por simplicidade, com o menor número de parâmetros possível. A dependência da temperatura é modelada com base em valores da literatura relacionada sob o pressuposto da existência de uma temperatura ótima para as taxas, piorando para temperaturas extremas. A dependência da precipitação, que mal é tratada em experimentos, é modelada sob o pressuposto de uma dependência monotônica descrita por uma lei de potência com valores estimados em ordens de magnitude a partir de dados da literatura. Por comparação com dados de campo de um indicador entomológico baseado no número de fêmeas Ae. aegypti capturadas por um programa de saúde pública na cidade de Caratinga (Minas Gerais, Brasil), o modelo apresentou correlação significativa (R = 0,75). O resultado mostra que a abordagem, se refinada, pode fornecer previsões para o tamanho da população.

Palavras-chave: Aedes aegypti; variáveis meteorológicas; sistemas dinâmicos

1 INTRODUCTION

The Aedes (Stegomyia) aegypti (hereafter referred to as Ae. aegypti ), a species of the Culicidae family of mosquitoes, is the main vector of arboviruses like dengue fever, yellow fever, zika, and chikungunya ²⁷ . The species is found mainly in tropical and subtropical regions of the world and its distribution is closely related to climate patterns determined by variables such as temperature ^{1),(11),(21),(30} , precipitation ¹⁴ , humidity ⁵ , and winds ^9),(22 . Those factors act in the productivity of breeding sites, metabolic processes, behaviour, and other characteristics ¹⁸ .

The development of Ae. aegypti takes place by complete metamorphosis ⁸ . The stages of development occur in two distinct environments: (a) aquatic, composed by the egg, larva and pupa stages, and (b) aerial with the winged stage of the vector ¹⁵ . To obtain nutrients for maturation and development of eggs, and energy reserves, adult gravid females perform blood feeding almost exclusively on humans ¹⁹ . Transmission of arboviruses to humans takes place when a female infected by the virus bites a susceptible human. Among the arboviruses, there are dengue fever, zika, chikungunya and yellow fever and others. According to the World Health Organization (WHO), just considering dengue fever, there occur 50 to 100 million new infections per year in more than 100 countries ²⁸ .

To avoid social impacts with epidemics of diseases transmitted by Ae. aegypti , the approach recommended by the WHO are the actions of entomological surveillance and control ²⁷ . The MIAEDES (Intelligent Aedes Monitoring), a system that make use of an adult mosquito trap called MosquiTRAP^®. The MosquiTRAP^® captures gravid females for monitoring the population of Ae. aegypti females.

In the surveilled municipalities, the traps are installed in the study area, spaced from each other by 200-250 meters, and monitored weekly. The surveillance data is summarised in terms of an entomological indicator MFAI (Mean Females Aedes Index), which is defined as the ratio of the total number of females caught during an epidemiological week to the number of positive traps inspected.

There is a long and well-known connection between meteorological and climatic factors and the rates of development stages of Ae. aegypti^7),(10 . Such a connection produces a seasonal population pattern for Ae. aegypti , and consequently, for the number of cases of the infections. High precipitations and high temperatures increase the reproduction and survival of Ae. aegypti , markedly increasing the mosquito population size, and, consequently, the risk of infections ^{3), (10} . The knowledge of the influence of the meteorological factors is important for the effectiveness of the environmental control. However, the adaptability of Ae. aegypti to urban environments and to human way of life, together with social conditions, allow the mosquito to survive throughout the year ^{3), (12} .

Mathematical and/or statistical models that consider the associations between meteorological variables and Ae. aegypti life cycle are developed with the purpose of describing the population dynamics of the vector, as well as the epidemics caused by that ^{1), (2), (7), (17), (24), (30} .

In this work, a mathematical model based on differential equations is used to describe the Ae. aegypti life cycle as a function of meteorological data (temperature and precipitation). The main goal is evaluate if a simple model under the effect of weather would be able to forecast field population. Hence it was presented a nonlinear dynamic system model that comprises: (1) four equations for the populations of the stages of development of the mosquito, designed for the ongoing surveillance program; (2) parametric dependencies of the rates of development on mean temperature and weekly accumulated precipitation.

The dependencies of the entomological rates on temperature and precipitation are modelled with aim of simplicity with the fewer number of parameters as possible. Temperature dependence is modelled based on values of the related literature with polynomials of degree 2, under the assumption of existence of a optimum temperature for the rates, and getting worse for extreme temperatures ^{1), (21), (30} . The dependence on precipitation, which is barely found in experiments, is modelled under the assumption of a monotonic dependence through a power law ²⁵ . The power law constants are estimated by orders of magnitude from the literature that refers to entomological rates, regardless if it refers to precipitation itself ^{1), (7), (17), (24), (29), (30} . Although the determination of those parameters from experiments is of limited realism, it is intend to verify the forecasting power of the model as a evaluation of that approach for eventual future refinement. For that purpose, the results of the computational model implementation will be compared to experimental data from Aedes females (MFAI) from the municipality of Caratinga, Minas Gerais, Brazil.

In Section 2, the entomological data of captures and meteorological data are characterised. In Section 3, the mathematical model for the life cycle of Aedes aegypti is detailed, considering the dependence with the meteorological variables. A description of the simulation and validation of the results is shown in Section 4. The final considerations are made in Section 5.

2 DATA BASE

The data used in this work refers to the epidemiological weeks 23-52 of 2009 and epidemiological weeks 1-51 of 2010 in the city of Caratinga (Minas Gerais, Brazil). The climate of the city is classified as tropical semi-humid. Average temperatures are around 22^◦C, with a mean average of 16 . 7^◦C and average maximum of 27 . 7^◦C (averages considering the period from 1961 to 2014). The average annual rainfall is 1 , 122 . 6mm ¹³ .

The weekly meteorological data, accumulated precipitation p (mm) and average temperature T (^◦C) for the study period were obtained through the MDBTR (Meteorological Data Bank for Teaching and Research) on the website of the National Meteorological Institute INMET (www.inmet.gov.br).

Experimental data from the Ae. aegypti females population were obtained by using MosquiTRAP® in a public heath surveillance from the study area. During the period of study of 96 epidemiological weeks, the MFAI was generated in the MI-Aedes program developed in the municipality by the biotechnology company Ecovec S. A. (Belo Horizonte, MG, Brazil).

3 MATHEMATICAL MODEL

To evaluate the viability of a predictive population model with dependencies on mean temperature and precipitation, it is presented a nonlinear dynamic system model that comprises: (1) four equations for the populations of the stages of development of the mosquito, designed for the ongoing surveillance program; (2) parametric dependencies of the rates of development on mean temperature and weekly accumulated precipitation. The population model (3.1) written as a system of differential equations is proposed to describe the dynamics of the population size of the four stages of development Ae. aegypti , considering the weekly variation of precipitation and average temperature. The population of eggs, aquatic forms (larvae and pupae), pre-bloodmeal and post-bloodmeal females are represented by E ( t ), A ( t ), F ₁( t ) e F ₂( t ), respectively (see Figure 1 ).

Figure 1:
Diagram of the model populations related to the life cycle of Aedes aegypti with rates of development as parametric dependent on temperature T and precipitation p .

$\left\{\begin{array}{l}\begin{array}{l}\frac{dE}{dt}=\varphi \left(p,T\right)\left(1-\frac{E\left(t\right)}{C\left(p,T\right)}\right)F_2\left(t\right)-{\alpha }_1\left(p,T\right)E\left(t\right)-{\mu }_E\left(p,T\right)E\left(t\right),\\ \frac{dA}{dt}={\alpha }_1\left(p,T\right)E\left(t\right)-{\alpha }_2\left(p,T\right)A\left(t\right)-{\mu }_A\left(p,T\right)A\left(t\right),\\ \frac{d{F}_1}{dt}={\alpha }_2\left(p,T\right)A\left(t\right)-{\alpha }_3\left(p,T\right)F_1\left(t\right)-{\mu }_{F_1}\left(p,T\right)F_1\left(t\right),\end{array}\\ \frac{d{F}_2}{dt}={\alpha }_3\left(p,T\right)F_1\left(t\right)--{\mu }_{F_2}\left(p,T\right)F_2\left(t\right),\\ C,\varphi ,{\alpha }_1,{\alpha }_2,{\alpha }_3,{\mu }_E,{\mu }_A,{\mu }_{F_1},{\mu }_{F_2}\ge 0,\forall p,T,t\in {\mathrm{ℝ}}_{+}.\end{array}\right.$ (3.1)

The coefficients of the model are parametric functions of the weekly precipitation index p and the mean temperature T . The oviposition rate carried out by post-bloodmeal females is symbolised by ϕ. The rate α ₁ corresponds to the development of eggs in an aquatic population. In turn, the rate of individuals of the aquatic population that develop to blood pre-bloodmeal females is represented by α ₂. Rate α ₃ represents the maturation rate of pre-bloodmeal females to blood post-bloodmeal females. The natural mortality rates of the populations involved are given by µ _E , µ _A , ${\mu }_{F_1}$ e ${\mu }_{F_2}$ , respectively, linear term, $\varphi \left(1-\frac{E}{C}\right)F_2$ includes itself the term $-\frac{\varphi E{F}_2}{C}$ , which stands to mitigate rate ϕ if the population of stage E becomes large enough compared to the value of C⁶ .

The relationship between the life cycle of Ae. aegypti and the temperature is strongly consolidated in the literature ^{1), (30} . In turn, the relationship between the population of Ae. aegypti and precipitation has not been sufficiently studied in laboratories or field experiments ²⁵ . The association between those dependencies is therefore ignored so far. Let $\mathit{\Pi }=\left(\varphi ,{\alpha }_1,{\alpha }_2,{\alpha }_3,{\mu }_E,{\mu }_A,{\mu }_{F_1},{\mu }_{F_2}\right)$ be a general representation of the model (3.1) rates, and let the parametric dependencies of the rates on precipitation p and temperature T be, respectively, $\mathit{\Pi }=\mathit{\Psi }\left(p\right)$ and $\mathit{\Gamma }=\mathit{\Gamma }\left(T\right)$ . Then, the design for the dependence of the rate Π on p and T could be constructed by supposing it a well behaved enough function (analytic or class C ⁿ ), such that $\mathit{\Pi }\left(\mathit{\Psi },\mathit{\Gamma }\right)=a_0+a_1\mathit{\Psi }+b_1\mathit{\Gamma }+O\left({\mathit{\Psi }}^2,{\mathit{\Gamma }}^2\right)$ . In this case, the association of the dependencies of the model rates on p and T would be considered in a general form. By simplicity, as a working hypothesis, and because the model (daily) rates Π are typically smaller than one, we retain the linear part to the purpose of the association:

$\mathit{\Pi }\left(p,T\right)=\mathit{\Psi }\left(p\right)+\mathit{\Gamma }\left(T\right).$ (3.2)

Hence, by such simplifying consideration, the dependencies on p and T turn to be additive, as if they were rates itself in parallel association. Furthermore, such approach helps in minimisation of the number of free parameters.

To provide a insight on the parametric dependency of rates on precipitation Ψ ( p ), let us focus in the Figure 2a and notice that the peaks pattern resembles the precipitation peaks pattern. This fact is also manifested in Figure 2b , where there are significant correlations ( p < 0 . 05) between precipitation and MFAI for lags ranging from 0 to 12. Because of that similarity, the dependency of the rates on precipitation is supposed to be a monotonously increasing function. Other premises adopted here are simplicity and the fewest possible parameters. For those reasons, the parametric function $\mathit{\Psi }=\mathit{\Psi }\left(p\right)$ was written in terms of a power law, as

$\mathit{\Psi }\left(p\right)={\mathit{\Psi }}_0+{\left(\frac{p-p_0}{p_1-p_0}\right)}^r\left({\mathit{\Psi }}_1-{\mathit{\Psi }}_0\right);\text{with}p\ge p_0.$ (3.3)

Figure 2:
(a) Time series of precipitation, temperature and MFAI data from Caratinga-MG, Brazil, and (b) cross correlation between precipitation and MFAI time-series. The crossing curve stands for the significance level. Lags between 1 to 12 are associated to significant correlations ( p < 0 . 05)

The values Ψ ₀ and Ψ ₀ are associated to p ₀ = 0 mm and p ₁ = 34 . 62 mm, that correspond in weekly base to the threshold of rainfall volume between subtropical and tropical rainforest (1800mm/year) ¹⁶ . This power law parametric equation can provide increasing, decreasing, constant functions which can be also concave up and down. In the present work, we consider r > 0, focusing on increasing dependencies.

More exact relations for the dependencies between rates and precipitation should be obtained by optimisation or more complex experiments. The values adopted in this work ( Table 1 ) belong to the spectrum of values from other studies, regardless if they are related to precipitation. Due to this lack of exact values the constants of Equation (3.3) are set here in terms of orders of magnitude.

Table 1:
Free parameters of the dependencies of the model rates on temperature (3.4) and precipitation (3.3) with references which based their choice.

Π( p,T )	Y0 (Reference)	Y1 (Reference)	a	b	T op
ϕ	1.060327 29	8.294997 29	-1.64×10-2	7.12	33.42
α1	0.073206 1	0.55 7)	-1.17×10-3	0.26	31.15
α2	0.036597 29)	0.116117 29	-5.51×10-4	0.14	28.23
α3	0.2 24	0.2 24	-1.92×10-3	0.30	27.50
μE	0.01 17	0.01 17	2.64×10-4	0.00664	27.34
μA	0.11 24	0.11 24	9.03×10-4	0,0157	23.74
${\mu }_{F_1}={\mu }_{F_2}={\mu }_F$	0.02709 30)	0.09169 30)	2.95×10-4	0,0305	23.58

Although there are many studies for the dependence of the rates on temperature ^{1), (21), (30} , in a first approach and as a working hypothesis, we adopt simple functions according to the following assumptions: (1) simple relation with a small number of free parameters; (2) the existence of an optimal temperature for the metabolic processes. As a consequence, there arises a maximum value for the rates of development ( Figures 3a , 3b and 3f ), or a minimum, for the mortality rates ( Figures 3c , 3d and 3e ). Once the temperature gets far from the optimum value, it is then expected that there decreases the fitness in the life cycle. Hence, for the case of rates of development ( α ₁, α ₂ and ϕ), as far as the temperature is from the optimum value, their values decrease until vanishing for very high or very low temperatures. In the case of mortalities ( µ _E , µ _A and µ _F ), their values increase as the temperature becomes far from the optimum value. To fulfill these conditions, the rates were modelled by simple polynomials of degree 2, concave up for mortalities, and down for development rates with vanishing values in the outer intervals from the roots to avoid negative values.

$\mathit{\Gamma }\left(T\right)=\mathit{a}{\left(T-T_{op}\right)}^2+\mathit{b}.$ (3.4)

Figure 3:
Parameter setting α₁, α₂, µ_E , µ_A and µ_F , considering experimental data available in ^{1), (21), (30} .

The values of a , b , and T _op were obtained by fitting the degree 2 polynomials to experimental data present in some works ^{1), (21), (30} . The idea is more to find the best polynomials of degree 2 according to data, than to provide a high goodness of fit. Figure 3 depicts the modelled parameters Γ( T ) together with the set of experimental data in related literature. The parameters of rate α ₃ were determined by imposing three conditions: the inverse of the mean lifetime of F ₁ of 5 days was set as the mean value of the function (3.4), and the two roots were chosen to be the temperatures 15 ^◦C and 40 ^◦C. Because outside of such interval the females do not bite, the values of T < 15 ^◦C and T > 40 ^◦C are set to be zero ^{4), (21} . In the cases of rates of development, outside the roots of the polynomials, the values in the model are fixed to zero, although it never happens in the present range of data.

The system (3.1) is non-autonomous because it has coefficients dependent on rainfall and temperature which, in turn, depend o time ²⁰ . However the temperature and precipitation data is updated only weekly and the time evolution is hourly. So the evolution of the system along the yer can be viewed as a succession of autonomous systems. Hence, in each week, the system parameters become constant in time. Under such assumption, by applying the condition of equilibrium $\frac{d\mathit{X}}{xt}=\mathbf{0}$ with $\mathit{X}=\left(E,A,F_1,F_2\right)$ in model (3.1), we identify a trivial equilibrium point P ₀ and a non-trivial equilibrium point P ₁, described by

$P_0=\left(E*,A*,F_1^{*},F_2^{*}\right)=\left(0,0,0,0\right)$ (3.5)

$P_1=\left(E**,A**,F_1^{**},F_2^{**}\right)=\left(\left(1-\frac{1}{Q_0}\right)C,\frac{{\alpha }_1}{{\alpha }_2+{\mu }_A}E**,\frac{{\alpha }_2}{{\alpha }_3+{\mu }_{F_1}}A**,\frac{{\alpha }_3}{{\mu }_{F_2}}{F}_1^{**}\right)$ (3.6)

in which Q ₀ is the vector reproduction number, given by

$Q_0=\left(\frac{{\alpha }_1}{{\alpha }_1+{\mu }_E}\right)\left(\frac{{\alpha }_2}{{\alpha }_2+{\mu }_A}\right)\left(\frac{{\alpha }_3}{{\alpha }_3+{\mu }_{F_1}}\right)\left(\frac{\varphi }{{\mu }_{F_2}}\right).$ (3.7)

The equilibrium points (see ^{Appendix A} ) are degenerated and trivial, are $P1=P0$ , if $Q_0=1$ . The point P ₀ is asymptotically stable and P ₁, unstable, if $0<Q_0<1$ , while P 1 turns to be asymptotically stable and P ₀ unstable, if $Q_0>1$ . The system shows transcritical bifurcation at $Q_0=1$ ²³ .

As the meteorological parameters vary weekly along the ranges $\left\{T_{min}\le T\le T_{max}\right\}\times \left\{0\le p\le p_{max}\right\}$ , the model parameters (3.2) vary. Therefore, Q ₀ has different values each week but, according the the parameters defined, it remains positive $Q_0\ge 3.9187$ along the range of ( p, T ) data. Consequently, the non trivial equilibrium point P ₁ varies, but remains asymptotically stable. In particular, the equilibrium point coordinate $F_2^{**}$ values may be represented by a surface according to p and T (see Fig. 4 ).

Figure 4:
Surface of the females coordinate $F_2^{**}=F_2^{**}\left(p,T\right)$ of the non trivial asymptotically stable equilibrium point, according to meteorological variables. The dependence on temperature reaches a maximum at a optimum temperature and increases monotonically with increasing precipitation.

4 NUMERICAL SOLUTION AND MODEL VERIFICATION

Model (3.1) was solved numerically with a 4th order Runge-Kutta Method using Matlab^® software with time step 0.001, corresponding to 7000 points per week of data. In the evaluation of the results, we are interested only in the evolution of the population F ₂( t ), since the data refer to the capture of pregnant females who seek the trap to make oviposition. The model is evaluated, minimising the mean square error given by:

$S^2=\frac{1}{2N}\sum _{k=1}^N{\left(I\left(k\right)-\lambda F_2\left(k-l\right)\right)}^2,$ (4.1)

where N is the number of epidemiological weeks of the data and I ( k ) is the value of MFAI data in the week k . Since the weekly monitoring of traps takes time, the lag between model and data is given by l . The scale of the model is arbitrarily related to the environmental capacity, so the λ is a scale factor to adjust the scale to the sample capture data (MFAI). Subject to condition $\frac{\partial S^2}{\partial \lambda }=0$ , one finds

${\lambda }_l=\frac{{\displaystyle \sum _{k=1}^N}I\left(k\right)F_2\left(k-l\right)}{{\displaystyle \sum _{k=1}^N}{\left(F_2\left(k-l\right)\right)}^2}.$ (4.2)

By setting l _max as the value for which the cross-correlation between MFAI series I ( k ) and F ₂ reaches its maximum, the scale factor between data and model can be fixed so that $\lambda ={\lambda }_{l_{max}}$ . The value of parameter r in equation (3.3) was determined empirically by a simple one-dimensional optimisation as the process is repeated for values or r such that $0<r\le 2$ , with step of 0 . 1 searching for the lowest value of S ². In the present case, it was obtained for r = 1 . 2, value used in the simulation. For simplicity, the carrying capacity C ( p, T ) of the medium was set to be constant, C = 1, defined as relative to the best environmental conditions, when all existing water containers are available, filled of water and larval nutrients.

To evaluate the result, it was used cross-correlation between the MFAI series I ( k ) and model female population corrected by the scale factor $f_2=\lambda F_2\left(k\right)$ (see Figure 5a ). The maximum of the correlation is 0.75 at lag l _max = 6. The cross correlation between data and model remains significant ( p < 0 . 05) for lags ranging from 0 to 14. Figure 5b depicts a comparison between MFAI and $f_2=\left(k-6\right)$ . The simple model reproduces with certain quality the field data of captures. More complex details like peaks are not fitted in general.

Figure 5:
(a) Cross-correlation between model population f ₂ and infestation index MFAI. The crossing line stands for significant level. Maximum correlation occurs for lag of l = 6 weeks (p¡0.05) (b) Plots of the infestation index MFAI together with the lagged, l = 6, and scaled f ₂ model population.

5 CONCLUSIONS

After the comparison between model and data of female captures, one can conclude the curves are similar, indicating that the proposed model is able to reproduce the data in an acceptable manner. y the population of F ₂( t ). In some intervals, the model population becomes larger than the MFAI indicator, and, in others, the opposite occurs. Although the modelled peaks are not as pronounced as the infestation data, it is possible to notice their existence. The provided infestation data is obtained from a a public health program and by this reason, it is subjected to entomological control with insecticides. However, with few simple hypotheses on the dependence of the entomological parameters with temperature and precipitation, the model was able to reproduce the baseline of the infestation. Although the determination of dependencies of model rates on temperature and precipitation were subjected to simplifying assumptions and, therefore, of limited realism, the capability of the model to reproduce and forecast the data according to meteorological parameters is good. The task to evaluate the approach of those dependencies on meteorological data was successful. Eventual future refinement on those parametric dependencies with more realistic temperature parameters could provide more accuracy. Since the modelled precipitation parameters are set only as order of magnitude, an approach via optimisation could search for candidates for the values of the dependence of the entomological rates of development with precipitation. The present parameters can be used as initial set of parameters for optimisation. The result of the optimised parameters could guide biological experiments. The model is able to predict infestation according to meteorological parameters in public heath control programs.

Acknowledgments

The authors express their gratitude to the financial support of the Minas Gerais State Research Support Foundation (FAPEMIG), to the Federal Centre for Technological Education of Minas Gerais (CEFET), and Ecovec S.A.

REFERENCES

• ¹ E.B. Beserra, C.R. Fernandes, S.A.D.O. Silva, L.A.D. Silva & J.W.D. Santos. Efeitos da temperatura no ciclo de vida, exigências térmicas e estimativas do número de gerações anuais de Aedes aegypti (Diptera, Culicidae). Iheringia. Série Zoologia, 99 (2009), 142-148.
• ² M. Burattini, M. Chen, A. Chow, F. Coutinho, K. Goh, L. Lopez, S. Ma & E. Massad. Modelling the control strategies against dengue in Singapore. Epidemiology & Infection, 136(3) (2008), 309-319.
• ³ D.D. Chadee, F.L. Williams & U.D. Kitron. Impact of vector control on a dengue fever outbreak in Trinidad, West Indies, in 1998. Tropical Medicine & International Health, 10(8) (2005), 748-754.
• ⁴ S.R. Christophers. “Aedes aegypti: the yellow fever mosquito”. CUP Archive (1960).
• ⁵ E.A.P.D.A. Costa, E.M.D.M. Santos, J.C. Correia & C.M.R.D. Albuquerque. Impact of small variations in temperature and humidity on the reproductive activity and survival of Aedes aegypti (Diptera, Culicidae). Revista Brasileira de Entomologia, 54(3) (2010), 488-493.
• ⁶ L. Esteva & H.M. Yang. Mathematical model to assess the control of Aedes aegypti mosquitoes by the sterile insect technique. Mathematical biosciences, 198(2) (2005), 132-147.
• ⁷ D.A. Focks, D.G. Haile, E. Daniels & G.A. Mount. Dynamic life table model for Aedes aegypti (Diptera: Culicidae): analysis of the literature and model development. Journal of medical entomology, 30(6) (1993), 1003-1017.
• ⁸ FUNASA. Programa Nacional de Controle da Dengue. Ministério da Saúde Vigilância Epidemiológica, (2002).
• ⁹ A.D.S. Gomes, C.J. Sciavico & Á .E. Eiras. Periodicity of oviposition of females of Aedes aegypti (Linnaeus, 1762)(Diptera: Culicidae) in laboratory and field. Revista da Sociedade Brasileira de Medicina Tropical, 39(4) (2006), 327-332.
• ¹⁰ D.J. Gubler. The global emergence/resurgence of arboviral diseases as public health problems. Archives of medical research, 33(4) (2002), 330-342.
• ¹¹ M.U.G. Kraemer, M.E. Sinka, K.A. Duda, A.Q. Mylne, F.M. Shearer, C.M. Barker, C.G. Moore, R.G. Carvalho, G.E. Coelho, W. Van Bortel et al. The global distribution of the arbovirus vectors Aedes aegypti and Ae. albopictus. elife, 4 (2015), e08347.
• ¹² G. Kuno. Review of the factors modulating dengue transmission. Epidemiologic reviews, 17(2) (1995), 321-335.
• ¹³ A.P.M. Moreira, V.A.L. C. & S.J. C. Tendências climáticas e anomalias de precipitação em CaratingaMG. In “Os Desafios da Geografia Física na Fronteira do Conhecimento” (2018), pp. 2000-2009.
• ¹⁴ V.S.G. Neto & J.M.M. Rebêlo. Aspectos epidemiológicos do dengue no Município de São Luís, Maranhão, Brasil, 1997-2002. Cad. Saúde Pública, 20(5) (2004), 1424-1431.
• ¹⁵ D.P. Neves. “Parasitologia humana”. Atheneu (2005).
• ¹⁶ S. Noguchi, A.R. Nik & M. Tani. Rainfall characteristics of tropical rainforest at Pasoh forest reserve, Negeri Sembilan, Peninsular Malaysia. In “Pasoh”. Springer (2003), pp. 51-58.
• ¹⁷ M. Otero, H.G. Solari & N. Schweigmann. A stochastic population dynamics model for Aedes aegypti: formulation and application to a city with temperate climate. Bulletin of mathematical biology, 68(8) (2006), 1945-1974.
• ¹⁸ A.F. Ribeiro, G.R.A.M. Marques, J.C. Voltolini & M.L.F. Condino. Associação entre incidência de dengue e variáveis climáticas. Rev Saúde Pública, 40(4) (2006), 671-6.
• ¹⁹ T.W. Scott, P.H. Amerasinghe, A.C. Morrison, L.H. Lorenz, G.G. Clark, D. Strickman, P. Kittayapong & J.D. Edman. Longitudinal studies of Aedes aegypti (Diptera: Culicidae) in Thailand and Puerto Rico: blood feeding frequency. Journal of medical entomology, 37(1) (2000), 89-101.
• ²⁰ L.D.S.B. Silva, R.T.N. Cardoso, J.L.A. Fernandes, C.A. Silva & A.E. Eiras. Modelo Entomológico Determinístico sob Efeito da Pluviosidade para o Aedes aegypti e o Aedes albopictus. TEMA (São Carlos), 19(2) (2018), 289-303.
• ²¹ G. Smith, D. Eliason, C. Moore & E. Ihenacho. Use of elevated temperatures to kill Aedes albopictus and Ae. aegypti. Journal of the American Mosquito Control Association, 4(4) (1988), 557-558.
• ²² R. Souza-Santos & M.S. Carvalho. Análise da distribuição espacial de larvas de Aedes aegypti na Ilha do Governador, Rio de Janeiro, Brasil. Cadernos de saúde Pública, 16(1) (2000), 31-42.
• ²³ S.H. Strogatz. “Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering”. CRC Press (2018).
• ²⁴ A. Tran, G. L’Ambert, G. Lacour, R. Benôıt, M. Demarchi, M. Cros, P. Cailly, M. Aubry-Kientz, T. Balenghien & P. Ezanno. A rainfall-and temperature-driven abundance model for Aedes albopictus populations. International journal of environmental research and public health, 10(5) (2013), 1698-1719.
• ²⁵ J. Waldock, N.L. Chandra, J. Lelieveld, Y. Proestos, E. Michael, G. Christophides & P.E. Parham. The role of environmental variables on Aedes albopictus biology and chikungunya epidemiology. Pathogens and global health, 107(5) (2013), 224-241.
• ²⁶ S. Wiggins. “Introduction to applied nonlinear dynamical systems and chaos”, volume 2. Springer Science & Business Media (2003).
• ²⁷ World Health Organization. Dengue and severe dengue. http://www.who.int/mediacentre/factsheets/fs117/en/ (2017). Acessado em 20/03/2018.
  » http://www.who.int/mediacentre/factsheets/fs117/en/
• ²⁸ World Health Organization. Neglected tropical diseases. http://www.searo.who.int/entity/vector_borne_tropical_diseases/data/data_factsheet/en/ (2018). Acessado em 10/09/2018.
  » http://www.searo.who.int/entity/vector_borne_tropical_diseases/data/data_factsheet/en/
• ²⁹ H.M. Yang. The transovarial transmission in the dynamics of dengue infection: Epidemiological implications and thresholds. Mathematical biosciences, 286 (2017), 1-15.
• ³⁰ H.M. Yang, M.L.G. Macoris, K.C. Galvani, M.T.M. Andrighetti & D.M.V. Wanderley. Assessing the effects of temperature on the population of Aedes aegypti, the vector of dengue. Epidemiology and infection, 137(08) (2009), 1188-1202.

Appendix A STABILITY OF THE EQUILIBRIUM POINTS

Let the system of equations (3.1) rewritten as

$\frac{d\mathit{X}}{dt}=\mathit{f}\left(\mathit{X}\right),\text{where}\mathit{X}=\left(E,A,F_1,F_2\right).$ (A.1)

The condition $\mathit{f}\left(\mathit{X}\right)=0$ revealed two critical points, one trivial, P ₀(Equation 3.5), and one nontrivial, P ₁(Equation 3.6). In the following it is performed the analysis of stability of the system.

After the expansion near the equilibrium points P _k ; k = 0;1 and computation of the Jacobian matrices $B_{P_k}={\left[\frac{\partial f}{\partial \mathit{X}}\right]}_{P_k}$ , k = 0 , 1, their characteristic polynomials are given, respectively, by

${\Phi }_{P_k}\left(\lambda \right)=a_{k4}{\lambda }^4+a_{k3}{\lambda }^3+a_{k2}{\lambda }^2+a_{k1}\lambda +a_{k0},$ (A.2)

with coefficients for P ₀,

$a_{04}=1,$

$a_{03}={\alpha }_1+{\mu }_E+{\alpha }_2+{\mu }_A+{\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2},$

$a_{02}=\left({\alpha }_1+{\mu }_E\right){\mu }_{F_2}+\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right)\left({\alpha }_1+{\mu }_E+{\mu }_{F_2}\right)+\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right),$

$a_{01}=\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right)+\left({\alpha }_1+{\alpha }_2+{\mu }_E+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2},$

$a_{00}=\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}\left(1-Q_0\right);$

and for P ₁,

$a_{14}=1,$

$a_{13}=({\alpha }_1+{\mu }_E)Q_0+{\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2},$

$a_{12}=({\alpha }_1+{\mu }_E)({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2})Q_0+({\alpha }_2+{\mu }_A)({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2})+({\alpha }_3+{\mu }_{F_1}){\mu }_{F_2},$

$a_{11}=({\alpha }_1+{\mu }_E)({\alpha }_2+{\mu }_A)({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2})Q_0+\left(\right({\alpha }_1+{\mu }_E)Q_0+{\alpha }_2+{\mu }_A)({\alpha }_3+{\mu }_{F_1}){\mu }_{F_2},$

$a_{10}=\varphi {\alpha }_1{\alpha }_2{\alpha }_3(1-\frac{1}{Q_0}).$

The stability of each equilibrium points P _k can be assessed from the sign of the real part of their respective eigenvalues. The polynomial ${\Phi }_{P_k}$ (A.2) with have roots with real part strictly less than zero if and only if all elements in the first column of the Routh Table (A.3) are nonzero and have the same sign ²⁶ .

$R^{\left(k\right)}=\left[\begin{array}{ccc}1& a_{k2}& a_{k0}\\ a_{k3}& a_{k1}& 0\\ \frac{a_{k2}{a}_{k3}-a_{k1}}{a_{k3}}& a_{k0}& 0\\ a_{k1}-\frac{a_{k3}^2{a}_{k0}}{a_{k2}{a}_{k3}-a_{k1}}& 0& 0\\ a_{k0}& 0& 0\end{array}\right]$ (A.3)

The elements $R_{11}^{\left(k\right)}=1>0$ , k = 0 , 1 are positive. Since all model rates remain positive by definition, $R_{21}^{\left(k\right)}>0$ , k = 0 , 1 is positive. As consequence, $R_{31}^{\left(k\right)}$ , k = 0 , 1, are positive if $a_{k2}{a}_{k3}-a_{k1}>0$ . For k = 0, it reads

$$\begin{gathered} a_{02}{a}_{03}-a_{01}=2({\alpha }_1+{\mu }_E)({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}){\mu }_{F_2}+2({\alpha }_2+{\mu }_A)({\alpha }_3+{\mu }_{F_1})({\alpha }_1+{\mu }_E+{\mu }_{F_2}) \\ +({\alpha }_1+{\mu }_E{)}^2({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2})+({\alpha }_2+{\mu }_A{)}^2({\alpha }_1+{\alpha }_3+{\mu }_E+{\mu }_{F_1}+{\mu }_{F_2}) \\ +({\alpha }_3+{\mu }_{F_1}{)}^2({\alpha }_1+{\alpha }_2+{\mu }_E+{\mu }_A+{\mu }_{F_2})+({\mu }_{F_2}{)}^2({\alpha }_1+{\alpha }_2+{\alpha }_3+{\mu }_E+{\mu }_A+{\mu }_{F_1})>0. \end{gathered}$$

In the case k = 1, it is expressed as

$$\begin{gathered} a_{12}{a}_{13}-a_{11}=({\alpha }_1+{\mu }_E{)}^2({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2})Q_0^2+2({\alpha }_2+{\mu }_A\left)\right({\alpha }_3+{\mu }_{F_1}){\mu }_{F_2} \\ +({\alpha }_1+{\mu }_E\left)\right({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2}\left)\right({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1})Q_0 \\ +({\alpha }_2+{\mu }_A{)}^2({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2})+{\mu }_{F_2}^2\left(\right({\alpha }_1+{\mu }_E)Q_0+{\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}) \\ +({\alpha }_3+{\mu }_{F_1}{)}^2({\alpha }_2+{\mu }_A+{\mu }_{F_2})+({\alpha }_1+{\mu }_E\left)\right({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}){\mu }_{F_2}{Q}_0>0. \end{gathered}$$

Hence, $R_{13}^{\left(k\right)}>0$ , k = 0 , 1. Since $a_{k2}{a}_{k3}-a_{k1}>0$ , k = 1 , 2, the analysis of $R_{14}^{\left(k\right)}$ restricts in the analysis of whether $a_{k1}{a}_{k2}{a}_{k3}-a_{k01}^2-a_{k3}^2{a}_{k0}>0$ , k = 0 , 1. For k = 0 the expression reads

$$\begin{gathered} a_{01}{a}_{02}{a}_{03}-a_{01}^2-a_{03}^2{a}_{00}={\left({\alpha }_1+{\mu }_E\right)}^3{\left({\alpha }_2+{\mu }_A\right)}^2\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right) \\ +{\left({\alpha }_1+{\mu }_E\right)}^3{\left({\alpha }_3+{\mu }_{F_1}\right)}^2\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right) \\ +{\left({\alpha }_1+{\mu }_E\right)}^3{\left({\mu }_{F_2}\right)}^2\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right) \\ +2{\left({\alpha }_1+{\mu }_E\right)}^3\left(a_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2} \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\left(a_2+{\mu }_A\right)}^3\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right) \\ +{\left({\alpha }_2+{\mu }_A\right)}^3{\left(a_3+{\mu }_{F_1}\right)}^2\left({\alpha }_1+{\mu }_E+{\mu }_{F_2}\right) \\ +{\left({\alpha }_2+{\mu }_A\right)}^3{\mu }_{F_2}^2\left({\alpha }_1+{\alpha }_3+{\mu }_E+{\mu }_{F_1}\right) \\ +2\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_2+{\mu }_A\right)}^3\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2} \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^3\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right) \\ +{\left({\alpha }_2+{\mu }_A\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^3\left({\alpha }_1+{\mu }_E+{\mu }_{F_2}\right) \\ +{\left({\alpha }_3+{\mu }_{F_1}\right)}^3{\mu }_{F_2}^2({\alpha }_1+{\alpha }_2+{\mu }_E+{\mu }_A \\ +2\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right){\left({\alpha }_3+{\mu }_{F_1}\right)}^3{\mu }_{F_2} \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\mu }_{F_2}^3\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right) \\ +{\left({\alpha }_2+{\mu }_A\right)}^2{\mu }_{F_2}^3\left({\alpha }_1+{\alpha }_3+{\mu }_E+{\mu }_{F_1}\right) \\ +{\left({\alpha }_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}^3\left({\alpha }_1+{\alpha }_2+{\mu }_E+{\mu }_A\right) \\ +2\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}^3 \\ +4{\left({\alpha }_1+{\mu }_E\right)}^2\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right) \\ +4\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_2+{\mu }_A\right)}^2\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right) \\ +4\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right){\left({\alpha }_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}\left({\alpha }_3+{\mu }_E+{\mu }_{F_2}\right) \\ +\varphi {\alpha }_1{\alpha }_2{\alpha }_3{\left({\alpha }_1+{\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_E+{\mu }_{F_1}+{\mu }_{F_2}\right)}^2>0. \end{gathered}$$

The expression is positive because all parameters of the model are positive. In the case k = 1, it reads

$$\begin{gathered} a_{11}{a}_{12}{a}_{13}-a_{11}^2-a_{10}{a}_{13}^2={\left({\alpha }_1+{\mu }_E\right)}^3{\left({\alpha }_2+{\mu }_A\right)}^2\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right)Q_0^3 \\ +{\left({\alpha }_1+{\mu }_E\right)}^3{\left({\alpha }_3+{\mu }_{F_1}\right)}^2\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right)Q_0^3 \\ +{\left({\alpha }_1+{\mu }_E\right)}^3{\mu }_{F_2}^2\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right)Q_0^3 \\ +{\left({\alpha }_1+{\mu }_E\right)}^3\left(a_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}{Q}_0^2\left(2Q_0+1\right) \\ +2{\left({\alpha }_1+{\mu }_E\right)}^2{\left(a_2+{\mu }_A\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^2{Q}_0^2 \\ +2{\left({\alpha }_1+{\mu }_E\right)}^2{\left(a_2+{\mu }_A\right)}^2{\mu }_{F_2}^2{Q}_0^2 \\ +2{\left({\alpha }_1+{\mu }_E\right)}^2{\left(a_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}^2{Q}_0^2 \\ +4{\left({\alpha }_1+{\mu }_E\right)}^2\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2}\right)Q_0^2 \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\left({\alpha }_2+{\mu }_A\right)}^3\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right)Q_0^2 \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^3\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right)Q_0^2 \\ +{\left({\alpha }_1+{\mu }_E\right)}^2{\mu }_{F_2}^3\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right)Q_0^2 \\ +\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_2+{\mu }_A\right)}^3{\left({\alpha }_3+{\mu }_{F_1}+{\mu }_{F_2}\right)}^2{Q}_0 \\ +\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_3+{\mu }_{F_1}\right)}^3{\left({\alpha }_2+{\mu }_A+{\mu }_{F_2}\right)}^2{Q}_0 \\ +\left({\alpha }_1+{\mu }_E\right){\mu }_{F_2}^3{\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}\right)}^2{Q}_0 \\ +2\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_2+{\mu }_A\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}\left(2Q_0+1\right) \\ +2\left({\alpha }_1+{\mu }_E\right){\left({\alpha }_2+{\mu }_A\right)}^2\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}^2\left(2Q_0+1\right) \\ +2\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right){\left({\alpha }_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}^2\left(2Q_0+1\right) \\ +{\left({\alpha }_2+{\mu }_A\right)}^3\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}\left({\alpha }_1+{\alpha }_3+{\mu }_E+{\mu }_{F_1}+{\mu }_{F_2}\right) \\ +\left({\alpha }_2+{\mu }_A\right){\left({\alpha }_3+{\mu }_{F_1}\right)}^3{\mu }_{F_2}\left({\alpha }_1+{\alpha }_2+{\mu }_E+{\mu }_A+{\mu }_{F_2}\right) \\ +\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right){\mu }_{F_2}^3\left({\alpha }_1+{\alpha }_2+{\alpha }_3+{\mu }_E+{\mu }_A+{\mu }_{F_1}\right) \\ +2{\left({\alpha }_2+{\mu }_A\right)}^2{\left({\alpha }_3+{\mu }_{F_1}\right)}^2{\mu }_{F_2}^2 \\ +2\varphi {\alpha }_1{\alpha }_2{\alpha }_3\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\alpha }_3+{\mu }_A+{\mu }_{F_1}+{\mu }_{F_2}\right)>0. \end{gathered}$$

It is also positive because all parameters of the model are positive.

Finally, for k = 0, $R_{15}^{\left(0\right)}=a_{00}=\left({\alpha }_1+{\mu }_E\right)\left({\alpha }_2+{\mu }_A\right)\left({\alpha }_3+{\mu }_{F_1}\right)\left({\mu }_{F_2}\right)\left(1-Q_0\right)$ and for k = 1, $R_{15}^{\left(1\right)}=a_{10}=\varphi {\alpha }_1{\alpha }_2{\alpha }_3\left(1-\frac{1}{Q_0}\right)$ . In both cases, the signal of the matrix element will depend on Q ₀to be $0<Q_0<1$ or $Q_0>1$ . Moreover, the factors $\left(1-Q_0\right)$ and $\left(1-\frac{1}{Q_0}\right)$ have opposite signals in those cases. Therefore, in the case $0<Q_0<1,R_{15}^{\left(0\right)}>0$ and the first column of the Routh Table for k = 0 will be strictly positive, causing P ₀to be stable. Yet, for k = 1, in the same case, $R_{15}^{\left(1\right)}<0$ , and there will have a change of sign in the last element of the first column of the Routh Table, causing P ₁to be unstable. Now, turning to the case $Q_0>1$ , we now have $R_{15}^{\left(0\right)}<0$ for k = 0 and $R_{15}^{\left(1\right)}>0$ for k > 1. Hence the change of sign in the first column of the Routh Tables switches and now P ₀becomes unstable and P ₁becomes stable. In the case Q ₀= 1, $R_{15}^{\left(k\right)}=0$ for k = 0; 1, P ₀and P ₁are coincident and trivial,the independent terms a _k0 of the characteristic polynomials ${\Phi }_{P_k}\left(\lambda \right)$ , k = 0;1 vanish and there take place a null eigenvalue, λ = 0. This fact together with the change of stability in ρ = 1 is sufficient to characterise the point as an transcritical bifurcation point.

Publication Dates

• Publication in this collection
  05 Apr 2021
• Date of issue
  Jan-Mar 2021

History

• Received
  12 Sept 2018
• Accepted
  15 Oct 2020