Basic Reproductive Number (R₀)

The basic reproductive number (R₀) is the average number of secondary cases generated by a primary case. If R₀ < 1, an epidemic can not be sustained. On the other hand, if R₀ > 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₀, 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, γ₂, δ, 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₀) 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₂). 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₀

Figure 1

Figure 1. . Histograms of the six distributed parameters appearing in equation 1 with sample size 10⁵. The transmission rate was assumed to be exponentially distributed with mean 0.25, our estimated transmission rate...

We carried out an uncertainty analysis on the basic reproductive number (R₀) to assess the variability in R₀ that results from the uncertainty in the model parameters. We used a Monte Carlo procedure (simple random sampling) to quantify the uncertainty of R₀ to model parameters when these parameters are distributed. Similar methods have been used before (12–14). Parameters (k, γ₂, δ, α) 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, γ₂, δ, α, l) 10⁵ times, holding q, p, and ρ fixed. We then computed R₀ from each set. A probability density function for R₀ is obtained and can be statistically characterized. Here, we characterize R₀ by its median and interquartile range.

Sensitivity Analysis for R₀

We performed a sensitivity analysis on R₀ to quantify the effect of changes in the model parameters on R₀. Hence, we rank model parameters according to the size of their effect on R₀. Partial rank correlation coefficients (12–15) were computed between each of the parameters (with the exception of p, q, and ρ, which were held fixed) and R₀ as samples were drawn from the distributions, thus quantifying the strength of the parameter's linear association with R₀. The larger the partial rank correlation coefficient, the larger the influence of the input parameter on the magnitude of R₀. Because the distribution of the parameter l (relative infectiousness after isolation) is not known, we also studied the sensitivity of R₀ 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₀

Figure 2

