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

1.   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 understanding of its transmission dynamics, in order to provide public health care authorities with evidence that could potentially help in the development of more effective prevention and control approaches. This article focuses on Zika virus (ZIKV), one of the most recent mosquito-borne pathogens that have emerged as a global public health concern, specifically in the 2016–2017 outbreak that took place in Costa Rica.

Originally isolated in 1947 from the Zika forest of Uganda ^[3] and later identified in humans in the 1950's ^[4], ZIKV remained as a relatively unknown arbovirus for nearly 60 years ^[5]. It was until recent outbreaks in most American countries ^[6,7] that the disease gained global attention, especially after its association with microcephaly and Guillain-Barré Syndrome (GBS) cases ^[8,9]. This flavivirus, closely related with the dengue virus, 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 albopictus ^[10]. It has also been documented that the ZIKV can be passed from a pregnant mother to her child during all trimesters and at time of delivery ^[11,12,13], through blood transfusion ^[14], laboratory accidents ^[15] and during sexual contact with an infected partner ^[16]. While several studies have evidenced that vector transmission remains the key driver of ZIKV transmission ^[17,18], sexual transmission has the potential to increase the risk of infection, epidemic size and prolong the outbreak ^[19,20]. Adding further complexity to the dynamics of the disease, it has been estimated that up to 80% of all ZIKV infections are asymptomatic ^[21], where even at low viremia levels, these cases have the potential to cause congenital Zika syndrome ^[22], which emphasizes the importance of silent ZIKV transmission.

To discuss potential control strategies of an infectious disease, such as ZIKV, the use of mathematical models can play an important role in the development of efficient control strategies ^[23]. In general, these models help decision makers focus on key aspects of the disease, and are now commonly used in the development of public health policies around the world ^[24,25]. An important threshold quantity that can be estimated with such models, is the basic reproductive number, $\mathcal{R}_0$, first presented in ^[26]. Its biological meaning describes the number of secondary infections caused by a typical infectious individual in a mostly susceptible population. Previous studies have estimated $\mathcal{R}_0$ to determine key parameters and region specific public health strategies ^{[27,28,29,30,31,32,33]}. This threshold quantity can be estimated using a variety of mathematical and statistical techniques, in this article we use an Approximate Bayesian Computation (ABC) approach in order to obtain parameter estimates and $\mathcal{R}_0$ distribution.

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. This data can provide public health authorities with information to assess viable public health strategies based on one of the key features 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 applied 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 (ABC). In Section 4, we provide the results and, in Section 5, we give our conclusions.

2.   Materials and methods

2.1.   Study area

Costa Rica is a country located in the Central America isthmus with a population of approximately 5,003,402 ^[34], of which an estimated 75% live in urban areas ^[35]. 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) ^[36]. The annual average temperature varies from 10 $^{\circ}$C to 27 $^{\circ}$C, with an annual rainfall that can go from 1000 mm per year, up to 4000 mm per year, depending on the region ^[37]. These 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 1980s the Ae. aegypti mosquito re-infested the country, after a brief period of eradication ^[38]. Since 1993, dengue virus has been endemic in Costa Rica, and three of the four serotypes have circulated the country ^[39]. In 2014, the chikungunya virus also emerged, causing outbreaks in different regions of the national territory ^[39].

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 ^[40]. Prior to this patient, a 55-year-old US tourist that had traveled to Nosara, Guanacaste from December 19-26, 2015, was diagnosed with Zika upon his return to Massachusetts ^[41]. 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 ^[42].

Since then, the virus spread rapidly throughout the seven Costa Rican provinces, mainly in areas near the coast, which are largely infested by the mosquito and have a high circulation of the other two arboviruses, dengue and chikungunya. No notified deaths associated with Zika or concomitant infections by Zika, dengue or chikungunya, have been identified in any of the provinces ^[43].

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.

2.2.   Data

In our model, we use the number of weekly clinically and laboratory-diagnosed cases of ZIKV infection reported during 2016-2017, to the Costa Rican Ministry of Health ^[39]. In Costa Rica, ZIKV infections are mainly identified during a clinical examination. Following national guidelines, laboratory testing is made only to those samples that meet the clinical and epidemiological criteria. Once endemic circulation in a region is determined, only samples from groups at risk, such as pregnant women, and newborns with microcephaly, are analyzed ^[44].

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–41) ^[39]. Through the first year of ZIKV introduction, 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 ^[45].

In 2017 ZIKV transmission prevailed, however, the number of reported cases were less than the previous year. Some factors that might have influenced this are: herd immunity and greater awareness by public health officials. A total of 2,414 cases were reported by the Ministry of Health ^[39], with a confirmation of 350 cases by the INCIENSA National Virology Reference Center ^[43].

2.3.   Single outbreak model

