key: cord-284060-6eonuc8x
authors: Siriprapaiwan, Supatcha; Moore, Elvin J.; Koonprasert, Sanoe
title: Generalized reproduction numbers, sensitivity analysis and critical immunity levels of an SEQIJR disease model with immunization and varying total population size
date: 2018-04-30
journal: Mathematics and Computers in Simulation
DOI: 10.1016/j.matcom.2017.10.006
sha: 
doc_id: 284060
cord_uid: 6eonuc8x

Abstract An SEQIJR model of epidemic disease transmission which includes immunization and a varying population size is studied. The model includes immunization of susceptible people (S), quarantine (Q) of exposed people (E), isolation (J) of infectious people (I), a recovered population (R), and variation in population size due to natural births and deaths and deaths of infected people. It is shown analytically that the model has a disease-free equilibrium state which always exists and an endemic equilibrium state which exists if and only if the disease-free state is unstable. A simple formula is obtained for a generalized reproduction number R g where, for any given initial population, R g < 1 means that the initial population is locally asymptotically stable and R g > 1 means that the initial population is unstable. As special cases, simple formulas are given for the basic reproduction number R 0 , a disease-free reproduction number R d f , and an endemic reproduction number R e n . Formulas are derived for the sensitivity indices for variations in model parameters of the disease-free reproduction number R d f and for the infected populations in the endemic equilibrium state. A simple formula in terms of the basic reproduction number R 0 is derived for the critical immunization level required to prevent the spread of disease in an initially disease-free population. Numerical simulations are carried out using the Matlab program for parameters corresponding to the outbreaks of severe acute respiratory syndrome (SARS) in Beijing, Hong Kong, Canada and Singapore in 2002 and 2003. From the sensitivity analyses for these four regions, the parameters are identified that are the most important for preventing the spread of disease in a disease-free population or for reducing infection in an infected population. The results support the importance of isolating infectious individuals in an epidemic and in maintaining a critical level of immunity in a population to prevent a disease from occurring.

The SEQIJR model (see, e.g., [6, 8, 9] ) is a generalization of the well-known SEIR model of an infectious disease in which a population at risk is separated into the four categories of susceptible (S), exposed (E), infectious (I ) and recovered (R) (see, e.g., [1, 12, 27, 29] ). A susceptible person is an uninfected person who can be infected through contact with an infectious or exposed person, an exposed person is someone who has come into contact with an infectious person but is asymptomatic, an infectious person is symptomatic, and a recovered person is someone who has recovered from the disease. It is assumed that a recovered individual cannot become infected again. In an SEQIJR model, two extra categories of quarantined (Q) and isolated (J ) are added to the SEIR model. A quarantined person is an exposed person who is removed from contact with the general population and an isolated person is an infectious person who is removed from contact with the general population, usually by being admitted to a hospital. Quarantine and isolation are important in controlling epidemics because they are effective methods of reducing the contact between infected and susceptible populations.

In this paper, we consider a generalization of the basic SEQIJR model developed by Gumel et al. [6] for the severe acute respiratory syndrome (SARS) outbreaks of November 2002-July 2003. The SEQIJR model of Gumel et al. for the SARS outbreaks is useful for study for several reasons. Firstly, Gumel et al. suggested parameter values which can be used to obtain physically reasonable models for the four SARS outbreaks in Beijing, Hong Kong, Singapore and Toronto. Secondly, considerable data are available on the SARS outbreaks and mathematical models and the effectiveness of control measures have been studied by many authors (see, e.g., [3, 7, 11, 13, [17] [18] [19] [20] [21] 26, 30] ).

In building their model, Gumel et al. assumed that there was no immunization (no vaccine was available for SARS), that susceptible people were born into the community at a constant rate and that the population was constant. They analyzed the effects of quarantine and isolation on the transmission of SARS and showed that an effective isolation policy was more important than quarantine or a combination of quarantine and a less effective isolation policy in reducing the transmission of the disease. A similar conclusion was obtained by Koonprasert et al. [9] in a numerical study of the effectiveness of quarantine and isolation. However, as noted by a number of authors (see, e.g., [4, 5, 14, 15, 24, 28] ), the effectiveness of isolation requires strict hospital hygiene to reduce the hospital-based (nosocomial) spread of a disease.

