Estimating the reproductive number, total outbreak size, and reporting rates for Zika epidemics in South and Central America
Deborah P. Shutt Carrie A. Manore Stephen Pankavich Aaron T. Porter Sara Y. Del Valle


[
    Adam M. Oberman
    December 30, 2023
=====================




As South and Central American countries prepare for increased birth defects from Zika virus outbreaks and plan for mitigation strategies to minimize ongoing and future outbreaks, understanding important characteristics of Zika outbreaks and how they vary across regions is a challenging and important problem. We developed a mathematical model for the 2015 Zika virus outbreak dynamics in Colombia, El Salvador, and Suriname. We fit the model to publicly available data provided by the Pan American Health Organization,  using Approximate Bayesian Computation to estimate parameter distributions and provide uncertainty quantification. An important model input is the at-risk susceptible population, which can vary with a number of factors including climate, elevation, population density, and socio-economic status. We informed this initial condition using the highest historically reported dengue incidence modified by the probable dengue reporting rates in the chosen countries. The model indicated that a country-level analysis was not appropriate for Colombia. We then estimated the basic reproduction number, or the expected number of new human infections arising from a single infected human, to range between 4 and 6 for El Salvador and Suriname with a median of 4.3 and 5.3, respectively. We estimated the reporting rate to be around 16% in El Salvador and 18% in Suriname with estimated total outbreak sizes of 73,395 and 21,647 people, respectively. The uncertainty in parameter estimates highlights a need for research and data collection that will better constrain parameter ranges.





§ INTRODUCTION


Mosquito-borne diseases contribute significantly to the overall morbidity and mortality caused by infectious diseases in Central and South America. Newly emergent pathogens, such as Zika virus in 2015, highlight the need for data and models to understand the public health impact of associated diseases and develop mitigation strategies to combat their spread.
In particular, since Zika virus is a newly emergent pathogen in the Americas, its impact on the naïve population is relatively unknown.

The disease was first discovered in isolation from a rhesus macaque in the Zika forest of Uganda in 1974 <cit.>. While infrequent human cases were confirmed in later years in both Africa and Southeast Asia, it was not until April 2007 that an outbreak outside of these traditional areas occurred on Yap Island in the North Pacific <cit.>, and this was followed by another outbreak occurring in French Polynesia beginning October of 2013 <cit.>.  However, the most significant Zika outbreak began within Central and South America in 2015 <cit.> and is currently ongoing. Thus far it has resulted in an estimated 714,636 infections, which includes both suspected and confirmed cases within Latin American and Non-Latin Caribbean countries (accessed Jan 10, 2017) <cit.>.
This study focuses on the behavior of the current burgeoning epidemic in Colombia, El Salvador, and Suriname.


 Zika is transmitted to humans primarily through bites from infected Aedes aegypti and Aedes albopictus mosquitoes. The transmission is in both directions, that is, infected mosquitoes infect humans and infected humans infect mosquitoes. Upon transmission of the virus from mosquito to human, an individual will become infectious within 3 to 12 days.

Symptoms of infection include fever, rash, joint pain, conjunctivitis, muscle pain and headache.  Recovery from Zika virus disease may require anywhere from 3 to 14 days after becoming infectious, but once contracted humans are immune from the virus for life.  Many people infected with Zika may be asymptomatic 
or will only display mild symptoms that do not require medical attention.  An estimated 80% of persons infected with Zika virus are asymptomatic <cit.>. Thus, there is a high occurrence of under-reporting for confirmed or suspected Zika cases.  In fact, the number of infected individuals who report their symptoms is estimated to be between 7% and 17% <cit.> of the total number infected by the virus.


Zika can also be transmitted vertically, as a mother can pass the virus to her child during pregnancy, and this can lead to a variety of developmental issues. Most notably, Zika is a cause of microcephaly and other severe fetal brain defects <cit.>.  Recent evidence further suggests an association between the virus and a higher incidence of Guillain-Barré syndrome, a disease in which the immune system damages nerve cells causing muscle weakness and sometimes paralysis <cit.>.  The potential for sexual transmission of Zika virus has also been confirmed <cit.>.  For instance, the Center for Disease Control (CDC) has determined that Zika can remain in semen longer than in other body fluids, including vaginal fluids, urine, and blood.  It should be noted, however, that while these latter means of transmission exist, the number of new human infections produced in this way is relatively low compared to mosquito-borne infections <cit.>, and therefore, in the interest of simplicity and reducing the number of parameters to fit, we will neglect them in formulating a model.


In general, mathematical modeling has been extensively used to understand disease dynamics and the impact of mitigation strategies <cit.>. A few recent papers have developed models that focus on the behavior of the Zika epidemic.  For example, Gao et al. <cit.>, proposed an SEIR/SEI model to understand the effects of sexually transmitted Zika in Brazil, Colombia, and El Salvador.  The authors did not, however, distinguish between countries as they
assume that the three nations of interest share common parameter values.  Additionally, Towers et al. <cit.> presented a model that incorporates spatial heterogeneity in populations at a granular level.  Their model focused specifically on the spread of Zika within Barranquilla, Colombia and included a sexual transmission term, but concluded that the effects are not significant enough to sustain the disease in the absence of mosquitoes. An approximation of the basic reproductive number was obtained using maximum likelihood methods.  Finally, a study of the impact of short term dispersal on the dynamics of Zika virus is analyzed in <cit.>.  The model formulated within <cit.> does distinguish between asymptomatic and symptomatic infected populations, and focuses on the estimation of the reproductive number between two close communities.

In contrast to previous studies, we adopt a more global approach to understanding the dynamics of the current Zika epidemic and present a model that can be used to study disease transmission at the country level.  We identify differences between countries in regards to parameter values and reproductive numbers.  We discuss the appropriateness of country level analysis for Colombia, El Salvador, and Suriname and quantify the uncertainty within the resulting biological parameter values.

In the following section, we develop a Susceptible-Exposed-Infectious-Recovered (SEIR) /SEI type model which distinguishes between the reported infected population, who are considered symptomatic, and the unreported infected population, who may be asymptomatic or experience symptoms that are not severe enough to seek medical attention.  The data available from the Pan American Health Organization (PAHO) serves to motivate the use of a split infectious population.  The number of reported cases of Zika, both suspected and confirmed, is reported from PAHO by country.  Hence, we develop a model for three countries of interest assuming that the dynamics of the disease within each may be associated with different parameter values. In Section 3, we address the additional complication that the total population of a country cannot be used as the initial susceptible population due to the biology and bionomics of the Aedes species and human contact with them, varying in regards to temperature, humidity <cit.>, sanitation, demographics and elevation <cit.>.  Not everyone within a country is equally susceptible to the Zika infection due to geographic diversity; hence, we calculate the unique at-risk population within each country to use as the initial susceptible population. Using the deterministic model at baseline parameter values, the at-risk population we compute yields realistic initial conditions that compare well with PAHO data. With the initial population sizes for the epidemic in each country known and fixed, we embed the deterministic system into a stochastic process to quantify the uncertainty of the parameter values. This method allows us to use biologically valid parameter ranges and Approximate Bayesian Computation (ABC) methods to obtain parameter distributions that are distinct to each country.  Section 4 provides a detailed explanation of our implementation of the ABC method along with the strengths and weaknesses discovered in applying this method to our particular model.  We summarize the results of the obtained posterior distributions in Section 5. The last section is dedicated to discussion and conclusions.























§ A DETERMINISTIC ZIKA MODEL




 §.§ Vector-borne SEIR Model with Case Reporting

line = [draw, -latex']


The spread of Zika relies primarily on interactions between humans and mosquitoes.  In Figure <ref>, members of the at-risk human population, S_h, are bitten by an infectious mosquito and become exposed, or infected  without yet being infectious, at a rate λ_h.  Humans within the exposed population, E_h,  progress from this state to the infectious compartment at a per capita rate, ν_h.  Note that there is a portion of infectious humans who are either asymptomatic or experience less severe symptoms, and therefore go unreported.  The parameter ψ denotes the proportion of humans who seek medical assistance and are reported as either suspected or confirmed Zika patients.  Thus, there are two categories for infectious humans: I_r_h, the reported infectious population, and I_h, the unreported infectious population.  Both infectious populations then recover at the same per capita rate of γ_h, where 1/γ_h is the average time spent as infectious.  Upon recovery,  humans acquire lifetime immunity.  Similarly, a member of the susceptible mosquito population, S_v,  becomes exposed at a rate λ_v when susceptible mosquitoes bite an infectious human, resulting in transmission.  The exposed mosquitoes, E_v, transition to the infectious compartment at a per capita rate ν_v, where 1/ν_v is the average extrinsic incubation period.  Once a mosquito is infectious it will remain so for the duration of its lifespan which can range between 8 - 35 days <cit.>.  Finally, because of the long duration of infection for mosquitoes relative to their lifespan, demography (i.e., births and deaths) is included within their population dynamics. In total, the system is described by the following eight coupled, nonlinear ordinary differential equations:



    . dS_h/dt   =- λ_h(t) S_h
    dE_h/dt   =  λ_h(t) S_h - ν_h E_h
    dI_r_h/dt   = ψν_h E_h - γ_h I_r_h
    dI_h/dt   =(1-ψ)ν_h E_h -γ_h I_h
    dR_h/dt   =γ_h I_r_h + γ_h I_h
    dS_v/dt   =μ_v N_v- λ_v(t) S_v - μ_v S_v
    dE_v/dt   =λ_v(t) S_v -ν_v E_v -μ_v E_v
    dI_v/dt   =ν_v E_v - μ_v I_v.
    }




The population state variables are described within Table <ref> and parameters are presented in Table <ref>.
The quantities N_h and N_v represent the total populations of humans and mosquitoes within the model, and remain constant.
We note that the equation for the evolution of R_h(t) decouples from the system, as it is determined merely by computing the remaining population values.











































 §.§ Basic reproductive number

We define the basic reproductive number, ℛ_0, as the expected number of secondary infections by a single infectious individual over the duration of the infectious period within a fully susceptible population <cit.>.  As there is more than one class of infectives involved, we utilize the next generation method to derive an explicit formula for ℛ_0, defined mathematically by the spectral radius of the next generation matrix <cit.>. We follow the process in <cit.> and define x = [E_h, I_r_h,I_h, E_v, I_v, S_h, R_h, S_v]^T; thus reordering the presentation of populations from the original system to ensure our calculations possess the correct biological representation.  Let ℱ_i(x) be the rate of appearance of new infections in compartment i.  We indicate the rate of transfer of individuals out of compartment i as 𝒱^-_i(x)  and  the rate of transfer of individuals into compartment i by all other means as 𝒱^+_i(x).  Thus, our system can be expressed in a condensed version as ẋ_i = ℱ_i(x) - 𝒱_i(x) where 𝒱_i(x) = 𝒱^-_i(x)-𝒱^+_i(x) for i = 1,..., 8.





Next, we compute F = [∂ℱ_i/∂ x_j(x_0)] and V = [∂𝒱_i/∂ x_j(x_0)] for the exposed and infected compartments, namely for  1 ≤ i,j ≤ 5, where x_0=[0,0,0,0,0,H_0,0,K_v] is the disease free equilibrium state with H_0 and K_v being the initial population sizes of humans and mosquitoes respectively, and obtain the following 5 × 5 matrices:


    𝐅 = (
    2[               0               0               0               0         β_hvσ_v;               0               0               0               0               0;               0               0               0               0               0;               0 β_vhσ_v K_v/H_0 β_vhσ_v K_v/H_0               0               0;               0               0               0               0               0 ])
            𝐕 = (
    2[        ν_h          0          0          0          0;     -ν_h ψ        γ_h          0          0          0; -ν_h(1- ψ)          0        γ_h          0          0;          0          0          0   μ_v +ν_v          0;          0          0          0       -ν_v        μ_v ])



Hence, we calculate the Reproductive Number as:

    ℛ_0 := ρ(FV^-1) = σ_v √(K_v β_hvβ_vhν_v)/√(H_0 γ_h μ_v(μ_v+ν_v))=√(R_hv R_vh)


where ρ(A) represents the spectral radius of the matrix A, and we have defined the quantities R_hv = (ν_v/μ_v + ν_v)(σ_v/μ_v)β_hv and R_vh = (K_v/H_0)(σ_v/γ_h)β_vh.


Here, R_hv is the expected number of secondary infections in a fully susceptible human population resulting from one newly introduced infected mosquito.  It is composed of the product of three terms.  The first term, ν_v/μ_v + ν_v, represents the probability that an exposed mosquito will survive the extrinsic incubation period.  The second term, σ_v/μ_v, is the number of human bites an infectious mosquito would make if humans were freely available.  The third term, β_hv, is the probability of transmission occurrence given that a human is bitten by an infected mosquito.  The number of secondary infections in a fully susceptible population of mosquitos resulting from one newly introduced infected human is represented by R_vh.  This value is also formed by the product of three terms.  The first, K_v/H_0, is the vector to host ratio.  The σ_v/γ_h term is the maximum number of bites an infectious human will experience before recovery without impediment to mosquito bites.  Finally, given that a susceptible mosquito bites an infectious human, β_vh is the probability of transmission from human to mosquito.  The type reproductive number, or expected number of secondary human cases resulting from one newly infectious human, is ℛ_0^T := (ℛ_0)^2 <cit.>.






§ INCIDENCE RATES AND PARAMETER ESTIMATION

Although the traditional approach when modeling disease dynamics is to assume the total population within a country is susceptible, we deviate from this convention. Specifically, we assume the susceptible population depends on the biology and binomics of the Aedes species as well as the country's geography and uniquely calculate the at-risk population (described in the following section). We then use the at-risk population as the size of the initial susceptible population within simulations.   While the deterministic model (<ref>) captures the dynamics of the epidemic fairly well, we follow the opinion of the authors in <cit.>, namely that more attention should be paid to how uncertainty in parameter estimates might affect model predictions.  Thus, in Section <ref>, we focus on incorporating and quantifying the uncertainty of the disease process and the parameter values by means of embedding the deterministic model into a stochastic process.



 §.§ Calculated At-Risk Population

The data utilized herein was reported by the Pan American Health Organization (PAHO) <cit.>, and is given by the number of cases, both confirmed and suspected, of Zika per week at the country level for Colombia, El Salvador and Suriname.  Simulations performed by using the entire country population as the number of initially susceptible humans mischaracterized the disease dynamics, leading to overestimates in the final size of an epidemic, see Appendix (Figure <ref>).  Since dengue and Zika occur in the same areas, share a common vector and have similar asymptomatic rates, we calculate the at-risk population size per country for a Zika outbreak based on historical data for dengue from 1995 to 2015 within Colombia, El Salvador, and Suriname <cit.>.  The year of highest incidence for dengue provides an approximation for the number of susceptible individuals in a fully naïve at-risk population, which coincides with the dynamics of a newly emerging pathogen, such as Zika, that would spread rapidly in a completely susceptible population.  The World Health Organization (WHO) released a report on dengue <cit.>, stating “...available results suggest that the actual number of cases of dengue may range from 3 to 27 times the reported [dengue] number."  Hence, we use a scalar multiple of the reported number of dengue cases during the highest incidence rate year to obtain a reasonable at-risk population count for the number of initial susceptible individuals with regards to a Zika epidemic.  Within each of the three countries, the at-risk population value is strictly less than the total country population.






Table <ref> indicates the year with the highest incidence rate and the reported number of cases for that year in each country of interest.  See the Appendix (Figure <ref>) for all historical data on dengue incidence rates for Colombia, El Salvador and Suriname.





We simulated Model (<ref>) using the ode45 solver in MATLAB with chosen baseline parameter values (Table <ref>) while assuming that a single infectious human exists at the start of the epidemic. The at-risk population size for Colombia is 2.75 times the number of reported dengue cases from 2010.  We estimated the multiplier for the at-risk population based on the best fit for the currently ongoing epidemic. The at-risk population size for El Salvador is 1.425 times the number of reported dengue cases from 2014, while the at-risk population size for Suriname is 7.75 times the number of reported dengue cases from 2005.  The initial size of the susceptible mosquito population is double the number of the at-risk population per country.   When starting our analysis with only one infected human, the simulated epidemic takes several weeks to ramp up to a detectable (reportable) level. If we begin simulations on the first day that a case is reported, the simulated peak occurs after the reported peak. Thus, we use shifted initial conditions (Table <ref>) that correspond to the population state sizes obtained approximately on day 38, 118, and 50 from the original simulation using Model (<ref>) of the epidemics in Colombia, El Salvador and Suriname, respectively.  We then repeat the process of simulating Model (<ref>) using the shifted initial conditions to obtain solutions for the number of reported cases whose peak reporting weeks align more closely in time with that of the data. The precise values of the shifted initial conditions for each country, found in Table <ref>, are of the form [S_h,E_h,I_r_h,I_h,R_h,S_v, E_v,I_v].







The number of reported cases from our simulations, using the calculated at-risk population and shifted initial conditions, are compared to the PAHO Zika data in Figure <ref>.  The solutions, which are similar to the data in peak and general shape, identify reasonable initial conditions.  Although the parameter ranges considered are biologically reasonable, it is unknown which values within these ranges are the most accurate.  To obtain an expected value for a given parameter and describe the associated uncertainty within this quantity, we embed the deterministic Model (<ref>) into a stochastic process enabling statistical inference.  This process is described in Section <ref>.

















































 §.§ An Embedded Stochastic Model



Previous literature has provided biologically relevant parameter ranges (see Table <ref>) for the model.  However, to obtain insight as to the distribution of these parameters across their ranges, we embed the deterministic system of ordinary differential equations (<ref>) within a discrete-time stochastic process to obtain a corresponding stochastic model (<ref>, <ref>) and perform an analysis in a Bayesian paradigm. This approach adds a probabilistic component to both the contact/infection process and the disease progression process by capturing uncertainty at each modeling stage, rather than just within the contact/infection process <cit.>.  Therefore, embedding deterministic models into stochastic processes serves to (a) create a stochastic model informed by population-level dynamics which includes uncertainty in the entire disease process rather than just due to data collection, (b) more adequately capture uncertainty present in the modeling framework, which is particularly important in small to medium sized epidemics, and (c) provide a convenient framework for fast computation of model parameters.  The process by which we embed Model (<ref>) into a stochastic process, <cit.>, is summarized in the following paragraphs.


Previously, we obtained specific rates for the transfer of both humans and mosquitos between their respective compartments, and these can be used to inform the stochastic analogues of these parameters.  Throughout, we impose that the new model be conservative, i.e. S_h + E_h + I_r_h + I_h+R_h = N_h and S_v + E_v +I_v = N_v where N_h and N_v are the total human and mosquito populations, respectively.  The assumption of conservation within the human population is plausible as the time span of the epidemic is much smaller than the average lifespan of a human.  We hold the mosquito population as constant for the convenience of calculations and analysis.

To approximate a continuous SEIR type model, we account for all events/transitions that may occur during time interval (i, i+h_i] to be assigned to the related compartment on day i.  Thus, the rate of change for a given population size can be approximated by the difference between the previous and current time steps.  The model now takes the form


    . 
    S_i+h_i   =S_i-E^*_i+h_i
    
    E_i+h_i   =E_i+E_i+h_i^*-I_i+h_i^*
    
    Ir_i+h_i   = Ir_i +ψ I_i+h_i^* - RIr_i+h_i^*  
    
    I_i+h_i   =I_i+(1-ψ)I_i+h_i^*-RI_i+h_i^*
    
    R_i+h_i   = R_i + RIr_i+h_i^* + RI_i+h_i^*
    
    Sv_i+h_i   =Sv_i +dEv_i+h_i^*+ dIv_i+h_i^*- Ev_i+h_i^*
    
    Ev_i+h_i   =Ev_i + Ev_i+h_i^* - Iv_i+h_i^* - dEv_i+h_i^*
    
    Iv_i+h_i   =Iv_i+Iv_i+h_i^*-dIv_i+h_i^*.
    }














