Generalized reproduction numbers, sensitivity analysis and critical immunity levels of an SEQIJR disease model with immunization and varying total population size.

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.

Introduction

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 (), exposed (), infectious () and recovered () (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 () and isolated () 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.

The SEQIJR model

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), (2), (3), (4), (5), (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 of the susceptible people from the to the class per unit time. We also assume that people in the class cannot be infected for a second time. where is an effective infectiousness rate of the population due to contact with asymptomatic exposed and quarantined populations () and symptomatic infectious and isolated populations (, ), and where are decay rate coefficients for , , , , and populations respectively.

Fig. 1

The flow chart of the SEQIJR model (generalized from: Gumel et al. [6]).

Table 1

Parameters for the SEQIJR model (rates are per day).

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

Although the basic SEQIJR model only contains 6 equations, one each for , we have found it convenient (see, e.g., [9]) to add 2 extra equations for easy computation of the total population size at time and the death rate due to disease at time . The equations of the Gumel et al. model can be obtained from Eqs. (1), (2), (3), (4), (5), (6) by setting and assuming that the total population is constant, i.e., that .

Equilibrium analysis

The system (1), (2), (3), (4), (5), (6) and (9) has two equilibrium states:

Disease-free:

Endemic: where

The equilibrium points are solutions of the system of nonlinear algebraic equations obtained by setting in the system (1), (2), (3), (4), (5), (6) and (9) .

We begin by solving the system of nonlinear algebraic equations for and as functions of to obtain the equations in (11) and (12). Then, on substituting these solutions into (2), we obtain Therefore, two possible values for are and the solution for in (11). The proof is complete.  □

For an endemic equilibrium point, all populations must be non-negative and the infectious population must be positive.

The endemic equilibrium exists with all populations non-negative and positive infected populations if and only if and .

We first prove that if then . From (11), we have Then since , we obtain Therefore, if , then and the endemic equilibrium does not exist.

We now prove that if and , then all populations are positive. If and , then . So we have Therefore, the endemic equilibrium exists and all populations are positive.

We next note that if , then from (11) and the endemic equilibrium does not exist.

Finally, if and The proof is complete.  □

Generalized reproduction numbers

We consider three cases for generalized reproduction numbers with the property that a given population is locally asymptotically stable if and unstable if . 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:

Basic reproduction number . This number is defined for a population that is 100% susceptible so that if and infected people enter the population then the number of secondary infections is less than the number of primary infections, whereas if then the number of secondary infections is greater.

Reproduction number for disease-free equilibrium.

Reproduction number for endemic equilibrium.

The basic reproduction number and disease-free reproduction number 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 can only be obtained by checking the eigenvalues of the Jacobian of the linearized system about the endemic equilibrium point.

The next-generation method

For the next-generation method, we write and regard as the next-generation “infected” population [23]. Then, can be regarded as the remaining populations. Note that, in the SEQIJR model in this paper, an infected person always enters the population first and then later enters the , , or populations. Then the model can be written in the usual next-generation form as , where Then, a disease-free reproduction number is defined as the largest eigenvalue of the matrix , where and are the Jacobian matrices of and at the disease-free equilibrium point . We have used the Maple software package to find an analytic formula for the largest eigenvalue of , i.e., for , and obtained the result: From the next-generation method, if the disease-free reproduction number , then the disease-free equilibrium point is locally asymptotically stable and if then it is unstable.

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

An alternative formula for the disease-free reproduction number is and the endemic equilibrium exists if and only if , i.e., if and only if the disease-free equilibrium point is unstable. The formula for the basic reproduction number is .

From (18) and using the definitions of , , and in (12), we have Then, the formula for the endemic infectious population in (11) can be written as: Therefore, if and only if , i.e., if and only if the disease-free equilibrium point is unstable.

To obtain the basic reproduction number , we consider the special case that the initial population is the 100% susceptible population. For this case, we have , which corresponds to or . Thus the basic reproduction number of the SEQIJR model is The proof is complete.  □

Generalized reproduction numbers and asymptotic stability of equilibrium points

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

In this section, we derive generalized reproduction numbers 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 and that the real part of at least one eigenvalue is positive for .

The linearized model of the system at an equilibrium point of (1), (2), (3), (4), (5), (6) can be written in the form: where and 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 , where is an eigenvalue of and 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 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 ()) must be negative and the product of the six eigenvalues (det ()) must be positive.

