Model Parameters and Outbreak Control for SARS

Tool for estimating basic reproductive number for the SARS outbreak suggests need for multiple methods of control.

Taiwan, 94% of SARS cases occurred through transmission in hospital wards (6), and similar effects occurred in Hong Kong and Singapore (7). Although the SARS epidemic was eventually controlled, the measures used to achieve that control varied greatly in scope from one place to another. Control of an outbreak relies partly on identifying what disease parameters are likely to lead to a reduction in the reproduction number R 0 . Here we calculate the dependence of R 0 on model parameters.

Methods
Two models of the SARS epidemic that incorporate the effects of quarantine and early detection of new cases but assume perfect isolation were recently introduced (8,9). A slightly different model was used to quantify the role that fast diagnosis and efficient isolation of patients played in Toronto's outbreak (10). This model predicted control in Toronto and showed that lack of immediate action would have been catastrophic (11). The model incorporates differences in the population's susceptibility (3) by dividing the population into classes S 1 (high risk) and S 2 (low risk). A low-risk group in the age range <19 years can be observed from the age-specific incidence in Hong Kong (3). The low-risk class (S 2 ) has a reduced susceptibility to SARS, measured by the parameter p (0 < p < 1). While p = 0 denotes no susceptibility to SARS, p = 1 indicates that both susceptible classes are equally susceptible to SARS. Initially, S 1 = ρN and S 2 = (1-ρ)N, where N is the total population size and ρ is the initial proportion of fully susceptible (S 1 ) persons. Susceptible persons exposed to SARS enter the exposed class (assumed to be asymptomatic) with a rate proportional to β and remain there for a mean incubation period of 1/k. The possibility of reduced transmission from the exposed class is included through the parameter q (0 < q < 1), a relative measure of infec-

Model Parameters and Outbreak
Control for SARS tiousness. Once symptomatic, exposed persons progress to the infectious class (illness not yet diagnosed), where they may recover at the rate γ 1 , die at rate δ, or enter the diagnosed class at rate α. Isolation mechanisms may be put in place in the diagnosed class to reduce their impact on transmission. The relative infectiousness after isolation has begun is measured by the parameter l (0 < l < 1) so that l = 0 denotes perfect isolation and l = 1 denotes ineffective isolation.

Basic Reproductive Number (R 0 )
The basic reproductive number (R 0 ) is the average number of secondary cases generated by a primary case. If R 0 < 1, an epidemic can not be sustained. On the other hand, if R 0 > 1, an epidemic typically occurs.
The basic reproductive number derived from our model (10) is given by the formula . This equation includes 10 parameters of which 2, the diagnostic rate (α) and the relative infectiousness during isolation (l), are widely recognized as being amenable to modification by medical intervention. The transmission rate (β) is defined as the number of persons infected per infectious person per day. This differs from R 0 , which is the average number of secondary cases that an infectious person generates when introduced into a susceptible population. Definitions for the remaining parameters are provided in Table 1.

Parameter Estimation
Baseline values for k, γ 2 , δ, and α are taken from the mean values estimated in reference 3. Because whether asymptomatic persons (exposed class) can transmit the disease is not known, we have fixed q = 0.1 (the relative infectiousness of exposed, asymptomatic persons) as in reference 10.
The model parameters Θ = (β, l) are fitted to Hong Kong data (2) by least squares fit to the cumulative number of cases C (t, Θ) (equation 1 in reference 10). All other parameters are fixed to their baseline values (Table 1). We used a computer program (Berkeley Madonna, R.I. Macey and G.F. Foster, Berkeley, CA) and appropriate initial conditions for the parameters for the optimization process, which was repeated 10 times (each time the program is fed with two different initial conditions for each parameter) before the "best fit" was chosen. The best fit gives β = 0.25 and l = 0.43. We also estimated the relative infectiousness after isolation (l) for the case of Singapore (l = 0.49) by following the least squares procedure described above. However, for the case of Toronto, not enough data were available on the initial growth of the outbreak. Hence, we only estimated l from Toronto data after control measures were put in place on March 26 (10,11), where l = 0.1. We used the transmission rate (β) obtained from Hong Kong data as the baseline value (Table 1).
We revised earlier estimates for ρ and p (10) (both affect R 0 ) using data from the age distribution of residents and the age-specific incidence of SARS in Hong Kong, as reported (3). The revised estimates are ρ = 0.77 (the initial proportion of the population at higher risk) and p = 1/3 (the measure of reduced susceptibility in S 2 ). The lower-risk subpopulation lies in the age range <19. It constitutes approximately 23% of Hong Kong's population (3). The fact that most of the SARS cases included in the epidemiologic studies of the Toronto outbreak (5) were transmitted in hospitals limits the use of general demographic data from Toronto in the estimation of ρ and p. Hence, we used the parameters estimated from the situation in Hong Kong. Baseline values for all the parameters are given in Table 1.

