Parameter estimates of the 2016-2017 Zika outbreak in Costa Rica: An Approximate Bayesian Computation (ABC) Approach

In Costa Rica, the first known cases of Zika were reported in 2016. We looked at the 2016-2017 Zika outbreak and explored the transmission dynamics using weekly reported data. A nonlinear differential equation single-outbreak model with sexual transmission, as well as host availability for vector-feeding was used to estimate key parameters, fit the data and compute the {\it basic reproductive number}, $\mathcal{R}_0$, distribution. Furthermore, a sensitivity and elasticity analysis was computed based on the $\mathcal{R}_0$ parameters.


Introduction
Over the past few decades, arthropod-borne viruses (arboviruses), particularly those transmitted by mosquitoes, have emerged and/or reemerged as significant public health problems worldwide [1]. Different conditions, such as climate change, globalization, human movement, and the lack of adequate vector control strategies, are among the complex factors, that have influenced the growing burden of arboviral infections around the globe [2]. The present threat that these viruses represent, highlights the importance of a better understating of its transmission dynamics, in order to provide public health care authorities with important evidence that will help in the development of more efficient prevention and control approaches. This research article, focuses on Zika virus (ZIKV), one of the most recent mosquito-borne pathogens that have emerged as a global public health concern [3], specifically in the 2016-2017 outbreak that took place in Costa Rica.
Originally isolated in 1947 from the Zika forest of Uganda [4] and later identified in humans in the 1950's [5], ZIKV remained as a relative obscure arbovirus for nearly 60 years, with fewer than 20 human infections confirmed from countries in Asia and Africa [6,7]. It was until 2007, that the disease started to gain global attention, after approximately three quarters of the residents of Yap Island, were infected by the virus [8]. This initial outbreak, was followed by a series of epidemics, first in the French Polynesia and other Pacific Islands in 2013-2014, later reaching South America in 2015 [9]. From there, the virus spread rapidly to Central and North America. By January 2018 a total of 223,447 confirmed ZIKV infections, and 3,720 confirmed congenital syndrome associated with ZIKV had been reported to the Pan-American Health Organization [10]. Countries in the Central American Isthmus have reported a total of 7,820 confirmed cases and 190 confirmed congenital syndrome associated with ZIKV [10]. In Costa Rica, during that same period, a total of 10,253 suspected case were reported by the Ministry of Health [11], of which 2,144 were confirmed by laboratory, six newborns and one patient with Guillain-Barré Syndrome have been identified as positive with the ZIKV [12,13]. However, the actual proportion of human population infected and the number of new cases are extremely difficult to verify [2].
In urban environments, ZIKV is mainly transmitted in a human-mosquito-human transmission cycle, with the female Aedes aegypti as a primary vector and to a lesser extent Aedes albopicuts [14,15]. However, the recent outbreaks, had led to the identification of multiple confirmed and probable secondary modes of transmission [16]. It has been documented that the ZIKV, can be passed from a pregnant mother to her child during all trimesters and at time of delivery [17,18,19], accidentally as a result of a needle stick injury within a laboratory [16] and during sexual contact with an infected partner [20], where the virus has been found in the semen of infected men [21,22], even when it is no longer detectable in blood, usually clearing from the semen after about 3 months [23]. The virus, has also been isolated from breast milk [24,25] and in blood transfusions [26]. Once a person is infected, the average incubation period is estimated to last between 2 to 7 days [27]. It has been reported that approximately 80 percent of people infected with the virus are asymptomatic [8]. When symptomatic, the clinical manifestations, are usually mild and self limiting [28]. The main complications caused by Zika are neurological, where cases of Guillain-Barré syndrome had appeared in infected adults [6]. In neonates infected with Zika, an increased number of nervous system malformations, in particular microcephaly, were reported [29]. The diagnosis of the virus, is based on clinical examination, however it is very challenging, because of the overlap of symptoms with other arboviral diseases [30], the cross-reaction between antibodies directed to ZIKV and other flaviviruses [31], and the limited reliable techniques to diagnose infection accurately in those who have symptoms [2]. To date, an accurate Zika diagnosis is best achieve by real-time RT-PCR [31] and there are no specific treatments available for ZIKV infections, therefore only symptoms such as pain and fever can be treated [27].
To discuss potential control strategies of an infectious disease, such as ZIKV, mathematical modelling plays an important role [32]. In general, these models can help decision makers to focus on key aspects of the disease, and are now commonly use in the development of public health policies around the world [33,34]. An important threshold quantity that can be estimated with such models, is the basic reproductive number, R 0 [35]. This threshold can be estimated using a variety of mathematical and statistical techniques, in this article we used an Approximate Bayesian approach (ABC) in order to obtain a parameter distributions and R 0 .
Given the data provided by the Ministry of Health of Costa Rica, weekly reported ZIKV cases, and using the statistical methodology described in Section 3, we estimated key model parameters, as well as the R 0 distribution. This, to assess viable public health strategies regarding the key feature of the model, host availability for mosquito feeding. The article is organized as follows: in Section 2, we provide details on the data and statistical methodology that was used to estimate parameters, as well as the description of the model used. In Section 3, we show the parameter estimation using the Approximate Bayesian Computation. In Section 4, we provide the results and, in Section 5, we give our conclusions.

