The 2017 plague outbreak in Madagascar: data descriptions and epidemic modelling

From August to November 2017, Madagascar endured an outbreak of plague. A total of 2417 cases of plague were confirmed, causing a death toll of 209. Public health intervention efforts were introduced and successfully stopped the epidemic at the end of November. The plague, however, is endemic in the region and occurs annually, posing the risk of future outbreaks. To understand on the plague transmission, we collected real-time data from official reports, described the outbreak’s characteristics, and estimated transmission parameters using statistical and mathematical models. The pneumonic plague epidemic curve exhibited multiple peaks, coinciding with sporadic introductions of new bubonic cases. Optimal climate conditions for rat flea to flourish were observed during the epidemic. Estimate of the plague basic reproduction number during the large wave of the epidemic was high, ranging from 5–7 depending on model assumptions. The incubation and infection period for bubonic and pneumonic plague were 4.3 and 3.4 days and 3.8 and 2.9, respectively. Parameter estimation suggested that even with a small fraction of the population exposed to infected rat fleas (1/10000) and a small probability of transition from a bubonic case to a secondary pneumonic case (3%), the high human-to-human transmission rate can still generate a large outbreak. Controlling rodent and fleas can prevent new index cases, but managing human-to-human transmission is key to prevent large-scale outbreaks.


Introduction
One of the deadliest natural disasters in human history was reported as the Black Death -attributed to the bacterium Yersinia pestis -killing about 50 to 200 million people in the 14th century [1]. Although plague was naturally widespread in ancient times, plague outbreaks occurred following the deliberate 5 use and propagation of this disease, serving as a bioweapon [2]. Nowadays, plague epidemics continue to pose a threat to humans, reporting continuous annual occurrence in five countries: Madagascar, Tanzania, Vietnam, China, and the USA [1,3]. This lethal bacterium can derive in several forms of plague maintaining its existence in a cycle involving rodents and their fleas [4]. While 10 sanitation and public health surveillance have greatly reduced the likelihood of a plague pandemic, isolated plague outbreaks are lethal threats to humankind.
The disease manifests in different clinical forms of plague: bubonic, pneumonic, and septicemic [1]. Human infection is primary driven by bubonic plague, as a result of being bitten by infected fleas. Additionally, direct contamination 15 with infective material can be an alternative transmission route [1]. Patients with bubonic plague can develop sudden onset of fever, headache, chills, tender and painful lymph nodes [5]. While plague can be successfully treated with antibiotics, if untreated, the bacteria can disseminate from the lymph nodes into the bloodstream causing secondary septicemic plague. In addition to the 20 symptoms presented in the bubonic plague, patients with septicemic plague undergo abdominal pain and possibly bleeding into the skin and other organs, at the same time skin and other tissues may turn black and die, especially on fingers, toes, and the nose [4]. However, the most fulminant form of the disease is driven by pneumonic plague that is the only form of plague that can spread 25 from person to person by infectious droplets. The incubation period of primary pneumonic plague is shorter than in the other forms of the disease with an average of 4 days [6]. The disease progresses rapidly and is nearly always fatal without prompt antibiotic treatments [7].
Patients with pneumonic plague often do not transmit the disease to anyone, 30 but in the right conditions one can infect many people and cause an outbreak [7,8]. This was observed in the last epidemic in Madagascar when an 31-yearold man travelled from the central highlands to the eastern city of Toamasina via the capital city, Antananarivo. He died in transit and dozens of his contacts subsequently became ill [9]. Since then cases of suspected plague have been 35 reported from many areas of Madagascar. On 13 September 2017, the Madagascar Ministry of Public Health notified WHO of an outbreak of pneumonic plague [5]. A total of 2417 cases of plague has been confirmed of which 77% were pneumonic plague, causing until now a death toll of 209 [9]. The Government of Madagascar with the supports of WHO and partners had focused their efforts 40 on strengthening the identification and treatment of patients and their contacts, increasing control of rodents and fleas, and practising safe and dignified burials [10]. These measures have been preventing new cases and deaths, however, the disease is endemic and occurs annually in the region and elsewhere [1,3], posing the risk of future outbreaks. This paper provides descriptive and numerical 45 analyses of the plague outbreak to facilitate further studies in evaluating the spread of the plague as well as targets for disease control and prevention.

