Modeling the transmission dynamics of Ebola virus disease in Liberia

Ebola virus disease (EVD) has erupted many times in some zones since it was first found in 1976. The 2014 EVD outbreak in West Africa is the largest ever, which has caused a large number of deaths and the most serious country is Liberia during the outbreak period. Based on the data released by World Health Organization and the actual transmission situations, we investigate the impact of different transmission routes on the EVD outbreak in Liberia and estimate the basic reproduction number R0 = 2.012 in the absence of effective control measures. Through sensitivity and uncertainty analysis, we reveal that the transmission coefficients of suspected and probable cases have stronger correlations on the basic reproduction number. Furthermore, we study the influence of control measures (isolation and safe burial measures) on EVD outbreak. It is found that if combined control measures are taken, the basic reproduction number will be less than one and thus EVD in Liberia may be well contained. The obtained results may provide new guidance to prevent and control the spread of disease.


Estimation of Basic Reproduction Number R 0 of EVD in Liberia.
In this part, we use the least-square fitting method to estimate parameter values in order to minimize the sum of squared errors between the actual data and solution of the equation (1). The accumulated number of infected cases with time N I (t) can be given by the following equation with N I (t) = I P(c) (t) + I S(c) (t), where I P(c) (t) and I S(c) (t) denote the cumulants of the I P (suspect cases) and I S (probable cases), respectively: The actual number of the ( ) N t I (the accumulated incidence) can be found in 24 . The estimation process is as follows: we construct a function = ∑ ( ) − ( ) t n I I 1 2 and find the suitable parameters value to make f to be least, where n is the number of actual data (In our model, n = 50). Biological meanings of parameters can be found in Table 1. By applying the real data in 24 and the equation (1), we can estimate the value of p = 0.1, β 1 = 0.1102, β 2 = 0.12 and r i = 0.0667. The fitting result for the accumulated incidence is given in Fig. 1.
The time range is from June 29, 2014 to October 7, 2014. As seen from Fig. 1, we can find that the slope of the fitting curve gradually increases which implies that epidemic is still aggravating. If control measures are not taken effectively, outbreak of EVD in Liberia is inevitable.   Sensitivity and Uncertainty Analysis of Basic Reproduction Number R 0 . It is well known that the basic reproduction number (R 0 ) determines whether the epidemic will persist or not. If R 0 > 1, the disease will be epidemic; otherwise, it will eventually vanish. As a result, it is meaningful to discuss the sensitivity and uncertainty analysis of R 0 . In our model, several crucial parameters (β 1 , β 2 , β H , β F , θ and η) determine the value of R 0 . For the sensitivity and uncertainty analysis, we adopt Latin hypercube sampling (LHS) to study the influence of parameters on R 0 . We randomly choose 1000 samples and the six parameters follow a normal distribution. On the basis of the 1000 samples, we can perform an analysis of R 0 by computing variable partial rank correlation coefficient (PRCC). Larger the absolute value of PRCC, denotes stronger correlation between the chosen parameters and R 0 . The values of PRCC are showed in Table 2, and it is obvious that the absolute value of PRCC for parameter β 2 is the largest, which indicates that β 2 is the most influential in Case-fatality ratio from infectious to death (δ 1 ) 0. 8 11 Case-fatality ratio from hospitalized to death (δ 2 ) 0. 4 10,11 Transmission coefficient with the suspected in the community (β 1 ) 0.