We index time by i, and the temporal offset, h_i, which may not be constant but is known.   All quantities represent counts and quantities denoted by an asterisk represent transition counts. In regards to the terms ψ I_i+h_i^* and (1-ψ)I_i+h_i^*, the parameter ψ may yield values which are not integer counts.  To be more liberal with the reporting size we calculated the smallest integer not less than the corresponding value of ψ I_i+h_i^*. To be more conservative with under-reporting size we calculated the largest integer not greater than the corresponding value of (1-ψ)I_i+h_i^*.   When individuals can transition into a compartment via multiple routes, (e.g., recovered individuals can recover either with or without being reported as infected, RIr_i+h_i^* and RI_i+h_i^*), two letters are used to denote the transition.  In these cases, the first letter denotes to which compartment the individual transitions and the second denotes the compartment from which the individual transitions.  The transition compartments which represent birth/death counts of the mosquito popultion on day i+h_i are denoted, dEv_i+h_i^* and dIv_i+h_i^*.  Note that these values are drawn based on the calculated population sizes at time i+h_i.  The compartments are labeled as follows: S - susceptible humans, E - latently infected humans, Ir - infectious reported humans, I - infectious unreported humans, R - recovered humans, Sv - susceptible mosquitoes, Ev - latently infected mosquitoes, Iv - infectious mosquitoes.

