Modeling Chikungunya control strategies and Mayaro potential outbreak in the city of Rio de Janeiro

Mosquito-borne diseases have become a significant health issue in many regions around the world. For tropical countries, diseases such as Dengue, Zika, and Chikungunya, became epidemic in the last decades. Health surveillance reports during this period were crucial in providing scientific-based information to guide decision making and resources allocation to control outbreaks. In this work, we perform data analysis of the last Chikungunya epidemics in the city of Rio de Janeiro by applying a compartmental mathematical model. Sensitivity analyses were performed in order to describe the contribution of each parameter to the outbreak incidence. We estimate the “basic reproduction number” for those outbreaks and predict the potential epidemic outbreak of the Mayaro virus. We also simulated several scenarios with different public interventions to decrease the number of infected people. Such scenarios should provide insights about possible strategies to control future outbreaks.


Introduction
In the last decades, Mosquito-borne diseases have become a significant health issue in many regions around the world. Projections indicate that around 2050, half of the population will be at risk of some arbovirus infection [1]. These arboviruses, which include diseases such as Dengue, Zika, and Chikungunya, are epidemic in most of the tropical countries. Besides temperature and humidity, human migrations and sanitation also contribute to the epidemic conditions in these places [2,3]. For example, around 300.000 people were infected by Dengue, Zika, or Chikungunya by the end of the 11 th week of 2019 in Brazil. This number represents almost three times the reported cases in 2018 for the same period [4]. These surveillance reports over time are essential in providing scientific-based information to guide decision making, resources allocation, and interventions [5]. The usage of mathematical models has demonstrated to be a powerful tool in contributing to these data analysis [6][7][8]. One of the most significant parameters a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 extracted from these analyses is the basic reproduction number R o . R o is defined as the number of secondary infections derived from one single infectious subject and is widely used as an epidemiological metric employed to describe the transmissibility of infectious agents [9].
Here we apply a compartmental mathematical model to investigate the dynamics of Chikungunya outbreaks in the city of Rio de Janeiro in Brazil. The model consists of ordinary differential equations that describe the transmission and the transition of the diseases in humans and vectors [10,11]. The model's parameters were extracted from the literature or obtained from the best fit from the data of Rio de Janeiro surveillance report for the years of 2016, 2018 and 2019 [12]. Based on these parameters, we estimate the basic reproduction number R o for Chikungunya outbreaks in those years. We also simulate a scenario predicting if the Mayaro virus could be a potential epidemic disease in Rio de Janeiro. Modifications in the standard model equations were implemented to introduce different possible interventions in order to decrease the number of infected people [13]. Those simulated interventions include actions such as killing adult mosquitoes by fogging, decreasing mosquitoes birth rate by removing places where the vector lays eggs, e.g., removing standing water and, decreasing the contact between an infected human with mosquitoes by stimulating repellent usage. A scenario containing all interventions was also performed for different intensities of those actions.

Materials and methods
In this work, we perform mathematical modeling of Chikungunya outbreaks in Rio de Janeiro for the years 2016, 2017 and 2019 [12]. The Chikungunya virus infects humans through mosquitoes as the disease vector. The model adopted here is a compartmental model known as SEIR (Susceptible, Exposed, Infected, and Recovered) [8,10,14]. The approaches using this class of models have been successful in modeling epidemic related to human vector dynamics [11,15]. Fig 1 presents a schematic description of this modeling.
The human disease flow is presented by the blue blocks where S is the susceptible proportion of humans, which become exposed to the virus, E, at a rate β h after the contact with infectious mosquitoes Z. After the latent period λ h , the exposed humans become infectious: either symptomatically I or asymptomatically I a ; the parameter ϕ determines the ratio between the infectious states. Finally, the infected humans recover reaching the state R at rate α.
In the case of the vectors disease flow, shown by the red blocks in Fig 1, the susceptible mosquitoes X become exposed Y at a rate β m after acquiring the virus from infectious humans. λ m defines the latent period for the exposed mosquito to transition to the infectious state Z. In this modeling, we assume that human mortality and birth rates are the same, keeping the human population constant. For the vectors, we set the parameter μ and μ o as the mortality and birth rate, respectively. The model is represented by the following set of differential equations: Table 1 shows the definition of each state in the model for both humans and mosquitoes. Those states will dynamically vary during the model simulation in which the parameter I related to the number of cases reported by the surveillance data will be the variable used in the fitting of the model simulation. Table 2 describes the parameters and ranges used in the model. Some data information comes from the literature, and for the parameter which we have no description, they will be obtained from the model best fit.
In this work we estimate the basic reproduction number R o by applying the next generation matrix method [21,22]. R o indicates the number of secondary infections derived from one single infectious subject and can be described by: Where ρ(K) is the spectral radius of the matrix K = PΓ −1 . P is the transmission matrix that contains the rates of humans to get infected by the vector and vice-versa: Γ is the transition matrix that takes into account the transitions from being exposed to become infectious: The mathematical solution of (1) gives an expression [11]: The Eq 4 has parameters in which there is no information available such as β h and β m . In order to estimate R o , these parameters will be obtained from the best fit of the model simulations using data from the surveillance reports [12].