Materials and Methods
Outbreak Data -Cumulative Cases. Data were manually inputted from separate reports of WHO [9], including the cumulative total numbers of clinical cases 50 (confirmed, probable, and suspected). The data can be found at the following link  systemsmedicine/plague2017/Cumulative Outbreak Data -by Disease Forms. Data were digitized from the figure reported from WHO [9], including the incidences classified by the three forms of the plague disease: pneumonic, bubonic, and septicemic. The data can be 55 found at the following link  systemsmedicine/plague2017/Classification and the digitized figure can also be found at the same repository.
Temperature and Precipitation Data. Data were requested from the National Centers for Environmental Information (Order #1133340 Custom GHCN-Daily CSV). The data can be found at  systemsmedicine/plague2017/Climate. 60 Descriptive analyses. With the aim of facilitating modelling works, we described dynamics and patterns of variables that have previously shown to be relevant to plague outbreaks, including temperature and precipitation [11].
Statistical estimate of the reproduction number. We estimated the reproduction number (R0) of Yersina pestitis using data of pneumonic cases (excluding the bubonic cases) during the second (large) wave of the epidemic, i.e., visually defined from 22/09/17 onwards. The serial interval of plague was assumed gamma distributed with shape and scale parameters are 5.4 and 0.9, respectively [12]. We reported the R0 estimates using several methods for comparison purposes: exponential growth (EG) [13], maximum likelihood (ML) [14], sequential 70 bayesian estimation (SB) [15], binomial assumption [16], and sub-exponential growth model [17,18]. Implementations in R of the last two methods are available at  systemsmedicine/plague2017/ are as follows: where S, E, I describe Susceptible, Exposed, and Infectious and the subscripts b and p denote bubonic and pneumonic form, respectively. The model assumes 75 that in a population of size N , only a small part S b = pN is exposed to infected rat fleas. The infected bubonic cases become infectious with a proportion ε progressing to pneumonic stage [8]. The overall transmission rates of the two modes flea-to-human and human-to-human are denoted by α and β, respectively.
It has been shown that climatic conditions favour the survival and reproduc-80 tion of fleas [19,20,21]. Fleas density fluctuates by season temperatures, and changes in their host (rats) also depend on the season which relates to the breeding patterns of rats [22]. To this end, we assumed that the density of infected rat fleas can be described as a sinusoidal function f irf = A + B sin(2π/12t) + C cos(2π/12t) [23] following the temperature fluctuation; here the average tem-85 perature of Madagascar in the period 1960-2008 [11] was used. As such, the overall flea-to-human transmission parameter α implicitly incorporates a scaling factor from temperature to rat fleas density.
We assumed interventions would reduce flea-to-human and human-to-human transmission rates via rodent and flea control, via active case finding and the 90 identification and prophylaxis of contacts. The interventions effect is assumed imperfect and has a logistic form as We fitted the PTM to the daily data of pneumonic and bubonic cases during the large wave of the epidemic curve which was visually defined from 22/09/17 100 onwards. Model parameters were estimated using the global optimisation algorithm Differential Evolution [24]. We derived also R0 using the next generation matrix [25]. Simulations and estimations were written in R using packages base [26], deSolve [27], and R0 [28] and in Python. Stochastic simulations were performed using a tau-leaping algorithm with a fixed time-step of ca. 15 minutes; 105 code and data are publicly available at  systemsmedicine/plague2017/.

Descriptive analyses
During August, bubonic cases appeared sporadically with almost no records of pneumonic form (Fig. 2). An increase in number of pneumonic cases was not 110 necessarily preceded by an increase in the number of bubonic cases (Fig. 2). It seemed to be the epidemic curves include several waves of incidence overlapping each other. Figure 3 shows that plague incidences emerged everyday in the weeks though in some weeks fewer cases were reported during the weekends. No distinctive 115 time lag was observed between the appearances of bubonic and pneumonic cases (Fig. 3). The incidences were negligible during the period when the Famadihana tradition was presumably practised. Precipitation measure exhibited no pattern before or during the outbreak but generally showed a dry climatic condition.
Average temperature appeared increasing and reached a higher level (above 23 120 degree Celsius) around the same time as the outbreak. August. The parameters suggest that only a small fraction exposed to infected rat fleas is enough to generate the observed epidemic. Table 1 shows that several statistical estimation methods gave a similar value for plague's basic reproduction number which is approximately 7. We also estimated R0 using the PTM. Evaluating the Jacobian of the PTM at the disease-free equilibrium yielded a threshold parameter R * such that R * ≥ 1 results in λ 6

Data
The outbreak period Fitted model with the estimate obtained from the next-generation-matrix (NGM) [25], for which we find the leading eigenvalue at t = 0 to be R0 = β δp = R * (0). In other words, at the beginning (S b + S p = N ) when there are no interventions  days and thus generate β × 1 δp = 6.54 secondary cases.

Mathematical models of infectious diseases have played a central role in under-
standing epidemics, providing an effective way of assessing disease transmission 135 as well as evaluating disease control and prevention strategies [33]. Mathematical modeling has proposed new vaccination strategies against influenza infection [34]; supported public health strategies for containing emerging influenza pandemics [35,36] and for the use of antiretroviral treatment for HIV-infected From modelling aspects, plague outbreaks can be more challenging because: (1) there is a continuous input of flea-to-human transmission (Fig. 2); this im-150 plies the observed epidemic curve can be a mixed of multiple waves of infected cases generated from different index cases. Thus, epidemic evaluations could risk to over-or under-estimating the consequences, e.g. the reproduction number or the end time of epidemics; (2) there are a known seasonal pattern of the plague epidemic [11] for which a direct measure of the rat flea population 155 does not exist; (3) The flea-to-human transmission as well as the transition from bubonic cases to pneumonic cases appeared stochastically driven (Fig. 2) and could be highly affected by interventions.
Here, the epidemic curves showed plague incidences appear sporadically during August (Fig. 2). But then a large increase in pneumonic cases was observed, 160 preceding by only a few bubonic cases. Estimates of the plague reproduction number showed that a high estimate is needed to capture this fast growing phase of the epidemic; the estimate doubled the previous estimates ranging 2.8-3.5 [12]. It can be speculated that superspearders were likely to exist in order to generate the larger number of cases in a short time as observed in the pneumonic 165 epidemic curve. This was a common pattern as described in plague outbreaks in Madagascar 2017 [9] and in the US 1919US , 1924US , and 1980. Considering the potential mixed of epidemic waves and outbreak locations in the used data, our estimate of R0 could be overestimated. However, the smaller estimates of R0 of this outbreak ranging 1.1-1.4 [40,41] could be due to the inclusion of Nonetheless, approaches adjusting for the propagated outbreak data are needed to further understand its effect on parameter estimates. For example, 175 how do the data of a non-synchronized and combined epidemic curve affect the reproduction number? Implementations and calculations using current methods raised also an issue that has been discussed elsewhere [31], that is estimating the reproduction number using serial interval could yield very high value during the first days of the epidemic growing phase as the denominator of the estimator is 180 extremely small during this period. Overestimate of R0 has also been observed in populations with heterogeneous contact pattern population [42].
It appeared that sporadic inputs of bubonic cases brought new index cases to the human-to-human transmission network constantly with an estimate of transition to secondary pneumonic is of 3% (Fig. 2). It followed that the pneu-185 monic epidemic curve exhibited multiple peaks and waves. As the transmission rate can be high, this observation suggest that vector control are key to prevent potential next waves of the epidemic. This might be practical as the fraction of population exposed to infected rat fleas was low ( human-to-human transmission but also provide a better chance for new bubonic cases to be treated early.
Experiments have shown that an optimal climate for rat fleas to flourish is a dry climate with temperatures of 20-25°C [21]. These conditions were observed during the outbreak: a generally dry weather with the optimal temperature 200 coincided with the period of high epidemic activity. As further shown in Fig. 4, these conditions would typically remain the same until next May, see further in Kreppel et al. [11]. This again stresses the role of vector control in preventing the next waves. It is worth noting that the global changes in climate conditions could pose a risk of irregular changes in vector biting rate and reproduction [43],

Supporting Information
The eigenvalues of the Jacobian of the PTM evaluated at the disease-free equilibrium (S b , S p , E b , E p , I b , I p )/N = (0,s p , 0, 0, 0, 0) have the form