INTRODUCTION
Middle East respiratory syndrome (MERS), is a respiratory illness caused by a novel coronavirus (CoV) [Reference de Groot1]. The disease was reported for the first time in Saudi Arabia in June 2012 and spread to several countries in Africa, Asia, Americas and Europe [2, Reference Zaki3]. The capability of human-to-human transmission has been observed in at least four hospital settings [4–Reference Drosten7]. Significantly, MERS-CoV shares certain similarities with the severe acute respiratory syndrome (SARS)-CoV that produced a global epidemic with more than 8000 human cases in 2002–2003 [Reference Marra8, Reference Peiris9]. First, a number of patients infected with both viruses developed an acute respiratory disease that in some cases resulted in death [Reference Drosten7, Reference Peiris9]. In this sense, MERS-CoV appears to be highly pathogenic with an estimated case-fatality rate of around 50%, although this might be an overestimation as many infected patients may not have sought hospital assistance [2]. Second, both MERS-CoV and SARS-CoV, belong to the genus Betacoronavirus and are closely related to coronaviruses isolated from bats [Reference Annan10–Reference Li12]. This strongly suggests that MERS-CoV and SARS-CoV may have been transmitted from bats to humans through intermediate species (e.g. camels for MERS-CoV and civet cats for SARS-CoV). Third, since they are new emerging viruses, there are no effective vaccines or antiviral treatments. The SARS epidemic is a clear example of how a networked health system can respond to a new threat to human health.
One of the best-characterized outbreaks during the SARS epidemic was in Hong Kong in 2003 where there were 1755 confirmed cases with 299 deaths (WHO; http://www.who.int/csr/sars/country/en/index.html). The outbreak began in mid-February caused by an infected person who travelled from Guangdong to Hong Kong [13]. An important fact in the generation of a model of the outbreak is that the Hong Kong health authorities quickly implemented contagion control procedures [Reference Riley14]. In general, two interventions were introduced to prevent the spread of SARS-CoV. The first was implementation of quarantine measures to isolate healthy people who had been in contact with infected people and therefore potentially in contact with the virus; isolating those that could be infected and asymptomatic during the incubation period and isolating and treating patients who had developed the disease. The other intervention was the application of protective measures by healthy people who were in contact with infected people to avoid becoming infected, such as respiratory protection for healthcare workers and daily disinfection of the environment of affected rooms [Reference Riley14]. These control interventions were implemented progressively in the Hong Kong Special Administrative Region from mid-March to late April [Reference Chau and Yip15]. The application of these procedures allowed the rapid control of the outbreak in the subsequent months. In this sense, it is estimated that the epidemic in Hong Kong ended in late June.
System dynamics has been proved to be a powerful instrument for analysing social, economic, ecological and biological systems [Reference Sterman16]. In addition, disease epidemiology has been studied using this approach, whereby system dynamics offers the practical application of concepts by computerized models that allow the systematic test of different scenarios and alternative policies [Reference Homer and Hirsch17–Reference Koopman, Jacquez and Chick19]. In this work, we performed a modelling of the Hong Kong SARS-CoV outbreak using system dynamics. The developed model contains five states, four flows, eight auxiliary variables and six parameters that interact through five differential and 12 algebraic equations. The parameters of the model were optimized following an iterative process of simulation to obtain a model that largely fits the data available to the epidemic. Moreover, the credibility of our model and its parameters are supported by both univariate and multivariate sensitivity analyses. The model reproduces how the implementation of control measures was effective in preventing the spread of infection to the rest of the population. Basically, these measures result in a sustained reduction in the frequency of contacts. At present, the application of similar measures for infection containment can help to prevent the spread of new emerging epidemics, such as the outbreak caused by MERS-CoV.
METHODS
Data sources and collection
Data on the total population of Hong Kong in 2003 was obtained from statistics of the Census and Statistic Department, The Government of the Hong Kong Special Administrative Region (http://www.censtatd.gov.hk/home/). Data on cumulative cases, deaths and recoveries during the SARS epidemics in Hong Kong was obtained from Global Alert and Response databases of the WHO (http://www.who.int/csr/sars/country/en/index.html).
Modelling
The model was developed following the four-step sequence proposed by system dynamics methodology [Reference Sterman16]. First, the real data from the Hong Kong SARS epidemics (Fig. 1) together with other evidences and our professional experience were used to create a mental modelling of the reality of the outbreak. Second, the model structure that is able to explain the evolution of the epidemics was represented as a Forrester diagram (Fig. 2). Third, the outbreak was mathematically modelled as a continuous dynamic process represented by a set of differential and algebraic equations (Tables 1–3). Finally, the model was optimized to fit the real data from Figure 1.
Fig. 1. Real data from the Hong Kong SARS outbreak. (a) Daily reported new cases, deaths, and recovered. (b) Cumulative cases, dead, and recovered.
Fig. 2. Graphical description of the model using Vensim software. Stock variables are represented inside boxes; flow variables are designed as arrows with a valve; the rest of the variables are auxiliary (inside circles) or parameters (in bold); single lines with arrow are used as connectors to specify that the destination variable is affected by the variable of origin.
Table 1. Differential equations for the five model stocks
* Hong Kong's population in 2003 was about 6·5 million people who were susceptible to infection with the SARS virus.
† It is considered that the Hong Kong outbreak originated from an infected person from the province of Guangdong [13].
Table 2. Algebraic equations for the four model flows
Table 3. Algebraic equations for the eight model auxiliary variables
A dynamic compartmental model provides a framework for the study of transport between different compartments of a system, i.e. well known epidemiological compartmental models [Reference Anderson and May20]. To explain how the protective measures taken by the government of Hong Kong allowed the rapid control of the epidemic we consider a dynamic model with five compartments (states) and four transitions (physical flows) between them. This model is based on two assumptions. First, the individuals are classified into five subgroups (susceptible, latent, infected, recovered, dead). Although, the last subgroup is not strictly needed, it is used to keep an account of those dead. Second, every day there is a different number of people: who are infected without symptoms of the disease (incidence); who develop signs and symptoms of the illness (sick per day); who recover from the disease (daily recovered), and who die (daily deaths) as consequence of the outbreak.
The model structure built with Vensim software (Ventana Systems Inc., USA) (Fig. 2) contains the five states and the four flows mentioned above and also eight auxiliary variables (inside circles) and six parameters (in bold). These elements are linked by the physical flows (double line with arrow) and by the information transmissions (single line with arrow) according to the mathematical model represented by the set of differential and algebraic equations (Tables 1–3).
The five differential equations in Table 1 establish the mass balance (inflows minus outflows) in the compartments, and as such, they describe the changes in the number of people (stocks) in the five subgroups. The four algebraic equations of Table 2 express that the physical transitions depend directly from the stock in the compartment of origin and indirectly from other stocks and parameters through the corresponding auxiliary variables. These equations involve three stocks (susceptible, latent, infected), three auxiliary variables (prevalence, contagion rate, recovery rate) and three parameters (incubation period, case fatality, disease duration). For instance, the incidence flow is proportional to the number of susceptible people, the prevalence and the contagion rate. The contagion rate or transmission coefficient theoretically depends on the number of contacts per unit time and the probability of effective contact, i.e. the probability that a contact between an infectious source and susceptible host results in a successful transfer whereby the susceptible host becomes infected.
The first three algebraic equations of Table 3 express the auxiliary variables (prevalence, contagion rate, recovery rate) used in the equations of Table 2. In turn, these depend on other auxiliary variables, stocks and parameters. For instance, prevalence is defined as the ratio between number of infected people and the total population [Table 3, equation (1)]. The recovery rate depends on the illness duration and case fatality [Table 3, equation (2)]. Equation (3) in Table 3 expresses the contagion rate which is directly dependent on the auxiliary variable ‘frequency of contacts’ and the parameter ‘infectivity’.
One of the key variables in the model is frequency of contacts because it tries to reproduce the control measures carried out by the Hong Kong Government to control the SARS outbreak. We assumed that: (1) the control measures gave, as consequence, a marked reduction of daily contacts in Hong Kong; (2) the control measures were based on the cumulative attack rate, measured as the ratio between cumulative cases and total population [Table 3, equation (6)]. Therefore, we decided to use the frequency of contacts as the daily contacts modulated by the repression Hill function (Fig. 3), in accordance with the equation (7) of Table 3. Despite the complexity level of biological systems several cases have been modelled using the Hill function, in order to simulate repressor activities of enzymatic reactions and the regulation mechanisms of several transcription factors [Reference Omony21]. In our model this repressor function allows the establishment of the relationship between frequency of contacts and the cumulative attack rate. Note that the cumulative attack rate is used as an active repressor, so half-maximal repression occurs when the cumulative attack rate is equal to the threshold, and almost total repression occurs when this cumulative attack rate is double the threshold.
Fig. 3. Normalized representation of repression Hill function. The solid line corresponds to a value of n = 3 [as used in equation (7) of Table 3]. The dotted line corresponds to the threshold repression function (when n tends to infinity). The representation is normalized to the threshold of cumulative attack rate on the abscissa axis and the value of the discontinuity (daily contacts) on the ordinate axis.
Equation (8) in Table 3 is used to report on the basic reproductive number (R 0), which is defined as the expected number of secondary infectious cases generated by an average infectious case in an entirely susceptible population [Reference Anderson and May20]. In our model this number is calculated as the ratio between the contagion rate and the recovery rate. If R 0 < 1, then the infected individuals in the total population fail to replace themselves, and the disease does not spread. However, if R 0 > 1, the number of cases generally increases over time and the disease spreads.
Parameterization
One of the most challenging issues in system dynamics modelling is to establish the value of the model parameters. The parameters can be estimated taking advantage of the known information available in the literature and can be optimized by means of an iterative process of simulation. Using this approach, we set values for the six parameters of our model which are summarized in Table 4.
Table 4. Values for the six model parameters
Infectivity expresses the ability of the pathogen to penetrate, survive and multiply in the host and it is measured through the secondary attack rate which is defined as the probability that infection occurs in susceptible persons within a reasonable incubation period following known contact with an infectious person or an infectious source. Epidemiological studies in Singapore showed that 80% of individuals infected with SARS-CoV did not cause secondary infections, suggesting low infectivity [22]. After the optimization process, the final value set for infectivity was 0·022, which is in agreement with a previous work examining the probability of transmission for SARS-CoV [Reference Pitzer, Leung and Lipsitch23].
The incubation period of the disease, during which time the individual is asymptomatic appears to be between 2 and 10 days with an average value of 5·3 days [Reference McBryde24]. Similarly, although the duration of the disease is variable depending on the case, the average of the most severe and the mildest cases is 26 days [Reference McBryde24]. During this time period the virus can be transmitted to other people.
Several studies have shown that the case fatality of SARS-CoV was also variable in the 2003 outbreak depending on different factors. It has been estimated that the mean case fatality for SARS in Hong Kong was around 0·17 [Reference Hirose25].
Social contact patterns have been shown to be highly relevant on the transmission dynamics of respiratory infections such as measles, rubella and influenza [Reference Fine and Clarkson26–Reference Anderson and May28]. Quantifying this parameter is critical for estimating the impact of such infections, for designing and targeting preventive interventions, and for modelling their impact. In our model the daily contacts parameter, indicating the mean number of contacts of a person in 1 day, was set to 16·8 after the optimization process, which is in agreement with previous studies quantifying these social mixing patterns [Reference Edmunds, O'Callaghan and Nokes29, Reference Mossong30].
The parameter ‘threshold of cumulative attack rate’ is critical in our model since it allows setting the value of the cumulative attack rate at which control measures were established by the Hong Kong authorities. This parameter was finally set at 7·8 × 10−5. Therefore, according to the Hill function of Figure 3, for any number of the cumulative attack rate below the 78 cases by million, the frequency of contacts will be greater than 8·4 (half of the daily contacts), and in the opposite case the frequency of contacts will be markedly reduced.
As explained before, our dynamic model was subjected to successive rounds of simulation and optimization in order to fine-tune the parameters of the model. In all these simulations we assumed that the Hong Kong epidemic originated from a traveller from Guangdong in China [13], and at the beginning of the outbreak the entire population of Hong Kong was susceptible to SARS-CoV infection. These assumptions are represented by the initial values in the five stocks (Table 1). Moreover, based on the reported data, the simulation period was set to 146 days with a time step of 1 day.
Sensitivity analysis
Parameters of system dynamics models are subject to uncertainty. Sensitivity analyses were conducted to provide insight into how uncertainty in the parameters affects the model outputs and which parameters tend to drive these variations. In this task it is essential to define the probabilistic distribution patterns of the model parameters, which are shown in Table 5, together with the references that support these patterns.
Table 5. Probabilistic distribution patterns of the model parameters
* The distribution values for normal distributions are mean and standard deviation, whereas for gamma distributions they are order/shift/stretch.
The most influential parameters were estimated by univariate analyses, in which changes in the model output were studied after disturbance in each parameter value independently. In complex models, univariate sensitivity analysis can be insufficient for a comprehensive study of the model. Simultaneous fluctuations in the value of more than one parameter may create an unexpected output change due to nonlinear relationships in different model components. Thus, to test the influence of simultaneous changes in the model parameters, the univariate analyses were followed by Monte-Carlo multivariate sensitivity analysis, in which the values of the six parameters were altered at the same time.
Software
Modelling process, simulations and sensitivity analyses were performed using Vensim DSS software v. 5.7a (Ventana Systems).
RESULTS
Figure 4 shows a graphic comparison between the simulation results using the parameters of Table 4 and the real data of the Hong Kong outbreak. We compared only the variables of the model for which records were found in the public databases, which are the same six variables shown in Figure 1.
Fig. 4. Simulation of the model. For simulation, the time period was set to 146 days with a time step of 1 day. The graphs represent the output of the simulation (solid lines) and the data reported for the SARS outbreak (dotted lines). (a) Sick per day, (b) cumulative cases, (c) recovered, (d) dead, (e) daily recovered, and (f) daily deaths.
The simulation output for the variable ‘sick per day’ fit the data reported by the Hong Kong authorities (Fig. 4a ), suggesting that the model was able to reproduce the epidemic curve. We observed that the number of new cases per day obtained in the simulation grew during the first 46 days. From this time, the number of new infections gradually fell to values <1 at later stages of the outbreak.
As a consequence, the auxiliary variable that stores the cumulative SARS cases showed a characteristic sigmoidal growth (Fig. 4b ), consistent with the real data. We observed that the number of SARS cases grew from one at the beginning of the epidemic to around 1800 at the final stage of the outbreak, similar to the 1755 cases reported by the authorities. Moreover, we observed a high fit between the model predictions and the real data.
As expected, the number of deaths and recovered people in the simulation also grew with a sigmoidal shape to reach values similar to those found in the public databases (Fig. 4c, d ). However, we observed a partial fit of these stock variables to the data reported during the outbreak. These slight mismatches were also observed in the flow variables for recovered and daily deaths (Fig. 4e, f ) that may be due to delays in reporting of cases by the authorities.
In short, looking the simulation results of Figure 3 we can conclude that our model is able to reproduce largely the most important indicators of the SARS epidemic that occurred in Hong Kong in 2003.
Focusing on the evolution of the basic reproduction number (Fig. 5), we note that during the early stages of the epidemic, R 0 >1, which is consistent with the observed disease spread. Moreover, after day 30, R 0 starts to drop to values <1, probably due to the implementation of containment measures by the Hong Kong authorities after the issuance of the first global alert against SARS on 12 March 2003. These results are consistent with a previous report showing the basic reproductive numbers for different SARS epidemic curves, which supports the notion that our model is able to largely replicate the disease outbreak in Hong Kong [Reference Wallinga and Teunis31].
Fig. 5. Basic reproductive number obtained from the simulation. R 0 was calculated as the ratio between the contagion rate and the recovery rate.
The results of the univariate sensitivity analysis are shown in Figure 6. We focused our attention on the epidemic wave although the analysis is possible in other variables, as shown in the multivariate analysis of Figure 7. Variations in case fatality, threshold of cumulative attack rate and disease duration induced little changes in the epidemic curve (Fig. 6–c ), while more extensive alterations in the epidemic wave were observed after changes in infectivity, daily contacts and incubation period (Fig. 6d–f ). Variations in the case-fatality parameter do not alter the output of the variable sick per day (Fig. 6a ), although other variables from the model such as recovered and dead are highly impacted (data not shown). Small changes in the shape and the maximum of the curve are observed after modification of the parameters ‘threshold of cumulative attack rate’ and ‘disease duration’ (Fig. 6b, c ). By contrast, changes in infectivity, daily contacts and to a lesser extent in incubation period significantly alter both the position and height of the maximum of the epidemic curve (Fig. 6d–f ).
Fig. 6. Sensitivity of the model to changes in parameters. Univariate analysis with 1000 simulations was performed for each parameter: (a) case fatality, (b) threshold of cumulative attack rate, (c) disease duration, (d) infectivity, (e) daily contacts, and (f) incubation period. The parameters were varied considering ranges and distributions shown in Table 5. The solid black line represents the simulation output and the grey area represents the 95% confidence bounds.
Fig. 7. Multivariate sensitivity analysis of the model. A Monte-Carlo sensitivity analysis with 1000 simulations was performed and the results for (a) sick per day, (b) infected, (c) recovered, and (d) dead are shown. The values of the six parameters were altered at the same time following the distributions described in Table 5. The solid black line represents the simulation output and the grey area represents the 95% confidence bounds.
The results of the multivariate sensitivity analysis are shown in Figure 7, we analysed the output of four variables: sick per day, infected, recovered, and dead (Fig. 7a–d , respectively). Variations in the model parameters clearly change the shape of the epidemic curve, altering both the position and the height of the maximum of this variable (Fig. 7a ). Similarly, the output of the variables infected and recovered is highly impacted by the changes in parameters carried out in the multivariate sensitivity analysis (Fig. 7b, c ). The alterations of variable outputs are clearly exemplified by the variable ‘dead’ (Fig. 7d ). Certain changes in the parameters of the model can explain an increase in the total number deaths during the epidemic rising from about 300 to about 700. The observed variations in model output when the value of the parameters is changed support the idea that this model might be able to simulate different scenarios and epidemic conditions.
DISCUSSION
System dynamics modelling has been successfully applied to study complex public health issues such as the design of optimal policies in healthcare [Reference Royston32], the impact of public health intervention in different situations [Reference Morilla33, Reference Hirsch and Miller34], and to study disease epidemiology [Reference Luginbuhl35–Reference Dangerfield, Fang and Roberts37]. In the latter, system dynamics technology has become a powerful tool to understand and predict the impact of infectious diseases. Epidemiological models can help health authorities to make recommendations regarding intervention to fight the spread of directly transmissible pathogens, especially when empirical data is limited. In this sense, mathematical models have been previously used to advise health policies against diseases such as pandemic influenza [Reference Longini38, Reference Halder, Kelso and Milne39] and SARS [Reference Wallinga and Teunis31, Reference Cauchemez40, Reference Pitman41].
Several models studying SARS transmission and interventions have been published. These are detailed hybrid stochastic and compartmental models that successfully explain the behaviour of the epidemic [Reference Riley14, Reference Wallinga and Teunis31, Reference Lipsitch42]. Here, trying to follow Ockham's razor principle, we have built a simpler deterministic model, which is also able to reproduce the behaviour of the epidemic based on its natural history and the intervention measures taken in Hong Kong. The use of a less complicated model could be helpful in understanding the disease epidemic and also facilitating its reuse under other conditions.
The epidemiological models depend on the consistency of the chosen parameters. Therefore, the accurate quantification of these parameters is critical to estimate the path of a disease, to predict the impact of possible interventions, and to inform planning and decision making. Here, we have combined reported information from the SARS epidemic with an iterative optimization process to set the final values for the model parameters. Under these conditions, the model output fits the epidemic curve observed in the Hong Kong SARS-CoV outbreak (Fig. 4).
Of the factors that influence the dynamics of infectious diseases, the person-to-person contact pattern has been shown to be essential in disease spread [Reference Grassly and Fraser43]. Our model takes this essential factor into account through the auxiliary variable ‘frequency of contacts’, which is dependent on the auxiliary variable ‘cumulative attack rate’ and the parameters ‘daily contacts’ and ‘threshold of cumulative attack rate’ (Fig. 2). A previous work showed that mixing patterns and contact characteristics were remarkably similar across the different European countries analysed in that study even though the average number of contacts recorded differed. Interestingly, the authors suggests that the results may well be applicable to other countries with similar social structures, and that the initial epidemic phase of an emerging infection in susceptible populations, such as SARS was in 2003, is likely to be very similar [Reference Mossong30]. During the SARS outbreak, health authorities, hospitals, and the overall population progressively implemented quarantine and protection policies to prevent the transmission of the disease [Reference Chau and Yip15]. With this in mind, we made the assumption that the intervention of the health authorities caused a decrease in the frequency of contacts, which in turn led a decrease in the rate of contagion. To mimic this event in our model, the number of daily contacts is regulated by means of a repression Hill function, which has been used to simulate repressor activities in complex biological systems [Reference Omony21].
The basic reproductive number is a key epidemiological variable that characterizes the potential of a disease to spread. Several works have estimated that prior to the first global alert the basic reproductive number for SARS was >1, correlating with an exponential increase in the number of cases. However, the implementation of effective control measures, such as quarantine, isolation, and strict hygiene practice in hospitals led a sudden decrease in R 0 [Reference Riley14, Reference Wallinga and Teunis31]. The fact that, in our model, R 0 drops to values <1 around the date of the global alert reinforces the idea that a decrease in the frequency of contacts is able to effectively simulate the effects of the control measures established in the first stages of the epidemic. It is important to note that estimations of R 0 in previous publications are considerably lower (around 2–4) than ours, which is almost 12 at the beginning of the outbreak [Reference Riley14, Reference Wallinga and Teunis31, Reference Lipsitch42]. Nevertheless, the estimated value of R 0 also differs in these works and the credible intervals surrounding these deterministic estimations were wide, reaching superior values of almost 8. This high variability can be explained in part by the superspreading events that occurred in SARS epidemics. Superspreading is an unusual situation, in which a single individual directly infects a large number of other people that has a large influence on the early course of the epidemic [Reference Lipsitch42]. Interestingly, the fact that our model does not explicitly account these events could partially explain the very high estimated value of R 0 at the beginning of the outbreak.
The reliability of the model parameters is supported by univariate and multivariate sensitivity analyses (Figs 6 and 7). Furthermore, sensitivity analysis is a powerful tool to analyse the influence of certain decision making in the evolution of the epidemic. Taken together, these findings strongly suggest that our model, together with other system dynamics models can be used by epidemiologists to investigate the likely consequences of future re-emergences of SARS-CoV based on analysis of the previous known epidemics. In addition, by adapting the key parameters of these models or with a little change in the model structures, they can be used to face emerging outbreaks of infectious diseases, such as the recent MERS-CoV epidemic. In this regard, there are several similarities and differences which should be taken into account when using this model. Both SARS-CoV and MERS-CoV may cause severe respiratory failure, extrapulmonary features such as diarrhoea and also mild or asymptomatic cases. In contrast with SARS, MERS has lower human-to-human transmission potential, affects predominantly older people with more comorbid illness and has a higher case-fatality rate [Reference Hui, Memish and Zumla44]. These factors would affect the output of the model variables (e.g. epidemic curve, cumulative cases, number of dead, etc.).
ACKNOWLEDGEMENTS
This research was supported by an intramural contract from Universidad Autónoma de Madrid awarded to E. Álvarez. The Institutional Grant awarded to the Centro de Biología Molecular ‘Severo Ochoa’ by the Fundación Ramón Areces is also acknowledged.
DECLARATION OF INTEREST
None.