Now, for the disease-free populations, we have , , and . Then, since and , and therefore For the endemic population, and then Now, from Eq. (16), we have and .

To evaluate the determinant of the Jacobian 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 () and the endemic equilibrium reproduction number () can be easily obtained from Eq. (31). The results are summarized in the following theorem.

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

Disease-free equilibrium.

Endemic equilibrium. Further, the endemic equilibrium exists if and only if .

For a disease-free population, we have 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 , it can be seen from Eq. (33) that when the endemic equilibrium exists and that when , 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 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

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 with respect to a parameter is defined by: The sensitivity indices can be calculated in at least three ways, namely, (1) by direct differentiation of formulas for , (2) by the method in Chitnis et al. [2] of linearizing the Eqs. (1), (2), (3), (4), (5), (6) about the equilibrium point to obtain a system of linear algebraic equations for the indices and then solving the set of equations or (3) by a Latin hypercube sampling technique (see, e.g., [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.

Sensitivity analysis of the disease-free reproduction number

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.

Table 2

Normalized sensitivity indices for disease-free reproduction number.

Param	Normalized sensitivity index
β	1
μ	νD_S−μD_E−μD_I+μβk_1D_QD_Jα_L[ε_ED_QD_J+ε_QD_J(γ_1−k_1α_Q)
	+ε_J(k_2(γ_1−k_1α_Q)−k_1D_Qα_J)].
ν	−νD_S
ε_E	βε_Eα_Eα_L
ε_Q	βε_Qα_Qα_L
ε_J	βε_Jα_Jα_L
γ_1	βα_QD_Jα_L[D_Jε_Q+k_2ε_J]−γ_1D_E
γ_2	βγ_2k_1D_QD_Jα_L[D_QD_Jε_E+γ_1D_Jε_Q+(k_1D_Q+k_2γ_1)ε_J]−γ_2D_I
σ_1	βσ_1k_1D_QD_Jα_L[D_QD_Jε_E+γ_1D_Jε_Q+k_2γ_1ε_J]−σ_1D_I
σ_2	−βσ_2α_Jε_JD_Jα_L
k_1	−βD_Jα_L[D_Jα_Eε_E+D_Jα_Qε_Q+k_2α_Qε_J]+γ_1+μD_E
k_2	−βk_2α_QD_QD_Jα_L(D_Jε_Q−με_J)
d_1	βd_1k_1D_QD_Jα_L[D_QD_Jε_E+γ_1D_Jε_Q+k_2γ_1ε_J]−d_1D_I
d_2	−βd_2D_Jα_Lα_Jε_J

The sensitivity indices for the basic reproduction number can be obtained from the results in Table 2 by setting , i.e., .

Sensitivity analysis of the endemic equilibrium point

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 . The sensitivity indices for the other endemic populations can then be easily obtained from the sensitivity indices for using the formulas in Eqs. (11) and (12).

The formulas for the sensitivity indices for can be computed from the formula in (11) as: where The formulas for the coefficients , , 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).

Par	ξ	ψ	ω
β	μα_L	μα_L	0
ν	−να_S	−να_N	0
ε_E	μβε_Eα_E	μβε_Eα_E	0
ε_Q	μβε_Qα_Q	μβε_Qα_Q	0
ε_J	μβε_Jα_J	μβε_Jα_J	0
γ_1	γ_1α_E(−D_S	γ_1α_E(−D_Sd_2k_2D_JD_Q.	γ_1α_E
	+μβD_Q(ε_Q+ε_Jk_2D_J))	+μβD_Q(ε_Q+ε_Jk_2D_J))
γ_2	γ_2k_1(−D_ED_S+μβε_Jk_1D_J	γ_2k_1(−D_Sd_2D_J(k_1+γ_1k_2D_Q)+μβε_Jk_1D_J	D_Eγ_2k_1
	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))
σ_1	σ_1k_1(−D_ED_S	σ_1k_1(−D_Sd_2γ_1k_2D_JD_Q	σ_1D_Ek_1
	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))
σ_2	−σ_2μβε_Jα_JD_J	σ_2α_JD_J(D_Sd_2−μβε_J)	0
k_1	−α_E(D_S(k_1−D_E)	−α_E(−D_Sd_2γ_1k_2D_JD_Q	(k_1−D_E)α_E
	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))
k_2	μβk_2α_QD_JD_Q(ε_J(D_Q−k_2)	k_2α_QD_JD_Q(−μβε_QD_J	0
	−ε_QD_J)	+(μβε_J−D_Sd_2)(D_Q−k_2))
d_1	d_1k_1(−D_ED_S	d_1k_1(−D_S(k_1+d_2γ_1k_2D_JD_Q)	d_1D_Ek_1
	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))	+μβ(ε_E+ε_Qγ_1D_Q+ε_Jγ_1k_2D_JD_Q))