Finally, the stochastic components of the model are given by

    . 
    E^*_i+h_i   ∼Bin (S_i,1-exp(-λ_h h_i) ),   
    I_i+h_i^*∼Bin (E_i,1-exp(-ν_h h_i)  )
    
    RIr_i+h_i^*   ∼Bin (Ir_i, 1-exp(-γ h_i)  ),   
    RI_i+h_i^*∼Bin (I_i,1-exp(-γ h_i)  )
    
    Ev_i+h_i^*   ∼Bin (Ev_i,1-exp(-λ_v h_i)),   
    Iv_i+h_i^*∼Bin (Iv_i,1-exp(-ν_v h_i))
    
    dEv_i+h_i^*   ∼Bin (Ev_i+h_i,1-exp(-μ_h h_i)),   
    dIv_i+h_i^*∼Bin (Iv_i+h_i,1-exp(-μ_h h_i)).
    }




















§ ABC ALGORITHM AND COMPUTATION


Previous investigations <cit.> have established biologically accepted ranges for the parameters Θ = [σ_v, β_hv, ν_h, ψ, β_vh, ν_v, μ_v] used within (<ref>), but the conversion of Model (<ref>) to Model (<ref>, <ref>) will incorporate uncertainty into the disease process, and thus the values of these parameters.  The stochastic model (<ref>, <ref>) also allows for Bayesian inference on parameter posteriors which have the form

    f(Θ | Y) ∝ f(Y | Θ)π(Θ)


