A metapopulation model for the 2018 Ebola outbreak in Equateur province in the Democratic Republic of the Congo

Ebola is a viral haemorrhagic fever with high mortality that has caused a number of severe outbreaks in Central and West Africa, the largest of which was in 2014-16 and resulted in 11,325 deaths. Although most previous outbreaks have been relatively small in comparison, the result of managing outbreaks places huge strains on already limited resources. Mathematical models matched to early case reporting data can be used to identify outbreaks that are at high risk of spreading. Here we consider the Ebola outbreak in Equateur Province in the Democratic Republic of the Congo, which was declared on 8 May 2018. We use a simple stochastic metapopulation model to capture the dynamics in the three affected health zones (Bikoro, Iboko and Wangata) and the capital city Kinshasa. We are able to rapidly simulative a large number of realisations and use a likelihood-free method to determine parameters by matching between reported and simulated cases. This likelihood-free matching has a number of advantages over more traditional likelihood-based methods as it is less sensitive to errors in the data and is a natural extension to the prediction framework. Using data from 8 to 24 May 2018 we are able to capture the exponential increases in the number of cases in three locations (Bikoro, Iboko and Wangata) and the probability of transmission to capital city Kinshasa, although our estimated basic reproductive ratio of 4.02 is higher than for previous outbreaks. Using data until 28 June 2018 we are able to infer the decrease in transmission due to public-health intervention measures, such that the reproductive ratio is predicted to drop below one around 16 May 2018 leading to decreasing numbers of cases. We believe this method of fitting models to data offers a generic approach that can deliver rapid results in real time during a range of future outbreaks. Author summary Ebola is an infectious disease that, if left untreated, is often fatal. In addition, the consequence of managing Ebola outbreaks places huge strains on already limited healthcare resources in affected countries. Mathematical models can be useful in identifying high-risk outbreaks, and deciding where to allocate resources. However, many existing models of Ebola cannot capture the spatial spread of the outbreak, require highly detailed data, or are too complicated to be used in real-time during an outbreak. In this paper we describe a framework that can capture the spatial spread of the 2018 Ebola outbreak in Equateur province in the Democratic Republic of the Congo, and is simple enough that it can be re-evaluated as new data become available. We used this framework to understand how transmission changed as a result of public-health interventions, and the risk of transmission to Kinshasa, which was a major public-health concern. This framework is highly flexible and can easily be adapted to new geographic regions, to use different data sources, or to answer pressing public-health questions. This research provides important insights into this Ebola outbreak, and has the potential to generate substantial impact on a range of pubic-health decision-making.


Introduction
shows a schematic representation of the model, while Table 1 summarises 204 the possible events and transition rates. Table 1. A summary of possible transitions and their rates for the compartmental model used to describe Ebola dynamics within a population of size N . There are ten free parameters in the model and five additional parameters that are calculated from the first ten. Model parameters are: β I community transmission rate; β H hospital transmission rate; β F funeral transmission rate; α −1 mean duration of incubation period; γ −1 IH mean time from symptom onset to hospitalisation; γ −1 IF mean time from symptom onset to death; γ −1 IR mean time from symptom onset to end of infectious period for survivors; γ −1 HF mean time from hospitalisation to death; γ −1 HR mean time from hospitalisation to end of infectious period for survivors; γ −1 F R mean time from death to traditional burial. θ 1 is calculated from θ such that θ% of infectious cases are hospitalised; δ 1 and δ 2 are calculated from δ such that the case fatality ratio is δ. Calculations are given in Appendix A.

206
To describe the spatial dynamics of Ebola infection we use a metapopulation model, 207 whereby the total population is split into K interacting sub-populations of sizes N i , i = 208 1, . . . , K. We define σ ij to be the proportion of epidemiologically relevant contacts that 209 individuals from population i have with individuals in population j, which we will sim-210 ply refer to as the coupling from population i to population j. We naturally have that 211 K j=1 σ ij = 1, and so the within population coupling (which we expect to be close to one) which susceptible individuals become infected, can then be written in terms of the coupling 214 parameters as: As such the transmission is assumed to be due to the movement of healthy susceptible 216 individuals visiting infected locations, such that the risk to individuals in population i is 217 related to the coupling terms σ ij .

218
The coupling is usually parameterised using mobility data [ We define v ij to be the the number of visits from population i to j. According to the 224 generalised gravity model this is proportional to N a i N b j /d c ij , where d ij is the distance be- introduce an additional scaling parameter A ∈ [0, 1] such that min i σ ii ≥ 0. Combining 232 each of these elements we define the coupling σ ij , i = j, to be: where a, b, c and A are parameters to be inferred from the epidemiological dynamics.  (Table 3). 253 We perform parameter inference using approximate Bayesian computation, a flexible 254 likelihood free approach that is straightforward to compute. A traditional likelihood based 255 methodology would require us to infer infection times and subsequent disease progression 256 (together with associated times) for each case, and is likely to be affected by the temporal 257 aggregation of observed cases as a result of detection and testing. Therefore for each reali-258 sation we calculate an error between the realised (C sim ) and observed (C obs ) cumulative 259 confirmed cases from 11 May 2018 until some later date T 1 . We define confirmed cases in 260 our model to be individuals who have moved from the infected class to either the hospi-261 talised or funeral class. We define the error as a weighted root mean square error, summed 262 over the four sub-populations: .
( The denominator in this expression is motivated by considering a Poisson distribution.

264
In a Poisson distribution the variance is equal to the mean, therefore we would normalise  such that: where T C defines the time at which control effects begin, and (1 − δ) determines the 292 reduction in transmission from all infectious classes.

293
These two new parameter values are inferred as detailed above, using the associated 294 errors. We start with the first 50 days of the outbreak (to 24 May 2018) and determine the 295 best 1,000 parameter sets (for all 7 parameters) for 10,000,000 parameter updates; these 296 are then fed into simulations for the first 55 days (to 29 May 2018) and we again iterate 297 the fitting procedure for a further 10,000,000 steps. This incremental fitting processes 298 mimics what would occur in real-time as more data becomes available on a daily basis and 299 parameters are refined. The re-fitting procedure is easily completed by an over-night run. 300 We perform incremental re-fitting until we reach 85 days (29 June 2018) by which point 301 we have had three weeks without any additional confirmed cases in the region.        There are ten free parameters in the model and five additional parameters that are calcu-710 lated from the first ten. The model parameters that can be derived from the free parameters 711 are: γ −1 HF , the mean time from hospitalisation to death; γ −1 HR , the mean time from hospital-712 isation to end of infectious period for survivors; θ 1 , the rate of transition from infectious to 713 hospitalised, which is calculated from θ such that θ% of infectious cases are hospitalised; 714 δ 1 and δ 2 , the effective case fatality ratios in the infected and hospitalised class, respec-715 tively, both which are calculated from δ such that the case fatality ratio is δ. These five 716 parameters are calculated as follows: