Potential Impact of Sexual Transmission on Ebola Virus Epidemiology: Sierra Leone as a Case Study

Background Sexual transmission of Ebola virus disease (EVD) 6 months after onset of symptoms has been recently documented, and Ebola virus RNA has been detected in semen of survivors up to 9 months after onset of symptoms. As countries affected by the 2013–2015 epidemic in West Africa, by far the largest to date, are declared free of Ebola virus disease (EVD), it remains unclear what threat is posed by rare sexual transmission events that could arise from survivors. Methodology/Principal Findings We devised a compartmental mathematical model that includes sexual transmission from convalescent survivors: a SEICR (susceptible-exposed-infectious-convalescent-recovered) transmission model. We fitted the model to weekly incidence of EVD cases from the 2014–2015 epidemic in Sierra Leone. Sensitivity analyses and Monte Carlo simulations showed that a 0.1% per sex act transmission probability and a 3-month convalescent period (the two key unknown parameters of sexual transmission) create very few additional cases, but would extend the epidemic by 83 days [95% CI: 68–98 days] (p < 0.0001) on average. Strikingly, a 6-month convalescent period extended the average epidemic by 540 days (95% CI: 508–572 days), doubling the current length, despite an insignificant rise in the number of new cases generated. Conclusions/Significance Our results show that reductions in the per sex act transmission probability via abstinence and condom use should reduce the number of sporadic sexual transmission events, but will not significantly reduce the epidemic size and may only minimally shorten the length of time the public health community must maintain response preparedness. While the number of infectious survivors is expected to greatly decline over the coming months, our results show that transmission events may still be expected for quite some time as each event results in a new potential cluster of non-sexual transmission. Precise measurement of the convalescent period is thus important for planning ongoing surveillance efforts.


Introduction
Recent reports suggesting the potential for sexual transmission of Ebola virus from convalescent survivors have raised a number of important questions regarding its impact on the final phase of the epidemic in West Africa [1,2]. Even once the worst hit countries of Guinea, Liberia, and Sierra Leone are declared free of Ebola virus disease (EVD), rare cases may still arise from the large number of remaining survivors. Importantly, sexual transmission is dependent on the frequency of infections rather than the density of available hosts, allowing chains of transmission to persist at low susceptible densities where non-sexual transmission would typically fail to occur [3]. Perhaps the most crucial element for bringing the epidemic to an end is maintaining vigilance in the community by preventing-or quickly responding to-new chains of transmission. Thus, it is important to investigate the potential impact of convalescent sexual transmission on the transmission dynamics in general, and on the tail of the epidemic in particular, to understand how long that vigilance might remain critical.
Follow-up studies on survivors of the 1995 outbreak in the Democratic Republic of Congo [4] and the 2000 [5] and 2007 [6] outbreaks in Uganda have raised awareness of what is now being termed "post-Ebola syndrome" (post-Ebola sequelae)-debilitating illnesses from myalgia to uveitis-which can persist for at least 21 months after the onset of symptoms. Though the virus is no longer detected in the blood after acute EVD symptoms disappear, active (replicating) virus has been documented in ocular fluid, rectal fluids, vaginal fluids, and semen [1,4,7,8]. Transmission to sexual partners was never confirmed in earlier outbreaks, but was suspected to have occurred in at least one instance [4]. Similarly, cases of sexual transmission of other hemorrhagic fever infections, notably by the closely related Marburg virus, have been suspected in the past [9,10]. Studies from the West African outbreak, showing viremia in semen 4-6 months after onset of symptoms in 65% of men tested (7-9 months in 26%) [1] and presenting molecular evidence of sexual transmission from a survivor 179 days after onset of symptoms [2], suggest that sexual transmission from convalescent men can and does occur.
Sexual transmission of Ebola virus from convalescent survivors is likely a rare event, but researchers have warned that it should be considered in epidemiological models that are used to predict the trajectory of an outbreak [11]. Without aiming to make a predictive model but rather to understand what aspects of the epidemic may be affected by inclusion of sexual transmission, we devised a novel formulation of the mathematical model for EVD transmission: SEICR (susceptible-exposed-infectious-convalescent-recovered), which includes a component for convalescent sexual transmission from convalescent survivors who maintain active Ebola virus replication. We illustrated the model by fitting it to weekly EVD incidence from Sierra Leone, the largest population of recovering survivors from the current West Africa epidemic. We performed sensitivity analysis to understand the influence of key unknown parameters, such as the duration of the convalescent period and the transmission probability per sexual contact. Considering the stochastic nature of such rare sexual transmission events, we also performed Monte Carlo simulations to explore the impact of sexual transmission on the epidemic tail in Sierra Leone.

