Effect of sexual transmission on the West Africa Ebola outbreak in 2014: a mathematical modelling study

The outbreak of the Ebola virus has resulted in significant morbidity and mortality in the affected areas, and Ebola virus RNA has been found in the semen of the survivors after 9 months of symptom onset. However, the role that sexual transmission played in the transmission is not very clear. In this paper, we developed a compartmental model for Ebola virus disease (EVD) dynamics, which includes three different infectious routes: contact with the infectious, contact with dead bodies, and transmission by sexual behaviour with convalescent survivors. We fitted the model to daily cumulative cases from the first reported infected case to October 25, 2014 for the epidemic in Sierra Leone, Liberia and Guinea. The basic reproduction numbers in these countries were estimated as 1.6726 (95%CI:1.5922–1.7573), 1.8162 (95%CI:1.7660–1.8329) and 1.4873 (95%CI:1.4770–1.4990), respectively. We calculated the contribution of sexual transmission to the basic reproduction number R0 as 0.1155 (6.9%), 0.0236 (2.8%) and 0.0546 (3.7%) in Sierra Leone, Liberia and Guinea, respectively. Sensitivity analysis shows that the transmission rates caused by contacts with alive patients and sexual activities with convalescent patients have stronger impacts on the R0. These results suggest that isolating the infectious individuals and advising the recovery men to avoid sexual intercourse are efficient ways for the eradication of endemic EVD.

A number of mathematical models have been developed to understand the transmission dynamics of Ebola and other infectious diseases outbreak from various aspects (see for instance 15,16, ). Furthermore, several mathematical studies for the transmission dynamics only focus on the direct transmission or only consider sexual transmission: Lewnard et al. constructed a mathematical model with individuals divided into susceptible, latently infected, infectious, recovered and infectious people who have died from the disease 33 . They fitted the model to the reported cumulative Ebola cases and deaths in Montserrado County, Liberia, and assessed various intervention strategies in controlling the spread of Ebola virus in the country. In another study, Drake et al. carried out a multi-type branching process model and used time-dependent parameters to reflect increasing intensity of intervention efforts, incorporate the rates of hospitalization, exposure of healthcare workers, and secure burial 34 . Similarly, Barbarossa et al. developed a compartmental model including virus transmission in the community, at hospitals, and at funerals, and they derived a final epidemic size by forecasting the total number of cases during the outbreak 27 . In a recent study, Abbate et al. developed an SEIR compartmental model that includes sexual transmission from convalescent survivors 16 . Using the model, they evaluated the potential effect of sexual transmission on Ebola virus epidemiologically. Similarly, Vinson et al. presented a mathematical model that comprised the direct contact and sexual transmission modes 15 . Mathematical analysis showed that reducing sexual transmission could play a significant role in eradication of the disease. Gao et al. 45 presented a mathematical model to investigate the impact of mosquito-borne and sexual transmission on spread and control of ZIKV and used the model to fit the ZIKV data in Brazil, Colombia, and El Salvador.
Inspired by these previous literatures, in this paper, without aiming to make a predictive model but rather to understand how the epidemic was affected by various transmission routes which include contact with infections, contact with dead bodies, and having sex with convalescent survivors, we develop a mathematical model which considers the three main features of the disease transmission and fit to the reported cumulative case data of infected cases in Sierra Leone, Liberia, and Guinea.
The article is organized as follows. In Section 2, we introduce the EVD model and present the expression of the basic reproduction number, and introduce the data sources. The model fitting and contribution of related component to the basic reproduction number R 0 in Sierra Leone, Liberia and Guinea are carried out in Section 3. In Section 4, we give a brief summary and some discussions.