The modifications introduced in our analysis include: (1) the existence of an immunization program, (2) changes in the total population due to births and deaths, (3) the derivation of generalized reproduction numbers for disease-free states and endemic equilibrium states, (4) a sensitivity analysis of the effect of variations in parameter values on the reproduction numbers for disease-free states and on the infected populations of the endemic equilibrium state, and (5) a derivation of a simple formula for the critical immunization level required to prevent spread of the disease in an initially disease-free population.

The organization of the paper is as follows. In Section 2, we define the mathematical model. In Section 3, we obtain analytical formulas for disease-free and endemic equilibrium points. Section 4 contains an analysis of the local stability of the equilibrium points and gives derivations of formulas for the generalized reproduction numbers. Formulas are derived for sensitivity indices in Section 5 and for critical immunization levels in Section 6. Section 7 and 8 contain results of numerical simulations using Matlab for the parameter values proposed by Gumel et al. [6] for four of the cities affected by the SARS outbreaks of 2002 and 2003, namely Beijing, Hong Kong, Singapore and Toronto. In Section 9, we give a discussion of results and conclusions.

A flow chart of the dynamics of our SEQIJR model is shown in Fig. 1 . The system of equations for the SEQIJR model that we use is given in Eqs. (1)- (6) . Definitions of the parameters in the model are given in Table 1 . For the immunization program, we assume that a vaccine is available and that a fraction ν of the susceptible population acquires immunity from vaccination per unit time. The effect of the immunization can then be modeled as a direct transfer of a fraction ν S of the susceptible people from the S to the R class per unit time. We also assume that people Table 1 Parameters for the SEQIJR model (rates are per day). Source: Adapted from Gumel et al. [6] Param Definition

Rate of inflow of susceptible individuals into a region or community through birth or migration. µ

The natural death rate for disease-free individuals ν

Rate of immunization of susceptible individuals β Infectiousness and contact rate between a susceptible and an infectious individual ε E Modification parameter associated with infection from an exposed asymptomatic individual ε Q Modification parameter associated with infection from a quarantined individual ε J Modification parameter associated with infection from an isolated individual γ 1 Rate of quarantine of exposed asymptomatic individuals γ 2

Rate of isolation of infectious symptomatic individuals σ 1 Rate of recovery of symptomatic individuals σ 2

Rate of recovery of isolated individuals k 1 Rate of development of symptoms in asymptomatic individuals k 2

Rate of development of symptoms in quarantined individuals d 1 Rate of disease-induced death for symptomatic individuals d 2

Rate of disease-induced death for isolated individuals in the R class cannot be infected for a second time.

where

is an effective infectiousness rate of the S population due to contact with asymptomatic exposed and quarantined populations (E, Q) and symptomatic infectious and isolated populations (I , J ), and where

are decay rate coefficients for S, E, Q, I , and J populations respectively.

Although the basic SEQIJR model only contains 6 equations, one each for S, E, Q, I, J, R, we have found it convenient (see, e.g., [9] ) to add 2 extra equations for easy computation of the total population size N (t) = S(t) + E(t) + Q(t) + I (t) + J (t) + R(t) at time t and the death rate D(t) due to disease at time t.

The equations of the Gumel et al. model can be obtained from Eqs. (1)-(6) by setting ν = 0 and assuming that the total population is constant, i.e., that Π = d 1 I + d 2 J + µN .

Theorem 1. The system (1)- (6) and (9) has two equilibrium states:

where

Proof. The equilibrium points are solutions of the system of nonlinear algebraic equations obtained by setting (6) and (9) . We begin by solving the system of nonlinear algebraic equations for S * , E * , Q * , J * , R * and N * as functions of I * to obtain the equations in (11) and (12) . Then, on substituting these solutions into (2), we obtain

Therefore, two possible values for I * are I * 1 = 0 and the solution for I * 2 in (11) . The proof is complete. □ For an endemic equilibrium point, all populations (S * 2 , E * 2 , Q * 2 , I * 2 , J * 2 , R * 2 , N * 2 ) must be non-negative and the infectious population I * 2 must be positive. Proof. We first prove that if I * 2 ̸ = 0 then α S > α N . From (11), we have

Then since µN * 2 = Π − α N I * 2 , we obtain

and therefore

Therefore, if α S ≤ α N , then I * 2 = 0 and the endemic equilibrium does not exist. We now prove that if α S > α N and µα L − D S α S > 0, then all populations are positive. If α S > α N and µα L − D S α S > 0, then µα L − D S α N > 0. So we have

Therefore, the endemic equilibrium exists and all populations are positive. We next note that if µα L − D S α S = 0, then from (11) I * 2 = 0 and the endemic equilibrium does not exist.