In general, mathematical models have been extensively used to understand the dynamics of vector-borne diseases ^{[46,47,48,49,50,51,52,53,54,55]}. Most recently, several papers have developed models specifically on estimating parameters for Zika ^[19,55]. For example, it has been estimated that $\mathcal{R}_0$ is most sensitive to the mosquito biting rate, mortality rate of mosquitoes ^[19] and to the human infectious period ^[56], where increased concentration of ZIKV in blood significantly increases the rate at which mosquitoes become infectious ^[57]. It has also been estimated that prior dengue virus infections could potentially protect individuals from symptomatic Zika ^[58] and that the effects of sexual transmission are not significant enough to sustain the disease in the absence of mosquitoes ^[33], but can play a role in sustaining the epidemic, in the epidemic size and in creating potential endemic scenarios ^{[19,20,59,60]}.

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. We included the element of host availability because adult Ae. aegypti mosquito prefers to rest indoors ^[61,62], disperses relatively short distances from their development site, ^[63] and feeds on humans during daylight hours whenever the host is available ^[61,64]. It has also been identified that the majority of meals taken by the mosquito, are from residents or neighbors of the house they emerged as adults ^[65] where individuals who spent more time at home are more likely to receive Ae. aegypti bites in their homes than other household residents ^[66]. In Costa Rica, the National Institute of Statistics and Census has identified that women represent the part of the population that tends to stay longer periods of time at home (inside) (as shown in 2), most of that time, doing domestic chores ^[67]. Hence, making women more vulnerable to mosquito bites. In a country largely infested by the vector, sexual transmission dynamics has the potential to spread ZIKV in non-endemic Aedes established regions ^[60]. Undermining the risk of sexual transmission, may lead to an underestimate of the potential for disease persistence ^[20].

The population is divided 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:

where $N_m = S_m+E_m+U_m+D_m+R_m$, $N_f = S_f+E_f+U_f+D_f+R_f$, $N = N_m+N_f$ and $N_z = S_z+E_z+I_z$. The transmission rate is given by $\lambda(t) = \frac{N_z}{N}\kappa p_{hv}$, where $\frac{N_z}{N}$ is the vector to human ratio, $\kappa$ 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. For sexual transmission within hosts the transition from susceptible to infected is measured by $\beta_{m, f} = \beta \times \sigma_{m, f}$, where $\beta$ is the transmission rate (weeks$^{-1}$) and $\sigma_{m, f}$ (dimensionless) is the sexual activity indicator for males and females, respectively. In Figure 3, we illustrate the model transitions.

We are interested in exploring the effects of host availability, as well as sexual transmission of ZIKV in Costa Rica. Here, host availability is modeled via the parameter $\xi_{m, f}$ (male and female), where $\xi_{m, f}\in [0, 1]$. This parameter serves as a reduction factor in the transmission dynamics from mosquito to host.

Moreover, using the next generation approach by ^[68], we computed the basic reproductive number, $\mathcal{R}_0$. However, the analytic expression is highly complex, therefore, we computed its value in Mathematica using the parameter estimates.

3.   Parameter estimation

We used weekly reported cases from ^[39] 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 ZIKV confirmed cases and the model solution after estimating parameters.

3.1.   Parameter estimation using ABC

We use the Approximate Bayesian Computation (ABC) parameter estimation method. 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, this method is 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 ^[69] for an exhaustive summary of the subject).

Some of the parameters were 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 explored in previous articles ^[54,70]. 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 ($\mu$, $\kappa$ and $\delta$) 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 (2.1). We used a genetic algorithm ^[71,72] and the quasi-Newton L-BFGS-B method ^[73] to obtain the minimum. With the first algorithm we make an initial search in the parameter space and then 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 ^[74] through the EasyABC package of R ^[75]. 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 through the optimization procedure.

4.   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 from the 2016-2017 ZIKV outbreak and model solutions based on the estimated parameters and initial conditions (see Table 2), and different model solutions based on $\xi_m$ and $\xi_f$ (host availability), as well as the respective $\mathcal{R}_0$ values. In Figure 4a, the model parameters are taken from Table 2 where the values of $\xi_m$ and $\xi_f$ are assumed to be $0.5$ and $0.8$, respectively. The model curve (solid line) is fitted to the weekly reported data (dots). Given the parameter estimates, the model suggests that $\mathcal{R}_0 = 1.519$, a sustainable outbreak, which coincides with the situation in 2016-2017 in Costa Rica. Previous studies conducted in Costa Rica with dengue and chikungunya ^[70] have shown that the introduction of a new infectious disease in the population proves difficult to contain even when strict public health policies are in place. Moreover, in Figure 4b, the dashed line represents a reduction in female host availability ($\xi_f = 0.5$), $\mathcal{R}_0$ is reduced to $1.434$. This presents a situation where even when reducing host availability the epidemic is still sustainable ($\mathcal{R}_0 > 1$), albeit with a reduction in the number of cases. In Figure 4c, a more extreme case where $\xi_m = \xi_f = 0.4$, host availability for males and females is reduced significantly, we observe an even lower number of ZIKV cases (dashed line) and $\mathcal{R}_0 = 1.368$. Still, the epidemic remains sustainable ($\mathcal{R}_0 > 1$) and suggests that in order to have $\mathcal{R}_0 < 1$ host availability needs to be very low, and hence, not a viable public health strategy on its own.

