Parameter fitting using time-scale analysis for vector-borne diseases with spatial dynamics

Vector-borne diseases are becoming increasingly widespread in a growing number of countries and it has the potential to invade new areas and habitats, either associated to changes in vectors habitats, human circulation or climate changes. From the dynamical point of view, the spatial-temporal interaction of models that try to adjust to such events are rich and challenging. The first challenges are to address the dynamics of the vectors (very fast and local) and the dynamics of humans (very heterogeneous and non-local). The objective of the present paper is to use the well-known Ross-Macdonald models, incorporating spatial movements, identifying different times scales and estimate in a suitable way the parameters. We will concentrate in a practical example, a simplified space model, and apply to Dengue’s spread in the state of Rio de Janeiro, Brazil.

Vector-borne diseases have now received increased attention because of their high 2 potential of dissemination [1,2]. Mosquitoes of the species Aedes Albopictus and Aedes 3 Aegypti are the most responsible for virus transmission, such as Dengue, Zika, 4 Chikungunya and Yellow Fever. Due to lack of vaccination, basic sanitation, climate 5 changes, and with increasing human mobility, such diseases are spreading and appearing 6 in new regions, where the climate favors the proliferation of vectors. For instance, we 7 may cite locations in Portugal, France and Italy where cases of Dengue and 8 Chikungunya have already been recorded, and the United States with cases of dengue 9 fever and Zika virus [3][4][5][6]. In addition, hosts may be infected in environments that are 10 different from their place of residence because mosquitoes do not travel long distances, 11 and this may lead to increased population heterogeneity and consequently in changes of 12 the disease dynamics.
bacteria have also been successfully tested, which prevents the vector from transmitting 22 the virus [11,12]. Some vaccines have been tested and others are in testing phase, but 23 the greatest difficulty is that such vaccines should be tetravalent, that is, they must be 24 effective against the 4 existing serotypes of the disease [13,14]. 25 Mathematical models used to describe indirectly transmitted infectious diseases have 26 the interesting characteristic of coupling the dynamics of hosts and vectors, whose 27 parameters have different time scales, the life cycle of mosquitoes is in days while 28 humans life cycle is in years. Dengue is an example of infectious disease in which this 29 occurs, in addition of having the particularity of different serotypes. Studies with 30 spatial networks, or meta-populations, provide a way to understand the interactions 31 between individuals in different scales, being a powerful tool to understand the 32 characteristics of transmission in communities, regions and countries incorporating 33 spatial heterogeneity [6,[15][16][17][18][19][20][21]. In addition, as in modelling the dynamics of several 34 vector-borne diseases, if the goal is to fit the model to real data, one has to deal with 35 the asymptomatic cases, reliable data, in particular for the mosquitoes population, 36 besides having to take into account the different time scales of vectors and hosts, which 37 makes it difficult to study and understand disease dynamics [22][23][24].

38
In this work, we consider a basic compartmental model that divides the human 39 population into susceptible S, infected I and recovered R (sometimes making it a little 40 simpler), coupled with susceptible S m and infected I m mosquitoes. This model depends 41 on parameters such as the mosquito mortality rate, the transmission rate and the total 42 vectors population, which are very difficult to measure. From this model we do the time 43 scale separation of humans and vectors. As a consequence, we have that the mosquitoes 44 equations do not appear explicitly in the model, and the new parameters depend only 45 on the total mosquitoes population. Finally, we incorporate spatial movements, 46 considering mobility between cities two to two. We will adjust this model to dengue 47 incidence data, whose cities were chosen based on [4]. Thereby, our main purpose is to 48 show that it is possible to have a good fit of this reduced model in the beginning of the 49 infection, as well as to estimate parameters related to mobility between two cities.

50
First of all, our approach will be deterministic, we will not take into account 51 stochastic effects or incorporate the element of chance in the models. Secondly, being 52 more precise, the goal is to consider the effects of the spatial dynamics into the 53 Ross-Mcdonald models and use it to fit to the real data. This can be done either by a 54 continuous space domain, which in turn will give us Partial Differential Equations, local 55 or non-local, or consider discrete networks in space, which will provide system of 56 ordinary differential equations ODE. There are advantages and disadvantages to both even if one proves that the second approach can be viewed as an approximation of the 60 first and that the dynamics must somehow converge. The second approach can more 61 easily be used to fit the real data, since they are always discrete in nature. Since, in this 62 work, we are interested in concrete data and fit the dynamics, we will concentrate in the 63 second model. For the continuous model, one can refer to [25][26][27].

64
The paper is organized as follows. First we set up the local dynamics, that will be 65 used to describe the dynamic in each city, we also identify the small parameter that will 66 be used. Next we set up the network dynamic, introducing a diffusion operator. With 67 this two ingredients, for completeness, we show a formal expansion that reflect the 68 general ODE singular perturbation results. Finally we can estimate our parameters 69 using a network found to represent the initial spread of the disease in the State of Rio  interaction dynamics between the compartments is described through a system of 78 ordinary differential equations (ODEs): In this model, a susceptible human becomes infected at a transmission rate β when in parameter µ h represents the birth/mortality rate of humans, γ is the recovery rate of 83 humans and µ m is the mosquitoes birth/mortality rate.

84
Assuming that birth and mortality rates are equal, we have that populations remain 85 constant over time, that is, and then work with the equivalent reduced system: Considering that the life expectancy of an adult female mosquito is about 10 which is equivalent to Making ε = 0 in the third equation, we have that I m (t) can be obtained as a 97 function of I(t) at any time t: Replacing this value of I m Eq (5) into the equations of S and I, it results in and then defining a = β ΩN m , b = Ω and c = µ m N h , we have a new equivalent system 100 without the variables related to mosquitoes: Mobility 102 The SIRS m I m model characterizes the dynamics of a disease within a population.

103
Thus, if we want to describe its transmission dynamics more realistically, we must 104 consider a mobility network that includes interaction between populations. Such a 105 network can be represented by a graph, where each node corresponds to a population 106 and an arrow that leaves the population r to s means existence of mobility of 107 individuals from population r to s [28].

108
Let N hr be the total human population that is registered in node r = 1, 2, ..., M , so 109 that the disease dynamics in each location is described by the S r I r R r S mr I mr model.

110
The parameter d rs ∈ [0, 1] corresponds to the mobility rate from the population r to s 111 per unit of time [18,28]. Accordingly, to include mobility in the Model (1), we must 112 consider the inflow and outflow of humans in each compartment, and since mosquitoes 113 do not move large distances, their respective compartments remain unchanged. Thus, 114 the system of equations representing mobility between cities r and s, r = s, is given by 115 November 18, 2019 4/12 with initial conditions S r (0), I r (0), R r (0), S mr (0), I mr (0). Here, we suppose that the 116 parameters are different for each location except the human birth/mortality rate µ h , 117 and the human recovery rate γ. The total human population for each patch is given by 118 N hr = S r + I r + R r , so Similarly, we include mobility in (7), resulting in the ODE system, named S r I r : Asymptotic expansion 121 Here we use a power series expansion to analyze the asymptotic behavior at ε = 0 of the 122 perturbed System (3) with mobility given by (8). In order to do that, let S, I and I m 123 be vectorial functions whose coordinates are denoted respectively by S r , I r and I mr for 124 r = 1, 2, ..., M . We consider the following singular perturbed system of ODEs: Now let us expand the solutions with respect to the small parameter ε, that is, let us 126 assume that vectorial functions S, I and I m given by (10) satisfy 127 S = S 0 + εS 1 + ε 2 S 2 + . . . I = I 0 + εI 1 + ε 2 I 2 + . . . and 128 I m = I m0 + εI m1 + ε 2 I m2 . . .

Thus, the time derivatives themselves set
November 18, 2019 5/12 which gives us from the right-hand side of (10) that Then, if we plug these expressions in the System (10), we obtain at ε = 0 that 132 Hence, we get as in (5) 133 I mr 0 = (Ω r /N hr )I r0 N mr (Ω r /N hr )I r0 + µ m and then, we deduce the reduced equations with initial condition S r0 (0) and I r0 (0). Thus, if we proceed as in (7) defining 135 a r = β r Ω r N mr , b r = Ω r and c r = µ m N hr , we can obtain the limit system (9) without 136 the mosquitoes equations. Hence, the solutions S and I of the System (10) can be 137 approximated to the solutions given by the limit equation (9). Indeed, under 138 appropriated assumptions, it has been shown in [29][Theorem 4.4] (see also [22,23]) that 139 the convergence is uniform in finite time with order O(ε) justifying our approach.
Parameter estimation 141 We have the reduced System (9) which does not depend directly on the mosquitoes 142 parameters. For the parameter's estimation, we employed the pomp package 143 implemented in language R [30], in order to obtain the values for: a r , b r , c r , d rs , γ and 144 the initial conditions of each location S r (0) and I r (0). The algorithm was applied to fit 145 the S r I r model to dengue data of Brazilian cities which presents evidence of mobility 146 according to the results obtained in [4].

147
The adjustment is done using the error of the least squares method. We set up a 148 function that will calculate the sum of the squared errors which consists of the 149 differences between the result obtained with the model and the data [30,32]. In our 150 work we do the fitting of two time series each time. The model with time scale 151 separation and considering mobility among two cities is given by: 12) and the sum of the squared errors is: where E 1 = I 1 − C 1 , E 2 = I 2 − C 2 , and C 1 and C 2 are the weekly incidences of City 1 155 and City 2, respectively.