The proof is complete. □

We consider three cases for generalized reproduction numbers R g with the property that a given population is locally asymptotically stable if R g < 1 and unstable if R g > 1. These generalized reproduction numbers are useful because they give a simple test for the local asymptotic stability of an initial population. Three useful cases are as follows:

1. Basic reproduction number R 0 . This number is defined for a population that is 100% susceptible so that if R 0 < 1 and infected people enter the population then the number of secondary infections is less than the number of primary infections, whereas if R 0 > 1 then the number of secondary infections is greater. 2. Reproduction number R d f for disease-free equilibrium. 3. Reproduction number R en for endemic equilibrium.

The basic reproduction number R 0 and disease-free reproduction number R d f of the SEQIJR model can be obtained from the next-generation method of van den Driessche and Watmough [23] or by using Lyapunov's first method, i.e., by checking the eigenvalues of the Jacobian of the linearized system about populations with all infected populations equal to zero (see, e.g., [16] ). However, as far as the authors are aware, the endemic equilibrium reproduction number R en can only be obtained by checking the eigenvalues of the Jacobian of the linearized system about the endemic equilibrium point.

For the next-generation method, we write X = (E, S, Q, I, J, R) T and regard E as the next-generation "infected" population [23] . Then, (S, Q, I, J, R) T can be regarded as the remaining populations. Note that, in the SEQIJR model in this paper, an infected person always enters the E population first and then later enters the Q, I , J or R populations. Then the model can be written in the usual next-generation form as d X

