Severe acute respiratory syndrome (SARS) is a highly contagious respiratory disease which is caused by the SARS Coronavirus. It is a serious form of pneumonia, resulting in acute respiratory distress and sometimes death. SARS emerged in China late 2002 and quickly spread to 32 countries causing more than 774 deaths and 8098 infections worldwide. One of the intriguing characteristics of the SARS epidemic was the occurrence of superspreading events. Super-spreading events for a specific infectious disease occur when certain infected individuals, called super-spreaders, infect more than the average number of secondary cases.
According to the U. S. Centers for Disease Control and Prevention, a person is a super-spreader if they cause more than 10 secondary infections. Such super-spreading individuals have been identified in the SARS outbreak and they are thought to have caused most of the secondary infections. For example, in Singapore, about 80% of infections have been attributed to only 5 superspreading individuals ( [7] ). One extreme super-spreading individual in Hong Kong caused more than 100 secondary infections ( [13] ). To contrast, in Singapore, most individuals caused 0 secondary infections ( [10] ). One differential equation model for super-spreading individuals was proposed before the recent SARS outbreak, by Kemper ([6] ). He presents a modified susceptibleinfected -removed (SIR) disease model to capture the effect of the super-spreading individuals in which the infected individuals are split into two different infected classes with different transmission rates. The superspreaders have a higher transmission rate, meaning more of their contacts with susceptible individuals result in a new infection than the regular infected individuals. This was the first model designed specifically to address the effect of super-spreading individuals on the course of infectious disease epidemics.
Due to the extreme influence of super-spreading individuals during the SARS outbreak, many mathematical epidemiologists developed models for the spread of SARS, which included the occurrence of the super-spreading events 1 . It is interesting to note how each incorporated the super-spreading events into their model. First are models where the super-spreading individuals have a higher transmission rate than the regular infected individuals, as in Kemper ([6] ). For example, Masuda, et al. ([10] ) developed a contact network model which had different transmission rates for regular infected individuals and for super-spreading individuals. Second are models in which parameters are taken from probability distributions, making the super-spreading individuals naturally appear as the right-hand tail of the distributions. Lloyd-Smith, et al. ( [9] ) developed a stochastic compartment model and Fang, et al. ([4] ) developed a spacial lattice combined with a deterministic compartment model wherein the individual reproduction numbers are drawn from a continuous probability distribution.
The goal of this paper is to capture the effect of the superspreading individuals using a modification of the classical SIR disease model. The modification is inspired by Li, et al. ([8] ) who determined that for SARS "the daily infection rate did not correlate with the daily total number of symptomatic cases but with the daily number of symptomatic cases who were not admitted to a hospital within 4 days of onset of symptoms." This means that the number of infected individuals is closely associated with the number of individuals who remain outside of isolation longer than most other infected individuals. These individuals remaining outside of isolation longer than normal are related to the disease severity. They have more contacts with susceptible individuals and thus more chances to spread the disease, in other words, they are the superspreaders. In our model, we split the infected individuals into two classes, with two different "removal" rates; these two rates determine how long an infected individual remains outside of isolation.
In this paper, we modify the classical SIR disease model to capture the effect of a super-spreading event, using the idea that super-spreading individuals stay out of isolation longer than individuals who are not super-spreading. A description of the SIP R model, the model equations, and some basic model properties are given in Section 2. We demonstrate the model using a small-scale and a largescale outbreak of SARS in Section 3. Our discussion follows in Section 4.
We conclude the introduction with a brief review of the classical SIR disease model. A more detailed review is given in Murray ([11] ), and in Ching, et al. ( [2] ) and the references given there.
The classical SIR disease model splits a fixed-size population, N , into three distinct classes: the susceptible individuals, S, those do not have the disease and can become infected; the infected individuals, I, those who have the disease and can infect susceptible individuals; and the removed individuals, R, those who have recovered, die, or moved into isolation. Individuals in the removed class gain permanent immunity and remain in the R class forever. Schematically, the SIR model is Figure 1 : SIP R model schematic, where S are the susceptible individuals, I are the regular infected individuals, P are the super-spreaders, R are the removed individuals, β is the transmission rate, b is the probability that a new infection will be a regular infected person, 1 − b is the probability that a new infection will be a super-spreading individual, ν 1 is the removal rate for a regular infected individual, and ν 2 is the removal rate for a super-spreading individual.
equations describing the SIR model is
where β is the transmission rate and ν is the removal 2 rate.
We describe the SIP R model, a modification of the classical SIR model, which captures the effect of superspreading individuals. Schematically, the model is given in Figure 1 . While the basic SIR model has one class of infected individuals, the SIP R model has a second class of infected individuals, the super-spreaders, denoted by the variable P . In the SIP R model we divide the population, of fixed size N , into four groups, namely the susceptible individuals, S, the regular infected individuals, I, the super-spreaders, P , and the removed individuals, R.
A susceptible individual can become infected through contact with either a regular infected individual or a super-spreading individual. Then with probability b the new infected individuals become regular infected individuals (move to class I) and with probability 1 − b the new infected individuals become super-spreading individuals (move to class P ). We assume most new infections are regular infected individuals, that is, b > 1 − b. We assume that the two infected classes, I and P , have the same transmission rate, β.
To capture the effect of the super-spreading individuals, we use the idea that individuals who are super-spreaders stay out of isolation longer than the regular infected individuals. The regular infected individuals stay out of isolation (the R class) for 1/ν 1 days and the superspreading individuals stay out of isolation for 1/ν 2 days, with 1/ν 1 < 1/ν 2 . Therefore, the model has two distinct removal rates namely, ν 1 and ν 2 , corresponding to regularly infected individuals, I, and super-spreading events, P , respectively. When a person is removed to the R class there is no possibility of becoming susceptible again, but rather they recover and gain permanent immunity, or die; in either case, they remain in the removed compartment forever.
Based on the previous descriptions and assumptions we formulate a system of four ordinary differential equations for the SIP R model.
subject to the following initial conditions
The SIP R model given by Equations (1) has the following properties.
• The SIP R model has a unique global solution.
• The components of the solution, S(t), I(t), P (t), and R(t), of the SIP R model are non-negative and bounded by N for all time, t ≥ 0.
• The SIP R model has equilibrium points (N, 0, 0, 0), (S * , 0, 0, R * ), for any 0 < S * < S 0 , and (0, 0, 0, N ).
• The individual reproduction number is
The properties of the SIP R are, in general, the natural modification of the corresponding properties of the classical SIR model. It is interesting to note that the R 0 for the SIP R model is the R 0 for the corresponding SIR model for each of the two infected classes, I and P , multiplied by the probability of a new infected individual becoming an I or a P , respectively.
The SIP R model can be used to analyze any infectious disease where super-spreading events have been identified. As a particular example, super-spreading events have been identified in one outbreak of measles as described in Paunio, et al. ( [12] ).
We use the SIP R model to study the spread of SARS on a small scale, in the Amoy Gardens apartment complex in Hong Kong, and on a large scale, in the entire Hong Kong Special Administrative Region. In both cases, we fit the model to the data using two parameter estimation methods. We begin with a description of the estimation procedures.
The general parameter estimation procedure was as follows. The SIP R system of equations (1) was solved repeatedly with parameter sets taken from the allowable range of possible parameters, and the least squares error between the cumulative number of cases of the solution and the actual data was computed. The least squares error was minimized. All computations were done with MATLAB using the function fminsearchbnd.
In parameter estimation Method 1, we used the parameters estimated with a fit of the classical SIR model to the data. For the data on the entire Hong Kong Special Administrative Region, we were unable to find SIR parameters. Therefore, we began with fitting the classical SIR model to the data. A summary of the SIR model parameters is given in Table 1 . We use the transmission rate of the SIR model as the transmission rate of the SIP R model. The transmission rate, the removal rate, and the initial number of susceptible individuals are used to determine the basic reproduction number R 0 . We estimate parameters b, ν 1 , and ν 2 so that together they have R 0 as calculated, and they satisfy b > 1 − b and ν 1 > ν 2 .
In parameter estimation Method 2, bounds for the parameters in the SIP R model were determined from the literature. Then, using these bounds, we estimate the parameters that give the best fit to the data. Details and references for these bounds are given in the next paragraph and in Subsections 3.1 and 3.2.
For estimation Method 2, we make some assumptions regarding the parameters that apply to both the small-and large-scale outbreaks. In situations with super-spreading individuals, it is assumed that most infected individuals are not super-spreaders, which is the case for the spread of SARS. SARS outbreak in Singapore 80% of infected individuals infected no one else. This leads us to set a lower bound on the probability of becoming a regular infected individual of 0.8, and we have 0.8 < b < 1. For SARS, after onset of symptoms regular infected individuals stayed out of isolation between 3 and 5 days ( [3] ); we assume the average and set ν 1 = 1/4. We assume that on average super-spreading individuals moved into isolation 10 days after they became infectious, that is, ν 2 = 1/10. Finally, many researcher have shown that the individual reproduction number, R 0 , for the SARS outbreak is between 1.5 and 4 ([1]). We assume this region for R 0 , which gives a constraint on the parameters β, b, ν 1 , and ν 2 .
We use the SIP R model to study the spread of SARS at the high rise apartment building Amoy Gardens, Block E, in Hong Kong. We assume that the cumulative number of confirmed cases of SARS in the Amoy Gardens, summarized in Table 2 Table 2 .
For parameter estimation Method 2, we survey literature to determine appropriate bounds on all of the SIP R model equations. We assume that the total number of infected individuals at t = 0 is 4, as in Ching, et al. ( [2] ). We conservatively assume 0 < β < 0.01. The bounds and constraints for each parameter as well as the estimated parameter values are given in Table 3 . The resulting number of confirmed cases of SARS are presented in Table 2 . We examine the first and third periods only, using the SIP R model. The first time period corresponds to quick spreading of the disease and we show this region has an R 0 value above 1. The third period corresponds to the end of the disease and correspondingly has an R 0 value below 1.
In 2003, the Hong Kong area had a total population of 6.803 million people. The number of susceptible individuals in Hong Kong during the SARS epidemic is not the entire population; the number of susceptible individuals must be approximated. Katriel formula to estimate the percent of the population that is susceptible during an epidemic. The number who are susceptible can be computed using the R 0 value and the percent of the population who became infected. Riley, et al. ([13] ) determine the R 0 value specifically for the SARS outbreak in Hong Kong is in the region 2.2 -3.7. Using this, we compute that the number of susceptible individuals falls in the range 1.8 million -3.1 million, where the smaller population correspond to larger R 0 values. We assume a large R 0 value, and, therefore we use the lower value, N = 1.8 million.
For parameter estimation Method 1, we use the parameters determined by fitting the data sets to the classical SIR model. The fixed and estimated parameters of the SIP R model are given in Tables 5 and 6 , corresponding to time period one and three, respectively. The resulting number of confirmed cases of SARS are presented in Table 4 .
For parameter estimation Method 2, we survey literature to determine appropriate bounds on all of the SIP R model equations. We assume that the total number of infected individuals at t = 0 is 1 for time period one, corresponding to the index case, and 470 for time period three 5 , corresponding to the number of confirmed cases on March 29, 2003. The bounds and constraints for each parameter as well as the estimated parameter values are given in Tables 5 and 6 . The resulting number of confirmed cases of SARS are presented in Table 4 .
We summarize the infections for each of the data sets and their fitted parameters in Table 7 . Specifically, the table contains the total number of individuals who became infected, the total number who became infected through contact with a regular infected (I) and super-spreading individuals (P ), the total number of regular infected (I) and super-spreading individuals (P ), and the individual reproduction numbers for the I and P classes (as SIR models).
We have presented a modification of the classical SIR disease model that captures the effect of super-spreading individuals on an infectious disease epidemic. Using an idea from the progression of the SARS outbreak, we distinguish the regular infected individuals from the superspreading individuals by how long they remain outside of isolation; the super-spreading individuals spend longer outside of isolation than most infected individuals. The model was fit to data from the SARS epidemic, using two different parameter estimation methods.
Parameter estimation Method 1, used parameters estimated for the classical SIR model. The resulting parameters show only a slight super-spreading behavior in all cases studied. For Amoy Gardens and the short-term outbreak in Hong Kong, the I and P classes are not distinguished by how long each stays outside of isolation. They are also not distinguished by their individual reproduction numbers. In each case, the two recovery rates are similar to eachother, and similar to the recovery rate Day   Feb 21 Mar 17  19  21  29  Apr 12  26  May 10  24  Jun 6  20  Confirmed  1  95  150  203  470  1108  1527  1674  1724  1750  1755  SIR  1 Table 7 : Summary of the outbreaks for each of the data sets and their fitted parameters: the total number of individuals who became infected; the total number who became infected through contact with a regular infected (I) and super-spreading individuals (P ); the total number of regular infected (I) and super-spreading individuals (P ); and the individual reproduction numbers for the I and P classes (as SIR models). of the SIR model estimated parameters. It is not surprising, therefore, that the final number of infected individuals, given by the final number of individuals in the removed class, in the two models are within just a few individuals.
The fitted parameters for the long-term outbreak in Hong Kong using parameter estimation Method 1 indicate that regular infected individuals move so quickly into isolation that we can disregard their influence on the disease spread. As evidence, over the entire course of the disease, they collectively only infect a total of 4.4 individuals. Due to the incredibly short time spent in the infected class, the individual reproductive number in this case is very low.
Though the outbreak in the entire Hong Kong Special Administrative Region was split into three time periods corresponding to three different disease dynamics, the fit parameters do not match the third time period data set as well as one might hope. This is noticed for the classical SIR model, and, since the SIR parameters are used for estimation Method 1, we also see this in the fit parameters for Method 1. One possible fix for this problem is to split the third time period into other periods in which the outbreak has common dynamics.
Parameter estimation Method 2, used research to set a priori bounds on the parameters of the SIP R model. In every case, the resulting parameters show superspreading behavior, for example, the two infected classes are distinguished by their individual reproduction rates. All three cases show that, on average, each regular infected individual infects less than 1 other individual, while each super-spreading individual infects more than 1 other individual.
Considering the specific case of the outbreak at the Amoy Gardens apartment complex, Method 2 predicts that, without any other intervention, almost all of the residents of Block E, will become infected. In this case, both the regular infected and super-spreading individuals have basic reproduction numbers larger than one. It is clear from the dire outcome predicted for the residents of the Amoy Gardens, Block E, super-spreading individuals must be brought into isolation as quickly as possible.
For the long-term outbreak in Hong Kong (in the third time period), the overall individual reproductive number is less than 1, which matches the notion that individuals were taking precautions to protect themselves, and the disease spread was slowing down. Using the SIP R model, we see that the regular infected individuals do have an individual reproductive number less than 1, however, the super-spreading individuals have a reproductive number larger than 1. Again, we see that it is imperative that the super-spreading individuals be brought into isolation as quickly as possible.
Finally, it is evident from the short-term outbreak in Hong Kong (in the first time period), that without any control measures the spread of SARS in Hong Kong would have been extreme. The model predicts that more than 1.7 million residents would have been infected over the course of the disease outbreak. Correspondingly, they would have experienced a large disease related mortality. The death rate for SARS is estimated to be around 10%, and so, Hong Kong would have lost an estimated 170,000 residents. Noting that the final death toll in Hong Kong was (though tragic) only 299 people, the control measures put in place by the government of Hong Kong saved thousands of lives.
The SIP R model is versatile; it can be used to examine an outbreak of any disease known to have super-spreading individuals, measles for example (see ( [12] )). On the other hand, the model was built using the idea that super-spreading individuals stay out of isolation longer than regular infected individuals, as during the SARS epidemic. (There are documented cases of SARS infected individuals violating strict isolation mandates 6 .) For diseases where the super-spreading behavior is a result of differing transmission rates, one should use the Kemper model ( [6] ).
Super-spreading individuals for infectious diseases pose a serious threat to public health. The SIP R model clearly demonstrates that infectious individuals must be removed from interactions with the susceptible individuals as quickly as possible to decrease the seriousness of an infectious disease epidemic.