where f(Y | Θ) is the stochastic data model, which for fixed values of Θ can be used to generate the random epidemic process, and π(Θ) is the prior distribution of the parameters.   Thus, given both the data model and the prior distribution, we are able to calculate the posterior distribution f(Θ | Y) up to a proportionality constant. This serves to update the distribution of the biologically accepted parameter ranges based on the actual epidemic data which we denote as Y.



The data model is constructed as a product of binomials with different sample sizes and probabilities across every time point. Determining a Maximum Likelihood Estimate (MLE) to provide a point estimate and standard deviation for parameters may provide estimates outside of the valid biological ranges.   In SEIR models, Markov Chain Monte Carlo (MCMC) methods produce parameter autocorrelations in chains which becomes problematic for tuning. We can avoid these obstacles by using Approximate Bayesian Computation (ABC).  This method was introduced by <cit.> to obtain an approximation of the true posterior distribution, f(Θ | Y).
ABC samples from the posterior by randomly selecting parameter values from the prior that could adequately generate the data.  In particular, random draws, Θ^*, from the prior distribution produce generated data sets, X, which are then compared to a given epidemic data set, Y, by means of a chosen distance metric, ρ(X,Y).  Those values of Θ^* that generate data sets which fit the given data will be accepted as valid draws from the posterior distribution, and this implicitly conditions the posterior on Y.  The algorithm of the computation is as follows, where N is the total number of accepted generated data sets X:




For j ≤ N


  * Draw Θ^*∼Unif(a_k,b_k), where (a_k,b_k) is the corresponding parameter range found in Table <ref> for all k = 1,...,7. 

  * Generate X, time series data of the number of reported infectives, from Model (<ref>, <ref>)

  * Calculate fitness of data using ρ(X,Y)

  * Set Θ_[j]←Θ^* if ρ(X,Y) ≤ϵ  and set j ← j+ 1 else return to Step 1.  




The metric, ρ(X,Y), is considered a distance between the generated data set, X, and the observed data set, Y, typically based on sufficient statistics of the parameter space Θ.  If ρ(X,Y) is small enough, i.e., ρ(X,Y) ≤ϵ for some fixed small value ϵ > 0, then

    f(Θ |ρ(X,Y) ≤ϵ ) ≈ f(Θ | Y).

One often uses sufficient statistics in defining ρ(X,Y). These are statistics for which the data distribution conditioned on the sufficient statistic is free of Θ (e.g., the sufficient statistic contains all information about Θ that the full data set contains).  By definition, the full data set is a sufficient statistic for Θ.  Since a sufficient statistic of lower dimension than the full data set cannot be readily computed for our model, we consider a pointwise envelope metric comparing point by point the L_1-norm at every time step.  Thus, if every data point of X and Y is close, then it directly follows that any statistics computed from X and Y, including sufficient statistics, will also be close.  For the ABC method, the exact posterior distribution of  Θ can be found by accepting the simulated data set in which X = Y.  This is computationally infeasible, however, so we instead consider X, such that  X_t/Y_t∈ (ϵ_1,ϵ_2) for all time steps t with a corresponding peak week and epidemic duration to Y, as an acceptable epidemic.  We call (ϵ_1,ϵ_2) the envelope of tolerance around the observed data set, Y, which is a common method for assessing a stochastic SEIR model fit <cit.>. In Figure <ref>, the observed data set is compared to randomly generated data sets without calculating a metric for validation of the drawn parameter values.  Figure <ref>  shows that many randomly drawn and biologically accepted values of Θ generate epidemic outcomes which have distinctly different characteristics than the observed epidemic.  Conversely, Figure <ref> demonstrates that epidemic outcomes of the accepted Θ^* values generate epidemics with similar characteristics to that of the observed epidemic in regards to the total number of new infections each day, peak week occurrence, and peak value.  Figure <ref>(b) and <ref>(c) show generated epidemics similar to the data in both peak and duration.  Note that this does not occur in Figure <ref>(a).  The accepted parameter values obtained using the ABC method yield outbreaks which vary greatly in peak occurrence, leading us to conclude that the dynamics of the country-level data are different than those in the deterministic model. Therefore the ABC method does not generate infectious curves similar to the Colombia data. Figure <ref>(a) was generated from acceptances using an envelope of (1/20,20).  A tighter envelope of (1/10,10) was computed; however, the method ran for approximately 35 days to obtain the same number of acceptances as the (1/20,20) envelope with no visible changes in the same plot as Figure <ref>(a).  Because acceptances for tighter envelopes are not possible, we can conclude the model is improperly specified to capture the dynamics of Colombia at the country level. Reasons for this outcome are discussed in Section <ref>.









