[estimated]
Transmission coefficient with the probable cases in the community (β 2 ) 0.12 [estimated] Transmission coefficient at the hospital (β H ) 0.062 11 Transmission coefficient during the funerals (β F ) 0.489 11 The initial number of susceptible sheep (S′ ) 3441700 11 The initial number of exposed individuals (E′ ) 20 24 The initial number of suspected individuals ( ′ ) I S 29 24 The initial number of probable individuals ( ′ ) I P 18 24 The initial number of hospitalized individuals (H′ ) 20 24 The initial number of cases dead but not yet buried (F′ ) 11 24 The initial number of R (R′ ) 23 24  determining the value of R 0 . Although the values of β H , β F , θ, η are less important for PRCC in contrast with β 2 , these parameters also have some impacts on the value of R 0 . Furthermore, we can find that the value of PRCC between β 1 or β 2 and R 0 is larger than those in β H or β F , indicating that contact transmission with infected cases (I S , I P ) in the community posses larger influence than contact transmission with the hospitalized cases or cases dead but not yet buried. Additionally, we find that there are positive correlations between β (β 1 , β 2 , β H , β F ) and R 0 , which suggests that the bigger the transmission coefficients, the larger value of R 0 . The negative correction between θ and R 0 indicates that increasing the patients' hospitalization rate is an effective method for controlling EVD in Liberia under the current situation.
Peak arrival time of EVD in Liberia. In this part, we show the relationship between the parameters (typically, we just show the results of parameters β 1 and β 2 ) and the peak arrival time of the EVD as well as the maximum value of I P (t).
In Fig. 2, it is obvious that β 1 (or β 2 ) and the peak arrival time of EVD is close to linear increment relationship. This figure illustrates that the peak arrival time arrives sooner when the β 1 (or β 2 ) is larger. In other words, we can take measures to decrease β 1 (or β 2 ) to cause the delay of peak arrival time and thus less people will be infected by EVD. In Fig. 3, we can conclude that the final scale of EVD outbreaks in Liberia will get larger when the β 1 (or β 2 ) increases which is consistent with the actual situation. If effective control measures are not taken, the epidemic will become more serious and there will be much more new infected cases in the future.
Control measures of EVD in Liberia. Now we consider the effects of existing control measures on EVD in Liberia. There are mainly two kinds of prevention strategies: isolation of the infected individuals (I S , I P , H) and safe burial of cases dead but not yet buried (F). For comparison, we need to examine the efficiency of isolation of I S , I P , H and safe burial of F. Firstly, we consider the model with isolation and safe burial measures as follows:   24 . Estimated basis reproduction number is 2.012, which is consistent with the real cases in Liberia. This figure indicates that EVD will spread as an endemic in the absence of the control measures.  In our results, we estimated that EVD in Liberia may arrival its peak value after 370 days from the day when the first infected case was detected. where d 1 , d 2 , d 3 are the isolation rates of suspected cases, probable cases and hospitalized cases in the community, and d 4 is the safe burial rate of cases dead but not yet buried in funerals. Assume that d 1 , d 2 , d 3 and d 4 are all non-negative and the remaining parameters are shown in Table 1.
In the following part, we consider the sensitivity analysis of isolation and safe burial on R c and its expression has the following form: where ( = , , …, ) Fig. 4(A). Although d 1 = 1, we have that R c = 1.2344 > 1, which means that only taking control measures on suspected cases (I S ) is not enough to control EVD. In the case with d 1 = 0, d 3 = 0, d 4 = 0, we find that R c may be less than one (see Fig. 4(B)). That is to say when isolation measure on probable cases is sufficient to take, EVD in Liberia will ultimately disappear.
In Fig. 4(C), we check the effect of isolation measure in the hospital on EVD spreading. It can be seen from this figure that R c = 1.5685 > 1 even if d 3 = 1, which means that only taking control measure in hospital is not enough to eliminate the EVD in Liberia. We also show the influence of safe burial measure on eradication of disease in Fig. 4(D). One concludes that only taking safe burial measure in funerals cannot induce the disappearance of disease due to that R c = 1.521 > 1 with d 4 = 1.
From the above analysis, we find that taking isolation measure on probable cases may be an effective method to control the prevalence of EVD in Liberia. However, the actual situation in Liberia is that it is nearly impossible to isolate all the probable cases. In that case, it is necessary to combine different control measures together. In Fig. 5, we show the influences of combined control measures on R c . We can see that it is possible to cause basic reproduction number R c to be less than one for a certain range in parameters space except for Fig. 5(F). Compared with single control measure, combined control strategies are more useful for EVD control in Liberia.