In Figure 6a, we use the estimated parameters (in Table 2) to analyze host availability as a function of the basic reproductive number, $\mathcal{R}_0$. Here, we confirm that in order to have $\mathcal{R}_0 < 1$, host availability for both males and females needs to be below $3%$. Furthermore, the viability of reducing host availability ($\xi_m$ and $\xi_f$) is not realistic. In turn, other more viable strategies, or a combination of, are required in order to have $\mathcal{R}_0 < 1$. In Figures 6b, 6c we use the same analysis using the mosquito biting rate ($\kappa$) and mosquito mortality rate ($\mu$) with different host availability values.

The posterior densities of the three active parameters and $\mathcal{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 through the coefficient of variation of the posterior distribution, attains their minimum values for $\mathcal{R}_0$ and $\mu$, as can be seen in the last column of Table 2. Therefore, the degree of accuracy with which we are estimating $\mathcal{R}_0$ is reliable.

We performed an elasticity and sensitivity analysis on the $\mathcal{R}_0$ parameters (see Table 3). Our analysis revealed that $\mathcal{R}_0$ is most sensitive to the probability of infection from mosquito to human ($p_{vh}$) and the exposed rate of mosquitoes ($\alpha_z$) followed by the host availability parameters $\xi_m$ (male) and $\xi_f$ (female).

The importance of the female host availability is in part due to a larger proportion of the female population out of the workforce ^[76] and therefore they stay home (inside) for longer periods of time ^[67]. This scenario gives the mosquito more feeding opportunities throughout the day. In contrast, the male population spends most of their time in workforce activities ^[67], in turn, males tend to be more mobile giving the mosquito less opportunities to feed.

These parameters increase $\mathcal{R}_0$, therefore creating opportunities to reduce incidence. Taking into account that vector control measures are the most common practice of prevention/control, we argue that targeted strategies to individuals that stay home long periods of time along with vector control can prove viable strategies in order to reduce ZIKV incidence in Costa Rica.

Currently, the most common techniques used to combat vector borne diseases in Costa Rica, include epidemiological and entomological surveillance, environmental management, public education and chemical control ^[78]. In spite of these efforts, the numerous larval habitats identified in households of one of the most affected area of the country, shows that people may not be taking all the actions necessary to eliminate mosquito larval habitats ^[79]. In a disease where human behaviour plays a significant role, control strategies must depend on the specific socioeconomic context and behavioral characteristics of the population ^[80]. Reducing the mosquito biting rate involve strategies surrounding repellents, clothing that minimizes skin exposure during daylight hours, as well as, the use of window and door screens ^[81]. Nonetheless, given the characteristics of the population and household fixtures of the most affected areas ^[82], these possibilities are likely to incur in an economic burden for the population at risk. Lastly, the role of sexual transmission was shown to be minimal. However, it is important that this transmission mechanism be considered in order to implement the appropriate intervention strategies in future ZIKV outbreaks.

5.   Conclusions

After the rapidly spread of ZIKV through various Latin American countries, the introduction of the virus to Costa Rica was expected ^[83]. 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 ^[44]. 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. In the report, it was established that since the introduction of ZIKV to the country, the cases of microcephaly almost doubled, especially in the coastal regions, which could indicate that the observed increase may be associated, in some level, with ZIKV ^[84].

It is evident that ZIKV presents a 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, makes the planning of effective prevention/control measures a unique challenge to public health officials ^[85,86]. In Costa Rica, a country greatly infested by the Aedes mosquito, a better understanding of the transmission dynamics and ecology of the vector, can provide a guide to introduce more efficient prevention and control interventions, that allows for a better use of the human and economic resources available for the control of vector borne-diseases.

We estimated the basic reproductive number distribution using the estimated parameters 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. These results go in hand with previous studies that evidenced that female adult Ae. aegypti are more likely to bite household residents that spend more time inside their homes ^[66], 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 ^[87]. This becomes more significant, based on the fact that women constitute the majority of the population that stays inside ^[67], 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 propagation 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 ^[88]. A bottom-up communicative approach at all levels of the public health system ^[89] could provide a greater interest in the general population to implement the mosquito control strategies widely recommended by health 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 ^[90]. 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 ^[91]. Because of the ongoing disease 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. Prevention and 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.

Acknowledgments

The author is solely responsible for the views and opinions expressed in this research. We would like to thank the Ministry of Health of Costa Rica for providing the data regarding reported Zika cases.

Conflict of interest

All authors declare no conflicts of interest in this paper.