Figure <ref> depicts the effect of generating accepted epidemics with smaller and smaller envelopes of tolerance on the posterior distributions.  We plot the kernel densities of both the reporting rate, ψ, and reproductive number, ℛ_0, distributions for El Salvador and Suriname.  While only two parameters are shown here, the same observations are found in the kernel density plots of the other parameters of interest.  We see in Figure <ref> that envelopes of various size produce similar but different kernel densities, with the tighter envelopes producing a more peaked density distribution.  A slight shift in the peak of the kernel densities is observed in Figures <ref>(a) and <ref>(d).  This would indicate the need to shrink the envelope further to obtain an estimation that is closer to the true posterior distribution.  However, attempts to generate data sets satisfying tighter envelopes could not produce acceptances, indicating that these envelope sizes are the best possible fits of the stochastic SEI_rIRS_vE_vI_v model to the observed country-level data.  Further analysis would require refining the dynamics of Model (1) and Model (2,3) or using data at a more refined spatial scale.





§ RESULTS



We assigned uniform distributions to the accepted biological ranges found in Table <ref>. After 10,000 acceptances from the ABC algorithm, the histograms and kernel densities for selected parameters from El Salvador are given in Figures <ref> and <ref>.  Selected parameters from Suriname are given in Figures <ref> and <ref>.  Figure <ref> displays the number of reported cases from the generated epidemics based on 10,000 accepted parameter values. The time series of the number of total cases during the generated epidemics of both reported and unreported infectives is found in Figure <ref>.  The histogram and kernel density plots of the total number of Zika virus cases, reported and unreported are shown in Figure <ref>.


































