d_2	−d_2μβε_Jα_JD_J	−d_2α_JD_J(μβε_J+D_S(D_J−d_2))	0

Critical immunization levels

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, is required to stop the disease spreading and , we have the critical immunization level where is the basic reproduction number. Further, the critical ratio of immune (recovered) population to total population is

Numerical results for Beijing

For the models given in Eqs. (1), (2), (3), (4), (5), (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 corresponding to a disease-free equilibrium, corresponding to an endemic equilibrium, and a critical value of corresponding to a value for the disease-free reproduction number of . Substituting the parameters for Beijing in Table 4 into the system of differential equations for disease transmission in Eqs. (1), (2), (3), (4), (5), (6) and (9), we obtain Eqs. (40).

Table 4

Parameter values for SARS outbreaks in Beijing, Canada, Hong Kong and Singapore [6].

Parameter	Value				Unit
	Beijing	Canada	HongKong	Singapore
Π	408	136	221	136	day−1
μ	0.000034	0.000034	0.000034	0.000034	day−1
β	0.23	0.2	0.15	0.21	day−1
ε_E	0	0	0	0	None
ε_Q	0	0	0	0	None
ε_J	0.82	0.36	0.84	0.2	None
γ_1	0.1	0.1	0.1	0.1	day−1
γ_2	0.5	0.5	0.5	0.5	day−1
σ_1	0.0413	0.0337	0.0337	0.0337	day−1
σ_2	0.0431	0.0386	0.0386	0.0386	day−1
k_1	0.1	0.1	0.1	0.1	day−1
k_2	0.125	0.125	0.125	0.125	day−1
d_1	0.0055	0.0079	0.0079	0.0079	day−1
d_2	0.0041	0.068	0.0068	0.0068	day−1

Disease-free equilibrium point and asymptotic stability

For the parameter values in Table 4 and an immunization rate , 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.

Dynamical solutions for disease-free equilibrium

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.

Table 5

Initial population values used for numerical integration.

Set	S(0)	E(0)	Q(0)	I(0)	J(0)	R(0)
1	12×10^6	10	0	3	0	0
2	13×10^6	1×10^5	1×10^5	1×10^4	4×10^3	2×10^5
3	14×10^6	2×10^5	2×10^5	2×10^4	5×10^3	3×10^5
4	15×10^6	5×10^5	4×10^5	4×10^4	1×10^4	6×10^5

Example of convergence to a disease-free equilibrium for four different sets of initial conditions. The , , , populations rapidly converge to zero and the and populations converge to equilibrium values and (Beijing).

(a) Susceptible population.

(b) Exposed population.

(c) Quarantined population.

(d) Infectious population.

(e) Isolated population.

(f) Recovered population.

Sensitivity analysis of disease-free reproduction number

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 for the parameter values of Table 4 and . For these parameter values, the sensitivity indices , , , , , and are positive and the remaining indices are negative.

Table 6

Normalized sensitivity indices and required % changes in parameter for 1% reduction in for disease-free equilibrium state (Beijing).

Param	Φ(R_df|h)	% change	Param	Φ(R_df|h)	% change
μ	+0.996	−1.004%	ν	−0.997	+1.003%
β	+1.000	−1.000%	ε_E	+0.000	+0.000%
ε_Q	+0.000	+0.000%	ε_J	+0.948	−1.055%
γ_1	−0.0049	+205.8%	γ_2	−0.0089	+112.2%
σ_1	−0.038	+26.2%	σ_2	−0.865	+1.156%
k_1	+0.0050	−198.9%	k_2	+0.0001	−7428%
d_1	−0.0051	+196.9%	d_2	−0.082	+12.2%

For these parameter values, and the solution converges to a disease-free equilibrium point. Then, a reduction in the value of 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 decrease in . For example, in order to get a decrease in the value of , it is necessary to decrease the values of , , , and by , , , and , respectively. Also, in order to get a decrease in the value of , it is necessary to increase the values of , , , , , and by , , , , , and , respectively. Since the natural death rate cannot be easily changed, the most effective methods of reducing 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 and increase the value of the isolation recovery rate .

Endemic equilibrium point and asymptotic stability

For the Beijing parameter values in Table 4 and , we obtain the endemic equilibrium point The value of the disease-free reproduction number is and the value of the endemic reproduction number is . Further, the eigenvalues of the Jacobian for the disease-free equilibrium and for the endemic equilibrium are:

Disease-free Endemic 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 .

Dynamical solutions for endemic equilibrium

Using Matlab, we integrated the system of Eqs. (1), (2), (3), (4), (5), (6) for 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 () and for the susceptible () and recovered () populations. The solutions for the , and 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.

Example of convergence to an endemic equilibrium for four different sets of initial conditions (Beijing). The graphs in the left column show the initial behavior for days. The graphs in the right column show the long-time behavior for days.

(a) Susceptible population.

(b) Infectious population.

(c) Recovered population.

In the sensitivity analysis we do not consider variations in , , and . 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 and the quarantined population 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 , we obtain the values shown in Table 7 for the normalized sensitivity indices of the endemic equilibrium populations. As noted above, the values of and 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 increases if the parameter value is increased and a negative index means that decreases if the parameter value is increased with the magnitude giving the rate of change.

Table 7

Sensitivity indices for endemic equilibrium state (Beijing).

Parameter	Sensitivity indices
	S2∗	E2∗	Q2∗	I2∗	J2∗	R2∗	N2∗
β	−1.0237	+0.3326	+0.3326	+0.3326	+0.3326	+0.3062	−0.0237
ε_J	−0.9703	+0.3152	+0.3152	+0.3152	+0.3152	+0.2902	−0.0225
ν	+0.0013	−0.0185	−0.0185	−0.0185	−0.0185	+0.0013	+0.0013
γ_1	+0.0055	−0.5017	+0.4983	−0.5017	+0.0206	−0.0009	+0.0006
γ_2	+0.0101	−0.0033	−0.0033	−0.9176	+0.0376	−0.0017	+0.0012
k_1	−0.0057	−0.4981	−0.4981	+0.5019	−0.0204	+0.0012	−0.0007
k_2	−0.0002	+0.0001	−0.9997	+0.0001	+0.0002	+0.0002	−0.00001
σ_1	+0.0418	−0.0136	−0.0136	−0.0891	−0.0497	−0.0089	+0.0037
σ_2	+0.9481	−0.3080	−0.3080	−0.3080	−1.2205	−0.2013	+0.0832
d_1	+0.0014	−0.0005	−0.0005	−0.0105	−0.0053	−0.0054	−0.0037
d_2	+0.0214	−0.0070	−0.0070	−0.0070	−0.0938	−0.0880	−0.0608

In order to reduce the spread of the disease, it is necessary to reduce the infectious population () and isolated population () 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 .

From the indices in Table 7, the most effective methods of reducing are, in order, to increase the isolation rate (), reduce the transfer rate from exposed to infectious population (), increase the quarantine rate () and reduce the contact rate between infectious and susceptible people (). The most effective methods of reducing are, in order, to increase the recovery rate of isolated individuals (), reduce the contact rate , and reduce the modification parameter . It can be seen from Table 7 that increasing the death rates and 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 , i.e., to reduce the death rate, are, in order, increase the recovery rate of infectious individuals (), reduce the contact rate and reduce the modification parameter , i.e., to reduce the level of nosocomial disease as discussed in the introduction.

Comparison of results for Beijing, Canada, Hong Kong and Singapore

Using the parametervalues 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.

Reproduction numbers and sensitivity indices

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.

Table 8

Reproduction numbers.

	Beijing	Canada	Hong Kong	Singapore
Disease-free	ν=0.01
R_0	4.031	1.708	2.804	1.082
R_df	0.01366	0.005787	0.009502	0.003668
Endemic	ν=0.000002
R_df	3.807	1.613	2.649	1.022
R_en	0.3756	0.6887	0.4863	0.9786

For each of the four regions, we have also computed the sensitivity indices of the disease-free reproduction number 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 are shown in Table 9. The table shows the top three indices for each region for a disease-free case (), and a critical immunization case .

Table 9

Ranking for best parameters to reduce .

Region	R_df<1	R_df≈1
Beijing	β,ν,(μ),ε_J	β,ε_J,σ_2
Canada	β,ν,(μ),ε_J	β,ε_J,σ_2
Hong Kong	β,ν,(μ),ε_J	β,ε_J,σ_2
Singapore	β,ν,(μ),ε_J	β,ε_J,σ_2

For , the rankings for all regions are the same. The most effective methods of reducing 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 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 . 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 , the most effective methods of reducing are, in order, to reduce the infectiousness and contact rate , reduce the isolation modification parameter and increase the recovery rate of isolated individuals .

Table 10 shows the rankings for the most effective methods for reducing the endemic infectious populations in the four regions for and for . For , the rankings for Beijing, Canada, Hong Kong are the same. The most effective methods of reducing are, in order, to reduce the contact and infectiousness rate , increase the recovery rate of isolated individuals and increase the value of the immunization rate . For Singapore, the most effective methods of reducing are, in order, to reduce the infectiousness and contact rate , increase the recovery rate of isolated individuals and increase the isolation rate .

Table 10

Ranking for best parameters to reduce .

Region	R_df≈1	ν value for R_df≈1	R_df>1	ν value for R_df>1
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

For , the most effective methods of reducing for Beijing are, in order, to increase the isolation rate , reduce the transfer rate from exposed to infectious and increase the quarantine rate . For Canada the most effective methods of reducing are, in order, to reduce the infectiousness and contact rate , increase the recovery rate of isolated individuals and increase the isolation rate . For Hong Kong the most effective methods of reducing are, in order, to increase the isolation rate , reduce the infectiousness and contact rate and reduce the transfer rate from exposed to infectious . For Singapore the most effective methods of reducing are the same as for the case.

Critical immunity levels

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 in the four regions. The difference in 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 () and the isolated population ().

Table 11

Critical immunization levels for prevention of disease.

	Beijing	Canada	Hong Kong	Singapore
Critical ν	1.03×10^−4	2.407×10^−5	6.135×10^−5	2.803×10^−6
Critical R/N ratio	0.7519	0.4415	0.6434	0.07616

Discussion and conclusions

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 for any given equilibrium state such that the state is locally asymptotically stable if and unstable if . As special cases, we have obtained simple formulas for the reproduction numbers for a 100% susceptible population (the basic reproduction number ) and for the disease-free equilibrium 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, can be interpreted as the rate of new infections from all infected individuals and 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 and formulas. We have also obtained simple formulas for a reproduction number for the endemic equilibrium state and shown that the endemic equilibrium is locally asymptotically stable if and does not exist if .

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 () reduce the overall infection rate (), increase the immunization rate, and reduce the modification factor for infection by isolated individuals (), (2) for critical immunization rates with , reduce , reduce and increase recovery rate of isolated individuals (), (3) for low immunization rates (endemic equilibrium) the most important for Beijing and Hong Kong is to increase the isolation rate of infectious individuals () 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 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.