Discussion
EVD is a lethal disease with a high mortality rate which has caused many deaths. Although the Liberian government has taken some control measures, the number of infectives of EVD increases continuously. Consequently, we consider the contribution of different settings for transmission of EVD in the estimation of R 0 with no effective control measures. Based on the parameters estimation and literature 11 , we obtain R 0 = 2.012, where the term of R 0 concerning the transmission during funeral is about R 0F = 0.736 and the contact transmission with suspected cases (I S ) in the community is about = . R 0 725 which implies that these two transmission routes play more important roles in EVD transmission in Liberia.
Our model for EVD transmission in Liberia is based on ordinary differential equations (ODEs) which can be analyzed by mathematical analysis and make prediction on the trend of EVD, and thus it has essential differences from the research on EVD by agent-based model or branching process model 25,26 . Moreover, different from the previous work [25][26][27][28][29] , we divide the infected individuals into two classes: suspected case (I S ) and probable case (I P ), which is more in line with the actual situations of EVD in Liberia. We estimate basic reproduction number R 0 = 2.012 of EVD, which confirms the results obtained by Castillo-Chavez et al. that the basic reproduction number of Ebola in Liberia is in the range of [1.9, 2.4] 30 .
In our results, we estimate that EVD in Liberia may outbreak after 370 days since the time when the first case was confirmed. Moreover, the final size of suspected infectives may achieve 22, 000 cases. That is to say, EVD in Liberia is not well controlled in the current situation. In this case, we need to find effective methods to curb the spread of EVD. Based on sensitivity analysis, we demonstrate that only taking single measure can not control the spread of EVD well, which is consistent with the conclusions posed by Khan et al. 27 . Furthermore, we find that taking several control strategies together may be effective for EVD control in Liberia, which highly matches the findings by Merler et al. that decrease of incidence at country and county level is attributable to the increasing availability of EVD treatment 25  The prevalence of the disease is not optimistic currently and the natural reservoir is still not identified 2,6-8 , and thus EVD may outbreak somewhere outside of Africa in the future. In the further study, we will try to define nature reservoir by using mathematical models 8,31 . At the same time, the good news is that there are some therapies for EVD 32 , which needs to be well checked on the effectiveness of EVD control. What is more, contact tracing is an effective method in controlling EVD 33,34 , and we will do some efforts to examine the influences of human behaviors in EVD control in details.

Method
Data. Time series of reported cases were collected from the World Health Organization and the Ministry of Health of Liberia. The data contains the I S , I P (The definition can be founded in 35 ), the deaths and the confirmed cases. Though the data does not contain the patient level information, they provide the best available data on the outbreak of EVD in Liberia. More details on data is available in ref. 24. Mathematical model. In order to make our model more reasonable, we must do some assumptions (transmission rules can be seen in Table 3): (a) Nearly all the population in Liberia was considered initially as the susceptible; (b) Assume no effective prevention measures before November, 7th, 2014; (c) If a suspected case goes to see a doctor the suspected case will be considered as a probable case; (d) Some misdiagnosed cases will return to be susceptible; (e) We only consider the spread in human beings.
As a result, we arrive at the following equations to model the transmission dynamics of EVD in Liberia without effective control measures (Transmission diagram can be seen from Fig. 6):    where S is number of susceptible individuals; E, number of exposed individuals; I S , number of the suspected individuals in the community; I P , number of probable individuals in the community; H, number of the hospitalized cases; F, number of cases who are dead but not yet buried; R, number of individuals removed from the chain of transmission 10 . Parameter β 1 is transmission coefficient with the suspected cases in the community; β 2 , transmission coefficient with the probable cases in the community; β H , transmission coefficient with the hospitalized cases; β F , transmission coefficient during the funerals; θ, proportion of suspected cases hospitalized; p, misdiagnosed proportion in the suspected cases; q, the proportion of suspected cases except the misdiagnosed; η, the proportion of exposed cases who enter the I P compartment; Case-fatality ratio from probable cases to death is δ 1 ; Case-fatality ratio from hospitalized to death is δ 2 ; , the mean duration of progression from exposed cases to the probable cases; m 1 2 , the mean duration of progression from exposed cases to suspected cases; the mean duration of progression from probable cases to hospitalized cases is The basic reproduction number. Through direct calculation, we obtain that the model (3) has a disease-free equilibrium = ( , , , , , , ) E S R 0 0 0 0 0 0 0 0 , and the formula for R 0 is the spectral radius of the next generation matrix. Following the method described in van den Driessche 23 , we consider the infected compartments satisfied by model (3), which has the following form: (2) ( , ) → ( − ,   In order to express simply, we do some marks: We have the following matrix:    Therefore, the basic reproduction number , R 0F are partial reproduction numbers induced by the suspected cases, probable cases, hospitalized cases and dead cases but not yet buried, respectively. Next we calculate the expressing of the reproduction number (R c ) of model (2). For simplicity, we do some marks: We have the following matrix: where  represents the rate of appearance of new infection and  denotes the rate of transfer of individuals. Furthermore, we can obtain: where: , where R cI S , R cI P , R cH , R cF are partial reproduction number induced by the suspected cases in the community, probable cases in the community, hospitalized cases and cases dead but not yet buried, respectively.