156
After set up the error function which is the objective function, we apply an 157 optimization algorithm in order to minimize its value. To each parameter and variable 158 can be given a lower and/or upper bound, then it is necessary to start with an initial 159 value which satisfies the constraints and from this point, the optimization algorithm will 160 search the parameter space for the value that minimizes the objective function. The 161 algorithm uses function values and gradients to build up a picture of the surface to be 162 optimized [30]. It is important to note that it is a ill-posed problem and a change in the 163 initial conditions may change the estimated value of the parameters.

164
Data 165 We will consider dengue data for our simulations, which were obtained from the 166 brazilian Notification Disease Information System (SINAN) [31].  [28], which used the ideas of [33,34] to estimate an effective network that explained 172 the epidemic in Rio de Janeiro. In this year, they detected that the disease started and 173 spread to such cities, before spreading to the whole state. The fixed parameters are the 174 total population of each city and the birth/mortality rate of humans (see Table 1), and 175 the other parameters need to be estimated.

177
Regarding the initial values of the parameters, we will consider that the initial 178 susceptible population is 90% of the total population of each city due to the fact that 179 dengue is endemic in Brazil, so S r (0) = 0.90 × N hr . The human recovery rate is initially 180 γ = 1 (week). Also the parameters of mobility start with the values d rs = 0.0001, being 181 able to assume values in the interval [0, 1], where the indices represent the proportion of 182 the population that leaves City r and goes to City s.