Then, a disease-free reproduction number is defined as the largest eigenvalue of the matrix FV −1 , where F and V are the Jacobian matrices of F(X ) and V (X ) at the disease-free equilibrium point (

. We have used the Maple software package to find an analytic formula for the largest eigenvalue of FV −1 , i.e., for R d f , and obtained the result:

From the next-generation method, if the disease-free reproduction number R d f < 1, then the disease-free equilibrium point is locally asymptotically stable and if R d f > 1 then it is unstable.

We will now prove that R d f can be rewritten in a much simpler form and that R d f > 1 also corresponds to the condition that the endemic equilibrium exists. 

Proof. From (18) and using the definitions of α S , α E , α Q and α J in (12), we have

Then, the formula for the endemic infectious population I * 2 in (11) can be written as:

Therefore, I * 2 > 0 if and only if R d f > 1, i.e., if and only if the disease-free equilibrium point is unstable. To obtain the basic reproduction number R 0 , we consider the special case that the initial population is the 100% susceptible population. For this case, we have S = N , which corresponds to ν = 0 or D S = µ. Thus the basic reproduction number R 0 of the SEQIJR model is

The proof is complete. □

In the previous section, we used the next-generation method of [23] to obtain formulas for the reproduction number R d f for the disease-free equilibrium point and for the basic reproduction number R 0 . We also showed that the endemic equilibrium point exists if and only if R d f > 1.

In this section, we derive generalized reproduction numbers R g by using Lyapunov's first method (see, e.g., [16] ) of linearizing the system about an equilibrium point and checking that the real parts of all eigenvalues of the Jacobian at the equilibrium point are negative for R g < 1 and that the real part of at least one eigenvalue is positive for R g > 1.

The linearized model of the system at an equilibrium point (S * , E * , Q * , I * , J * , R * ) of (1)- (6) can be written in the form:

where x = (S, E, Q, I, J, R) T − (S * , E * , Q * , I * , J * , R * ) T and J * is the Jacobian matrix at an equilibrium point given by:

where (we omit the * on the populations in the matrix to save space)

As usual for linear systems, we can assume that Eq. (22) has a solution of the form x(t) = e λt v, where λ is an eigenvalue of J * and v is the corresponding eigenvector.

For the 6 × 6 Jacobian matrix in (23) it is not possible to obtain useful analytical expressions for all six eigenvalues either by a direct analysis or by solving the characteristic equation

where I is a 6 × 6 identity matrix. It is also not possible to obtain useful analytical tests for the real parts of eigenvalues to be negative from the six Routh-Hurwitz conditions (see, e.g., [16] ). However, it is possible to check two of the six Routh-Hurwitz conditions for the real parts of the six eigenvalues to be negative. These conditions are usually stated in terms of determinants derived from the coefficients of the characteristic polynomial. However, we will use the two equivalent necessary conditions for the real parts of the six eigenvalues to be negative that the sum of the six eigenvalues (trace (J * )) must be negative and the product of the six eigenvalues (det (J * )) must be positive. Now, for the disease-free populations, we have L *

since βε E α E < α L and µα L < D S α S , and therefore

For the endemic population, L * 2 > 0 and then

Now, from Eq. (16), we have

and trace(J * 2 ) < 0.

To evaluate the determinant of the Jacobian J * in (23), we first used the symbolic algebra toolbox in Matlab to find the following analytic expression:

Then, after considerable straightforward algebra, we obtained the following equivalent expression for the determinant.

where a generalized reproduction number can be defined as:

Formulas for the disease-free reproduction number (R d f ) and the endemic equilibrium reproduction number (R en ) can be easily obtained from Eq. (31). The results are summarized in the following theorem.

Theorem 4. Necessary conditions for stability of the two equilibrium points of the SEQIJR model are:

1. Disease-free equilibrium.

2. Endemic equilibrium.

Further, the endemic equilibrium exists if and only if R en < 1.

Proof. For a disease-free population, we have L * 1 = 0 and

, and therefore from (31) the reproduction numbers for the disease-free cases are as stated in part 1 of the theorem. These results are in agreement with the results from the next-generation method.

A simple formula for the endemic equilibrium number can be obtained by substituting the expressions in Eqs. (11) , (15) and (16) for the endemic equilibrium populations into Eq. (31). We then obtain the result given in part 2 of the theorem.

Since the endemic equilibrium exists only if µα L − D S α N > 0, it can be seen from Eq. (33) that R en < 1 when the endemic equilibrium exists and that R en → 1 when I * 2 → 0, i.e., at the transition between the endemic and disease-free equilibrium points. □ Note: In this example, it is possible to obtain the reproduction numbers from the trace and determinant of the Jacobian at an equilibrium point because the change in stability occurs due to a single real eigenvalue changing value from negative to positive. A change in value of the determinant cannot be used to detect a change in stability if the change is due to a complex conjugate pair of eigenvalues crossing the imaginary axis, e.g., if an Andronov-Hopf bifurcation occurs (see, e.g., [10] ). Therefore, a negative trace and a positive determinant for the Jacobian of a system of 6 equations are necessary but not sufficient conditions for the real parts of all eigenvalues to be negative and for an equilibrium point to be locally stable. However, for any given set of parameter values it is very easy to check if R g < 1 corresponds to asymptotic stability of an equilibrium point by numerically computing all eigenvalues of the Jacobian matrix with any of the standard mathematical software packages such as Matlab, Mathematica, Maple etc.

Sensitivity analysis is important as it can be used for determining the parameters which are of most importance in reducing the level of a disease (see, e.g., [2, 22] ). The normalized sensitivity index for a quantity Q with respect to a parameter h is defined by:

The sensitivity indices can be calculated in at least three ways, namely, (1) by direct differentiation of formulas for Q, [25] ). In this section, we will use the method of direct differentiation as it gives analytical expressions for the indices. However, as a check on our formulas we have compared numerical values from our formulas with numerical results we obtained using the method of Chitnis et al.

The normalized sensitivity indices can be computed from the formula in (19) as

We used the symbolic algebra package in Matlab to derive normalized sensitivity indices for the fourteen parameters of the model from Eq. (35) and then simplified the expressions to obtain the list shown in Table 2 .

The sensitivity indices for the basic reproduction number R 0 can be obtained from the results in Table 2 by setting ν = 0, i.e., D S = µ.

In this section, we derive formulas for the sensitivity indices for the endemic equilibrium populations. In particular, we derive indices for the endemic infectious population I * 2 . The sensitivity indices for the other endemic populations can then be easily obtained from the sensitivity indices for I * 2 using the formulas in Eqs. (11) and (12). The formulas for the sensitivity indices for I * 2 can be computed from the formula in (11) as:

where

The formulas for the coefficients ξ h , ψ h , ω h of the sensitivity indices in (37) are shown in Table 3 . The formulas for the natural birth rate Π and natural death rate µ are not included in Table 3 because these parameters are very difficult to change for human populations. As stated previously, the indices can also be computed using the method of Chitnis et al. [2] . As a check on our results, we have compared numerical results obtained from the formulas in Table 3 with results obtained from the method of Chitnis et al. Table 3 Formulas for endemic sensitivity coefficients in Eq. (37). 

From the disease-free reproduction number, we can easily find a formula for a critical immunization level required to prevent the spread of the SEQIJR disease. Since, R d f = µα L D S α S < 1 is required to stop the disease spreading and D S = µ + ν, we have the critical immunization level