§ CONCLUSIONS AND DISCUSSION



Within the present study, a new stochastic model was formulated to describe the spread of the Zika virus within Colombia, Suriname, and El Salvador. The variability in per capita susceptibility within each country was introduced by uniquely calculating the at-risk population based on historical data for dengue virus, whose epidemic characteristics are similar to that of Zika.  The initial population state sizes were fit to this data in order to create similar values of weekly reported cases of Zika, both suspected and confirmed, to those reported by the PAHO using baseline values of parameters to evaluate Model (ode).  Once the at-risk population was estimated, the initial population state sizes for these nations were fixed in order to estimate parameter values and obtain more informative distributions than the accepted biological parameter ranges.  The deterministic system (ode) was then embedded into a stochastic process to obtain a more general stochastic model (<ref>, <ref>).  For each of the three nations, an ABC algorithm was implemented using (<ref>, <ref>) to compute approximate posterior distributions of the parameters conditioned on the data.  By obtaining these posterior distributions, the uncertainty in parameter values for each country can be quantified by way of informative statistics, such as the mean, mode, median and variance, to more accurately describe rates within the system.


Properties of the disease within El Salavador and Suriname were accurately described by the model.  El Salvador is estimated to have a mean reporting rate of 16.5% which is near the upper bound of the previous estimate of <cit.> with a credible interval of [12.5%, 22%].  The mean values of the forcing terms, λ_h and λ_v, are approximately 0.28 and 0.31 with credible intervals of [0.0191, 0.7782] and [0.0119, 0.9244], respectively.  Suriname has a mean reporting rate greater than the predicted interval of <cit.> at 18.8% with a credible interval of [13%, 27%]. The forcing terms of Suriname, λ_h and λ_v, had mean values of 0.17 and 0.43 with credible intervals of [0.0165,0.4054] and [0.0429,1.1241], respectively.

The basic reproductive number, ℛ_0, was defined such that ℛ_0^2 yields the number of secondary human infections within a fully susceptible population arising from a single new human infection.  In Suriname, the mean value of ℛ_0^2 was 5.31, while the mean reproductive number was 4.35 in El Salvador.  These quantities are similar to other predictions for the mean of ℛ_0^2 found in <cit.>.

Though this analysis has provided additional insight into the spread of the disease, our methods were unable to accurately estimate the aforementioned statistics for Colombia (see Section <ref>) as the ABC method revealed a poor fit for the data obtained from this nation.  One possible reason for this poor fit could be the appearance of a second peak in the epidemic, see Figure <ref>(a).  While there is a distinct and large peak during EW 26, a smaller second peak occurs during EW 34.  This second peak may result from reporting error as this data was in no way cleaned once obtained from the PAHO website.  Another possible explanation is a second outbreak in a disjoint location from the first outbreak.  In particular, the extreme topographical variations within Colombia could lead to a delay in the spread of the disease across the entire country.  Thus, we conclude that the epidemiological characteristics of Zika in Colombia must be studied at a more granular spatial level, differentiating amongst regions or even counties and cities.