Methods and Materials
Model Formulation. To study the transmission of Ebola in West Africa we developed a EVD model, where the population is divided into six groups: the susceptible individuals S, exposed individuals E, symptomatic and infectious individuals I, deceased infected individuals D, convalescent survivors who maintain active Ebola virus replication R 2 , and fully recovered and immune people R 1 , such that the total population is N = S + E + I + R 1 + R 2 + D. According to the natural history of EVD, some people have an innate immunity to Ebola and the diseases may not be fatal to all 46,47 , so the exposed individuals may go directly to the recovery class R 1 . Because we considered a relatively short time frame (<1 year), we assumed for simplicity that the total population is constant and hence we omit population births and deaths when modeling the outbreak. Additionally, we assumed that people who recover from the Ebola infection develop antibodies that last for at least 10 years, so recovered individuals will not move to compartment S in a short time. A flow of individuals from one class to another is depicted in Fig. 1. The transmission process of Ebola virus is then described by the following a system of six ordinary differential equations:  We assumed susceptible individuals S could be infected via three paths: person-to-person transmission β 1 , sex with convalescent survivors β 2 and deceased individuals β 3 . 1/γ and 1/η represent the average durations of incubation and infection, respectively. δ represents the proportion of exposed E progress to the infectious I, θ is the case fatality rate. Deceased individuals can be buried directly during funerals at rate ρ. The average duration for convalescent patients to recover is assumed as 1/p. Basic reproduction number. The basic reproductive number R 0 is the average number of secondary cases caused by an infectious individual during his/her entire infectious period. To calculate R 0 , we used the next generation matrix method 48,49 (2) The reproduction number is given by ρ(FV −1 ), and In order to interpret R 0 in terms of the dynamics of Ebola spread, the basic reproduction number for the system breaks down to three components: the new infection induced by contacts with infectious individuals R I , sexual activities with convalescent survivors R S , and contacts with dead bodies R D .
Data and Parameter Estimation. We obtained data of the daily cumulative cases which is publicly available online from the Centers for Disease Control and Prevention 50 , which allows the user to obtain information about the daily number of Ebola cases (confirmed, probable and deaths) in Guinea, Sierra Leone, Liberia. We used the start date for each country as when the first infectious case was reported. According to the available data provided by the CDC, the first cases reported were on May 27th, 2014 for Sierra Leone, March 27th, 2014 for Liberia and March 25th, 2014 for Guinea. As of October 25, 2014, the cumulative cases and deaths in each country are shown in Fig. 2. We use the cumulative cases up to October 25, 2014 to fit model (1).
The values of parameters for model (1) are summarized in Table 1. To parameterize our model, we adopted the epidemiological parameters from literatures. For the total population size, N, we used a total population of 7,079,162 for Sierra Leone, 4,390,737 for Liberia, and 11,805,509 for Guinea, as provided by the World Bank 51 . In view of the reported cases from the EVD outbreak in West Africa in 2014, we chose the mean incubation period 1/γ as 9.31 days and the mean infectious duration 1/η as 7.41 days. The average duration from death to burial 1/ρ is assumed as 2 days. The disease-related death rations, θ, according to previous studies 52 , were the case fatality rates 0.690 for Sierra Leone, 0.723 for Liberia, and 0.707 for Guinea. We assume that Sierra Leone, Liberia and Guinea share the same values of γ, η, and ρ [see Table 1]. Although some parameters are not available, we can still make some reasonable assumptions. We used the value 87.35 for 1/p based on published values from the literature, and for all three countries.
We model cumulative cases as a Poisson-distributed random variable because the Poisson distribution describes the number of observed events in an interval of time. We calibrate the model by sampling from the posterior distribution of parameter vector θ|y = {β 1 , β 2 , β 3 , δ}|y, where vector y is derived from and Y(t) denotes the reported cumulative cases. We conduct sampling via Markov Chain Monte Carlo using a Metropolis-Hastings acceptance rule. The posterior density is The prior density f Θ (θ) is the joint probability of four univariate priors. We consider that β 1 , β 2 , β 3 , and δ are distributed according to (0, 1). The program was implemented in R version 3.3.1. We sampled from 30,000 MCMC iterations and discarded the first 10,000 samples as a burn-in period. On the basis of these 20,000 samples, the point estimates and 95% confidence intervals for the transmission coefficients were calculated, and the three kinds contribution to R 0 can be calculated based on existing parameter values. The results are shown in Table 1. In order to check the influence of different prior distributions for estimated parameters, we choose (0, 1), .
.  Fig. 3). This shows that the prior distributions have little influence on parameter β 1 . Similar results can be seen for parameters β 2 , β 3 , and δ.