where R 0 is the basic reproduction number. Further, the critical ratio of immune (recovered) population to total population is R *

For the models given in Eqs. (1)- (6) and (9), We have computed the equilibrium points, generalized reproduction numbers, sensitivity indices, dynamical behavior and the critical immunization levels by using the parameter values estimated by Gumel et al. [6] for the SARS epidemics in Beijing, Canada, Hong Kong and Singapore. In this section, we will show the results for Beijing. We will then briefly discuss any differences between the Beijing results and the results for the other three regions in Section 8.

The meanings of the parameters in the SEQIJR model are defined in Table 1 and the parameter values that we use are shown in Table 4 . To examine the effect of immunization, we will use values of ν = 0.01 corresponding to a disease-free equilibrium, ν = 0.000002 corresponding to an endemic equilibrium, and a critical value of ν corresponding to a value for the disease-free reproduction number of R d f ≈ 1. Substituting the parameters for Beijing in Table 4 into the system of differential equations for disease transmission in Eqs. (1)- (6) and (9), we obtain Eqs. (40).

d J dt = 0.5I + 0.125Q − 0.047234J, Table 5 Initial population values used for numerical integration. 

For the parameter values in Table 4 and an immunization rate ν = 0.01, the disease-free equilibrium point is

and the endemic equilibrium point does not exist as some populations are negative. The disease-free and basic reproduction numbers are:

That is, the disease-free equilibrium is locally asymptotically stable and the population with 100% susceptible population will be unstable. These conclusions can also be checked by using Matlab to compute the eigenvalues of the Jacobian (23) for the disease-free equilibrium and the 100% susceptible population. These eigenvalues are: Since the real parts of all eigenvalues for the disease-free equilibrium are negative and the real part of one eigenvalue of the 100% susceptible population is positive, the disease-free equilibrium state will be locally asymptotically stable and the 100% susceptible population will be unstable.

Using the Matlab program we have numerically integrated the system of Eqs. (40) for the four sets of initial conditions given in Table 5 to obtain the dynamical solutions shown in Fig. 2 . The graphs show that, for the four different sets of initial conditions, the four infected populations first decrease and then increase before finally converging to zero, i.e., the initial points are not locally asymptotically stable even though there is long-term convergence to the disease-free equilibrium point. However, it can also be seen that at the longer times, where the population values are close to the disease-free equilibrium point, the solutions converge asymptotically to this equilibrium state. This behavior for the longer times in Fig. 2 is in agreement with Eq. (41) and the values of the reproduction numbers and eigenvalues given above. 

Using the analytical results given in Table 2 , we have computed the normalized sensitivity indices shown in Table 6 of the disease-free reproduction number R d f = 0.01366 < 1 for the parameter values of Table 4 and ν = 0.01. For these parameter values, the sensitivity indices

and Φ(R d f |k 2 ) are positive and the remaining indices are negative. Table 6 Normalized sensitivity indices Φ(R d f |h) and required % changes in parameter for 1% reduction in R d f for disease-free equilibrium state (Beijing).

For these parameter values, R d f = 0.01366 < 1 and the solution converges to a disease-free equilibrium point. Then, a reduction in the value of R d f corresponds to a faster approach to equilibrium, i.e., to a faster disappearance of the disease. The values in column 2 of Table 6 are the sensitivity indices and the values in column 3 are % changes in parameter values required to give a 1% decrease in R d f . For example, in order to get a 1% decrease in the value of R d f , it is necessary to decrease the values of β, µ, ε J , k 1 and k 2 by 1.000%, 1.004%, 1.055%, 198.8% and 7428%, respectively. Also, in order to get a 1% decrease in the value of R d f , it is necessary to increase the values of ν, σ 2 , d 2 , σ 1 , γ 2 , d 1 and γ 1 by 1.003%, 1.156%, 12.2%, 26.2%, 112.2%, 196.9% and 205.8%, respectively. Since the natural death rate µ cannot be easily changed, the most effective methods of reducing R d f are to reduce the value of the infectiousness and contact rate β, increase the value of the immunization rate ν, reduce the value of the isolation modification parameter ε J and increase the value of the isolation recovery rate σ 2 .

For the Beijing parameter values in Table 4 and ν = 0.000002, we obtain the endemic equilibrium point Therefore, the disease-free equilibrium point is unstable and the endemic equilibrium point is locally asymptotically stable. As shown in Gumel et al. [6] (see also [9] ), an endemic equilibrium point is also obtained for ν = 0.