Study area
Costa Rica, is a country located in the Central America isthmus, bordered by Nicaragua to the North, Panama to the south, the Caribbean Sea to the east and the Pacific Ocean to the south-west. It has a territorial extension of 51,100 square kilometers, politically divided in seven provinces, 82 municipalities and 479 districts. The total population, is approximately 5,003,402 people [36], of which an estimated 75% live in urban areas [37]. The country has a natural environment rich in biodiversity, and tropical weather with two well defined seasons, a dry season (December-April) and a rainy season (May-November) [38]. The annual average temperature varies from 27 • C to 10 • C, with an annual rainfall that can go from 1500 liters per squares meter, up to 5000 liters per squares meter, depending on the region [39]. All tropical conditions that greatly influence Aedes mosquito population and behaviour.
The Ministry of Health, is the main entity responsible for the surveillance of arboviral diseases, and has a long standing experience on its management. In the 1980's the Ae. aegypti mosquito re-infested the country, after a brief period of eradication [40]. Since 1993, dengue virus has been endemic in Costa Rica, and three of the four serotypes have circulated the country [11]. In 2014, chikungunya virus also emerged, causing outbreaks in different regions of the national territory [11].
In Costa Rica, Zika was first documented in January 2016, as an imported case from a 25 year old Costa Rican male who contracted the virus after a trip to Colombia [41]. Prior to this patient, a 55 year old US tourist that had travelled to Nosara, Guanacaste from 19 to 26 December 2015, was diagnosed with Zika upon his return to Massachusetts [42]. Subsequently, in February 2016, the first two autochthonous cases were confirmed by the Ministry of Health, one from a 24 year old pregnant woman, and another from a 32 year old woman, both of them residents of Nicoya, Guanacaste [43]. Since then, the virus spread rapidly throughout the seven provinces, mainly in areas near the coasts, which are largely infested by the mosquito and have high circulation of the other two arboviruses, dengue and chikungunya. In 2016, the most affected province was Puntarenas, with an incidence rate of 791 cases per 100000, followed by Guanacaste (513 cases per 100000), both of them in the Pacific coast, and Limón (210 cases per 100000) in the Caribbean coast [11]. The municipalities that reported the highest incidence rates were: Garabito (3159 cases per 100000 population), Orotina (2900 cases per 100000) and Esparza (2450 cases per 100000), all of them located in the province of Puntarenas. In 2017, transmission continued, however, the number of cases were less than the previous year. Contrary to 2016, the most affected province was Limón (367 cases per 100000), followed by Puntarenas (56 cases per 100000) and Guanacaste (45 cases per 10000) [11]. The municipalities with the highest incidence rate were: Siquirres (690 cases per 100000 population), Guácimo (526 cases per 100000 population) and Pococí (395 cases per 100000 population). No notified deaths associated with Zika or concomitant infections by Zika, dengue or chikungunya, have been identified in any of the seven provinces [13].
(b) Zika incidence in Costa Rica, 2017. Figure 1: Zika incidence per municipality, Costa Rica, 2016-2017. In 2016, the municipality with higher incidence was Garabito, a municipality located in the Pacific coast of Costa Rica [11]. In 2017, the municipality with higher incidence was Guácimo, a municipality located in the Caribbean coast [11].
In Figure 1, we illustrate the concentration of ZIKV cases in 2016 and 2017 in Costa Rica. Note that the vast majority of cases were reported in coastal regions where temperature is ideal for mosquito prevalence.