Transmission model
We extended a SEIR (susceptible-exposed-infected-recovered) modeling framework, which has been extensively used to describe EVD transmission [12][13][14], by adding a component for sexual transmission from convalescent survivors who maintain active Ebola virus replication (Fig 1). The resulting SEICR model has five states: susceptible, S, exposed, E, symptomatic and infectious, I, convalescent, C, fully recovered and immune, R, and dead, D. The model is represented by the following set of ordinary differential equations (ODEs): where N = S + E + I + C + R denotes the total population size. We assumed the non-sexual transmission rate, β(t), to be initially constant (β 0 ) before it starts to decay exponentially due to the effect of control interventions and behavior change after timeτ: β(t) = β 1 + (β 0 -β 1 )e -k(t-τ) [12]. The sexual transmission parameter, β s , can be described as the product of two parameters (β s = ηq) that we will consider separately: η is the per sex act transmission probability of Ebola virus from convalescent men, and q is the daily rate at which they engage in sexual intercourse. The number of convalescent individuals are scaled by p, which is the proportion of convalescent survivors who are sexually active men. 1/σ and 1/γ represent the average durations of incubation and symptomatic infection, respectively. f is the case fatality rate. The average duration after which convalescent patients recover completely and shed no further replicating Ebola virus from their body is given by 1/α. We assumed that sexual transmission is frequencydependent [3,15,16], i.e., the probability that the sexual partner of a convalescent man is susceptible is given by S/N. The basic reproductive number, R 0 , for the SEICR model can be calculated using the nextgeneration matrix method [17,18] and is given by where S 0 is the initial number of susceptible individuals (see S1 Appendix). When α goes to infinity or either β s = 0 or p = 0, the equation reduces to the basic reproductive number in absence of sexual transmission: R 0;N ¼ bS 0 g . Thus, the second term represents the contribution of sexual transmission from convalescent patients to the overall R 0 :