Uncertainty Analysis for R 0
We carried out an uncertainty analysis on the basic reproductive number (R 0 ) to assess the variability in R 0 that results from the uncertainty in the model parameters. We Baseline values for k, γ 2 , α, ρ, p and δ have been taken from reference 3. b β = 0.25 is our estimated transmission rate in Hong Kong. c l = 0 means perfect isolation, while l = 1 means no isolation.
used a Monte Carlo procedure (simple random sampling) to quantify the uncertainty of R 0 to model parameters when these parameters are distributed. Similar methods have been used before (12)(13)(14). Parameters (k, γ 2 , δ, α) were assigned a different probability density function (PDF) (Figure 1), which is taken from reference 3. The relative measure of infectiousness of persons after isolation procedures are put in place (l) was assumed to be uniformly distributed in the interval (0 < l < 1). The observed heterogeneity in transmission rates during the SARS epidemic is modeled here by assuming that β is distributed exponentially with mean 0.25 person -1 day -1 (our estimate of the transmission rate in Hong Kong). Parameters q, p, and ρ are fixed to their baseline values (Table 1). We sampled the set of six parameters (β, k, γ 2 , δ, α, l) 10 5 times, holding q, p, and ρ fixed. We then computed R 0 from each set. A probability density function for R 0 is obtained and can be statistically characterized. Here, we characterize R 0 by its median and interquartile range.

Sensitivity Analysis for R 0
We performed a sensitivity analysis on R 0 to quantify the effect of changes in the model parameters on R 0 . Hence, we rank model parameters according to the size of their effect on R 0 . Partial rank correlation coefficients (12)(13)(14)(15) were computed between each of the parameters (with the exception of p, q, and ρ, which were held fixed) and R 0 as samples were drawn from the distributions, thus quantifying the strength of the parameter's linear association with R 0 . The larger the partial rank correlation coefficient, the larger the influence of the input parameter on the magnitude of R 0 . Because the distribution of the parameter l (relative infectiousness after isolation) is not known, we also studied the sensitivity of R 0 to various distributions of l. Distributions of l used for the Monte Carlo calculation of the partial rank correlation coefficients are: a) l ∼ β (a = 2, b = 2) where β is used to denote a beta distribution. Here, the likelihood of l is bell-shaped with mean 0.5 and variance 0.05; b) l ∼ β (a = 1, b = 2), the likelihood of l decreases linearly in the [0,1] interval; and c) l ∼ β (a = 2, b = 1), the likelihood of l increases linearly in the [0,1] interval.