183
Recalling that by definition: a r = β r Ω r N mr , b r = Ω r and c r = µ m N hr . As a 184 consequence of considering the System (11) with time scale separation, we have that the 185 mosquitoes equations do not appear in the system, however the values a r , b r and c r , 186 r = 1, 2 depend on the mosquito's parameter, being necessary to make an initial 187 assumption for its values. Lets consider that initially, µ m = 7/10 (weeks), so we use an 188 initial approximation for β r as β r = 2γ, Ω r as Ω r = 1µ m and N mr = mN hr . Actually, 189 by defining the parameter ε = µ h /µ m , it follows that Thus, we have that the total mosquito population N mr , is the only parameter of the 191 mosquitoes being used, which interferes with the initial value of the parameter a r .

192
Therefore, we will analyze the fit of the model to the data considering: N mr = N hr , 193 N mr = 2N hr , N mr = 3N hr and N mr = 4N hr , and also by varying the initial amount of 194 infected individuals I r = 1, 2, 3, 4. We will present the results with the value of N mr 195 that best fits the data, that is, the one that results in smaller quadratic error.

196
Simulations with the mosquitoes population larger than 4 times the human population 197 were not successful.

198
The results of the estimated parameters are presented in Table 2   1.750001e+03 0.000e+00 2.179996e+02 6.746 5.1194e-02 716631 3.8e+00 In the first column is shown the pairs of cities City 1 and City 2, and in the following columns, the results of the parameters obtained for each city. The first parameter related to mobility of each pair of cities is d 12 and the second is d 21 , corresponding to the proportion of the population leaving City 1 and going to City 2 and vice versa.  Table 2.  Table 2.  Table 2.  Table 2.   Table 2.

204
From a SIRS m I m model with humans and vectors, we made a separation of time scales 205 and reduced the system to an equivalent system independent of the mosquitoes 206 equations, and then we add mobility in this new system. So, the new parameters depend 207 only on the respective parameters for humans and the total mosquitoes population. 208 We made some considerations about the initial value of the parameters, and then we 209 applied an algorithm to fit the model to dengue data of some cities in Rio de Janeiro 210 state. The results showed that the reduced model was able to successfully adjust the 211 beginning of the dengue outbreak for all pairs of cities, also obtaining values for the 212 parameters related to mobility. This gives us an indication that human mobility 213 actually has influence on the spread of dengue.

214
In addition, we have that the model with time-scale separation S r I r depends only on 215 one parameter related to mosquitoes, whereas the SIRS m I m model contains three 216 parameters of the vectors, which makes it more difficult to fit the model to the data due 217 to the lack of available information about mosquitoes, since the time series of dengue 218 provide only the amount of humans infected weekly. This gives us perspectives for a 219 future study involving more complex mobility networks applied to analyse vector-borne 220 diseases.