Using Matlab, we integrated the system of Eqs. (1)-(6) for ν = 0.000002 for the four sets of initial conditions given in Table 5 . Since the dynamical behavior for the four infected populations are similar, we will only show the solutions for the infectious population (I ) and for the susceptible (S) and recovered (R) populations. The solutions for the S, I and R populations are shown in Fig. 3 first for a short integration time to show the initial behavior and then for a long integration time to show the convergence to an endemic equilibrium. It can be seen that the convergence to the endemic equilibrium is very slow and that many disease outbreaks occur during this period. This type of behavior was very common in "childhood" diseases such as measles, mumps, poliomyelitis, rubella, whooping cough etc. before vaccines became available for these diseases. Although it is difficult to see in Fig. 3 , we found that epidemics can recur when the fraction of the susceptible population exceeds a certain critical level due to births or if the population of the immune recovered population falls below a certain critical level due to deaths. Similar results for the dynamical behavior of populations in the endemic region have also been reported by Koonprasert et al. [9] in a study on the effectiveness of quarantine and isolation in an SEQIJR model.

In the sensitivity analysis we do not consider variations in Π , µ, ε E and ε Q . The parameter Π represents the natural birth rate and the parameter µ represents the natural death rate for the population. These parameters cannot be changed easily for human populations. Also, we omit the modification parameters for infection from the exposed population ε E and the quarantined population ε Q since these parameters have been assumed to be zero. We therefore only consider sensitivity indices for the remaining eleven parameters.

Using the analytical results given in Table 3 and Eq. (36), and the parameter values in Table 4 with ν = 0.000002, we obtain the values shown in Table 7 for the normalized sensitivity indices of the endemic equilibrium populations. As noted above, the values of R d f = 3.807 > 1 and R en = 0.3756 < 1 show that the endemic equilibrium is locally asymptotically stable. Since the sensitivity indices are functions of parameter values, the values in Table 7 will change if parameter values are changed. The interpretation of these indices is that a positive index means that I * 2 increases if the parameter value is increased and a negative index means that I * 2 decreases if the parameter value is increased with the magnitude giving the rate of change.

In order to reduce the spread of the disease, it is necessary to reduce the infectious population (I * 2 ) and isolated population (J * 2 ) since these populations are capable of direct transmission of the disease to uninfected people. It is also desirable to reduce the number of deaths from the disease, i.e., to increase the total surviving population N * 2 . From the indices in Table 7 , the most effective methods of reducing I * 2 are, in order, to increase the isolation rate (γ 2 ), reduce the transfer rate from exposed to infectious population (k 1 ), increase the quarantine rate (γ 1 ) and reduce the contact rate between infectious and susceptible people (β). The most effective methods of reducing J * 2 are, in order, to increase the recovery rate of isolated individuals (σ 2 ), reduce the contact rate β, and reduce the modification parameter ε J . It can be seen from Table 7 that increasing the death rates d 1 and d 2 of the infectious and isolated individuals would also reduce the infectious and isolated populations. This method is clearly not an acceptable control method for a human disease. However, this method was adopted as an effective control method in the case of the spread of avian flu among chickens in China and "mad-cow disease" in England, and it has also been used as an effective method of controlling some other diseases among animals and plants.

The most effective methods for increasing the total surviving population N * 2 , i.e., to reduce the death rate, are, in order, increase the recovery rate of infectious individuals (σ 2 ), reduce the contact rate β and reduce the modification parameter ε J , i.e., to reduce the level of nosocomial disease as discussed in the introduction. 

Using the parameter values given in Table 4 , we have computed and compared the disease-free and endemic results for the four regions. The overall behavior of the results for all four regions is similar, but with some differences in detail. We will compare the results for the reproduction numbers, the sensitivity indices and the critical immunity levels. Table 9 Ranking for best parameters to reduce R d f .

Region

Ranking for best parameters to reduce I * 2 .

Region

Beijing β, σ 2 , ν 0.0001 γ 2 , k 1 , γ 1 0.000002 Canada β, σ 2 , ν 0.000024 β, σ 2 , γ 2 0.000002 Hong Kong β, σ 2 , ν 0.000061 γ 2 , β, k 1 0.000002 Singapore β, σ 2 , γ 2 0.0000028 β, σ 2 , γ 2 0.000002