Uncertainty Analysis for R 0
The resulting R 0 distribution lies in the interquartile range 0.43-2.41, with a median of 1.10. Moreover, the probability that R 0 > 1 is 0.53. The same Monte Carlo procedure, but with fixed values of l = 0.1 and α = 1/3 day -1 for Toronto (i.e., after implementing control measures on March 26), give a median and interquartile range for the distribution of R 0 = 0.58 (0.24-1.18) ( Table 2). Similarly, a lower rate of diagnosis α = 1/4.85 day -1 and the relative infectiousness after isolation in Hong Kong (l = 0.43) and Singapore (l = 0.49) gives R 0 = 1.10 (0.44-2.29) and 1.17 (0.47-2.47), respectively ( Figure 2). Perfect isolation (l = 0) gives R 0 = 0.49 (0.19-1.08). Especially noteworthy is that even in cases when eventual control of an outbreak is achieved (Toronto and a hypothetical case of perfect isolation), 25% of the weight of the distribution of R 0 lies at R 0 > 1. Furthermore, the median and interquartile range of R 0 are larger when p = 1, as has been assumed (8). In Figure 3 we show the (β, l) parameter space when R 0 < 1 obtained from our uncertainty analysis (14). Figure 1. Histograms of the six distributed parameters appearing in equation 1 with sample size 10 5 . The transmission rate was assumed to be exponentially distributed with mean 0.25, our estimated transmission rate in Hong Kong. Here l is assumed to have a beta distribution (l~ β [1,2]). Alternative distributions for l were also used as described in the text. All other distributions were taken from reference 3.
The transmission rate β and the relative infectivity during isolation (l) are the most influential parameters in determining R 0 . The systematic decline in R 0 with increasing l in the range [0,1] is illustrated in Figure 4. Furthermore, our results do not change if we assume the three distributions mentioned in the Methods section (sensitivity analysis) for the parameter l. Table 3 shows the partial rank correlation coefficients for the other three possible distributions of l. The transmission rate is ranked first independent of the distribution of l. The relative infectiousness after isolation is ranked second when l comes from distributions (a) and (b) and ranked third when it comes from distribution (c) (see Methods). Our sensitivity analysis is corroborated by computing local derivatives on R 0 (see online Appendix at http://www.cdc.gov/ncidod/ EID/vol10no7/03-0647_app.htm). Because bounds exist on how much a given parameter can change in practice, achieving control (i.e., R 0 < 1) can require changing parameters other than those with the highest partial rank correlation coefficient. For example, reference 10 showed that control of the outbreak in Toronto relied on both a reduction in l and 1/α, even though α is ranked fairly low by the partial rank correlation coefficient.

Conclusion
We have estimated R 0 for the cases of Toronto, Hong Kong, and Singapore (Table 2) through an uncertainty analysis shown in equation 1. Our estimates for R 0 agree with the empirical R 0 observed from the data of the first week of the SARS outbreak in Singapore (8). A stretched exponential distribution fits the resulting distributions of R 0 for the different locations ( Figure 2). Even though the median of R 0 is <1 when perfect patient isolation is assumed (l = 0), we find that 25% of our R 0 distribution lies at R 0 > 1. That is, implementing a single method for control may not be sufficient to contain a SARS outbreak. Control may require modifying more than one parameter amenable to intervention. In our model, these parameters include the diagnostic rate (α), the relative infectiousness after isolation has begun (l), and the per capita transmission rate (β). The fact that α and l are not independent, but are tightly coupled, favors control.
Our expression for R 0 incorporates the effects of diagnosis-isolation strategies. Moreover, our approach includes differential susceptibility (p) and effective population size (ρ). Most models take p = 1, even though data from Hong Kong show that a low-risk subpopulation lies in the age range <19, approximately 23% of Hong Kong's population (3). The assumption p = 1 thus overestimates R 0 .
Our sensitivity analysis shows that the transmission rate (β) and the relative infectiousness after isolation in   hospitals (l) have the largest effect on R 0 . With the exception of a few measures, such as closing schools, no clear policies would modify β directly. This means that a substantial effort must be (and has been) made by the medical community to modify other parameters, such as the diagnostic rate. Hence, the strong sensitivity of R 0 to the transmission rate β indicates that efforts in finding intervention strategies that manage to systematically lower the contact rate of persons of all age groups promise an effective means for lowering R 0 . Such strategies may include using face masks (the probability of transmission per contact may be reduced), washing hands, and avoiding large crowds (large public events). Associated with the role of screening, diagnosis, and the effective isolation of patients is the issue of cost. We cannot ignore or minimize the value of stringent quarantine measures and the probability of compliance combined with the economic effect of lost wages (thousands were quarantined in Taiwan, Hong Kong, and Singapore [17]), the costs associated with screening at airports and hospitals, the cost associated with closing hospitals; and the costs associated with isolating SARS patients and exposed persons (see online Appendix for a brief discussion).  (Table  1). l = 0 denotes perfect isolation; l = 1 denotes no isolation. Figure 4. Boxplot of the sensitivity of R 0 estimates to varying values of l, the relative infectiousness after isolation has begun. l = 0 denotes perfect isolation while l = 1 denotes no isolation. The boxplot shows the median and the interquartile range of R 0 obtained from Monte Carlo sampling of size 10 5 . Table 3. Partial rank correlation coefficients (PRCCs) between each of the input parameters and R 0 from Monte Carlo sampling of size 10 5 for different distributions of the relative infectiousness after isolation (l) as described in the text Probability distribution Input parameters in order of decreasing PRCC (shown in parenthesis) β (a = 2, b = 2) β (0.92), l (0.57), δ (0.53), γ 2 (0.35), α (0.32), k (0.13) β (a = 1, b = 2) β (0.90), l (0.60), δ (0.44), α (0.39), γ 2 (0.26), k (0.12) β (a = 2, b = 1) β (0.92), δ (0.60), l (0.51), γ 2 (0.40), α (0.22), k (0.11)