Model parameters
Since the number of sexual transmission events was likely to be small and little information is currently available on the nature of each transmission event, we fitted only the non-sexual (SEIR) deterministic EVD transmission model to weekly incidence of confirmed and probable cases in Sierra Leone as reported in the WHO patient database [19] (S1 Fig). The data set was extended with weekly incidence from the situation report for the most recent weeks when no data was available in the patient database. In order to account for variability in the accuracy of reporting, we assumed that the number of reported cases follows a negative binomial distribution with mean as predicted by the model and dispersion parameter ϕ [20]. We derived maximum likelihood estimates (MLE) of the following model parameters [14,21]: the baseline transmission rate β 0 , the time τ at which transmission starts to drop, the rate k at which transmission decays, and the dispersion parameter ϕ. For the fitting procedure, we assumed that there were no sexual transmission events, i.e., we set β S to zero. The basic reproductive number in absence of sexual transmission is R 0,N = β 0 N 0 /γ, and the reproductive number in presence of partially effective control interventions is R 1,N = β 1 N 0 /γ, with N 0 being the population size of Sierra Leone. We explored value ranges for sexual transmission parameters (Table 1) based on information from the current epidemic [22] and studies of human immunodeficiency virus [23,24]. Remaining parameters were based on published values from the literature ( Table 1).  [12][13][14]. The red elements (convalescent individuals and additional transmission term) were added to account for sexual transmission. Deterministic model and sensitivity analysis We solved the system of ODEs numerically using the function 'ode' from the 'deSolve' package in the R software environment for statistical computing [29]. We compared the following response variables of the model: the epidemic peak number of exposed, E, acute, I, and convalescents, C, cases; the cumulative number of EVD cases, deaths, and recoveries; the date at the epidemic peak; the daily and cumulative incidence of sexual transmission; and the date at which the last symptomatic case either died or entered into convalescence ("day of last case"). We defined the day of last case as the time when the number of symptomatic and infectious individuals, I, dropped below 0.5. We considered the following parameters for the sensitivity analysis: the per sex act transmission probability of Ebola virus from convalescent men (η), the proportion of convalescent survivors who are sexually active men (p), the rate at which they engage in sexual intercourse (q), and the rate at which convalescent patients recover completely and shed no further replicating Ebola virus from their body (α). The sensitivity of the response variables to changes in R 0 was explored simultaneously as a comparison. We generated 1000 parameter combinations from the uniform ranges, log-transformed [0.5x -2x] for the parameter values for η, p, q, and α, given in Table 1 via Latin hypercube sampling using the Huntington and Lyrintzis correlation correction method (function 'lhs' from R package 'pse') [30]. We then calculated partial rank correlation coefficients (PRCCs) using 50 bootstrap replicates [31].

Monte Carlo simulations
We performed stochastic simulations of the SEICR model with and without sexual transmission using Gillespie's algorithm [32]. We specifically investigated the following response variables from the simulations: the cumulative number of EVD cases, the size and date of the epidemic's peak incidence (daily number of new symptomatic infections), and the date of last case (last day that symptomatic infections, I, fell below 1). Summary statistics were based on the results of 1000 simulation runs for each transmission scenario. We calculated the average of the peak and total cumulative number of EVD cases by including all simulations runs, i.e., also the simulations that rapidly go extinct. In contrast, the average dates of the epidemic peak and last case were based on the simulated epidemic trajectories over which 50 or more cases were accumulated.

Contribution of sexual transmission to overall R 0
Assuming a conservative baseline scenario (η = 0.001 and 1/ α = 3 months; Table 1), the reproductive number of a convalescent infection, R 0,C , is 0.0024. This corresponds to only 0.12% of the overall R 0 of 2.0224. Increasing the convalescent period from 3 to 6 months, the contribution of R 0,C (0.0051) to the overall R 0 rises to 0.25%. The equation for R 0,C (see Methods) illustrates that doubling the per sex act transmission probability has the same impact as doubling the convalescent period. It is important to note that the relative contribution of sexual transmission to the overall reproductive number rises as the effective reproductive number drops during the epidemic due to the effects of control interventions and decreasing density of susceptible hosts (see S2 Fig).