The results for the basic reproduction numbers and the disease-free and endemic reproduction numbers are shown in Table 8 . For each case, we have numerically integrated the differential equations and checked the asymptotic stability of equilibrium points by computing the value of the reproduction numbers and by computing the eigenvalues of the Jacobian matrices. In all cases we have shown that the numerical results agree with the theoretical analysis.

For each of the four regions, we have also computed the sensitivity indices of the disease-free reproduction number R d f and the endemic equilibrium populations for all of the parameters in the model.

From these calculations, we can rank the indices in order of effectiveness in reducing a disease. The results for the disease-free equilibrium R d f are shown in Table 9 . The table shows the top three indices for each region for a disease-free case (R d f < 1), and a critical immunization case R d f ≈ 1.

For R d f < 1, the rankings for all regions are the same. The most effective methods of reducing R d f are, in order, to reduce the infectiousness and contact rate β between infectious and susceptible people, increase the value of the immunization rate ν and reduce the modification parameter ε J for the infectiousness and contact rate of isolated people. It can be seen from Table 9 that reducing (index is positive) the natural death rate µ is the third most effective method of reducing the value of R d f . However, since the natural death rate cannot be varied easily for human populations, this is not an acceptable control method for a human disease. For R d f ≈ 1, the most effective methods of reducing R d f are, in order, to reduce the infectiousness and contact rate β, reduce the isolation modification parameter ε J and increase the recovery rate of isolated individuals σ 2 . Table 10 shows the rankings for the most effective methods for reducing the endemic infectious populations I * 2 in the four regions for R d f ≈ 1 and for R d f > 1. For R d f ≈ 1, the rankings for Beijing, Canada, Hong Kong are the same. The most effective methods of reducing I * 2 are, in order, to reduce the contact and infectiousness rate β, increase the recovery rate of isolated individuals σ 2 and increase the value of the immunization rate ν. For Singapore, the most effective methods of reducing I * 2 are, in order, to reduce the infectiousness and contact rate β, increase the recovery rate of isolated individuals σ 2 and increase the isolation rate γ 2 . For R d f > 1, the most effective methods of reducing I * 2 for Beijing are, in order, to increase the isolation rate γ 2 , reduce the transfer rate from exposed to infectious k 1 and increase the quarantine rate γ 1 . For Canada the most effective methods of reducing I * 2 are, in order, to reduce the infectiousness and contact rate β, increase the recovery rate of isolated individuals σ 2 and increase the isolation rate γ 2 . For Hong Kong the most effective methods of reducing I * 2 are, in order, to increase the isolation rate γ 2 , reduce the infectiousness and contact rate β and reduce the transfer rate from exposed to infectious k 1 . For Singapore the most effective methods of reducing I * 2 are the same as for the R d f ≈ 1 case.

The critical values of the immunization rate ν and the critical ratio of recovered to total population can be calculated from Eqs. (38) and (39) and parameter values in Table 4 for the SARS epidemics in Beijing, Canada, Hong Kong and Singapore in 2002 and 2003. The results for the four regions are shown in Table 11 . The results show that there is a big difference between the critical immunization levels and critical immunized ratio in the four regions. The difference is due to the different values of the basic reproduction number R 0 = α L α S in the four regions. The difference in R 0 values in turn is mainly due to differences in the values of the infection factor β and the modification factors for infection from the quarantined population (ε Q ) and the isolated population (ε J ).

We have studied an SEQIJR model for infectious disease transmission which generalizes the model developed by Gumel et al. [6] for the SARS epidemics of 2002 and 2003 by including immunization and a varying total populations size.

We have shown analytically that this model has a disease-free equilibrium state which always exists and an endemic equilibrium state which exists if and only if the disease-free equilibrium is unstable. We have defined generalized reproduction numbers R g for any given equilibrium state such that the state is locally asymptotically stable if R g < 1 and unstable if R g > 1. As special cases, we have obtained simple formulas for the reproduction numbers for a 100% susceptible population (the basic reproduction number R 0 = α L α S ) and for the disease-free equilibrium R d f = µα L (µ+ν)α S by both the next-generation method [23] and from the eigenvalues of the Jacobian of the linearized model about the disease-free equilibrium. In these formulas, α L can be interpreted as the rate of new infections from all infected individuals and α S as the rate of reduction of the susceptible population due to transfers to other population groups including to the recovered immune group. The effect of the immunization rate (ν) in reducing the reproduction number and increasing the stability of disease-free equilibrium states can be clearly seen by comparing the R 0 and R d f formulas. We have also obtained simple formulas for a reproduction number R en for the endemic equilibrium state and shown that the endemic equilibrium is locally asymptotically stable if R en < 1 and does not exist if R en > 1. We have also derived formulas for sensitivity indices of the disease-free and basic reproduction numbers, and the endemic equilibrium populations for changes in the values of model parameters.