Result
Model Fitting. Model fitting was performed using MCMC method. In numerical experiments, the daily number of Ebola cases for each country is used to fit the cumulative number of infected cases. The model fits the cumulative reported cases in Sierra Leone, Liberia, and Guinea well generally (see Fig. 4). Figure 4 demonstrates that the fitting result of Sierra Leone is the best of these three countries, as Liberia and Guinea reported data with fluctuations, especially Guinea. In any case, the cumulative numbers of disease in October continues to increase in three countries. It should be noted that our model does not a consider any interventions, since, intervention efforts were escalated significantly after October 1, 2014 20 .   Fig. 5). The simulated results indicate that Liberia is the most seriously affected area in three countries at that time. In contrast, our R 0 agrees with previous works summarized in Table 2  In order to understand the contribution of related components to R 0 more clearly, we calculated the basic reproduction numbers of either contact with alive infectious transmission or sexual transmission or contact with dead bodies transmission (see Table 1). More precisely, in Sierra Leone, the basic reproduction numbers of infectious component R I accounted for 0.9331 (55.9%) of R 0 , the sexual transmission component R S accounted for 0.1155(6.9%) and the death component R D accounted for 0.6218 (37.2%). For Liberia, we obtained 1.1053 (62.1%) for R I , 0.0236 (2.8%) for R S and 0.6250 (35.1%) for R D . As for Guinea, above three components account for 1.0890 (72.9%), 0.0546 (3.7%), and 0.3493 (23.4%), respectively (see Fig. 6).

Contribution of related component to
We found that the infectious component basic reproduction number was around 1 in the three countries, suggesting the most important factor for the spread of the epidemic is the virus transmission via contact with alive infectious individuals. In addition, we found that the epidemic was spread through contacting dead bodies in the funeral. Thus, we conjecture that the spread of Ebola can be effectively controlled by isolation and proper handling of dead bodies. The value of R S is less than 10% in three countries, and not enough to cause an outbreak of disease. However, sexual transmission can prolong the outbreak of Ebola.

Sensitivity Analysis and Disease Control.
It is well known that the basic reproduction number (R 0 ) is a very important parameter in the infectious disease model, which determines whether the could spread. In our model, R 0 is determined by the parameters of γ, θ, η, ρ, p, β 1 , β 2 , β 3 and δ. In order to identify the impacts of thses parameters on Ebola transmission and prevalence, we used the Latin hypercube sampling method and partial rank correlation coefficient (PRCC) 39 . Using model (1), 2000 samples are randomly generated by assuming a uniform distribution for each parameter based on values from Table 1.
We choose parameters of interest as the input variables, and the value of R 0 as the output variable. The PRCC values of eight parameters and corresponding p-values are computed. Given the common ways of controlling measures in three countries, and to avoid repeated work, here, we only presant Sierra Leone's sensitivity analysis. The results of Sierra Leone are illustrated in Table 3 and shown in Fig. 7. The larger the absolute value of PRCC, the stronger the correlation between the input parameters and R 0 . From Table 3, it is evident that η, ρ and p have negative impact upon R 0 , while γ,θ, β 1 , β 2 , β 3 and δ have positive impact. Furthermore, Fig. 7 shows that contact with alive infectious transmission rate β 1 , infectious class's removal rate η, sexual transmission rate β 2 and contact with dead bodies transmission rate β 3 are the most sensitive parameters on R 0 .
Through the result of sensitivity analysis, we focus on quantifying the effects of β 1 , η, β 2 , and β 3 on the basic reproduction number. In the numerical simulations, we choose γ = 1/9.31, ρ = 1/2, p = 1/87.35, θ = 0.69, δ = 0.9511, and vary the four controllable parameters from 0.001 to 0.3 for β 1 , from 0.001 to 0.01 for β 2 , from 0.1 to 0.7 for η and from 0.0001 to 0.1 for β 3 . Fixing values of η and β 3 , β 2 and β 3 , β 1 and η simultaneously, we plot the contour figure of R 0 in terms of β 1 and β 2 , β 1 and η, β 2 and β 3 are displayed in Fig. 8, to explain the specific value of controllable parameters when the basic reproduction number can be brought below 1.
From Fig. 8(a,b), we find that it is possible to reduce the basic reproduction number to be below 1 with a certain range of the parameters, which means reducing the rate of contacting with infectious individuals may be an effective method to control the outbreak of Ebola. In Fig. 8(c), we show the influences of β 2 and β 3 on the basic reproduction number. We can see that even if we can implement proper handling of dead bodies to reduce β 2 , the disease is hard to be controled with a high sexual transmission outbreaks.

