MVSE: an R-package that estimates a climate-driven mosquito-borne viral suitability index

Background Viruses such as dengue, Zika, yellow fever and chikungunya depend on mosquitoes for transmission. Their epidemics typically present periodic patterns, linked to the underlying mosquito population dynamics, which are known to be driven by natural climate fluctuations. Understanding how climate dictates the timing and potential of viral transmission is essential for preparedness of public health systems and design of control strategies. Results We developed the Mosquito-borne Viral Suitability Estimator (MVSE) software package for the R programming environment. The package estimates a suitability index based on a climate-driven mathematical expression for the basic reproductive number (R0) of a well established mathematical model for the transmission dynamics of mosquito-borne viruses. By accounting for local humidity and temperature, as well as viral, vector and human priors, suitability can be estimated for specific host and viral species, as well as different regions of the globe. Conclusions Here, we describe the background theory and biological interpretation of the new suitability index, as well as the implementation, basic functionality, research and educational potentials of the MVSE R-package. The package is freely available under the GPL 3.0 license.


32
Recently, we have developed a climate-driven, mathematical model of mosquito-borne viral 33 transmission based on ordinary differential equations, which has been used to study the 34 entomological and epidemiological determinants of the 2012 dengue outbreak in the island of 35 Madeira (Portugal) (Lourenço and Recker, 2014), the 2014 dengue outbreak in Rio de Janeiro 36 (Brazil) , and the 2015-2017 ZIKV outbreak in Feira de Santana (Brazil) 37 . In these case studies, the model was fit to notified epidemic curves using 38 Markov chain Monte Carlo, allowing the estimation of epidemiological parameters such as the basic 39 reproduction number (R0) and the effective reproduction number (Re). The success of this 40 framework stems from its data-driven approach, including temporal climatic series used to 41 parameterise both vector and viral variables under mathematical relationships derived from 42 experimental studies. The framework includes estimation of the effect of local climatic variables on 43 viral and vector parameters, which enables it to be generalized to regions with different climates.

44
Furthermore, by allowing specific vector, viral and host parametrization, this framework is also 45 general enough to be applied to different host-pathogen systems.

46
In this article, we translate our experience with this modelling framework into the MVSE 47 R-package, which is capable of estimating and analysing a mosquito-borne viral suitability index.

48
Given local climatic variables, this index can be used as a proxy for the timing and scale of 49 transmission potential, and is applicable to different viruses, vectors and regions. The main goal of 50 this package is to offer a software tool free of the mathematical complexity of epidemiological 51 models and of the need to depend on epidemiological time series (incidence) of the virus under 52 study. Here, we describe in detail the package's functionalities using real world examples and argue 53 for its public health and academic potential. Recently, we have shown that a suitability measure called index P can be derived from the formula 56 of R0 without the need to fit epidemic curves, as long as the relationships of some of R0's 57 parameters with climate are known. The index P was shown to be correlated with the timing and 58 scale of both local mosquito infestation measures and dengue notifications in Myanmar (Pn et al., 59 2017). Here, we describe the rationale behind the index P in detail, and demonstrate its potential in 60 other epidemiological contexts.
There are a total of eight parameters in the expression of R0, four of which are climate-independent 73 (the human life-span 1/µ h , the transmission probability from infected human to mosquito per bite 74 φ h→v , the human infectious period 1/σ h , and human incubation period 1/γ h ), and four are 75 climate-dependent (the life-span of adult mosquitoes 1/µ v (u,t) , the extrinsic incubation period 1/γ v (t) , 76 the daily biting rate a v (u) and the probability of transmission from infected mosquito to human per 77 bite φ v→h (t) ). Climate-dependent parameters are defined as functions dependent on humidity (u) and 78 temperature (t), which have previously been determined in experimental studies through laboratory 79 estimates of entomological data under various climate conditions. A list of these functions, sources 80 and details can be found in the Supporting Information S2 Text. partially driven by M, and are also determined by P (u,t) (equation 1). In our mathematical 87 framework, P (u,t) 's contribution to this oscillatory behaviour stems from the four climate-driven 88 3/16 parameters: the life-span of adult mosquitoes, the extrinsic incubation period, the daily biting rate, 89 and the probability of transmission from infected mosquito to human per bite.

90
Theoretically, the potential for outbreaks is determined by the epidemic thresholds of R0>1 or 91 Re>1. P (u,t) is a positive number that can be interpreted as the absolute potential of an adult 92 female mosquito. Thus, if at least one female mosquito exists per human (M >= 1) and P (u,t) > 1, 93 then R0 = M P (u,t) > 1 and epidemic growth is possible. However, interpreting P (u,t) on its own 94 means that no direct assessment on the classic epidemic thresholds can be made. We thus do not 95 make particular interpretations of P (u,t) > 1, but instead argue and show that its absolute value is 96 informative for transmission potential when assessed locally in time or between regions. 97 P (u,t) can be estimated for any region for which humidity and temperature are available, and can be 98 parameterized for any species of virus, host or vector. We thus define P (u,t) as the mosquito-borne 99 viral suitability index P, using a complementary terminology to existing vector suitability indices 100 (Gardner et al., 2018;Kraemer et al., 2015a,b). A key difference from vector suitability indices is 101 that the numerical scale of the index P has direct biological interpretation: P is the reproductive 102 (transmission) potential of an adult female mosquito.