Effect of sexual transmission parameters on epidemic dynamics
The two key unknown parameters of sexual transmission are the per sex act transmission probability, η, and the rate at which convalescent survivors fully recover, α. Both parameters were found to have very small effects on the peak number of infected or exposed patients ( Fig  3A). The total number of recovered individuals is reached more slowly the longer the convalescent period (Fig 2B), which is not an effect caused by η (Fig 3b). While the convalescent periods (1/ α = [3-9 months]) and the values of η (η = [0.0005-0.002]) we explored create very few extra cases (Figs 2B and 3B), sensitivity analyses revealed that a higher per sex act transmission probability, η, a higher proportion of sexually active convalescent individuals, p, or a higher frequency of sex acts, q, have larger impacts on the total number of cases than would proportional increases in the convalescent period (see S3A Fig). Sensitivity analyses also revealed that these sexual transmission parameters could produce a small delay in the epidemic peak, more so than would changes in the convalescent period (see S3B Fig).
The number of sexual transmission events expected from the baseline scenario (η = 0.001 and 1/ α = 3 months) is 31.2, the majority of which will occur around the peak of the epidemic (Figs 2C and 3D) and thus likely go undetected. Doubling either η or 1/ α results in nearly equal increases in the incidence and cumulative number of sexual transmission events (Figs 2C, 2D, 3C and 3D), with either leading to roughly double the number of sexually transmitted cases over the course of the whole epidemic (> 60 cases). It should be noted that the total number of cases increases more than by simply the number of sexual transmission events, because each sexual transmission event results in a new potential cluster of non-sexual transmission. The day of last case is affected more by the convalescent period than the per sex act transmission probability (represented by vertical lines in Figs 2A and 3A), a result confirmed by the sensitivity analysis (see S3A Fig). The tail of the epidemic will depend on a small number of events that are likely to be affected by stochastic processes, thus we used Monte Carlo simulations to explore this behaviour.

Impact of sexual transmission on the epidemic tail
We performed stochastic simulations of the EVD transmission model to investigate the epidemic dynamics when the number of new cases becomes small, i.e, during the tail of the epidemic. Comparing model simulations while assuming a convalescent period of 3 months to those without sexual transmission confirmed the deterministic results that sexual transmission from convalescent survivors does not lead to a significant increase in the cumulative number of  (Fig 4A, 4B and 4C). This conservative period of potential sexual transmission, which has recently been shown to extend well beyond 3 months in at least 65% of patients [1], lengthened the average date on which the last active case could be detected by nearly three months (non-STI: 548 +/-4.0 days; STI:  Fig 4D, 4E, 4G and 4H). The width of the tail (s.d. = 151 days) was such that 23.4% of the 529 simulated epidemics that accumulated at least 50 cases still experienced symptomatic individuals 730 days (two years) after the start of the epidemic (Fig 4H). Strikingly, when the convalescent period was extended from 3 months to 6 months, the projected length of the epidemic increased to a mean of 1088 days (+/-15.5), with 84.0% of the 538 sustained epidemics taking over two years to end (Fig 4F and 4I). However, the average number of new cases produced remained small (11,869 +/-663 cases; W = 482790, p = 0.18). Importantly, there is greater variance in the tail of the epidemic when sexual transmission is considered, and this uncertainty grows with the length of the convalescent period (Fig 4G, 4H and 4I). To understand the differential impacts of the convalescent period (1/ α) and the per sex act transmission probability (η) on the epidemic tail, the effect of two-fold reductions in η on the average duration of the epidemic under both 3 and 6 month convalescent periods were compared (S1 Table). Cutting the per sex act transmission probability in half led to statistically significant reductions in the length of the epidemic, but this was 3-fold less effective than equivalent changes in the convalescent period (S4 Fig and S1 Table), corroborating the results of the deterministic sensitivity analyses (above). We also note that reducing η did not greatly reduce the enormous variance observed with the longer convalescent period.