We also estimated the total number of infected people in El Salvador to be 72,721 (with an estimated 12,107 reported cases) and the total number of people infected in Suriname to be 21,390 (with an estimated 4,132 reported cases). So, about 95% of the at-risk (high-risk) populations were infected by the end of the outbreak. The predicted reported country-level incidence for El Salvador is 0.0019 and for Suriname is 0.0073. However, the true country-level incidence for El Salvador is predicted to be 0.0119 while for Suriname it is 0.0396 assuming total populations of 6,117,145 and 540,612, respectively. From a public health perspective, our model indicates that about 6 times as many people were infected than were reported in El Salvador and 5 times as many were infected than were reported in Suriname. Depending on the percent of the at-risk population that was pregnant during the outbreak, our model suggests that a larger number of birth defects than indicated by the reported number of cases could be expected in these countries. Interestingly, although the values for the mosquito extrinsic incubation period and mosquito lifespan had relatively wide ranges even among mean, median, and mode values, the probability of a mosquito surviving the incubation period (Tables 12 and 13) was quite similar across statistics and countries at about 60%. This corresponds to the extrinsic incubation period lasting around 2/3 of the average mosquito lifespan. A major difference between predicted parameter values for El Salvador and Suriname occurs in the transmission probabilities. While the values for β_hv and β_vh for El Salvador were similar (both estimated to be about 0.42), in Suriname, the mosquito-to-human probability of transmission, β_hv was consistently less than half that of the human-to-mosquito transmission, β_vh, with a median of 0.21 and 0.56, respectively. Since these terms capture many intrinsic uncertainties in the transmission process, it is hard to interpret the meaning of this difference. It could be an artifact of the model or could indicate a reduced efficiency of the mosquitoes in Suriname in passing on the virus.

In conclusion, our research provides important parameter estimates for the spread of Zika in El Salvador and Suriname, along with uncertainty quantification and credible intervals for those parameters. We estimated the basic and type reproduction numbers and the total number of people infected - quantities needed to inform assessments of economic cost and risk, among other factors. We found that the type reproduction number is higher in Suriname than in El Salvador, indicating a higher risk in Suriname. This could be explained by differences in climate between the two countries or in other socio-economic or geographic factors affecting mosquito-borne disease transmission. Additionally, our methods could be applied to other countries or regions experiencing outbreaks to estimate region-specific parameters and provide decision makers with important information about surveillance and control both at present and in the future. Another advantage of this method is that it can indicate what scale is appropriate for these calculations. For example, we found that a country-level analysis of the Colombia data was not appropriate. In the future, it would be interesting to apply the model to a regional Colombia data set.


In future studies, the initial population state sizes obtained in the calculations of the respective at-risk populations could be considered parameters themselves. Hence, one could perform a statistical analysis on acceptable ranges for such initial values to quantify the uncertainty of these populations. For instance, this could be done by generating a posterior distribution for the initial susceptible population for each country, as well as the parameters in the system.  Additionally, the mosquito population for each country was assumed to be constant for convenience.  However, utilizing a more realistic model of the total mosquito population which changes in time (similar to the methods in <cit.>) may yield different results, and this would also be a suitable direction for future research. Finally, the recovery period, 1/γ_h, was held constant in the current study due to the assumption that its mean value had been medically established.  Still, another investigation using the methods developed herein, but considering a uniformly-distributed prior for the recovery period, may provide better insight into the distribution and expected value of this period. We conclude that additional studies are needed to fully understand Zika virus transmission dynamics. However, our research suggests that the reporting rates in El Salvador and Suriname are quite low and thus, additional surveillance systems may be needed to measure the true burden of Zika in these countries.











§ ACKNOWLEDGMENTS

This work was supported by NSF SEES grant CHE - 1314029, NSF RAPID (DEB 1641130), and NSF EDT grant DMS-1551229. SP was partially supported by NSF under grant DMS-1614586.  SD was partially supported by NIH/NIGMS/MIDAS under grant U01-GM097658-01. LANL is operated by Los Alamos National Security, LLC for the Department of Energy under contract DE-AC52-06NA25396. Approved for public release: LA-UR-17-20963. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.



§ APPENDIX



 §.§ Results of Model <ref> using total country populations as the initial size susceptible populations









 §.§ Other parameter distributions










































    

    
        










    




 §.§ Other Reproductive Number distributions






acm