104
In previous studies on ZIKV and DENV, we estimated R0 and Re by fitting a dynamic 105 transmission model to notified epidemic curves within a Bayesian Markov chain Monte Carlo 106 (bMCMC) framework, whereby unobserved parameters were estimated by their resulting posterior 107 distributions Lourenço et al., 2017a;Lourenço and Recker, 2014). Here, 108 however, our goal is to estimate suitability for transmission, through P (u,t) , independently of the 109 availability of epidemiological data. (e,f) Estimated index P per day with mean (red) and 95% confidence interval (grey). Priors were assumed to be the same for Recife and São Paulo. See Supplementary Information S2 Text for prior distributions.
almost every year (Faria et al., 2017a,b). In the notified dengue data that we obtained for  Table), the total incidence per 100k individuals in Recife was approximately 11-fold 130 higher than it was in São Paulo (2184.3 notifications versus 195.2).

131
The input (climate time series, user-defined priors) and index P output for the two cities is 132 summarized in Figure 1 (for a more complete set of possible outputs see S1 Text). Here, for 133 simplicity, and for a lack of evidence of significant differences between the two cities, we assume all 134 informed priors to be the same (Figures 1 b, Paulo shows pronounced oscillations, with P presenting substantial troughs.

142
These key differences in index P are in accordance to what is known about the transmission 143 potential of mosquito-borne viruses in the two cities. For instance, the stabler seasonality patterns 144 and higher index P in Recife suggests year-round potential for single adult mosquitoes to contribute 145 to epidemic expansion; and can thus help explain why dengue incidence in Recife can be one order 146 of magnitude higher than in São Paulo. Concurrently, the substantial troughs of index P estimated 147 in São Paulo suggest very low transmission potential off-season, which may partially explain why 148 persistence of mosquito-borne viral lineages between seasons is rarely observed in the city.

5/16
Index P seasonality in Recife and São Paulo 150 We next present two visual outputs from MVSE, comparing some of the seasonal differences in the 151 estimated index P between Recife and São Paulo (for a more complete set of possible outputs see S1 152 Text). Figure 2 shows the timing of highest index P and sensitivity of P to each of the climatic the week with highest suitability is identified across all estimated index P solutions of that year (solutions' mean and 95% CI in Figure 1). The frequency of each peak week shown across the 1000 simulations. (c,d) Two dimensional sensitivity of mean index P per humidity (x) and temperature (y) observed time point (day). Each point is a combination of observed climate variables, colored according to the mean index P estimated (color scale on the right). The white dots (over the black link) mark the mean humidity and temperature of each month over the period of the data (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016); the floating circles with numbers identify each month's number.
For each year and city, the week with highest suitability is identified across all estimated index P 156 solutions (simulations, of which the mean and 95% CI is in Figures 1 e,f). As seen in Figures 2 a,b, 157 the two cities vary in their detected peak weeks. In Recife, suitability tends to peak between March 158 and May (regional autumn), with the occasional occurrence of higher suitability in June (in 159 particular in the years 2012 and 2013). In contrast, suitability peaks earlier in São Paulo, between 160 late December and early February (regional summer).

161
MVSE also offers visual output which highlights the regional sensitivity of the index P to climatic 162 changes in time (Figures 2 c,

7/16
The estimated mean suitability presents high variation across Brazil (Figure 3 a). As expected, the 180 southern regions present the lowest transmission potential (colder colors), while the highest potential 181 is estimated in the center and along the northern coast (warmer colors). The elevated regions of the 182 country, running in parallel to the eastern coast, present a pattern of intermediate suitability. 183 We further identify the month of peak suitability for each region (Figure 3 b). This discretization 184 highlights a seasonal gradient in the month of maximal suitability across Brazil. In the southern 185 regions, peak suitability is during the early summer months, and moving north along the eastern 186 and then northern coasts shifts peak suitability further into regional autumn. Notably, no region of 187 Brazil is seen to experience peak suitability during the core winter months (July -September, 188 darker blue).

189
We also map the mean monthly index P for a subset of months (Figure 3 c). In these maps, the correlation (ρ) shows that the index P is highly and positively correlated with dengue notifications 202 in Brazilian cities. These results are similar to those we have previously demonstrated for Myanmar 203 at the country and city level (Pn et al., 2017).

205
In this manuscript, we have introduced in detail a mosquito-borne viral suitability measure which 206 we named index P. We apply and exemplify the capacities of index P to estimate transmission 207 potential and characterize it in space and time using real world examples. We also introduce and 208 provide MVSE, an R-package that implements various functionalities allowing to estimate and 209 analyse the index P in regions for which humidity and temperature time series are available. 210 We believe that both MVSE and the index P are unique in the current literature. Other useful 211 suitability indices exist, but these are based on complex approaches dependent on many variables 212 from various sources, with no direct biological interpretation on the measured scale, and for which 213 no freely available estimation tools exist. In contrast, index P is simply based on the R0 expression 214 of a mosquito-borne dynamic model and two climate variables, has a direct biological interpretation, 215 and is easily calculated through MVSE, a complete and freely available software tool.

216
With its aforementioned mathematical simplicity, dependency on only two climatic variables and 217 freely available R-package MVSE, we foresee the potential of the index P not only for direct 218 research purposes, but also for education (e.g. epidemiological courses on arboviruses) and field 219 ento-epidemiology (e.g. entomologists and epidemiologists collecting real-time data). New ideas and 220 methodological updates will be implemented in the MVSE package as the user base will keep 221 growing. Future expansions of the package could include functionalities regarding the total number 222 of female mosquitoes per human (M), or the inclusion of other local variables which may be related 223 to the risk of mosquito-borne viruses (e.g. socio-economic variables, rainfall, mosquito control 224 initiatives, etc.). Any future updates on the methodology and R-package will be made available in 225 the package's online repository (see S1 Text).   Table).