Figure 2. . Empiric (dots) and stretched exponential estimated probability density function Prob(R₀) = a exp[–(R₀/b)^c] (solid line) (16) of R₀ for the cases of Toronto (a = 0.186, b =...

Figure 3

Figure 3. . (β, l) parameter space when R₀ < 1 obtained from the uncertainty analysis (black dots). The deterministic (β, l) level curve when R₀ = 1 is shown in by the...

The resulting R₀ distribution lies in the interquartile range 0.43–2.41, with a median of 1.10. Moreover, the probability that R₀ > 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.58 (0.24–1.18) (Table 2). Similarly, a lower rate of diagnosis a = 1/4.85 day^–1 and the relative infectiousness after isolation in Hong Kong (l = 0.43) and Singapore (l = 0.49) gives R₀ = 1.10 (0.44–2.29) and 1.17 (0.47–2.47), respectively (Figure 2). Perfect isolation (l = 0) gives R₀ = 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₀ lies at R₀ > 1. Furthermore, the median and interquartile range of R₀ are larger when p = 1, as has been assumed (8). In Figure 3 we show the (β, l) parameter space when R₀ < 1 obtained from our uncertainty analysis (14).

Sensitivity Analysis for R₀

Figure 4

Figure 4. . Boxplot of the sensitivity of R₀ estimates to varying values of l, the relative infectiousness after isolation has begun. l = 0 denotes perfect isolation while l = 1 denotes...

The transmission rate β and the relative infectivity during isolation (l) are the most influential parameters in determining R₀. The systematic decline in R₀ 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₀ (see Appendix). Because bounds exist on how much a given parameter can change in practice, achieving control (i.e., R₀ < 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.

Local Sensitivity Analysis of the Basic Reproductive Number

The sensitivity analysis approach through an exhaustive sampling of the parameter space provides a global measure of the sensitivity of model parameters. Another approach is to compute the sensitivity indices of the model parameters through local derivatives (18). This approach only provides a local measure as the sensitivity indices can change when the parameter values change. Here, we use local sensitivity analysis to corroborate our global sensitivity analysis results and discuss how this approach can be applied in the analysis of cost as part of a policy of outbreak control.

Let λ represent any of the 10 nonnegative parameters, β, ρ, p, q, k, γ₁, γ ₂, δ, α, and l, that define the basic reproductive number of our model (19)

(1)

If a "small" perturbation δλ is made to the parameter λ, a corresponding change will occur in R₀ as δR₀, where

The normalized sensitivity index Ψ_λ is the ratio of the corresponding normalized changes and is defined as

(2)

An approximation of the perturbed value of R₀, in terms of the sensitivity index is

where the 10 normalized sensitivity indexes are

with η = p(1–ρ)+ρ and γ₂ = a γ₁/( α –γ₁). For the values of the parameters used in this model, the sensitivity indices Ψ_β, Ψ_ρ, Ψ_p, Ψ_q, and Ψ_l are positive, Ψ_k = –Ψ_q, and the remaining indexes are negative. Furthermore, since all of the indexes (except Ψ_β) are functions of the parameters, the sensitivity indexes will change as the parameter values change.

For our specific case where β = 0.25, q = 0.1, p = 1/3, k = 0.15707, α = 0.2061, γ₁ = 0.035285, γ₂ = 0.0426, δ = 0.0279, and ρ = 0.77 and Toronto (l = 0.1) or Hong Kong (l = 0.43), the normalized sensitivity indices are computed. The sensitivity indices and the associated percentage changes needed to affect a 1% decrease in R₀ are given in Tables A1 and A2. Since the effective rate of patient isolation and the average rate of diagnosis are two feasible intervention strategies, we examine how changes to the parameters l and α affect R₀.

Let us first consider the outbreak in Hong Kong. The value α = 0.2061 means that the mean time to diagnose an infected person's illness is approximately 4.85 days. The sensitivity index Ψ_α = –0.1933 means that a 5.2% increase in α, which in turn requires a decrease of 5.7 hours of mean time to diagnosis, would result in a decrease of approximately 1% in R₀. Similarly, the sensitivity index Ψ_l = 0.5183 suggests that a 1.9% decrease in the value of l, that is going from 0.43 to 0.42 isolation effectiveness,2 results in a 1% decrease in R₀. In other words, individually a 5.2% increase in a or a 1.9% decrease in l both result in approximately a 1% decrease in R₀. For the particular values of the parameters chosen for Hong Kong, the most effective way to reduce R₀ is to decrease the transmission rate β and the parameter l (improve the effective isolation rate). In the case of Toronto, Ψ_α = –0.4758 means that a 2.1% increase in α, results in a 1% decrease in R₀, whereas Ψ_l = 0.2001 means a 5% decrease in l also results in a 1% decrease in R₀.

As can be seen from these two examples, the importance or ranking of the sensitivity indices can change as the values of the parameters change. Specifically, the sensitivity indices Ψ_l and Ψ_α satisfy the relationship

(3)

Figure A1

Figure A1. . Level curve of (l,α), where l(α+2γ₁+2δ) = δ+γ₂.

For the particular values of the parameters given above, the Figure A1 shows the level curve for the pair (l,α), where l(α+2γ₁+2δ) = δ+γ₂. The particular parameter values are either for Toronto (l,α) = (0.1, 0.2064) or for Hong Kong (l,α) = (0.43, 0.2064). Choosing the parameter values (l,α) below the level curve means that ||Ψ_l|| < ||Ψ_α||and the converse is true if (l,α) is chosen above the curve. Along the level curve, the magnitude of the sensitivities is equal. Notice that the level curve divides the parameter space into two regions, each of area A_above and A_below, respectively. Since A_above >> A_below, Ψ_l will be the dominant sensitivity index for randomly chosen (l,α).

One aspect of implementing an efficient intervention policy is the fact that limited resources are available. If one assumes, for example, that the strategies of isolation and diagnosis have associated 1% incremental costs in implementation of δC_I and δC_D, respectively, then a mixed strategy should be formulated that maximizes the effectiveness of a combined intervention. Specifically, if x denotes the magnitude of percentage decrease in l, and y denotes the % increase in a and assuming that there is a maximum amount of total additional resources available (δC_T), then the total additional cost of a new mixed isolation and diagnosis intervention policy must satisfy the inequality δC_Ix+δC_Dy < δC_T. Since the objective is to maximize the decrease in R₀, this means we want to maximize the objective function P = ||Ψ_l||x+||Ψ_α||y under the appropriate constraints. In a more general setting, additional nonlinear constraints could be involved, which case would require one to solve a nonlinear optimization problem. The situation in which the cost of diagnosis of infected persons may be much greater than the cost of isolation or vice versa is certainly of interest.