Data
Weekly reported data of Zika virus cases from 2016-2017 was obtained from the Ministry of Health [11]. Between January and December 2016, a total of 7,820 cases were reported from hospitals and clinics around the country. The peak of patients was observed during the months of August and September (epidemiological week [31][32][33][34][35][36][37][38][39][40][41], as shown in Figure 2 [11]. That first year, the laboratory responsible for coordinating the virological surveillance of arbovirus at a national level, the INCIENSA National Virology Reference Center, analyzed 6,297 samples, of which 1,794 were positive [12].

Single outbreak model
There is an ample history of mathematical models for vector-borne diseases [44,45,46,47,48,49,50,51,52,53], among others. Moreover, other models specifically on estimating parameters for dengue, chikungunya and Zika include: [53,54,55,56,57,58]. We introduce a nonlinear differential equation single-outbreak epidemic model that describes the Zika dynamics with sexual transmission and host availability for mosquito feeding in Costa Rica. Essentially, individuals who stay home (inside) have a higher chance of being bitten by a mosquito [59]. In the case of Costa Rica, females conform the smaller proportion of the workforce, hence a larger proportion of females tend to stay at home (inside). In Figure 3, we illustrate the model transitions. The population is di- Sexual activity probability for males and females, respectively α m,f Per-capita exposed rate of males and females, respectively δ m,f Per-capita diagnosis rate after infection γ m,f Per-capita recovery rate of males and females, respectively β z Transmission rate (mosquito-human) ξ m,f Host availability probability of a mosquito feeding on a host µ Per-capita mortality rate of vectors κ Mosquito biting rate p vh Probability of infection from mosquito to human vided as follows: S m,f -susceptible males and females, respectively, E m,f -exposed males and females, respectively (infected but not infectious), U m,f -infected (undiagnosed and infectious) males and females, respectively, D m,f -infected (diagnosed) males and females, respectively and R m,f -recovered (immune) males and females, respectively. For the vector we have: S z (susceptible mosquitoes), E z (latent mosquitoes) and I z (infected mosquitoes). The following is the system of nonlinear differential equations: (1) And the transmission rate is given by λ(t) = Nz Nm+N f κp hv , where κ is the mosquito biting rate and p vh is the probability of infection from mosquito to human. Model parameters and description are shown in Table 1.
Primarily, we are interested in analyzing the effects of host availability, as well as sexual transmission. Here, the host availability was modeled via the parameter ξ m,f (male and females), where ξ m,f ∈ [0, 1]. This parameter serves as a reduction factor in the transmission dynamics from mosquito to host. In the case of sexual transmission, σ m,f (males and females), where σ m,f ∈ [0, 1], measures the sexual activity in the host population.
Furthermore, using the next generation approach by [60], we computed the basic reproductive number, R 0 . However, the analytic expression is highly complex, therefore, we computed its value in Mathematica using the parameter estimates.

Parameter estimation
We used weekly reported cases from [11] and the Approximate Bayesian Computation (ABC) to fit the model, estimate parameters and the basic reproductive number distribution from the 2016 − 2017 Zika outbreak in Costa Rica.
In Figure 3, we show the number of Zika confirmed cases and the model solution after estimating parameters.

Parameter estimation using ABC
The parameter estimation method that we use in this article is the Approximate Bayesian Computation (ABC). This method seeks to approximate the posterior distribution of the parameters through algorithms where the evaluation and optimization of the likelihood is not performed. On the other hand, many of these methods are based on sequential algorithms that allow the approximation of the posterior density using sampling schemes such as rejection sampling, MCMC or sequential Monte Carlo sampling (see [61] for an exhaustive summary of the subject).
Some of the parameters were kept fixed in order to follow the ecology constraints of the vector, as well as some conditions of the transmission dynamics of the disease that have been studied in previous articles ( [52], [58]). The total set of parameters (fixed and estimated), as well as the initial conditions of the system are presented in Table 2.
The initial conditions of the susceptible populations and the initial values of the three remaining parameters (µ, κ and δ) are obtained through the minimization of the sum of squared differences between the observed cases and the total number of cases per week according to the model (1). We used a genetic algorithm ( [62], [63]) and the quasi-Newton L-BFGS-B method ( [64]) to obtain the minimum. With the first algorithm we make an initial search in the parameter space and with the second algorithm the search for the first step is improved using as initial value the result of the genetic algorithm.
Once we obtained the initial values of the three parameters, we used a rejection sampling scheme for the ABC ( [65]) through the EasyABC package of R ( [66]). The prior densities for the three parameters were truncated normal to assure that the parameters are strictly positive, and their corresponding means are the initial values obtained above.

Results
A summary of the main estimation results is shown in Table 2. It contains the point estimates and the Bayesian 95% prediction intervals of the active parameters and the assumed values in case of fixed quantities.
In Figure 4 we show the data and model solutions based on the estimated parameters and initial conditions (see Table 2) and different model solutions based on ξ m and ξ f , as well as the respective R 0 values. In Figure 4a, the value of ξ m = 0.5 and ξ f = 0.8. Furthermore, in Figure 4b, we have ξ m = 0.5 and ξ f = 0.5, with a reduction in female host availability for mosquito feeding. This leads to a significant reduction in ZIKV cases and an R 0 = 1.434. In Figure 4c, a more extreme case where ξ m = ξ f = 0.4, host availability for males and females is reduced, we observe an even lower number of ZIKV cases and an R 0 = 1.368.
The posterior densities of the three active parameters and the R 0 are shown in Figure 5. We obtained the posterior densities using the ABC approach under a rejection scheme with 100000 samples, from which we selected the 1000 samples that guaranteed the smallest mean square errors between the number of observed and theoretical diagnosed cases. It is important to note that the level of precision of the parameter estimates, measured though the coefficient of variation of the posterior distribution, attains their minimum values for R 0 and µ, as can be seen in the last column of Table 2. Therefore, the degree of accuracy with which we are estimating R 0 is largely reliable. Table 2: Model parameters, initial conditions and R 0 . We performed an elasticity and sensitivity analysis on the R 0 parameters (see Table  3). Our analysis suggests that host availability, specifically female host availability, is the most sensitive parameter in the model. This suggests viable opportunities for prevention and control strategies. This phenomena is in part due to a large proportion of the female population is out of the workforce. The Costa Rican census [59] indicates that females old enough to work, that is age 15 or older, approximately 26% are not in the workforce. This scenario gives the mosquito more feeding opportunities throughout the day assuming they stay home (inside). In contrast, approximately 13% of the male population is out of the workforce [59], in turn, males tend to be more mobile giving the mosquito less opportunities to feed.

Fixed parameters
We observe that the most sensitive parameters are host availability (ξ m,f ) followed by the transmission rate from mosquito to human (β z ) and biting rate (κ). These parameters   increase R 0 , therefore creating opportunities to reduce incidence based on strategies in accordance with these parameters may be viable solutions to reduce ZIKV incidence in the future. However, it has been very difficult in the past to significantly reduce the adult mosquito population just by fumigating and other common strategies. In many cases, fumigation is solely done in residential surroundings and not inside homes where female mosquitoes tend to live, close to their food source. Reducing the biting rate could entail strategies surrounding repellents, and in some vulnerable areas this possibility incurs in an economic burden for the at-risk population. Hence, not an optimal or viable solution.

Conclusions
After the rapidly spread of ZIKV through different Latin American countries, the introduction of the virus to Costa Rica was expected [68]. Despite efforts made from public health authorities, the virus has circulated the seven provinces and has caused important consequences in the health and well-being of the Costa Rican population, in particular, due to neurological complications. According to data provided by the Costa Rican Social Security Fund (CCSS), from January 2016-December 2017 there were a total of 21 hospitalizations due to Zika, 13 of which were from women in a reproductive age. During that same period of time, 323 pregnant women and six newborns with microcephaly were confirmed to have the virus [69]. However, the actual number of newborns affected by ZIKV could be higher, as evidenced in a report published by The Center for Congenital Diseases in Costa Rica, where they established that since the introduction of ZIKV to the country, the cases of microcephaly almost doubled, and by 2017, the number exceeded almost four times, when compared to the period used as the base line by the Center (2011-2015), affecting all seven provinces in different proportions, being the most affected Guanacaste, Puntarenas and Limón, which could indicate that the observed increase may be associated, in some level, with ZIKV [70]. It is evident that ZIKV possesses an unique challenge to public health authorities around the world. The multiple routes of transmission, the large number of asymptomatic cases, as well as, its link with microcephaly and other neurological disorders, forever changed the public health approach to this virus [71,72]. In Costa Rica, a country greatly infested by the Aedes mosquito, a better understanding of the transmission dynamics, can provide a guide to introduce more efficient prevention and control interventions, that allows a better use of the human and economic resources available for the control of vector borne diseases.
Despite not having an analytic expression for the basic reproductive number, we estimated its distribution using the estimated parameters from Table 2 by means of a simulation-based Bayesian approach (ABC). This method has the advantage that it does not require strong assumptions on the data which are necessary if we use classical estimation methods, such as maximum likelihood or least squares. Furthermore, the Bayesian nature of ABC allows us to assume that some parameters of the model are random with the additional advantage that the modeling of those parameters includes the understanding of their variation.
Using data from the Ministry of Health's surveillance system, our study showed that the risk of ZIKV transmission in Costa Rica is greater in the population that spends most of its time inside. This results go in hand with previous studies that evidenced that females adult Ae. aegypti disperse relatively short distances [73], and highlights the need for more effective community engagement strategies, in order to enable residents of different areas in the communities to make informed health decisions that will influence their overall well-being [74]. This becomes more significant, based on the fact that women constitute the majority of the population that stays inside, which gives an extra layer of relevance to the findings. The study also evidenced that for Costa Rica, the sexual transmission route for the virus plays a secondary role in the dispersal of the disease, however, public health officials need to remain vigilant and provide information to the general population about the risk of sexual transmission of ZIKV.
Although participatory approaches from the communities have been used for many years in the control of mosquito borne diseases, its proper implementation has been difficult to achieve [75]. A bottom-up communicative approach at all levels of the public health system [76] could provide a greater interest in general population to implement the mosquito control strategies widely recommended by the public officials around the country. A better understanding of the specific needs of each region, each one with different social, cultural and economic characteristics, makes it fundamental to plan mosquito control strategies according to its specific needs [77]. Not taking into account the great heterogeneity of houses and neighborhoods where the mosquito completes its life cycle, along with the limited resources and personnel trained in the control of the vectors, are at least in part, attributed to the failure of previous prevention and control actions [78]. Because of the ongoing transmission, and the risk of recurrence of an outbreak, health care authorities around the country need to remain vigilant and establish a comprehensive arboviral disease surveillance. An education and prevention control interventions adjusted to each specific region, with a more active involvement of the communities members is recommended to achieve a more efficient control of the ZIKV and other mosquito-borne diseases.