Discussion
Our study shows that the length of the convalescent period will determine whether or not sexual transmission of Ebola virus from recovering patients will have a profound effect on the length of time it will take for the epidemic to completely fade. For Sierra Leone, we found that an average convalescent period of 3 months, and a per sex act transmission probability of 0.1%, could extend the EVD epidemic in Sierra Leone by an average of 83 days (95% CI: 68-98 days). Such a scenario would be consistent with the occurrence of a small number of sexual transmission events during the end-phase of the epidemic. However, assuming an average convalescent period of 6 months led to simulated epidemics whose tails were much more variable, and much longer, despite a lack of significant increase in the total number of cases. So far, the reported cases of sexual transmission of EVD remain rare [1,2]. Hence, the per sex act transmission probability of Ebola virus from male convalescent survivors is unlikely to be higher than 0.1%, and might well be below this value. Our sensitivity analysis indicates that the duration of the EVD epidemic is heavily influenced by the period during which convalescent men can transmit sexually, calling for a better understanding of the persistence and duration of infectivity of Ebola virus RNA in convalescent patients.
We extended an accepted modeling framework that has been widely used to describe the epidemic trajectories of EVD outbreaks and epidemics [12][13][14]. To our knowledge, this is the first study using mathematical modeling to assess the potential impact of sexual transmission of Ebola virus on the epidemic in West Africa. In addition, given the generality of the model, this is also the first model that investigates the impact of including a secondary transmission route from convalescent individuals. Similar compartmental models have been formalized to account for more realistic infectious periods, including both infectious relapse [33,34] and progression through classes of varying stages of infectivity [35,36]. However, none of these models included a change in transmission mode between the primary and subsequent infectious classes. This model, then, may also have implications for other pathogens with this kind of secondary transmission route (e.g. some adenoviruses [37]; see [38] for examples across mammal species) including other neglected tropical diseases, such as African sleeping sickness [39], other hemorrhagic fevers that display pathologies similar to EVD [9,10], and the most recent emergence of Zika virus [40,41].
In the absence of a better understanding of sexual transmission of EVD, mathematical modeling currently remains the only tool to explore its potential impact on the epidemic trajectory. There are several extensions to this work that future models should consider including in order to make accurate predictions. These include: transmission by other bodily fluids (e.g., vaginal secretions [4,7]), asymptomatic infection [42,43], considering spatial aspects of both social and sexual transmission [44], and heterogeneity in sexual behavior [3,45]. Sexual behaviors, for instance, are often specific to a given culture, and may change drastically in response to such a devastating epidemic that can destroy entire communities and create stigma, disrupting existing social and sexual contact networks [46,47]; the lack of these details particularly in developing countries presents a challenge for parameterizing a more complex model. We assumed the duration of convalescence to be exponentially distributed. Eggo et al. [48] fitted a series of unimodal distributions to the data on Ebola virus RNA detection in semen recently reported by Deen et al. [1] and found that the convalescence period could be best described by a gamma distribution. Furthermore, Deen et al. [1] found that the cycle-threshold values decreased over time, indicating that the Ebola virus load in semen and the viral infectivity might also decrease during the convalescence period. It is also critical to measure the length of time viral particles persisting in the seminal fluid remain infectious. Molecular techniques to detect intermediate (replicative) positive-sense RNA stages of the virus, infection of human cell lines in tissue culture, or tests in animal models are some typical methods. Retrospective studies using phylodynamics could also prove helpful for estimating this type of parameter [25,49]. Uncertainty in the data is not limited to sexual transmission; we fitted our model to weekly incidence of confirmed and probable cases in Sierra Leone, but did not take into account potential underreporting as others have done recently [36]. In addition, the incidence data to which we fit our model will, for the most part, be driven by direct transmission of the virus and thus, to better parameterize and estimate the risk of sexual transmission, we would need data with greater resolution (e.g. knowledge of which cases were caused by sexual transmission events). Another caveat is that EVD is known to exhibit superspreading characteristics [50,51], and superspreading events could lead to explosive regrowth of the epidemic after the occurrence of a new case through sexual transmission [50]. And finally, like other negative-sense single-stranded RNA ((-)ssRNA) viruses [52], the species currently circulating in West Africa has been estimated to have high substitution rates [26,53,54]. This rapid evolution detected throughout the current outbreak zone suggests that within-or between-host adaptation of the virus leading to pro-longed persistence in the seminal fluids is possible. However, without significant attenuation of EVD's high mortality and morbidity virulence, evolution of sexual transmission becoming the primary route of spread is highly unlikely, as the subsequent infections that arise from a sexual transmission event will be caused primarily through transmission during the acute non-sexual transmission phase of the infection.
Awareness of the potential for sexual transmission has led to WHO issuing recommendations that ask convalescent men to abstain from sexual activity as much as possible and to use condoms for up to 6 months after the onset of symptoms [28]. Condom use and social awareness of the risks of sexual transmission during convalescence should reduce the per sex act transmission probability (η) and the frequency of sex acts (q), respectively. Our results suggest that while such interventions should reduce the number of sporadic sexual transmission events, they will not necessarily reduce the overall number of cases nor the length of time during which the public health community must stay vigilant in responding to these sporadic cases because they will not affect the time during which convalescent survivors can shed infectious virus (1/α). This is especially poignant since adherence to these recommendations will never be 100%, particularly after the threat from symptomatic individuals passes. Thus, our results suggest that the current requirement for declaring a region free from EVD (42 days following either death or a second negative RT-PCR test of the blood from the last known patient), may be too short. Sierra Leone was first declared free from EVD on 7 November 2015 [55], but the case of a young woman who died from Ebola in January [56] highlights the need for the 90-day period of enhanced surveillance after the declaration has been made. The WHO report that 10 such "flare-ups", or cases with no apparent link to the original acute symptomatic transmission chains, have been identified throughout the region, and are suspected to have resulted from contact with infectious survivors [57]. However, none of these events have caused a major resurgence of new cases. This is likely primarily due to the continuing vigilance, awareness, and resources provided by public health infrastructures, and is reflected by the basic reproduction number in presence of control interventions, R 1,N < 1 (Table 1). A relaxation of current response and surveillance efforts could allow a rare sexual transmission event to propagate a new epidemic.
As more data about the convalescent survivors of EVD becomes available, this and future mathematical modeling studies will help to better understand the potential epidemiological consequences of sexual transmission on EVD epidemics. Precise estimates of key parameters are important for providing convalescent survivors with sound advice that balances protection of the community with the harm that could come from unnecessary stigmatization [46,58,59].
Supporting Information S1 Fig. Dynamics of Ebola virus disease (EVD) epidemic in Sierra Leone. Model fit to weekly incidence of confirmed and probable cases are shown together with data from the patient database as reported by WHO (circles) [19]. The shaded area corresponds to the 95% prediction interval, assuming that the number of reported cases follows a negative binomial distribution. (PDF) The effects of twofold changes in convalescent period (1/ α) and in per sex act sexual transmission probability (η) on the average duration of Monte Carlo simulated Ebola virus epidemics. Box and whisker plots show the mean and variance for the length of those epidemics (number of days symptomatic cases remained in the population) which sustained 50 or more total cases (out of 1000 simulations). Statistical significance of comparisons were corrected for multiple tests following the Games-Howell method, and all are reported in Supporting Information S1 Table. Asterisks denote p-values of ( Ã ) < 0.05 and ( ÃÃÃ ) <0.0001. (PDF) S1 Table. Sexual transmission parameter value effects on epidemic length. Pairwise effects of two-fold changes in convalescent period (1/ α) and in per sex act sexual transmission probability (η) on the average duration of Monte Carlo simulated Ebola virus epidemics. Out of 1000 simulations for each set of parameters (α and η), the mean length (number of days) of the epidemics having at least 50 total cases (n = number of simulations) is given in the diagonal. Statistical significance of each pairwise comparison is given above the diagonal (t statistic (df) and p-value) and were corrected for multiple tests using the Games-Howell method. The relative reduction in the length of the epidemic (% fewer days) is given below the diagonal. The red and blue values are those referenced in the results of the manuscript. (PDF) S1 Appendix. Calculating R 0 . (PDF) S2 Appendix. R codes and data file for modeling sexual transmission of Ebola virus disease. Sierra Leone incidence data. As reported by the World Health Organization on their Ebola outbreak data web page. Parameter Estimation (SEIR). Implementing deterministic EVD models (SEIR, SEICR). Implementing sensitivity analysis. Implementing Monte Carlo simulations of EVD models. Summary statistics of Monte Carlo simulations. (ZIP)