Stability, Bifurcation and Chaos Analysis of Vector-Borne Disease Model with Application to Rift Valley Fever

This paper investigates a RVF epidemic model by qualitative analysis and numerical simulations. Qualitative analysis have been used to explore the stability dynamics of the equilibrium points while visualization techniques such as bifurcation diagrams, Poincaré maps, maxima return maps and largest Lyapunov exponents are numerically computed to confirm further complexity of these dynamics induced by the seasonal forcing on the mosquitoes oviposition rates. The obtained results show that ordinary differential equation models with external forcing can have rich dynamic behaviour, ranging from bifurcation to strange attractors which may explain the observed fluctuations found in RVF empiric outbreak data, as well as the non deterministic nature of RVF inter-epidemic activities. Furthermore, the coexistence of the endemic equilibrium is subjected to existence of certain number of infected Aedes mosquitoes, suggesting that Aedes have potential to initiate RVF epidemics through transovarial transmission and to sustain low levels of the disease during post epidemic periods. Therefore we argue that locations that may serve as RVF virus reservoirs should be eliminated or kept under control to prevent multi-periodic outbreaks and consequent chains of infections. The epidemiological significance of this study is: (1) low levels of birth rate (in both Aedes and Culex) can trigger unpredictable outbreaks; (2) Aedes mosquitoes are more likely capable of inducing unpredictable behaviour compared to the Culex; (3) higher oviposition rates on mosquitoes do not in general imply manifestation of irregular behaviour on the dynamics of the disease. Finally, our model with external seasonal forcing on vector oviposition rates is able to mimic the linear increase in livestock seroprevalence during inter-epidemic period showing a constant exposure and presence of active transmission foci. This suggests that RVF outbreaks partly build upon RVF inter-epidemic activities. Therefore, active RVF surveillance in livestock is recommended.


