Rift Valley fever: An open-source transmission dynamics simulation model

Rift Valley fever (RVF) is one of the major viral zoonoses in Africa, affecting humans and several domestic animal species. The epidemics in eastern Africa occur in a 5-15 year cycle coinciding with abnormally high rainfall generally associated to the warm phase of the El Niño event. However, recently, evidence has been gathered of inter-epidemic transmission. An open-source, easily applicable, accessible and modifiable model was built to simulate the transmission dynamics of RVF. The model was calibrated using data collected in the Kilombero Valley in Tanzania with people and cattle as host species and Ædes mcintoshi, Æ. ægypti and two Culex species as vectors. Simulations were run over a period of 27 years using standard parameter values derived from two previous studies in this region. Our model predicts low-level transmission of RVF, which is in line with epidemiological studies in this area. Emphasis in our simulation was put on both the dynamics and composition of vector populations in three ecological zones, in order to elucidate the respective roles played by different vector species: the model output did indicate the necessity of Culex involvement and also indicated that vertical transmission in Ædes mcintoshi may be underestimated. This model, being built with open-source software and with an easy-to-use interface, can be adapted by researchers and control program managers to their specific needs by plugging in new parameters relevant to their situation and locality.


Introduction
Rift Valley fever (RVF) is caused by the Rift Valley fever virus (RVFv), which belongs to the genus Phlebovirus in the family Bunyaviridae. RVF is one of the major viral zoonoses in Africa, affecting man and several domestic animal species [1,2].
A syndrome compatible with RVF was first described in the Rift Valley of Kenya in the early 1900s and the virus was isolated in the 1930s [3]. The known range of RVFv is shown in Fig 1. RVF was confined to eastern and southern Africa until about 1975. Since then it has expanded its range first to Egypt (1977), then to western Africa (ca. 1980) and finally to the PLOS ONE | https://doi.org/10.1371/journal.pone.0209929 January 9, 2019 1 / 27 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 linked to other sources of flooding, e.g. the construction of a hydroelectric dam along the Senegal river [12,13]. In the past, the above was the traditional view of the epidemiology of RVF, but recently there is more and more evidence of so-called inter-epidemic transmission: previously unnoticed low-level viral transmission in all species involved [12,[14][15][16][17][18]. In Tanzania, human involvement in RVF inter-epidemic transmission has been reported in the past [19,20]. During the 2006/07 RVF epidemic in eastern Africa, livestock and people in the Kilombero valley in Tanzania were affected [21]. Two serological surveys in this region since this last epidemic, one in livestock and one in people, effectively showed the presence of inter-epidemic transmission in the area [17,22].
RVF is transmitted to humans and other mammalian hosts, both livestock and wild ruminants (e.g. cattle, buffalo, sheep, goats and camels) through mosquito (e.g. Culex spp., AEdes spp. and Mansonia spp.) and other arthropod vector bites [1,2,16,23]. AEdine mosquitoes are capable of transovarial (= vertical) transmission of RVFv to the eggs, which can survive long droughts (several years) and hatch when new water arrives during e.g. the ENSO phenomenon, resulting in infected larvae and adult mosquitoes [2]. The highest risk for humans to become infected is through direct and indirect contact with infectious animal materials (blood, body fluids or tissues of viraemic animals). AErosol formation during e.g. milking or consumption of raw milk, meat or blood form another risk for transmission [13,[24][25][26][27][28]. An established treatment method or a vaccine for humans currently does not exist. Control of the disease needs to be done through vaccination of livestock and preventive measures by humans [29,30].
Clinical manifestation in humans can go from only mild illness, including fever, muscle pain, joint pain and headache to severe forms with ocular disease, meningo-encephalitis or haemorrhagic fever [29,31]. The disease manifests itself in livestock through morbidity and mortality in newborns and abortions during all stages of the pregnancy. This has devastating effects on livestock populations and has severe economic repercussions for livestock keepers [2,26,32,33].
Quantitative analysis and simulation modelling of RVFv dynamics have been undertaken on several occasions. Note that the list that follows cites only typical examples and that many more publications exist dealing with RVF modelling. The analytical models use environmental characteristics and range from post-hoc predictions of where outbreaks were to be expected during the 2006-2007 epidemic in East Africa [10] over statistical modelling in order to identify landscape features related to RVFv transmission [34] to the identification of ranges of potential vectors [35]. Simulation models include temporal models using differential equations [36] with extensions to spatial components [37]. Risk analysis of introduction into new territory (in casu The Netherlands) [38] has also been carried out. An overview of compartmental models, applied to the simulation of RVF dynamics, is provided by Danzetta and colleagues [39].
The existing models all suffer from being closed, inaccessible and specialised. The combination of R/RStudio 1 with the libraries shiny and deSolve offers the possibility to develop open-source, easily applicable, accessible and modifiable models that can, on the one hand, be adapted to a specific situation with minimal programming effort and, on the other hand, be perused by the epidemiological researcher to study different scenarios and/or the effects of different parameter settings. The model presented here has been developed for the specific situation in East Africa, but as explained above, it can easily be adapted to other areas/situations, mostly by switching on or off certain parameters or parameter groups or by the inclusion of extensions with minimal new coding. The model presented in this paper is thus to be considered a research tool, allowing the user to study the effect(s) of different scenarios in order to better understand RVFv transmission dynamics and the mammalian hosts and arthropod vectors involved, and ultimately to assist in the formulation of new research questions. The model is not a predictive tool, as too much uncertainty still exists with regards to the actual dynamics of inter-epidemic transmission of the virus.

Model-General description
The model describes the RVFv transmission dynamics in six species (human population, domestic animal population and four vectors) in three different areas. The model attempts to offer maximal flexibility, whilst remaining manageable. The model allows for migration of the various species between the different areas. The different compartments in the model are presented in Table 1 and a simplified schematic representation of the model is shown in Fig 2. Each human and animal population consists of a susceptible S, exposed E, infective I and removed R (= recovered/immune) compartment. There is a flow back from the removed to the susceptible compartment in both populations, i.e. immunity is not lifelong. All individuals are born susceptible and a proportion of the pregnant infected animals abort. Vectors A and B allow for vertical transmission: infected females (I compartment) transmit infection to their eggs (Q compartment), where the virus survives until the larvae hatch and the resulting adults are infective. Vector A furthermore has the possibility of long-term dormancy in the egg stage (both infected and non-infected).
A challenge lies in the correct modelling of the vector dynamics. More specifically, a point of attention is the distribution of feeding individuals over the different host populations (both species-wise and zone-wise). Vectors can feed on the two modelled host species (human and domestic animal), but they can also use alternative hosts (especially so in the forest zone). The latter means there is no increased mortality in case the two main hosts are not available, but this of course also influences infection prevalence in the vector population. The vector populations are furthermore limited by a density-dependent oviposition rate. The approach currently taken uses the following basic parameters (see Vector feeding and infection rates for details): • ε: proportion of vector X feeding on host Λ in zone i; it is the user's responsibility to ensure that the sum of the various ε per species per zone does not exceed one  • k j X : maximum number of vector X individuals in zone j ('carrying capacity') El Niño events are currently modelled to occur every ten years. Additionally, the user is given the opportunity to include annual overall climate variability through the choice of a random series of 'dry' or 'wet' years and a seasonal within-year variation in egg eclosion to model seasonal effects on vector population size. Finally, there is the possibility of including a 'fixed' annual domestic animal movements between zones 1 and 2, simulating seasonal transhumance of (e.g.) cattle between the plateau and the floodplain. Details are to be found in Seasonality and El Niño effect.  Table 1 for a list of all compartments. It is understood that the necessary calculations for these density-dependent ovipositionfeeding, climatic variability and transhumance processes slow down the model considerably. It was therefore decided to rewrite part of the code, doing the preparatory computations before calling the deSolve routines (using the classical Runge-Kutta 4 th order method), in C++ (making use of the RCCP library). This speeds up execution by a factor of about sixty, but of course means a lower accessibility of the code. Therefore, a slower version, entirely written in R is also offered. Full details on how to install and run the model are given in the accompanying user's manual S1 Appendix. The R and C++ code is provided in S2 Appendix.

Model-Differential equations
For every zone i(i = 1, 2, 3), we compute the differential equations of each compartment of the human, the animal and the vector populations.

Human population
Eq 5 describes the rate of change in the susceptible animal host compartment: refers to the newborn individuals, respectively born from uninfected and infected individuals and corrected for population density to simulate removal (sales) in function of herd size, other two zones and individuals losing their immunity and losses through natural mortality, animals becoming infected and emigration out of Zone i. Eq 6 describes the rate of change in the animal host exposed (incubating) compartment in Zone i: b i M M i S refers to the animals becoming infected, to the losses through natural mortality, changing from incubation to the infective stage and emigration from Zone i. Eq 7 describes the rate of change in the animal infective compartment in Zone i: x M M i E refers to the individuals becoming infective, to the losses through natural mortality, disease-specific mortality, recovery and emigration from Zone i. Eq 8 describes the rate of change in the recovered (immune) animal compartment in Zone i: a M M i I refers to the animals having recovered (gained immunity), to losses through natural mortality, loss of immunity and emigration from Zone i.
Eq 9 describes the rate of change in the infected-egg compartment of Vector A in Zone i: to the production of infected eggs (product of total biting rate, egg production rate, density-dependent correction and vertical transmission rate) while ðm i A Q þ t s t i A ÞA i Q refers to losses through mortality and hatching (in function of El Niño and seasonal flooding through τ s ). Eq 10 describes the rate of change in the uninfected-egg compartment of Vector A in Zone i: refers to the densitydependence corrected production of uninfected eggs both by infected adult vectors (absence of vertical transmission) and uninfected adult vectors while ðm i A P þ t s t i A ÞA i P refers to losses through mortality and hatching (in function of El Niño and seasonal flooding through τ s ). Eq 11 describes the rate of change in the uninfected-adult-vector compartment in Zone i: t s t i A A i P refers to the newly 'hatched' adults (note that stages intervening between egg and adult are omitted, requiring adjustment of hatching and mortality rates), The differential equations describing the dynamics of Vector B are identical as those for Vector A, the only difference being the possible presence of dormant eggs in the latter and not in the former.

Vector C
Vector D is identical to Vector C.

Population totals
Vector feeding and infection rates  are the basic parameters used to compute carrying capacity etc. of a zone vis-à-vis its resident vectors. The present approach is to compare the total number of bites (successful feedings, . . .-for sake of brevity referred to as 'bites' from now on) the vectors can inflict upon the hosts per time unit with the total number of number of vector bites the host populations can sustain (given their resistance, evasive behaviour, . . .). The minimum value of these two is used to compute the actual number of bites given per vector and/or the number of bites suffered per host. It is understood that this approach may introduce a number of parameters whose values are only vaguely known at best, but an attempt was made to avoid unrealistic numbers of vectors interacting with a single host, i.e. host numbers determine vector numbers. At the same time, the possibility is offered to include so-called alternative hosts, which can be used by the vectors when the hosts included in the model are insufficient, in order to avoid vectors disappearing when host population levels are too low.
Parameters 36 and 37 are computed from the simulation output: The potential maximum number of vector bites (all vector species) on whole host population Λ j is computed as: This is compared with the maximum number of bites the same host population can 'sustain' (see above for more details): The 'availability' of host population Λ j (i.e. the proportion of the potential bites actual inflicted on the host population in question) is the ratio of parameter 39 over parameter 38 with a maximum of unity: The actual number of bites by vector X k on the whole host population Λ j is thus: The individual biting rate of vector X k on host Λ j per time unit becomes: The total individual biting rate of vector X k on all host populations per time unit therefore is the sum of the respective ω kj : The biting rate of vector X k on alternative hosts (with O alt = number of alternative hosts) is defined as: The proportion of infection in vector X k feeding on all modelled hosts species is computed as (the reference to the zone is left out, I L j being the number of infective individuals of host Λ j ; β wl refers to the infection picked up from game animals and it is added only in the case of Zone-3-dwelling vectors): The infection rate of host Λ j being subjected to the actual number of bites by the various vectors and/or interacting with other infectious hosts is calculated as (φ j 0 , j refers to the number of transmitting hosts [domestic animal] met by one receiving host [a person] per time unit; The second and third terms of the logarithm function of Eq 46 are currently implemented only for animal-to-human direct transmission.

Seasonality and El Niño effect
Simulating an annual (seasonal) animal transhumance between Zone 1 and Zone 2 is possible: animals move to Zone 1 on day d 1 and move back to Zone 2 on day d 2 . This is achieved through the generation of 0/1 indicators, which are to be multiplied with the movement rate: Hatching of dormant eggs of Vector A can be regulated on a seasonal basis as well as periodically through El Niño events in Zone 1 (d 3 and d 4 are respectively the start and end of the annual flooding, π φ is the proportion proportion of Zone 1 that is seasonally flooded; d 5 and d 6 are respectively the start and end of the El Niño event): Annual variation (e.g. because of wet and dry years) and seasonal variation in vector egg eclosion (τ S ) in all three zones can be included in the model: the current approach is by penalising hatching rates during dry years (hatching rate becomes a fraction -π δ -of normal rates) and by allowing hatching rates in normal and dry years to vary seasonally according to a cosine curve (see the accompanying user's manual S1 Appendix for examples on different parameter settings). The different possible combinations are as follows in Table 2:

Model-Calibration
The model is calibrated using data that were extracted from two studies in the Kilombero Valley in Tanzania (Morogoro region, [17,22]: the principal findings of these studies were the presence of inter-epidemic RVFv circulation in human and domestic animal populations and the location of so-called infection 'hot-spots' away from the floodplain and in fact closer to forested areas on the plateau. The Kilombero Valley region consists of a seasonally inundated floodplain between the densely forested escarpments of the Udzungwa mountains to the northwest and the grass covered Mahenge mountains to the southeast. The valley receives an average annual rainfall of 1200-1800 mm and the average monthly temperature ranges between 25℃ and 32℃. The valley has a diverse ecology and demography with villages consisting largely of numerous distinct groups of houses located on the margins of the floodplain where rice cultivation is the predominant economic activity. Other land use types include hunting, fishing, forestry, pastoral livestock rearing and cultivation of other crops. Several mosquito species inhabit the valley, including known vectors of RVFv, such as Culex spp., AEdes spp. and Mansonia spp. [17,22,40]. The zones, the two mammalian hosts and the four vector populations modelled are in this case: • Areas   Parameter values (ranges) for this scenario are given in Tables 3, 4 and 5. The model was run for 27 years, thereby modelling three El Niño events (years 1, 11 and 21) allowing the model to reach quasi-equilibrium conditions and generating output six years after the last ENSO, which could be compared with the observations made during the field studies [17,22].

Results
The graphical output (showing results for the years 20-27) for the simulations over a period of 27 years, using the standard parameter values as shown in Tables 3-5 Table 6.

Discussion
A model on RVFv transmission in the Kilombero valley in Tanzania was run for 27 years to include three El Niño events (and thus three RVF epidemics), to allow the model to reach a state of 'equilibrium' and to allow model output during a period of 4-7 years after the epidemic to coincide with published observations [17,22]. The model is a complex interaction of density-dependent birth, death and transmission processes and as such very sensitive to certain parameter values. The model was explored by means of scenarios and no attempt was made to include a sensitivity analysis.
Most parameters could be kept at values within the ranges found in the literature, by adjusting the values of other parameters to acceptable values, based on expert opinion. In this respect, a major influence is exerted by ν, the maximum number of bites 'supported' by an individual host. The value itself directly determines the (e.g.) seroprevalence levels, but this parameter also introduces a competition between the various vector species, as at present it is assumed that the 'available' bites are distributed proportionally between the different vectors.
The effect can be seen in Table 6, when comparing lines one and (e.g.) nine: Culex on its own, being a more efficient vector, yields higher seroprevalence values than the standard setting, where it must share the biting opportunities with AEdes.
The exception to the above was the vertical transmission rate (trans-ovarial transmission rate) for AE. mcintoshi. The range found in [50] (0-8.5%) is not sufficient to carry the virus from one epidemic to another in the absence of other vectors to ensure inter-epidemic  [42,43] π MB p_mb Probability to transmit infection from bovine to AE. aegypti 0.89 (77-100%) [42,43] π MC p_mc Probability to transmit infection from bovine to Culex sp1 0.81 (78-84%) [42,43] π MD p_md Probability to transmit infection from bovine to Culex sp2 0.81 (78-84%) [42,43] π MH p_mh00 Probability to transmit infection from bovine to people 0.001 user-defined    [42,43] π DM p_dm Probability to transmit infection to bovine upon Culex sp2 bite 0.07 [42,43] https://doi.org/10.1371/journal.pone.0209929.t005  (Table 6). Using the standard values, predicted seroprevalence levels in humans and cattle at different times after the El Niño event were comparable to those observed. Seroprevalence is estimated to be 13.2% in people and 12.3% in cattle, six years after an El Niño event. The field studies found similar overall seroprevalence levels of 11.7% in people and 11.3% in cattle, five to six years after the 2006/07 RVF epidemic in the area [17,22]. The results are also in line with previous studies across Africa with evidence of inter-epidemic transmission of RVF [1,15,16]. The dynamics of levels of seroprevalence are of course in the first place dependent on the value employed for the loss-of-serotitre rate: currently a daily value of 1/900 is used, based on a single, rather vague  reference [41]. Inclusion of a wildlife reservoir ( Table 6, second line) did not have a significant effect on the predicted levels of seroprevalence.
The simulated seroprevalence levels in Table 6 in both the human and livestock populations show a gradual decline during the years after an epidemic event (El Niño), which seems to imply low numbers of infective bites during inter-epidemic periods, reflecting the generally low numbers of mosquitoes in the absence of heavy rainfall associated with the El Niño events. People and cattle transiting in the forest (zone 3, Figs 5 and 8) are exposed to infectious bites every year from the AE. aegypti and Culex sp.2 populations (Figs 11 and 14): the mosquitoes are constantly infected from the wildlife reservoir [56]. People and cattle remaining in the villages   Figs 3 and 6) are minimally exposed on an annual basis with high exposure rates occurring only every ten years (Figs 9, 10, 12 and 13). Infection thus principally spreads to the villages and floodplains by humans and cattle temporarily residing in the forest zone.
The AE. mcintoshi population in the floodplains (Fig 9) is the one maintaining the infection inside the dormant eggs. Adult mosquitoes do not survive the drier period following the El Niño event and only some eggs hatch every year during the partial seasonal flooding of the plain. Substantial hatching occurs during flooding related to the El Niño event in the East African region, releasing the infection and starting the epidemics. The infection is picked up by Culex sp.1 present in this area. The human population acquires the infection first, followed by the cattle population. From there on, the epidemic spreads to the village and the forest with migrating cattle and people.
As indicated by lines three and four of Table 6 (with the current standard parameter settings), AE. mcintoshi on its own is not able to explain the high seroprevalence found in both humans and cattle [17,22], not even when including annual partial flooding of zone 1 accompanied by eclosion of part of the dormant eggs. The same can be said for AE. aegypti, despite it being resident in the village and forest zones, although it must be understood that in this case the low values for vertical transmission were maintained.  Table 6 examine different scenarios with an efficient Culex vector in the village and forest zones. Introduction of infection, either by means of a wildlife reservoir (line seven) or through the introduction of an infective animal, allows for maintenance of the infection within the host and vector populations. Because of the interaction between the different vectors for host-feeding opportunities, the more efficient Culex vector on its own (without competition from Aedes species) results in higher infection transmission and higher seroprevalence levels. Again, a lot more detailed observations are required to properly quantify this aspect of the transmission dynamics.
Mosquito species in the forested environment (AE. aegypti and Culex sp.2) (Figs 11 and 14) had high annual infection rates. On the other hand, mosquitos in the residential area (AE. aegypti and Culex sp.2) and in the floodplain (AE. mcintoshi and Culex sp.1) have low infection rates (Figs 9, 10, 12 and 13) with peak rates occurring only during or immediately after an El Niño event and subsequent RVF epidemics in the East African region [57].
The model presented here needs further calibrating with datasets from other regions where there are similar or dissimilar ecologies compared to our study area in order to extend and/or improve usability of the model in different geographical, climatic settings. This model, being built with open-source software and with an easy to use interface, can be adapted by researchers and program managers to their specific needs by plugging in new parameters relevant to their situation and locality. Its use can be further expanded by including disease prevention and control interventions to model potential impact of these veterinary and public health measures on disease in people and domestic animals, for example vaccination, quarantining and vector control programs.