For the numerical study of the model, we have used values of model parameters estimated by Gumel et al. [6] for the SARS epidemics in 2002-2003 in Beijing, Canada, Hong Kong and Singapore. For each of the four regions we have considered parameter values corresponding to a disease-free equilibrium and an endemic equilibrium. From the computed sensitivity indices for each of the regions, we have shown that the most important parameters for reducing infection are, in order, (1) for high immunization rates (ν ≈ 0.01) reduce the overall infection rate (β), increase the immunization rate, and reduce the modification factor for infection by isolated individuals (ε J ), (2) for critical immunization rates with R d f ≈ 1, reduce β, reduce ε J and increase recovery rate of isolated individuals (σ 2 ), (3) for low immunization rates (endemic equilibrium) the most important for Beijing and Hong Kong is to increase the isolation rate of infectious individuals (γ 2 ) and the most important for Canada and Singapore is to reduce β.

The numerical results for the reproduction numbers and the sensitivity indices clearly support the importance of isolation and to a lesser extent quarantine in controlling the spread of a disease.

We have obtained very simple formulas in terms of the basic reproduction number R 0 for critical values of the immunization rate and the ratio of the recovered (immune) and total populations that are required to stop the spread of a disease. Our numerical results for the four regions considered in this paper clearly support the importance that is now widely accepted, except by anti-vax campaigners, of maintaining a critical immunization level of "herd immunity" to stop the spread of a disease. The differences in the critical values of the immunization rate in the four regions clearly show the importance of factors such as reducing infection rates, and increasing isolation and quarantine, in stopping an epidemic that can occur when the immunization ratio falls below a critical level.

Unfortunately, in 2002 and 2003 there was no known immunization available for SARS. Therefore, the immunization results in this paper are not relevant for the SARS outbreaks of 2002-2003. However, the inclusion of immunization in the study of an SEQIJR model is important as the SEQIJR model is a very general model which can be applied to a number of other diseases that are also mainly transmitted by person to person contact and for which immunization is available, e.g., influenza and "childhood diseases" such as measles, mumps, whooping cough etc. A similar model with the addition of a bird reservoir of the disease could also be used to analyze avian influenza. SEIR models have, of course, been studied by many authors and there is now an extensive literature on them. However, we would note that the formulas for generalized reproduction numbers, sensitivity indices and critical immunization levels for the SEQIJR model could easily be applied to SEIR models by removing the quarantined and isolated compartments and suitably redefining some of the parameters.

Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model

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

SARS spreads new outlook on quarantine models

Modeling the Invasion and Spread of Contagious Diseases in Heterogeneous Populations

Modelling strategies for controlling SARS outbreaks

SARS outbreak

Using a Mathematical Model to Study the Effect of Quarantine and Isolation for Controlling SARS Outbreaks

An analysis of quarantine and isolation in an SEQIJR model

A major outbreak of severe acute respiratory syndrome in Hong Kong

Global dynamics of a SEIR model with varying total population size

Transmission dynamics and control of severe acute respiratory syndrome

Characterizing transmission and control of the SARS epidemic: Novel stochastic spatio-temporal models

Curtailing transmission of severe acute respiratory syndrome within a community and its hospital

Introduction to Dynamic Systems: Theory, Models and Applications

Efficiency of quarantine during an epidemic of severe acute respiratory syndrome-Beijing

Evaluation of control measures implemented in the severe acute respiratory syndrome outbreak in Beijing

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

Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza

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

A double epidemic model for the SARS propagation

Viral dynamics of an HIV model with latent infection incorporating antiretroviral therapy

Simulating the SARS outbreak in Beijing with limited data

Dynamical analysis and perturbation solution of an SEIR epidemic model

Critical role of nosocomial transmission in the toronto SARS outbreak

Analysis and control of an SEIR epidemic system with nonlinear transmission rate

A discrete epidemic model for SARS transmission and control in China

This research is supported by the Centre of Excellence in Mathematics, the Commission on Higher Education, Thailand. This research was also supported by the Science and Technology Research Institute of King Mongkut's University of Technology North Bangkok. We would also like to thank colleagues in the Department of Mathematics, Faculty of Applied Science, King Mongkut's University of Technology North Bangkok for encouragement and critical discussions.