Introduction
Rift Valley fever (RVF) virus, a member of the genus phlebovirus and family Bunyaviridae, which has been isolated from at least 40 mosquito species in the field [1], infects both wild and domestic animals and humans. The RVF epizootics and epidemics are closely linked to the occurrence of the warm phase of the El Nino/Southern Oscillation (ENSO) phenomenon [2]. This phenomenon is characterized by elevated Indian Ocean temperatures which lead to heavy rainfall and flooding of habitats suitable for the production of immature Aedes and Culex mosquitoes that serve as the primary RVF virus (RVFV) vectors in East Africa [3,4]. Studies have shown that the life cycle of RVFV has distinct endemic and epidemic cycles. During the endemic cycle the virus persists during dry season/inter-epizootic periods through vertical transmission in Aedes mosquito eggs [3]. Aedes eggs need to be dry for several days before they can mature. After maturing, they hatch during the next flooding event large enough to cover them with water [5,6]. The eggs have high desiccation resistance and can survive dry conditions in a dormant form for months to years. At the beginning of the rainy season, Aedes mosquitoes quickly multiply into large numbers before declining due to the need for dry conditions for egg maturation [9]. There can be a second peak in mosquito densities at the end of the rainy season if there is a gap in rainfall for several days [5]. When these mosquitoes lay their eggs in flooded areas (including dambos), transovarially infected adults may emerge and transmit RVFV to nearby domestic animals, including sheep, goats, cattle, and camels. High viremias in these animals may then lead to the infection of secondary arthropod vector species including various Culex species [7].
Epizootic/epidemic cycles are driven by the subsequent elevation of various Culex mosquito populations, which serve as excellent secondary vectors if immature mosquito habitats remain flooded for a long enough period [4]. Their eggs require water to mature and hatch and the mosquitoes survive the dry season in adult form and during the rainy season, the population of Culex mosquitoes reaches a maximum towards the end of the season [9]. The propagation of these secondary vectors may spread the virus to additional infection in animal and human, causing an outbreak.
The disease is known to occur in outbreaks that come in cycles of 5-15 years in the Eastern Africa region and the Horn of Africa [10].
We observe that RVF outbreaks are highly linked to seasonal variations on rainfall, which is in turn reflected through seasonal fluctuations in mosquito population densities. Aedes eggs require water to hatch and dry condition for maturation, and at the beginning of the rainy season quickly grow to large numbers while Culex eggs require water to mature and hatch, and survive dry season in adult form and during the rainy season reach maximum numbers towards the end of the season. Thus, fluctuations in both seasons (wet and dry) favour the complex dynamics of both mosquito species. Hence the complexity observed on the dynamics of RVF virus transmission and maintenance.
The interplay between the internal nonlinear dynamic of ecological systems and various external factors that affect them, makes understanding of population fluctuation a unique problem [11].
Mathematical models have been developed in order to provide a better understanding of the nature and dynamics of the transmission and persistence of the disease, as well as predict outbreaks and simulate the impact of control strategies [9,12,17,18]. Most of these models considered constant mosquito oviposition rates, ignoring effects of seasonal fluctuations in the mosquito population size. Furthermore, some have ignored the effects of vertical transmission and secondary vectors [18] and some only considered Aedes species [9]. Temperature, rainfall and humidity have great influence in all stages of mosquito development from the emergence and viability of eggs, to the size and longevity of adults [19,20]. Recently, Mpeshe et al. [21] modified their previous study [18] to include vertical transmission in Aedes species and climate-driven parameters. These models provide important insights but do not investigate the stability dynamics and attractors structures of the model when there are external forces in the density of vector populations.
The most common manifestation of external forcing is through seasonality including both natural (e.g. the occurrence of the warm phase of the El Nino/Southern oscillation phenomenon) and induced (e.g human deforestation or human pollution).
Studies for understanding dynamical consequences of regular and stochastic external forcing are still ongoing but poorly understood [22][23][24][25]. To the best of our knowledge, no systematic investigation of stability and attractor structures of a realistic RVF model comprising two populations of mosquitoes (Aedes and Culex) and one livestock host population with two infected classes (asymptomatic and symptomatic) and seasonal variation on mosquito oviposition rates has been carried out. Based on the model proposed by Gaff et al. [12], we investigate a two vector and one host epidemic model, to capture the dynamical behaviour of both the disease free and endemic equilibria, the effects of seasonality on mosquito oviposition rates (b 1 ,b 3 ), parametrized by d 1 , d 3 and effects of asymptomatic class in livestock (parametrized by 1{h 2 ). We prove existence and global stability of both the disease-free and the endemic equilibria in the absence of secondary vectors (I 3~0 ), as well as the existence and local stability of both disease free and endemic equilibrium points of the overall model. We then investigate the structures of model attractors through bifurcation analysis, taking as bifurcation parameters d 1 and d 3 the strengths of seasonality of mosquito oviposition rates. The bifurcation diagrams with simultaneous variation of seasonal forcing on the oviposition rates of the two mosquito species reveal the complexity induced by their interactions. The understanding of possible state space scenarios through bifurcation analysis is helpful for understanding RVF epidemiological data with its seasonality aspects. To obtain robust analysis we then compute the largest Lyapunov exponents, Poincaré maps and maxima return maps.
The section methods gives a detailed description of the model and its parameters. In section results the model is used to study the dynamic behaviour of the disease stability and bifurcation analysis. Simulations are performed to investigate model dependence on initial condition and attractors structures of the model applying an external forcing on mosquito's oviposition rates.

Methods
Gaff et al. [12] proposed a one host and two vectors population model for RVF with vertical transmission in Aedes vectors to study the transmission of RVF and the impact of vertical transmission on the persistence of the disease. Chitnis et al. [9] analysed a RVF model with vertical transmission for Aedes mosquitoes and included asymptomatic class for livestock and removed one population of mosquitoes.
The model presented in this paper adopts a similar structure as in Gaff et al. [12]. We introduce an asymptomatic class for livestock [9], because for many species of livestock, RVF virus infection are frequently subclinical [26,27]. As the main purpose of this study is to study the dynamic behaviour of the disease, influenced by changes in climate and oscillation of rainfall, we include seasonal variation in the oviposition rates of both Aedes and Culex mosquitoes.
We divide the livestock population into four classes: susceptible, S 2 , asymptomatic, A 2 , infectious, I 2 , and recovered (immune), R 2 . Livestock enter the susceptible class through birth (at a constant rate). Birth rates are important because after an outbreak, herd immunity can reach 80% and the proportion of susceptible livestock must be renewed through birth or movement before another outbreak can occur [28]. When an infectious mosquito bites a susceptible animal, there is a finite probability that the animal becomes infected. Since the duration of the latent period in cattle is small relative to their life span, we do not model the exposed stage. Many adult cattle do not exhibit clinical signs apart from abortion of foetuses [6,26], thus, include an asymptomatic class for infectious animals that transmit the virus at a lower rate than those with acute clinical symptoms. After being successfully infected by an infectious Aedes and/or Culex mosquito, livestock move from the susceptible class S 2 to either the infected symptomatic I 2 or asymptomatic A 2 class. After some time, the symptomatic and asymptomatic livestock recover and move to the recovered class, R 2 . The recovered livestock have immunity to the disease for life. Cattle leave the population through a per capita natural death rate and through a per capita disease-induced death rate only for symptomatic livestock. The size of the livestock population is given by We divide the Aedes and Culex mosquitoes population into three classes: susceptible, S a , exposed, E a , and infectious, I a . The subscripts a~1 and a~3 represent Aedes and Culex mosquitoes, respectively. Female mosquitoes (we do not include male mosquitoes in our model because only female mosquitoes bite animals for blood meals) enter the susceptible class through birth. The virus enters a susceptible mosquito, S a , with finite probability, when the mosquito bites an infectious animal and the mosquito moves to the exposed class, E a . After some period of time, depending on the ambient temperature and humidity [29], the mosquito moves from the exposed class to the infectious class, I a .
To reflect the vertical transmission in the Aedes species, compartments for uninfected P 1 and infected U 1 eggs are included. As the Culex species cannot transmit RVF vertically, only uninfected eggs P 3 are included. Mosquitoes once infected remain infectious during their lifespan. Mosquitoes leave the population through a per capita natural death rate. The size of each adult mosquito population is N 1~S1 zE 1 zI 1 for adult Aedes mosquitoes and N 3~S3 zE 3 zI 3 for adult Culex mosquitoes. The three populations are modelled with carrying capacity K 1 ,K 2 ,K 3 , for Aedes, livestock and Culex respectively. While in [12], the total number of mosquito bites on cattle depends on the number of mosquitoes, in our model, the total number of bites varies with both the cattle and mosquito population sizes. This allows a more realistic modelling of situations where there is a high ratio of mosquitoes to cattle, and where cattle availability to mosquitoes is reduced through control interventions [9].

Mathematical Model
The state variables in Table 1 and parameters in Table 2 for the RVF model ( Figure 1) satisfy the following system of equations: Aedes Livestock Culex where from the model flowchart in Fig.1 Following the approach in [9], s a , where a~1 for Aedes and a~3 for Culex is the rate at which a mosquito would like to bite livestock (related to the gonotrophic cycle length), and s 2 is the maximum number of bites that an animal can support per unit time (through physical availability and any intervention measures on livestock taken by humans). Then, s a N a is the total number of bites that the mosquitoes would like to achieve per unit time and s 2 N 2 is the availability of livestock. Thus, the total number of mosquito-livestock contacts is half the harmonic mean of s a N a and s 2 N 2 , In addition to having the correct limits at zero and infinity, this form also meets the necessary criteria that b bƒmin(s a N a ,s 2 N 2 ) where b b is the total number of bites per unit time. The total number of mosquito-livestock contacts depends on the populations of both species. We define We defined the force of infection from mosquitoes to livestock, l a 2 (t), as the product of the number of mosquito bites that one animal has per unit time, b 2 , the probability of disease transmission from the mosquito to the animal, b 2a , and the probability that the mosquito is infectious, I a =N a . We define the force of infection from livestock to mosquitoes, l 2 a (t), as the force of infection from infectious (symptomatic and asymptomatic) livestock. This is expressed as the number of livestock bites one mosquito has per unit time, b b a ; the probability of disease transmission from an infected (asymptomatic) animal to the mosquito, b a2 (b b a2 ); and the probability that the animal is infectious, I 2 =N 2 (A 2 =N 2 ). Therefore the forces of infection are given by: The model system (1,2,3) is biologically relevant (solutions are positive) in the set V~( P 1 ,U 1 ,S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 ,R 2 ,P 3 ,S 3 ,E 3 ,I 3 ) [ R 13 z : P 1 ,U 1 , S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 ,R 2 ,P 3 ,S 3 , Lemma 1. The model system (1,2,3) is well-posed in V which is invariant and attracting.

Basic Reproduction Number
For epidemiology models, a quantity, R 0 is derived to assess the stability of the disease free equilibrium [12]. R 0 represents the the number of individuals infected by a single infected individual during his or her entire infectious period, in a population which is  [30]. When R 0 v1, if a disease is introduced, there are insufficient new cases per case, and the disease cannot invade the population. When R 0 w1, the disease may become endemic; the greater R 0 is above 1, the less likely stochastic fade out of the disease can occur. To compute this threshold we use the next generation operator approach, as described by Diekmann et al. [31] and van den Driessche and Watmough [32] as well as to describe the conditions for which the disease-free equilibrium points lose stability.
Since the model incorporates both vertical and horizontal transmission, R 0 for the system is the sum of the R 0 values for each mode of transmission determined separately [33], To compute each component of R 0 , the model equations in vector form are the difference between the rate of new infection in compartment i, F i and the rate of transfer between compartment i and all other compartments due to other processes, V i [32], (see Appendix S1). Then, R 0 is given by where R 0,V~b 1 q 1 m 1 and

Basic Reproduction Number for periodic environment
In periodic environment, the basic reproduction number is the generalization of the R 0 in non periodic environment. It is known as the transmissibility number R R 0 , which is defined as the average number of secondary cases arising from the introduction of a single infectious individual into a completely susceptible population at a random time of the year [34]. Thus, R R 0 is defined through the spectral radius of a linear integral operator on a space of periodic functions, given by the integral operator G j (see Appendix S1), As proposed by Bacaer [36], the transmissibility number R R 0 is given by where Re(.) is the real part of (.). G 0 is the basic reproduction number for the non-seasonal model, obtained when d i~0 .
The size of R R 0 is reduced compared to R 0 when oviposition rates are constant, and this makes it slightly difficult for the virus to invade the population with such fluctuations on the transmission rates [36].
From G 0 the following sub-reproduction numbers R 21 ,R 12 ,R 23 ,R 32 can be obtained: R 21 is the number of new infections in livestock from one infected Aedes mosquito and is given by representing the product of the probability that the Aedes mosquito survives the exposed stage c 1 c 1 zb 1 , the number of bites on livestock per mosquito l 0 1 N 0 2 , the probability of transmission per bite b 21 , and the infectious lifespan of Aedes mosquito 1=b 1 .
R 12 is the number of new infections in Aedes mosquitoes from one infected (asymptomatic or symptomatic) animal, and is given by the weighted sum of new infections resulting from asymptomatic and symptomatic livestock This is the product of the number of bites an animal receives l 0 1 N 0 1 , the probability of transmission per bite (b b 12 for an asymptomatic animal and b 12 for symptomatic animal), and the duration of the infective period ( 1 e e 2 zb 2 for an asymptomatic animal and 1 e 2 zb 2 zm 2 for symptomatic animal) weighted by the probability that an animal either becomes asymptomatic or symptomatic upon infection. R 23 is the number of new infections in livestock from one infected Culex mosquito and is given by This is the product of the probability that the Culex mosquito survives the exposed stage c 3 c 3 zb 3 , the number of bites on livestock per mosquito l 0 3 N 0 2 , the probability of transmission per bite b 23 , and the infectious lifespan of Culex mosquito 1=b 3 .
R 32 is the number of new infections in from an infected (asymptomatic or symptomatic) animal and is given by the weighted sum of new infections resulting from asymptomatic and symptomatic livestock This is the product of the number of bites one animal receives l 0 3 N 0 3 , the probability of transmission per bite (b b 32 for an asymptomatic animal and b 32 for symptomatic animal), and the duration of the infective period ( 1 e e2zb2 for an asymptomatic animal and Culex mosquitoes ~~ that an animal either becomes asymptomatic or symptomatic upon infection.
If q 1 w0, R 0 increases because vertical transmission directly increases the number of infectious mosquitoes and indirectly increases the transmission from livestock to mosquitoes and back to livestock.
Evaluating J at the disease-free equilibrium and using basic properties of matrix algebra, it is evident from the characteristic polynomial of J that the following eigenvalues l 1~{ m 1 ,l 2~{ h 1 ,l 3~{ m 2 ,l 4~{ h 3 ,l 5~{ m 3 have negative real part and the remaining reduced matrix is The stability of a disease-free equilibria should be established from the eigenvalues of the reduced Jacobian matrix (21). To simplify the computations, we perform the following operations on matrix (21): first we add the first row to the third one and take the resultant as the new third row; second we multiply the second row by c 1 =(c 1 zm 1 ) and add it to the new third row, then take the resultant as the new third row and at last we multiply the sixth row by c 3 =(c 3 zm 3 ) and add it to the last row and maintaining the rest as it is, we obtain the following matrix From the basic properties of matrix algebra, it is evident from the characteristic polynomial of J 2 that the following eigenvalues l 1~{ h 1 ,l 2~{ (c 1 zm 1 ) and l 3~{ (c 3 zm 3 ) have negativẽ 0.5 Stability analysis of the model (1,2,3) without Culex species In the absence of Culex species, I Ã 3~0 , equation (18) can be written as Equation (24) has two possible solutions implies an existence of a disease-free equilibria and the case I Ã 1 =0 implies an existence of an endemic equilibria. Let us now derive conditions under which positive endemic equilibria exist. For I Ã 1 =0, we get I Ã 1 is epidemiologically meaningful, that is, which can be written in the form for R 1 0 ƒ1 and exactly one endemic equilibrium point (EE), The result in Theorem 1 indicates the impossibility of backward bifurcation in the RVF model system (1,2,3) without Culex species since it has no endemic equilibrium when R 1 0 v1. This explains that the model (1,2,3) without Culex species has a globally asymptotically stable disease-free equilibrium whenever R 1 0 ƒ1. In its simplest form, backward bifurcation in epidemic models usually implies the existence of two subcritical endemic equilibria when the basic reproductive number for R 1 0 v1, and a unique supercritical endemic equilibrium for R 1 0 w1 [37]. Thus, a unique real part and the remaining reduced matrix is positive endemic equilibrium exists only when R 1 0 w1. We note that the increase in complexity of an epidemic model (by adding more infected classes, for example) can lead to backward bifurcation and even more complicated phenomena associated with endemic equilibria [37]. However, increase in complexity of the proposed RVF model does not appear to give rise to more complex behaviour with regard to endemic equilibria. 0.5.1 Local stability of DFE, X 0 1 . In the absence of secondary vector (Culex species) that serve as RVF outbreak amplifiers the Jacobian matrixJ J(X 0 ) in (23) reduces to J(X 0 1 )~b The characteristic equation corresponding to the above Jacobian matrix is where Here Aw0 for b 1 q 1 m 1 v1, Bw0^Cw0 for b 1 q 1 m 1 v1^R 1 0 v1. Thus the equation (27) has no root which is positive or zero (Descartes' rule of sign). The equation (27) will only have negative roots or complex roots with negative real part if AB{Cw0 (according to Routh-Hurwitz criteria), that is, b 1 q 1 m 1 v1^R 1 0 v1. Thus, the system (1,2,3) without Culex species is stable about the interior equilibrium X 0 1 and the following result holds: Theorem 2. For R 1 0 v1 the model system (1,2,3) without Culex mosquitoes has a unique DFE point which is locally asymptotically stable in V 1 . 0.5.2 Global asymptotic stability of DFE, X 0 1 . To ensure that the disease elimination is independent of the initial sizes of the populations, we need to show that the disease-free equilibrium X Ã 1 is globally asymptotically stable (GAS). This is established using the approach proposed in Castillo-Chavez et al. [38]. There are two conditions that if met guarantee the global asymptotic stability of the disease-free state. First, system (1,2,3) without Culex mosquitoes must be written in the form: where X R m denotes (its components) the number of uninfected individuals and Z R n denotes (its components) the number of infected individuals including latent and infectious. U 0~( x 0 ,0) denotes the disease-free equilibrium of this system.
(H1) For dX dt~F (X ,0), X 0 is globally asymptotic stable (H2) G(X ,Z)~AZ{Ĝ G(X ,Z),Ĝ G(X ,Z) §0 for (X ,Z) V 1 where~D G Z (X 0 ,0) (see [31] for more details) is an M-matrix A are nonnegative) and V 1 is the If the system (28) satisfies the above two conditions then the following Theorem holds.
Then, X Ã 1 is globally asymptotic stable in V 1 0 5V 1 . Proof 3. Global stability of the EE is explored via the construction of a suitable Lyapunov function. Let us consider the following function: V (P 1 ,U 1 ,S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 ) where e i w0 for i~1,2, e and e Further simplification yields S 2 ) 2 m 2 S 2 zF (P 1 ,U 1 ,S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 ) where By theorems hypothesis, where strict equalities holds only when, Furthermore, for all I 1 ,S 2 ,A 2 ,I 2 §0, because the arithmetic mean is greater than or equal to the geometric mean. Thus, F ƒ0 for P 1 ,U 1 ,S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 w0.
Hence, _ V V ƒ0 for all P 1 ,U 1 ,S 1 ,E 1 ,I 1 ,S 2 ,A 2 ,I 2 w0 and is equal to zero for P 1~P and X Ã 1 is the only equilibrium state of the system on this plane. Therefore, the largest compact invariant set in V 1 0 such that _ V V ƒ0 is the singleton X Ã 1 which is the endemic equilibrium point. LaSalle's invariant principle [40] guarantees that X Ã 1 is globally asymptotically stable (GAS) in V 1 0 , the interior of V 1 . 0.6 Stability analysis of the overall model (1,2,3) The overall model system (1,2,3) describes the epidemiological and ecological complexity involved on RVF dynamics. Theorem 2 in van den Driesche and Watmough [32] states that the local stability of the disease-free equilibrium of the model can be determined by its basic reproduction number, R 0 . However, in host-vector models where multiple transmission cycle are observed to occur as in the case of our model (vertical transmission, host to Aedes infection, Aedes to host infection, host to Culex infection and Culex to host infection) the basic reproductive number obtained via next-generation method does not give the number of host infected by a single host if there an intermediate vector, but rather the geometric mean of the number of infections per generation [41]. Therefore, in our case the local stability of the disease -free equilibrium, X 0 , (10) of the model is established through the Routh-Hurtwitz criteria [42,43], and the following result holds.
Theorem 5. The model system (1,2,3) always has the diseasefree equilibrium X 0 . If b1q1 Proof 4. To prove the stability of the equilibrium point X 0 we use the Jacobian matrix (23) of the linearised system, which yield the following characteristic polynomial: where n 1~m3 zẽ e 2 zm 2 ze 2 zm 2 zm 2 zm 1 {b 1 q 1 , n 2~m3ẽ e 2 zm 2 ð Þ 1{c 1 ð Þzm 3 e 2 zm 2 zm 2 ð and n 4 w0 for R 1 0 v1^R 3 0 v1^R 0 v1. Thus the equation (37) has no root which is positive or zero (Descartes' rule of sign). Therefore equation (37) will only have negative roots or complex roots with negative real part if n 3 (n 2 n 1 {n 3 ){n 2 1 n 4 w0 (according to Routh-Hurwitz criteria), that is, b 1 q 1 m 1 v1^R 1 0 v1^R 3 0 v1R 0 v1. Thus, the system (1,2,3) is locally asymptotically stable about the interior equilibrium X 0 . 0.6.1 Existence and uniqueness of endemic equilibrium, X Ã . The existence of the endemic equilibrium in V 1 , is determined by equation (18). Taking A~g 3 l 7 ,B~g 4 l 7 , C~b 1 N 1 g 3 l 6 ,D~m 1 m 2 l 5 and E~b 1 N 1 g 4 l 6 , equation (18) can be written as Solving equation (38) for fI Ã 1 ,I Ã 3 g we get . The existence of positive I Ã 3 is given by the following inequalities: 1 N 1 l 6 l 7 and C{D A~b , we get that the meaningful inequality is Since I Ã 1 w0, then C{D should be positive. C{D is the expression on the numerator of equation (25), which was verified to be positive whenever R 1 0 w1 and b1q1 v1. This gives the  threshold for the endemic persistence. Therefore the following result holds: Theorem 6. The RVF model (1,2,3) has a unique endemic equilibrium point X Ã whenever R 1 0 w1 and The result in Theorem (6) indicates that depending on vertical transmission efficiency, if the Aedes basic reproduction number R 1 0 w1 and I Ã 1 satisfy the inequality it is sufficient to cause an outbreak, since secondary vectors (Culex species) co-exist and serve as disease amplifiers. Figure 2 shows the region where I Ã 3 is strictly positive when varying both I Ã 1 and R 1 0 . That is, in region II both infected Aedes and Culex co-exist while in region I only infected Aedes exist. This confirm the analytical results obtained above. The existence of infected Culex at endemic equilibrium depend on the existence infected Aedes and initial spread of the disease R 1 0 . Thus, Aedes species have the potential to initiate the epidemic through transovarial transmission and the potential to sustain low levels of the disease during post epidemic periods.

Bifurcation and chaos investigation on the RVF model
To provide some numerical evidence for the qualitative dynamic behaviour of the model (1,2,3), time series with both transient and permanent regimes, phase portraits, Poincaré maps, bifurcation diagrams, Lyapunov exponents have been used to assess model sensitive dependence on initial conditions and return maps are used to illustrate the above analytical results and for determining new dynamics as the parameters vary. We start by introducing a simple case of seasonality on time dependent oviposition rates of mosquito populations (Aedes and Culex): where b 1 and b 3 are the baseline parameters of the oviposition rates of Aedes and Culex mosquitoes respectively, T~1 year, d 1 and d 3 are the external forcing amplitudes for the two species of mosquitoes respectively, which represent the strength of seasonality that controls the magnitude of the fluctuations. When d 1~d3 :0, the model reduces to a non-seasonal model and the system possesses two types of equilibria: disease free and endemic equilibria. When the magnitude of the external forcing parameters d 1 ,d 3 is sufficiently small, d 1 ,d 3 (0,1) the system responds with oscillations of the same annual period as external forces (see Figs.3 (a) and (b)). However with larger values (for instance d 1~7 0,d 3~1 :1) the system shows other modes of oscillations (see Figs.3 (c) and (d)) with period 5 as confirmed by Poicaré maps Fig.4. In all this section, the system is integrated numerically with the fifth order Runge-Kutta algorithm [44]. The initials conditions h 2 0.6 [6,9] Probability of an infected host moving to the symptomatic stage, dimensionless (1{h 2 ) 0.4 [6,9] Probability of an infected host moving to the asymptomatic stage, dimensionless s 1 ,s 3 0.33 [5,9] Number of times one Aedes, Culex mosquito would want to bite a host per Day, if it were freely available. This is a function of the mosquito's gonotrophic cycle (the amount of time a mosquito requires to produce eggs) and its preference for livestock blood, Day 21 The maximum number of mosquito bites a host can sustain per Day. This is a function of the host's exposed surface area, the efforts it takes to prevent mosquito bites (such as switching its tail), and any vector control interventions in place to kill mosquitoes encountering hosts or preventing bites, Day 21 b 2a 0.21 [6,9] Probability of transmission of infection from an infectious mosquito to a susceptible host given that a contact between the two occurs, dimensionless, where a~1 and a~3 b a2 0.7,0.15 [6,9] Probability of transmission of infection from an infectious host to a susceptible mosquito given that a contact between the two occurs, dimensionless, where a~1 and a~3 b b a2 0.30 [6,9] Probability of transmission of infection from an asymptomatic host to a susceptible mosquito given that a contact between the two occurs, dimensionless 1=c a 6 [12,15] is the average duration of the mosquitoes latent period, Days, where a~1 and a~3 1=e 2 4 [5,9,16] is the average duration of the infectious period I 2 , Days ∈ and other values are P 1 (0)~1000, U 1 (0)~999, E 1 (0)~0, I 1 (0)~1, S 2 (0)~1000, A 2 (0)~0, I 2 (0)~0, R 2 (0)~0, P 3 (0)~1000, S 3 (0)~5000, E 3 (0)~0, I 3 (0)~0, K 1~1 0000, K 2~2 000 and K 3~1 0000. The parameter values are shown in Table 2. 0.7.1 Time series simulations. Figure 3 depicts the time evolution of the sum of infectious Aedes and Culex mosquitoes, I 1 zI 3 and sum of infectious asymptomatic and symptomatic livestock for different values of d 1~0 :6,d 3~0 :6; d 1~7 0,d 3~1 :1 and d 1~2 4:7,d 3~1 :1. In (a) the number of infectious mosquitoes oscillates yearly reaching the same maximum. In (c) the quantity I 1 zI 3 also oscillates with first peak of above 500 around the second year. In (c) we notice a long lasting peak of about 500 infectious mosquitoes in the interval 18-25 months, which is likely to cause an inter-epidemic outbreak. Fig.3(b) shows a constant low oscillation, high peaks around second and fifth year in (d) and high peaks around second and fourth year in (f). Note that the internal figures describes the permanent regime which represent the dynamics where the system is expected to adapt to the external forcing. The time series for d 1 ,d 3 (0,1), also show that the total of infected vectors I 1 zI 3 and infected livestock A 2 zI 2 stay quite away from zero, avoiding the chance of extinction in stochastic system with reasonable size (see Figs.3 (a) and (b)). This is due to the fact that for d 1 , d 3 (0,1) vector oviposition continues In the region d 1 w1,d 2 w1, Figs.3 (c)-(f), we observe fluctuations in the total number of infected from reasonable small peaks (describing RVF post-epidemic activities) to very low values, which in this case drive almost surely the system to extinction.

0.7.2
Phase portrait diagrams and Poincaré maps. Instead of studying the entire complicated trajectories, important information is encoded in the phase plane. This approach allow us to analyse geometrically the total dynamics of the system. Varying d 1 ,d 3 the state space plots show a rich dynamical behaviour with bifurcations from limit cycles, multiperiodic oscillation to completely irregular behaviour which is usually the fingerprint of chaos (see Fig.4). Poincaré map is a useful tool for analysing the dynamics of a nonlinear system. It allows good insight for global dynamics of the system by displaying the types of attractors of the system [45]. The successive iterations of the map are defined as: seasons. This is not the case of East African region, where we have two rainy seasons (long and short) and a dry season, where intervals of interepidemic periods.
under this former we expect stochastic extinction during some throughout the year, albeit at lower rates during unfavourable The attractor is generated by sampling the system stroboscopically at time corresponding to the multiple of the period T~2p=V. We have used 100,000 points and a period of one year. Figures 5 (a) and (b) with (d 1~0 :6, d 3~0 :6) show that the system is attracted by a limit cycle, because of the presence of a single dot. In this case the system is periodic. In (c) and (d) with (d 1~7 0, d 3~1 :1) we notice a presence of a few dots, thus, the system is multi-periodic and in (e) and (f) with (d 1~2 4:7, d 3~1 :1) we notice a strange attractor which is usually a sign of a chaotic system. 0.7.3 Maxima return maps of I 1 zI 3 ,A 2 zI 2 for state phase plots. We have used maxima return maps in order to get supplementary classification of different dynamics for parameters d 1 and d 3 . For a time selected as t max , at which I 1 zI 3 2 zI 2 have a local maximum, we have plotted the number of infected mosquitoes and livestock respectively at time t max and at the next local maximum t returnmax . Figures 6 (a) and (b) show that all consecutive maxima coincide with themselves as shown by a single dot. In (c) and (d), we notice that consecutive maxima are few and different as a sign of irregularity, and in (e) and (f), we observe that a dot rarely comes back to the same point. Thus, the fingerprint of chaotic attractor is clearly visible now with the maxima return maps analysis.

Lyapunov exponents and bifurcation
diagrams. The largest Lyapunov Exponent (LE) is quantitatively characterized by the average rate of separation of infinitesimally close trajectories in the phase space for a dynamic system. It can be used to determine how sensitive a dynamical system is to initial conditions [46]. In general for a N-dimensional dynamical system described by a set of equations dX i dt~F i (X,t), the LEs are defined by [35]: where l i is the i th LE and EdX i t E is the distance between the trajectories of the i th component of the vector field F at time t. Recall that exponential divergence in the phase space is given by the LEs. If the largest LE is less than or equal to zero, then the system may be regarded as periodic or quasi-periodic. Otherwise, if the largest LE is positive the system may have an irregular or chaotic behaviour. Another important fact to be mentioned is that negative LE does not, in general, indicate stability, and that positive largest LE does not, in general indicate chaos [47].
In Figs.7 (a)-(d) we have computed the bifurcation diagrams with respect to d 1 , the external forcing amplitude on the response of the RVF model. Figures (e) and (f) show the maximal LE after infinitesimal perturbation of 10 {10 in the initial conditions. In Fig.7(e), the maximal LE is positive for d 1 > 60 and around 50 and 25. In Fig.7(f), the maximal LE is positive for 15 = d 1 = 34 and for d 1 > 85. Figure 7 shows the bifurcation diagrams of the local maxima of infectious mosquitoes and livestock undergoing forward forking bifurcation from period-1 to period-6 oscillatory type behaviour. In Fig.7(a), local maxima extrema I 1 of infectious Aedes species undergo irregular behaviour for d 1 > 65, which is the fingerprint of chaos. Fig.7(b) shows irregular behaviour for 15 = d 1 = 34 and d 1 > 85, with large number of periods. In Figs.7 (c) and (d), we observe almost the same qualitative behaviour with the same parameters, but with notable difference in the value of the local maxima of the overall infectious mosquitoes fuelled by the elevation of several secondary vectors which serve as disease amplifiers. When d 3~1 :1 the local extrema A 2 zI 2 undergoes irregular behaviour for 15 = d 1 = 34 and d 1 > 85, with large number of periods Fig.7(h).
We observe from Fig.7(e) that for a fixed d 3~0 :1 and varying d 1 (0ƒd 1 = 62) the largest Lyapunov exponent is fairly negative indicating stable limit cycles and multi-periodicity with some shift to positive values as the system bifurcates through period doubling routes to chaos. Above d 1~6 2 a positive Lyapunov exponent clearly moves away from zero, indicating deterministically chaotic attractors. For a fixed d 3~1 :1 and varying d 1 Fig.7(f) the largest Lyaponov exponent fairly confirms the behaviour seen through bifurcation diagrams with positive values on the chaotic regions. 0.7.5 Interaction between Culex and Aedes oviposition rates. In the preceding section we have fixed the value of d 3 , while investigating the bifurcation behaviour when d 1 is varying. In Figure 8 we have computed the maximal LE when those two parameters are varying. For 20 = d 1 ,d 3 = 100, the maximal LE is negative, then the system is sensitive to initial conditions. For low values of d 3 and 18 = d 1 ,d 3 = 45, the maximal LE is positive. Another remarkable fact is observed when d 1 is around 10 no matter the value of d 3 , the maximal LE will be positive. This shows i fact that, Aedes oviposition rate is predominant in leading Recall that in certain Aedes species of the subgenera Neomelaniconion and Aedimorphus, the female mosquitoes transmit RVF virus vertically to their eggs [3]. When these mosquitoes lay their eggs in flooded areas, transovarially infected adults may emerge and transmit RVF virus to nearby domestic animals which may then lead to the infection of secondary arthropod vectors species including various Culex [48]. Thus, there is an initial quantity of primary infected vectors required to trigger an outbreak. Fig.8 shows that if the control magnitude of fluctuations in Aedes oviposition rate is around 10, and the number of newly transovarially infected mosquitoes is amplified by nearby domestic animals, then, the number of infected (in both host and vector) will be sufficiently enough to cause subsequent elevation of secondary vectors, including Culex species, and consequently trigger an outbreak.

Discussion and Conclusion
The proposed model accounts for the population dynamics of both livestock and mosquitoes (Aedes and Culex) and seasonal changes in weather that heavily affects the vector population size. Mosquito density varies over seasons, and the contact rates and vector oviposition rates vary dynamically based upon both host and vector densities since female mosquitoes need blood for oviposition. Qualitative analysis of the model showed that there exists a domain where the model is epidemiologically and mathematically well-posed. We then analysed the existence and stability of both disease free and endemic equilibria.
Dynamical analysis shows that when R 0 v1, then the disease dies out and when R 0 w1 the disease become endemic. A suitably constructed Lypunov function is used to determine the global stability of the endemic equilibrium of the model without Culex species and the existence of the endemic equilibrium of the overall model is seen to exist whenever meaning that the co-existence of the infectious host, Aedes and Culex mosquitoes is subjected to the number of infected Aedes mosquitoes.
We have used visualisation techniques to study the behaviour of RVF epidemic model under external forcing in the mosquito oviposition rates. The bifurcation diagrams show the emergence with increase in external forcing parameters d 1 ,d 3 of Hopf and pitchfork modes of bifurcation. That they have much larger amplification of infection levels that can take place if the system is encouraged to switch to multi-periodic mode. In transition, further amplification can occur if the multi-periodic mode becomes unstable and the system moves into chaotic state before finding an alternative stable periodic mode (e.g. Fig.7).
On the bifurcation diagrams the highest maximum number of infectious Aedes mosquitoes is only observed for values of d 1 (d 1 v10) with different values of d 3 , meaning that for the disease to trigger an inter-epidemic a certain number of infectious Aedes mosquitoes is necessary. This confirm the analytical results obtain in section 0.6, as well as results obtained in [9] which showed that when mosquito populations follow seasonal patterns with large amplitudes, vertical transmission could play a significant role in long-term persistence of a pathogen. Another important conclusion is that even with a low maximum number of infectious individuals, the bifurcation diagrams show that if for fixed d 3~1 :1 and varying d 1 the system becomes chaotic in the interval 15 = d 1 = 35, meaning that unpredictable and possibly uncontrolled low levels of inter-epidemic activities may occur, leading to higher morbidity in livestock. Hence observed fluctuations in RVF outbreak data and non deterministic nature of RVF interepidemic activities could now be better understood considering fluctuations on both rainy and dry season as significant factor.
A sero-survey study done in livestock approximately four years after the 2006/07 RVF outbreak in Tanzania, showed a linear increase in seroprevalence in the post-epidemic annual cohorts implying a constant exposure and presence of active foci transmission [10]. Figure 3 (d) and (f) demonstrate this behaviour which is shown to come in cycles of 5 to 7 years approximately, as well as fluctuations in the total number of infected from reasonable small peaks (describing RVF post-epidemic activities) to very low values. During these periods of low troughs for the total number of infected, the virus survive through vertical transmission in Aedes species and among wild animals as reservoirs [49]. Note that, this recurrent low level RVF virus activity during inter-epidemic periods, in East African region in particular, infects 1{3% of livestock herds annually [50]. Generally, these infections pass undetected where there is no regular active surveillance in the livestock and human populations [10]. This suggests that RVF outbreaks partly result from build up RVF inter-epidemic activities for it has been observed that optimum climatic conditions (temperature and rainfall) only and presence of mosquitoes can not completely explain the RVF outbreaks [51].
Simulation of the interaction between the two populations densities of Aedes and Culex by varying the magnitudes of external forcing d 1 and d 3 of the oviposition rates b 1 3 have opened a new window of research on the potential of Aedes species to initiate RVF outbreaks and sustain low endemic levels of the disease during inter-epidemic periods. This result concurs with the Chitnis et al. [9] suggestion that vertical transmission is required for interepidemic persistence.
One of the main objectives of this study was to investigate the possibility of prediction of RVF outbreaks with the aim of controlling RVF incidence. We have shown that seasonality may induce irregular behaviour on the disease dynamics. It has been shown that the interaction between oviposition rates of Aedes and Culex mosquitoes makes prediction more complex. In fact, higher irregularity are naturally expected in the higher seasonality forcing. However, our proposed model has shown that the complexity occurs even for a relatively low level of the magnitude of seasonal forces. We have also found that seasonal Aedes oviposition rate is most likely to generate uncontrollable behaviour than Culex seasonal oviposition rate. This study is of great epidemiological significance as it highlights a high uncertainty in RVF outbreak prediction by a simple theoretical mathematical model including seasonal influence in mosquito populations. In addition, the model including external seasonal forcing on mosquito oviposition rates shows ability to mimic the linear increase in livestock seroprevalence as reported in Sumaye et al. [10], with first post-epidemic peak around the second year, a following peak larger than the previous one around the fifth year (see Fig.3 (d) and (f)).
Currently, two types of RVF vaccine for animals exist: a live vaccine and inactivated vaccine. However, the current live vaccine can not be used for prevention and prevention using the inactivated vaccine is almost impossible to sustain in RVF affected countries for economic reasons [6,21,52]. Then, the possible alternative of controlling RVF transmission remains in keeping the vector population at the lowest levels. Therefore we argue that locations that may serve as RVF virus reservoirs should be eliminated or kept under control to prevent multi-periodic outbreaks and consequent chains of infections. We also recommend a systematic surveillance in the livestock or human population in order to monitor inter-epidemic RVF activities.
This study is not exhaustive and can be extend to include humans not just as dead ends [18] but also as disease amplifiers since it has been demonstrated that humans have potential to transmit the virus, particularly to Aedes mosquito species [8]. Also, including ticks on the model may help to explain and gain more insights on the understanding of disease dynamics and enhance control strategies, since ticks have been reported to play a role on disease transmission [51]. For mathematical convenience and tractability of the model, we made several assumptions, thus our results are driven by the model formulation and structure. A step toward a more quantitative and qualitative study is viable by relaxing some of the assumptions made and incorporating more epidemiological features of the disease as well as the use of a double periodic function and inclusion of stochasticity in order to capture the dynamic of the two rainfall seasons in East Africa (long and short rainy seasons), where the disease is likely to be more predominant. Further studies are needed to enhance the understanding of RVF epidemic and inter-epidemic activities in order to provide further insights in assessing the current and future control strategies.

Supporting Information
Appendix S1 A1. Computation of the basic reproduction number. A2. Computation of the basic reproduction number in periodic environments. (PDF)