key: cord-0698391-eb8ysg2b
authors: Wang, Wendi; Ruan, Shigui
title: Simulating the SARS outbreak in Beijing with limited data
date: 2004-04-07
journal: J Theor Biol
DOI: 10.1016/j.jtbi.2003.11.014
sha: 5e85a5bfca539ebee61755f48923854277fb2034
doc_id: 698391
cord_uid: eb8ysg2b

We propose a mathematical model to simulate the SARS outbreak in Beijing. The model consists of six subpopulations, namely susceptible, exposed, quarantined, suspect, probable and removed, as China started to report SARS cases as suspect and probable separately from April 27 and cases transferred from suspect class to probable class from May 2. By simplifying the model to a two-compartment suspect-probable model and a single-compartment probable model and using limited data, we are able to simulate the SARS outbreak in Beijing. We estimate that the reproduction number varies from 1.0698 to 3.2524 and obtain certain important epidemiological parameters.

Atypical pneumonia, called fei dian in Chinese, was first reported in Guangdong Province, China in November of 2002 (WHO, 2003a . As of February 11, 2003, there were 305 atypical pneumonia cases (105 were health workers) and five related deaths in Guandong Province. 3 Despite the urgent announcements and warnings issued by the Health Department of the Guandong Province in late January, 4 very quickly the emerging disease spread to some other regions in Mainland China as well as Hong Kong, Singapore, Vietnam, Canada, etc. In March of 2003, the World Health Organization, for the first time in its history, issued a globally warning about the disease (WHO, 2003a) and officially defined the disease as severe acute respiratory syndrome (SARS) (WHO, 2003d) . In late June of 2003, the disease was under control globally, but it had spread to 32 countries and regions causing about 800 deaths and more than 8000 infections (WHO, 2003e) . 5 The first SARS cases in Beijing were reported in early March, 2003. 6 Compared with the outbreaks in Hong Kong, Singapore, Vietnam and Toronto, in mid April the WHO Beijing team estimated that Beijing may have had as many as 200 cases of SARS, rather than the 37 officially reported cases (WHO, 2003c) and recommended that Beijing improve its reporting system. On April 20, the Health Minister of China and the Mayor of Beijing City were surprisingly removed from their posts. 7 On April 21, China's State Council Information Office reported that Beijing had 339 confirmed (probable) and 402 suspected SARS cases. 8 Since then China has reported public information about SARS cases on a daily and more accurate basis (Ministry of Health of P. R. China, 2003) .

In late April and May, Beijing became the hardest hit city in the world with the most probable SARS cases. After a vigorous campaign to combat the spread of SARS, such as closing schools, theaters, Internet cafes and other recreational centers, creating the first SARSonly hospital, Beijing Xiaotangshan Hospital, 9 and quarantining more than 30,000 exposed individuals, 10 the outbreak of SARS in Beijing was under control in June. On June 24, Beijing was removed from the WHO SARS list. By then there had been 2521 probable cases and 191 deaths (WHO, 2003e; Ministry of Health of P. R. China, 2003) .

It has been recognized that a coronavirus is responsible for SARS (Drosten et al., 2003; Ksiazek et al., 2003; Peiris et al., 2003; WHO, 2003b) . Although its transmission modes are not completely known, droplets and contacts appear to be the predominant transmission modes. The cases at universities and hospitals in Beijing suggest that airborne transmission is also possible. Since there are no valid medicines or vaccine for SARS, measures to control the spread of SARS have taken two major forms: isolation of symptomatic individuals and quarantine and close observation of asymptomatic individuals (WHO, 2003f) .

Mathematical models have recently been used to examine the transmission dynamics and model the control of SARS in Hong Kong (Chowell et al., 2003; Riley et al., 2003) , Singapore (Chowell et al., 2003; Lipsitch et al., 2003) and Toronto (Chowell et al., 2003) . The objective of this paper is to study the transmission dynamics of the spread of SARS in Beijing at the population level. As China started to report SARS cases as suspect and probable separately from April 27 (Ministry of Health of P. R. China, 2003) following the WHO definitions (WHO, 2003d), we classify the population of Beijing into six classes, namely susceptible ðSðtÞÞ; exposed ðEðtÞÞ; quarantined ðQðtÞ), suspect infective (I s ðtÞ), probable infective (I p ðtÞÞ and removed ðRðtÞÞ: It is assumed that individuals in the E class are exposed to SARS virus but not isolated from the community and individuals in the Q class are isolated from the community.

We first propose a deterministic model consisting of the six components, which is a system of six ordinary differential equations with various parameters modelling the feature of SARS transmission. However, unlike the SARS outbreaks in Hong Kong, Singapore and Toronto, which have been well studied and for which detailed data are available in the literature (see Donnelly et al., 2003; Lee et al., 2003; Tsang et al., 2003 for Hong Kong; Lipsitch et al., 2003 for Singapore; Booth et al., 2003; Poutanen et al., 2003 for Toronto) so that parameter values can be estimated from these reports, so far there have been few published reports about the SARS outbreak in Beijing. Moreover, the data about Beijing are very limited (Ministry of Health of P. R. China, 2003; WHO, 2003e) , with only the daily suspect cases, probable cases (including number of cases transferred from previously reported suspect cases), cleared suspect cases, number of patients released from hospitals, and deaths given (see Table 1 ). Because of the limitation of the data, we simplify the general model to two simpler models. The first is a two-component model consisting of only the suspect infectives ðI s ðtÞÞ and the probable infectives ðI p ðtÞÞ: The daily suspect and probable cases will be regarded as the inputs, the cleared suspect cases will be the output for I s ðtÞ; the number of SARS patients released from hospitals plus the deaths will be the output for I p ðtÞ: The model also include a term representing the cases transferred from I s ðtÞ to I p ðtÞ:

The second is a one-compartment model for the probable variable I p ðtÞ only. We will use the data in Table 1 to model the input, output and transmission terms. Then we use the mathematical models to make predictions of future incidence and to simulate the impact of additional control strategies.

Let f rs ðI s ðtÞÞ denote the removal rate of the suspected infective individuals at time t; f sp ðI s ðtÞÞ the transition rate from suspect class to probable class at time t; f rd ðt; I p ðtÞÞ the sum of the removal rate and the death rate of the probable infective individuals at time t: Statistics show that some individuals who are infectious are not quarantined at Beijing partially because some infective individuals in hospitals failed to report some of their close contacts for various reasons. These individuals are the sources for the spread of SARS outside the quarantined sites. Further, the data (Ministry of Health of P. R. China, 2003) show that the valid contact rate of SARS in Beijing depends on the numbers of infected individuals. One of the reasons is that the people in Beijing have been taking various protection measures such as reducing social activities, using mouth-masks, strengthening immune functions, etc.

Let lðI p Þ denote the valid contact rate of SARS at Beijing. Then the transition rate from the S compartment to the E compartment is f se :¼ lðI p ÞðlI p =ðS þ E þ R þ lI p ÞÞS; where l is the fraction of probable infectives who can transmit the disease. The quarantine rate for the E class is q (see the flow diagram in Fig. 1 ). k 01 and k 02 are fractions of the transition rates from E to I s and from E to I p ; respectively. k 1 and k 2 are fractions of the transition rates from Q to I s and from Q to I p ; respectively. k 3 is the fraction of the transition rate from Q to S: Using Fig. 1 , we have the model (1) 38 (4) 27 (2) 20 (3) 50 (4) 54 (2) 83 (2) 

In order to use valid data, we simplify (2.1) into two submodels. The first is a two-compartment model for I s and I p and the second is a scalar equation for I p only.

Our purpose is to consider the dynamics of suspect and probable cases. Note that the first two terms in the right-hand side of the fourth equation (i.e. the I s equation) of (2.1) represent the net flow rate from the E and Q compartments to the I s compartment (because of the I p ). Let us denote this input rate by f 0 : To fit this function and others, we set April 26, 2003 as time zero and take days as time units. Thus, April 27, 2003 is the first day, April 28, 2003 is the second day, etc. Then, we use the curve-fitting package (cftool) from Matlab 6.5, in which polynomial, exponential, power, rational, Weibull, Gaussian and Custom equations can be chosen and certain techniques can be used to fit data according to the mean-square error (MSE). By analysing the SARS data for Beijing from April 27 to June 10, 2003 (see Table 1 ), we can see that it can be approximated by a function aðI p ÞI p ; where (see Fig. 2 Let D denote the first two terms of the right-hand side in the fifth equation (i.e. the I p equation) of (2.1). Then D represents the net flow rate at which individuals enter the I p compartment directly from the E and Q compartments (which depends on I p ). By the curvefitting method, we obtain (see Fig. 5 For the removal rate of the probable class (recovered plus death), using the data from April 27 to June 10 in Table 1 , we can see that it can be approximated by a function cðtÞI p ; where (see Fig. 6 ) cðtÞ ¼ 0:00367ðt=10Þ 5 À 0:03791ðt=10Þ 4 þ 0:1384ðt=10Þ 3 À 0:2063 ðt=10Þ 2 þ 0:01156 t Fig. 5 . The direct transition rate from the E and Q compartments to I p compartment versus I p ; the solid line is made from DðI p Þ and the dots are plotted using the data in Table 1. Summarizing the above discussions, we obtain the simplified two-compartment suspect-probable model dI s dt ¼ aðI p ÞI p À bðtÞI s À eðtÞI s ; dI p dt ¼ DðI p Þ þ bðtÞI s À cðtÞI p : ð3:6Þ

First, note that the first three terms in the right-hand side of the I p equation in (2.1) represent the net flow rate from the E; Q and I s compartments to the I p compartment (which caused by I p ). If we denote this input rate by u; by analysing the data in where April 27 is specified as t ¼ 1 (see Fig. 7 ). The transmission rate b reflects the protection measures at two stages. At the first stage, people took protection measures such as reducing social activities, Table 1 . Fig. 7 . The per capita input rate of probable individuals to I p compartment versus I p ; the solid lines represent bðI p Þ and the dots are plotted using the data in Table 1. wearing mouth-masks, strengthening immune system, etc. The rate is a linear decreasing function of the number of infected individuals, because the more infectious individuals there are the more frightened the society is, so that more severe measures are taken. The second stage coincides with the period of keeping maximal protection measures. In this period, all the main protection measures have been performed, and its intensity stays invariant so that the rate mainly depends upon the number of infectious individuals. Furthermore, from Fig. 7 we can see that the values of bðI p Þ is quite low at the second stage, which implies that the disease can be eradicated.

The rate of removal and death for the probable class is same as cðtÞI p (see Eq. (3.5) and Fig. 6) . Hence, the I p equation in (2.1) is reduced to dI p dt ¼ bðI p ÞI p À cðtÞI p ; 1pt: ð3:8Þ

Using the suspect-probable model (3.6) and choosing April 27 as an initial point ðt ¼ 1Þ; the computer simulations are given in Fig. 8 for suspect I s :

Let us now change the initial data from I s ð1Þ ¼ 1191; I p ð1Þ ¼ 980 to I s ð1Þ ¼ 0; I p ð1Þ ¼ 1; that is, we assume that Beijing was completely susceptible and only one probable member was introduced into on April 27. Then by keeping the parameter values in (3.6) unchanged, numerical calculations from (3.6) show that there will be 16 suspected individuals on the 17th day.

Furthermore, if we let aðI p Þ ¼ À0:1257I p ðtÞ=980 þ 0:2712; 1pt; DðI p Þ ¼ ðÀ0:06881I p ðtÞ=980 þ 0:1461ÞI p ðtÞ; 1pt; that is, we keep the per capita suspect input rate and the direct transition rate from the E and Q compartments to I p compartment after the 17th day as the same as those before the 17th day. Basically, this means that intervention measures are not maintained. Now keep cðtÞ; bðtÞ and eðtÞ unchanged and take I s ð1Þ ¼ 0; I p ð1Þ ¼ 1: Then there will be 418 suspect individuals when t ¼ 43: This indicates the importance of maintaining intervention measures.

The probable model (3.8) fits the data well (see Fig. 9 ), where we fixed cðtÞ ¼ 0:1491 for tX45; which was obtained by fitting the data after 45 days.

We now consider the reproduction number for the SARS outbreak in Beijing. Since the exact infection period is not clear at this moment, which is crucial for the computation of the reproduction number, we look for an estimation. We claim that the mean infection duration of probable individuals varies from 5 to 16. The lower bound 5 is obtained from the study by Donnelly et al. (2003) in which it was reported that the mean time from the onset of clinical symptoms to the admission to hospitals is about 5 days. If there is no risk to transmit the disease by a probable individual quarantined in a hospital and the transmission occurs only in the time interval from the onset of clinical symptoms to the admission to a hospital, the five days can be used as a lower bound for the infection duration. The upper bound 16 is derived from Zhong (2003) in which it was reported that ''few probable individuals are infectious at the onset of symptoms, the infectiousness of probable individuals increases as the infected process develops, the maximal infection force occurs with highest temperature and heavy coughing in the second week and then falls quickly as the temperature decreases''. Now assume that the infection force is equally distributed in the time interval from the onset of symptoms to the 16th day. From the expression of cðtÞ; we can see that the expectation that an infective stays in I p compartment is 4.9436 during the first 5 days and 15.0298 during the first 16 days. By the standard method for the calculation of the reproduction number (see, for example van den Driessche and Watmough, 2002), we can see that a reasonable estimation for the reproductive number is 1:0698pR 0 p3:2524:

Suppose that bðI p Þ defined in (3.7) is valid for all tX1: If we take cðtÞ as defined in (3.5) for 1ptp45 and cðtÞ ¼ 0:1491 for t > 45; which is the mean value of the per capita removal rate of probable individuals from April 23 to June 21, and choose an initial value as I p ð1Þ ¼ 1 in (3.8), then we have the dynamical behavior shown in the left picture of Fig. 10 . From the figure, we can see that there will be 41 infective individuals when maximal control measures are taken if we introduce one infectious individual into a completely susceptible population in Beijing. Moreover, that figure shows that the number of infective individuals returns to 1 after 61 days and dies out as time evolves thereafter. Now suppose that the infection force from the first stage always takes effect, i.e. we take bðI p Þ ¼ À0:09589I p ðtÞ=980 þ 0:2164 for tX1: This means that the maximal control measures were not maintained. If we keep all others unchanged as the case in previous paragraph, then the disease can spread quickly (see the right graph in Fig. 10 ). In this case, the disease will be persistent at a level of 688 infective individuals and there will be 1000 infective individuals on the 43th day. Furthermore, if we replace the transmission rate by bðI p Þ ¼ À0:0958kI p ðtÞ=980 þ 0:2164 for tX1; where k represents the intensity of protection measures, numerical calculations show that the disease can be persistent at a level of 14 infectives and there will be 35 infective individuals on the 36th day even though k ¼ 50: This indicates the necessity of maintaining the maximal control policies for a suitable period of time in order to eradicate the disease.

To learn more about the mechanism of controling the disease successfully in Beijing from the point of view of modelling, we suppose that the SARS disease can activate for a long time and take cðtÞ as its mean value 0.03349 in first the 45 days. Then, Eq. (3.8) can be partitioned into two equations dI p dt ¼ I p ðÀ0:09589I p ðtÞ=980 þ 0:2164Þ À 0:03349I p ð4:1Þ and dI p dt ¼ 0:01029ðI p ðtÞ=1991Þ 2:842 I p À 0:03349I p : ð4:2Þ

Eq. (4.1) has an endemic equilibrium I p ¼ 1869:35 which is globally stable. For Eq. (4.2), there is an endemic equilibrium I Ã p ¼ 3015:8 which is unstable. Furthermore, I p ðtÞ decreases to zero as t evolves if I p ð0ÞoI Ã p ; and increases as t evolves if I p ð0Þ > I Ã p : If protection measures are merely used and the intervention policies are not adopted, then only (4.1) takes effect. It follows that the disease can be persistent at a level of 1869 probable individuals when t is large. Furthermore, if the intervention measures are implemented so late that the probable numbers are greater than 3016, then the above analysis shows that the disease will be persistent. This again shows that it is important to perform strong intervention measures before the number of infectives reaches I Ã p :

There have been a few published reports on studying the transmission dynamics of SARS in Hong Kong, Singapore and Toronto (Chowell et al., 2003; Lipsitch et al., 2003; Riley et al., 2003) , but so far there are few published scientific reports about the transmission dynamics in China Zhou et al., in press), where the virus was originally found. Beijing, the hardest hit city in the world with the most SARS cases, had already been removed from the WHO SARS list in June 2003. However, questions and mysteries about spread, reporting and treatment of SARS in Beijing remain. One of the fundamental questions is that how many of the reported suspect SARS cases were indeed probable cases and how the suspect individuals affect the transmission of the disease.

In this study, we proposed a mathematical model to simulate the SARS outbreak in Beijing. Differing from the SARS models in Chowell et al. (2003) , Lipsitch et al. (2003) and Riley et al. (2003) , our model consists of six subpopulations, namely susceptible, exposed, quarantined, suspect infectives, probable infectives and removed. By simplifying the model to a two-compartment suspect-probable model and a single-compartment probable model and using limited data, we are able to simulate the SARS outbreak in Beijing. We estimated that the reproduction number varies from 1.098 to 3.2524 and obtained certain important epidemiological parameters. The basic reproduction number R 0 is defined as the expected number of secondary cases produced in a completely susceptible population by a typical infective individual. This quantity determines the potential for an infectious agent to start an outbreak (Anderson and May, 1991; Brauer and Castillo-Chavez, 2000; Diekmann and Heesterbeek, 2000) . We found that the incidence rate has the characteristic of two stages. The first stage is the process of taking protection measures and quarantine policy and the second stage coincides with the process of maintaining control measures. Our study showed that it is necessary to implement the maximal control measures in the second stage for a certain period in order to eradicate the disease. Furthermore, the control measures in the second stage should be implemented before a threshold for the number of probable cases is reached.

The Chinese health professionals have substantial experience with a variety of anti-viruses, antibiotics, Chinese herbal medicines, and corticosteroids, and with using assisted ventilation in the treatment of SARS patients (Breiman et al., 2003; Zhong and Zeng, 2003) .

The SARS outbreak in Beijing shows that better and careful surveillance of suspected cases, rapid and accurate reporting of infected cases, transparent communication with the public about the disease, rigorous application of knowledge through public health measures are effective in fighting emerging infectious disease. Delaying reporting and implementation can result in disastrous public health consequences.

As more clinic and epidemiological data about the SARS outbreak in Beijing become available (Jiang et al., 2003; Liu et al., 2003; Xie et al., 2003) , we expect that we would be able to estimate the parameter values in the general model (2.1) and would have a better understanding of the SARS transmission dynamics in Beijing by analysing the general system (2.1). Also, in modelling real diseases such as SARS, stochastic processes should be considered. We leave these for future study.

Infectious Diseases of Humans

Mathematical Models in Population Biology and Epidemiology

Role of China in the quest to define and control severe acute respiratory syndrome

Clinical features and short-term outcomes of 144 patients with SARS in the Greater Toronto area

SARS outbreaks in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism

Mathematical Epidemiology of Infective Diseases: Model Building, Analysis and Interpretation

Epidemiological determinants of spread of causal agent of severe acute respiratory syndrome in Hong Kong

Identification of a novel coronavirus in patients with severe acute respiratory syndrome

Analysis of severe acute respiratory syndrome in Beijing

A novel coronavirus associated with severe acute respiratory syndrome

A major outbreak of severe acute respiratory syndrome in Hong Kong

Transmission dynamics and control of severe acute respiratory syndrome

Clinical features and therapy of 106 cases of severe acute respiratory syndrome

Bulletin of SARS epidemic situation

Coronavirus as a possible cause of severe acute respiratory syndrome

Identification of severe acute respiratory syndrome in Canada

Transmission dynamics of etiological agent of SARS in Hong Kong: impact of public health interventions

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Empirical analysis and forecasting for SARS epidemic situation. Beijing Daxue Xuebao

Case definitions for surveillance of severe acute respiratory syndrome (SARS) (www

Consensus document on the epidemiology of severe acute respiratory syndrome (SARS

Analyses on one case of severe acute respiratory syndrome 'super transmitter' and chain of transmission. Zhonghua Liu Xing Bing Xue Za Zhi (Chin

Clinical treatment references for SARS

Our strategies for fighting severe acute respiratory syndrome (SARS)

A discrete epidemic model for SARS transmission and control in China

We would like to thank the three anonymous referees for their criticisms, comments and suggestions.