Results and discussion
The usage of the SEIR model to investigate diseases epidemics provides a tool to quantify different parameters in outbreaks. The basic reproduction number R o is the most important quantity, and it is defined as the number of secondary infections caused by an infected individual [3,23]. It estimates the potential of an outbreak to occur in the case of R o > 1 [10,24]. The knowledge of R o also gives insights into the understanding of the epidemiology of a particular disease and its spreading changes over time and geography [25].
The number of secondary infections in humans from an infected human, defined as R T , the type reproduction number [26], can be obtained by R o squared (R T = R o 2 ) [11]. The information of R T can be used to estimate the number of people that need to be isolated or vaccinated (Q) to contain the epidemic using the relation Q = 1 − 1/R T . We present the application of the SEIR model in the data of Chikungunya in Rio de Janeiro-Brazil in different years (2016, 2018 and 2019). We also provide an estimation of the potential outbreak of Mayaro virus in Rio de Janeiro. Some additions on the model are also proposed in a way to simulated possible interventions in the epidemic control [13,27].

Basic reproduction number-R o
The data which contains the weekly number of infected people of Chikungunya outbreak in Rio de Janeiro was obtained from the surveillance report for the years 2016, 2018, and 2019, which are publicly available [12]. The ratio between the number of reported cases and the total population is presented in red circles Fig 2. The total number of cases reported in 2019 is 22896 until the 26 th week when the data was collected. This number is almost two times higher in comparison with 2016 and 2018 in 52 weeks, 14203 and 10700, respectively. In 2017 the total number of cases was 1870, which will not be used in this study. Then, the SEIR model was applied to fit the incidence data where the upper and the lower bound of the model parameters were set to vary in a range described in the literature. The transmissions coefficients λ h , λ m , α, β h and β m values are obtained from the model best-fit since there are no values reports in literature about these parameters for Rio de Janeiro [11]. The simulated number of cases from the best-fit are presented as bars in Fig 2A, 2B Table 3. Table 3 also shows in the last row, the estimation of R o using Eq (4)  Considering other epidemic diseases in Rio de Janeiro as Dengue and Zika, which are transmitted by the same vector, the parameters estimated here are similar to other registered outbreaks studies [25,28]. It is worthwhile to mention that these parameters give insights from a city as a whole, which is invariant on how heterogeneous the sanitary conditions could be at different neighborhoods [7,25,29,30].
The last column in Table 3 presents the estimated parameters for Mayaro virus using the data for the Chikungunya outbreak in 2018, which is the most recent complete data available.
The assumption is based on the similarities between these two alphavirus in which they can also be transmitted by the same mosquito vector: Aedes aegypti [6,31,[34][35][36]. The estimated R MAYV o for Mayaro presents values between 1.18 and 3.51 for the lower and upper limits. Even the lower bound R MAYV o is greater than 1, suggesting that Mayaro has the potential to be an epidemic disease as recent reports are signaling for different locations [37][38][39].

Sensitivity analysis
The usage of the SEIR model allows us to estimate the importance of each parameter in the characterization of an outbreak by performing a parameter sensitivity analysis (Fig 3). It was carried out a hundred thousand Monte Carlo simulations sampling around the ±5% range from the best fit value of the Chikungunya epidemic data of 2018. From these simulations, the time-shift on the peak and the total number of infected people were obtained. Additionally,  the Spearman correlation for each parameter with the time-shift and the total number of infected people was calculated. [11]. In Fig 3 the transmission parameters (β h and β m ) in both: peak shift and total population infected, exhibit as an important contribution to the outbreak characterization [40]. Fig 3A shows that the rate of recovery plays a significant role in the time scale of the epidemic. The parameter presented a negative correlation, which means that positive variations in this parameter will modify the shift by delaying the peak reaching point; the same is observed about the mosquito birth-death rate. On the other hand, Fig 3B shows that by decreasing the magnitude of the transmission coefficients or increasing the value of γ and μ may reduce the number to total cases of infected people. These parameters modification can be related by triggering public interventions, which will be discussed in the next section. Intuitively, a fast recovery would decrease the probability of a mosquito to bit an infected person, while the increase on the birth-death rate of the mosquito will create similar effect: an infected mosquito would die quickly, therefore it may not be able to transmit the disease. Sensitivity analyses based on the outbreak time dependence are presented in the Supporting Information.

Interventions
In this section, we will discuss the outcome of different possible intervention strategies to control the epidemic disease spreading [13,23]. The simulations were carried out using the Chikungunya epidemic outbreak data from Rio de Janeiro in 2018. The first approach simulates the action of killing adult mosquitoes, which is related to the use of insecticide as fogging. In the model, this strategy appears as an increase in the mosquitoes mortality rate μ presented in Eq 5: where θ(. . .) is the unitary step function, μ c is the natural rate of birth/death of the mosquito, C is the cumulative number of infected people and, ω is the parameter related to the intensity of the fogging action reflected in the mosquito death rate μ. The fogging action is triggered when the cumulative number of infected people C, described in Eq 6, reaches the value C p which is 30% of the total number of cases from the real data. For all the interventions discussed in this study, the trigger event will be the same as the one present here in the fogging action.
Fig 4 presents the distribution of the number of cases and the cumulative number of cases as a function time for different fogging intensities, A and B, respectively. In Fig 4A, once the fogging action starts, the number of cases per week stop to grow and starts to decrease over time. The strength of the parameter ω dictates how fast these curves decay. The cumulative number of cases also reflects the fogging action for different intensities, as presented in Fig 4B. The total number of case drops to 70.0% when ω = 0.25 presented in dashed red line and drops to 54.0% when ω = 0.5 shown in the dotted yellow line when compared with the real data without fogging action ω = 0.0.
The second simulated intervention is the action of reducing the birth rate of the vector. This approach can be associated with population orientation or better social sanitary conditions. These actions may produce a decrease in the number of places where the mosquitoes lay the eggs such as standing water, for example. Here, this intervention appears in the model by decreasing the mosquito birth rate μ o , as shown in Eq 7: where θ(. . .) is the unitary step function, μ c is the natural rate of birth/death of the mosquito, C is the cumulative number of infected people, C p represents the cumulative amount of infected people needed to trigger the action and, ω o is the parameter related to how efficient are the population and government actions in preventing the vector from laying the eggs, which led to decrease the mosquito birth rate μ o . Fig 5 shows the distribution of the number of cases and the cumulative number of cases as a function time for different intensities of the mosquito birth rate reduction. Similar behavior as the fogging intervention is observed here. In Fig 5A,   intervention ω o = 0.0 presented in Fig 5B. Although the behavior is similar to the fogging action, the response of decreasing the mosquito birth rate to the total number of cases is less efficient during the outbreak. It is worthwhile to mention that this kind of action, different from the fogging, can be inherited and passed to the following years and avoid new outbreaks to occurs in the future.
The third and the last studied intervention acts as the reduction of the rate in which infected humans transmit the disease to the mosquitoes. This effect can be associated as a quarantine action, isolating infected people or, more realistic, the usage of repellents by the infected human. Both strategies go in the direction of decreasing the contact between infected humans and the vector which is simulated using the Eq 8.
where θ(. . .) is the unitary step function, β c is the natural rate at which humans infect mosquitoes, C is the cumulative number of infected people, C p represents the cumulative amount of infected people needed to trigger the action and, � is the parameter that modulates how intense is the decrease in the rate at which humans infect mosquitoes β m . Fig 6 presents the results for this last intervention. The curves in Fig 6A and 6B show a similar pattern, as observed in the other two previous actions. The number of cases shows a decay after the intervention starts for different intensities of �. In Fig 6B the total number of cases curves present results closer to the birth control intervention than the fogging action. For � = 0.25, the total number of cases decreased to 82.7% of the initial value, meanwhile, for � = 0.5 this number drops to 64.1%.
A combined simulation applying all the three interventions was carried out, and the results are presented in Fig 7. The combined intervention presents, as expected, the most effective strategy to decrease the number of infected people. The distribution curve of the number of cases per week shows more intense decay in Fig 7A. The total number of cases drops to 54.3% from the initial value when ω, ω o , � = 0.25 and reduce to 43.6% for the parameters ω, ω o , � = 0.5 in Fig 7B.

Conclusion
The last three Chikungunya outbreaks in the city of Rio de Janeiro, Brazil, were modeled using the SEIR model and estimates the Basic Reproduction Number R o for the years 2016, 2018, and 2019. The simulations results register values greater than 1 for all of them, and 2019 is the most severe, even though the data was limited for the first six months. The calculation of R o gives a global overview of the impact and scale of the outbreak. Sensitivity analyses were performed to indicate, quantitatively, the importance of each parameter to the epidemic profile in different stages of the outbreak. A more detailed approach could take into account the number of infected people in each neighborhood with different sanitary conditions, and such details are not explored in this work. This study was expanded to include the Mayaro virus, which was reported as an emerging disease in South America [37,39,41]. Based on the assumption  that Mayaro and Chikungunya viruses have a similar spreading mechanism [37], since both viruses have the same vector [31,36,42,43], we used parameters fitted from the Chikungunya outbreak from 2018 to estimate the R MAYV o from Mayaro. The results indicate that Mayaro has the potential to be an epidemic disease in Rio de Janeiro with R MAYV o values in a range of 1.18 and 3.51. Also, to possibly stop or at least decrease the intensity of an outbreak, three interventions strategies were proposed by modifying the basic equations of the SEIR model. These interventions are associated to the increase of the vector mortality rate by fogging techniques, the decrease of mosquito birth rate by decreasing the amount of places where the mosquito lay the eggs and, the decrease of the rate in which humans transmit the disease to mosquitoes by the isolation of infected people or the usage of repellent. Although those simulations do not retract real data, they can contribute to discussions about public and government policies directions.