Conclusion and Discussion
Ebola has a high mortality rate and it has caused in many West Africa countries, especially in Sierra Leone, Liberia and Guinea in 2014, which led to many deaths. In this paper, we developed a mathematical model to explore the transmission dynamics of Ebola virus. As a comprehensive consideration of the transmission of the virus, the model includes the effect of contact with infections, contact with dead bodies, and sexual contact with convalescent survivors. Furthermore, we qualitatively analyzed the impact of sexual transmission on the prevention and control of the Ebola. By fitting the model to the cumulative numbers of infected cases (Fig. 4), we obtained some estimates of the parameters (Table 1). Based on reasonable parameter values and the existing literature, we calculated the basic reproduction number R 0 and 95% confidence interval for Sierra Leone, Liberia and Guinea as 1.67 (95%CI:1.59-1.76), 1.82  Our results seem to be in accordance with previous work carried out in this field (Table 2).
To be more precise, we consider the contribution of different settings for transmission of Ebola to R 0 with no effective control measures. In our results, the ranking proportional values of R I , R S , and R D accounting for the basic reproduction number was the same in the three countries, the biggest proportional value among the three countries is R I , then R D , and the smallest is R S . The results imply that R I and R D transmission routes play more crucial roles in the Ebola transmission in Sierra Leone, Liberia and Guinea. Similar components of R 0 were identified by Barbarossa 27 , but in our study, the value of R I is the largest, which may due to the fact that R I includes both the infections generated in the community R C and that in the hospital R H in their model. However, contribution of sexual transmission to overall R 0 is 0.12% in 16 .
Sensitivity analysis is vital to identify key parameters and find effective control strategies for combatting the spread of the disease. The result of sensitivity analysis indicates that contact with alive infectious transmission   rate β 1 , infectious class's removal rate η, sexual transmission rate β 2 and contact with dead bodies transmission rate β 3 are the most sensitive parameters to R 0 than any other parameters. The simulation results show β 1 imacts the most of the value of R 0 , which indicates that protective among the community and medical workers should be strengthened, and measures to isolate infected people in turn reduce the rate of infection in the population. Remarkably, the sexual transmission component R S accounted for 6.9%, 2.8% and 3.7% in Sierra Leone, Liberia and Guinea respectively. However, the contribution of sexual transmission to overall R 0 is 0.12% in 16 . Our results show that the potential risk of sexual transmission will cause greater difficulties in terms of controlling the spread of EVD.   The knowledge of the potential for sexual transmission has led to the WHO to make recommendations that require recovered men to avoid sexual intercourse and use condoms as much as possible 53 . Moreover, the public health community must be vigilant in responding to these cases. Most of the population in the affected area are still susceptible, and as long as the virus is not eradicated, the risk of an outbreak could still be high.