key: cord-1050708-ge871ypq
authors: Safi, Mohammad A.; Gumel, Abba B.
title: Qualitative study of a quarantine/isolation model with multiple disease stages
date: 2011-11-01
journal: Appl Math Comput
DOI: 10.1016/j.amc.2011.07.007
sha: b56d89c897e5a628dee42e58d9121832e9371ff2
doc_id: 1050708
cord_uid: ge871ypq

Recent studies suggest that, for disease transmission models with latent and infectious periods, the use of gamma distribution assumption seems to provide a better fit for the associated epidemiological data in comparison to the use of exponential distribution assumption. The objective of this study is to carry out a rigorous mathematical analysis of a communicable disease transmission model with quarantine (of latent cases) and isolation (of symptomatic cases), in which the waiting periods in the infected classes are assumed to have gamma distributions. Rigorous analysis of the model reveals that it has a globally-asymptotically stable disease-free equilibrium whenever its associated reproduction number is less than unity. The model has a unique endemic equilibrium when the threshold quantity exceeds unity. The endemic equilibrium is shown to be locally and globally-asymptotically stable for special cases. Numerical simulations, using data related to the 2003 SARS outbreaks, show that the cumulative number of disease-related mortality increases with increasing number of disease stages. Furthermore, the cumulative number of new cases is higher if the asymptomatic period is distributed such that most of the period is spent in the early stages of the asymptomatic compartments in comparison to the cases where the average time period is equally distributed among the associated stages or if most of the time period is spent in the later (final) stages of the asymptomatic compartments. Finally, it is shown that distributing the average sojourn time in the infectious (asymptomatic) classes equally or unequally does not effect the cumulative number of new cases.

Recent studies suggest that, for disease transmission models with latent and infectious periods, the use of gamma distribution assumption seems to provide a better fit for the associated epidemiological data in comparison to the use of exponential distribution assumption. The objective of this study is to carry out a rigorous mathematical analysis of a communicable disease transmission model with quarantine (of latent cases) and isolation (of symptomatic cases), in which the waiting periods in the infected classes are assumed to have gamma distributions. Rigorous analysis of the model reveals that it has a globally-asymptotically stable disease-free equilibrium whenever its associated reproduction number is less than unity. The model has a unique endemic equilibrium when the threshold quantity exceeds unity. The endemic equilibrium is shown to be locally and globally-asymptotically stable for special cases. Numerical simulations, using data related to the 2003 SARS outbreaks, show that the cumulative number of disease-related mortality increases with increasing number of disease stages. Furthermore, the cumulative number of new cases is higher if the asymptomatic period is distributed such that most of the period is spent in the early stages of the asymptomatic compartments in comparison to the cases where the average time period is equally distributed among the associated stages or if most of the time period is spent in the later (final) stages of the asymptomatic compartments. Finally, it is shown that distributing the average sojourn time in the infectious (asymptomatic) classes equally or unequally does not effect the cumulative number of new cases.

Ó 2011 Elsevier Inc. All rights reserved.

Since the pioneering works of Sir Ronald Ross, Kermack and McKendrick (see, for instance, [15, 16, 23] ), numerous mathematical models have been designed and used to gain insight into transmission dynamics of emerging and re-emerging diseases of public health interest. The models, typically of the forms of deterministic or stochastic systems of non-linear differential equations, are used to evaluate various control strategies such as: vaccination, the use of antibiotics or antivirals, quarantine, isolation, etc. Of the aforementioned control strategies, the use of quarantine (of individuals suspected of being exposed to the disease) and isolation (of those with clinical symptoms of the disease) are the most commonly used (since the beginning of recorded human history). These measures have been used in the control of numerous diseases such as leprosy, plague, cholera, typhus, yellow fever, smallpox, diphtheria, tuberculosis, measles, ebola, pandemic influenza and, more recently, severe acute respiratory syndrome (SARS) [3, 11, [18] [19] [20] 22, 29, 31, 32] . Furthermore, quarantine and isolation are popularly applied to combat the spread of animal diseases such as bovine tuberculosis, rinderpest, foot-and-mouth, psittacosis, 0096-3003/$ -see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.07.007

The total population at time t, denoted by N(t), is sub-divided into six disjoint classes of susceptible (S(t)), exposed (E(t); with m exposed stages), quarantined (Q(t); with m quarantined stages), infectious (I(t); with n infectious stages), hospitalized (H(t); with n hospitalized stages) and recovered (R(t)) individuals, so that

H j ðtÞ þ RðtÞ:

In this paper, unlike in [18] , it is assumed that the fraction of infected contacts that can be traced and quarantined at the time of infection is very small. Furthermore, it is assumed that the total population is large in comparison to the size of the infected individuals (N ) E + I + Q + H + R). Consequently, the quarantine of susceptible individuals (feared exposed to the disease) is unlikely to have a significant impact on the disease transmission dynamics. Hence, the quarantine of susceptible individuals is not considered in this study (see also [7] ). In other words, in this study, quarantine refers to the isolation of exposed (latently-infected) individuals only.

The susceptible population is increased by the recruitment of individuals into the community (assumed susceptible), at a rate P. Susceptible individuals may acquire infection, following effective contact with infectious individuals (in any of the n infectious stages) at a rate k, where

It is assumed that infected individuals in the classes E i , Q i (with i = 1,2,. . . , m) and H j (with j = 1,2,. . . , n) do not transmit infection (i.e., it is assumed that exposed individuals do not transmit infection, and that quarantine and isolation measures are implemented in a perfect manner). Although some of these assumptions may not be entirely realistic in some epidemiological settings, such as in the transmission dynamics of influenza (where transmission by infected individuals without disease symptoms occurs), they help in making the mathematical analysis of the resulting large system of non-linear differential equations more tractable. Further, in (1), b is the effective contact rate (contact capable of leading to infection).

The population of susceptible individuals is further decreased by natural death (at a rate l), and increased when recovered individuals lose their infection-acquired immunity (at a rate w). Thus, the rate of change of the susceptible population is given by dS dt ¼ P þ wR À kS À lS:

The population of exposed individuals in stage 1 (E 1 ) is generated by the infection of susceptible individuals (at the rate k). This population is decreased by progression to the next exposed stage (E 2 ; at a rate a 1 a), quarantine (at a rate r 1 ) and natural death (at the rate l), so that

The population of exposed individuals in stage i (with 2 6 i 6 m) is generated by the progression of individuals in stage E iÀ1 into the stage i (at a rate a iÀ1 a). It is decreased by progression to the next exposed stage (at a rate a i a), quarantine (at a rate r i ) and natural death (at the rate l), so that

The population of infectious individuals in stage 1 is generated when exposed individuals in the final (m) stage develop symptoms (at the rate a m a). It is decreased by progression to the next infectious stage (I 2 ; at a rate d 1 j), hospitalization (at a rate / 1 ), natural death (at the rate l) and disease-induced death (at a rate d 1 ). This gives

The population of infectious individuals in stage j (with 2 6 j 6 n) is generated by progression of individuals in stage j À 1 (at a rate d jÀ1 j). It is decreased by progression to the next infectious stage (at a rate d j j), hospitalization (at a rate / j ), natural death (at the rate l) and disease-induced death (at a rate d j ). Individuals in the final (n) stage of infectiousness recover (at a rate c 1 = d n j). Thus,

Exposed individuals in stage 1 are quarantined at the rate r 1 . The population of quarantined individuals in stage 1 is decreased by progression to the next quarantined stage (at a rate b 1 a) and natural death (at the rate l). Thus,

Similarly, the population of quarantined individuals in stage i (with 2 6 i 6 m À 1) is generated by the quarantine of exposed individuals in stage E i (at the rate r i ) and the progression of quarantined individuals in stage Q iÀ1 into the stage Q i (at a rate b iÀ1 a). It is decreased by progression to the next quarantined stage (at a rate b i a) and natural death (at the rate l). Thus,

It should be mentioned that the parameters r i (i = 1,2,. . . , m) can be used to model progressive refinement of quarantine measures in the population, by assuming smaller values of r i at the beginning and higher rates for later stages (e.g., for m = 3, we can assume smaller values for r 1 and r 2 , but a higher value for r 3 ; i.e., r 1 < r 2 < r 3 ).

The population of hospitalized individuals in stage 1 is generated by the hospitalization of quarantined individuals in the final stage (m; at the rate b m a) and infectious individuals in stage 1 (at the rate / 1 ). It is decreased by progression to the next hospitalized stage (at a rate c 1 j), natural death (at the rate l), and disease-induced death (at a rate d n+1 ). Thus,

The population of hospitalized individuals in stage j (with 2 6 j 6 n) is generated by the hospitalization of infectious individuals in stage j (I j ) (at the rate / j ) and the progression of hospitalized individuals in stage j À 1 (H jÀ1 ) into the H j class (at a rate c jÀ1 j). It is decreased by the progression to the next hospitalized stage (at a rate c j j), natural death (at the rate l) and disease-induced death (at a rate d n+j ). Individuals in the final n stage of hospitalized recover (at a rate c 2 = c n j). Thus, dH j dt

. . . ; n À 1; and, dH n dt ¼ / n I n þ c nÀ1 jH nÀ1 À ðc 2 þ l þ d 2n ÞH n :

As in the case for the of quarantine measures discussed above, the parameters / i (i = 1,. . . , n) can also be used to model the progressive refinement of isolation (in hospital; so that, for n = 3, we can have / 1 < / 2 < / 3 ). Finally, the population of recovered individuals is generated by the recovery of non-hospitalized and hospitalized infectious individuals in the final n stage (at the rates c 1 and c 2 , respectively). It is decreased by the loss of natural immunity (at the rate w) and natural death (at the rate l), so that dR dt ¼ c 1 I n þ c 2 H n À ðw þ lÞR:

It should be stated that, in the above formulation, a i , b i , c j , d j (i = 1,2,. . . , m; j = 1,2,. . . , n) are constants. Furthermore, it is assumed that the distributions of exposed, quarantined, infectious and hospitalized periods are exponential, given by

In (2), T E i ¼ 1=a i a, T I j ¼ 1=d j j, T Q i ¼ 1=b i a and T H j ¼ 1=c j j are the mean exposed, quarantined, infectious and hospitalized periods, respectively. The relations in (2) are such that:

That is, the respective mean time spent in a given infected compartment (e.g., 1/j for the hospitalized compartment, H) is shared among the various stages in that compartment. In other words, the time period 1/j is distributed equally (if c 1 = c 2 = Á Á Á = c n = n) or unequally (if c 1 -c 2 -Á Á Ác n -n) between all the H j (j = 1,2,. . . , n) stages. Hence, this formulation extends the formulation in [7] , where these periods are equally distributed among the relevant stages (for all the infected compartments, E, Q, I, H), by allowing for equal or unequal distribution of the sojourn times in asymptomatic (1/a) and symptomatic (1/j) compartments. In line with [7] , it is assumed that the mean exposed and quarantined periods are the same (1/ a) and the mean infectious and hospitalized periods are the same (1/j). Let,

It follows from (2) and (4) where the associated exposed, infectious, quarantined and hospitalized periods are given, respectively, by (see also [7, 33] )

The above formulation ( (3) and (4)) reduces to that in [7] by setting a i = b i = m (for i = 1,. . . , m) and c j = d j = n (for j = 1,. . . , n).

In other words, it should be emphasized that the main distinction between the formulation in the current study and that in [7] is that, here, it is assumed that the sojourn periods in each of the four compartments, E, I, Q, and H, given by 1/a, 1/j, 1/a and 1/j, respectively, are distributed (not necessarily equally) among the various sub stages (whereas, these periods are distributed equally at each related stage in [7] ). Eichner et al. [6] considered 9 latent and 19 infectious stages to model the transmission dynamics of pandemic influenza. It is worth stating that although the sums defined in (4) are gamma distributed, the actual (true) total number of infected individuals, E true , I true , Q true and H true , given, respectively, by

are not necessarily gamma distributed. However, the different sums in (4) have the same means, with their respective sums given in (5), but different variances. Thus, putting all these formulations and assumptions together, it follows that the model for the transmission dynamics of an infectious disease in the presence of exposed, quarantine, infectious and isolation periods, subject to gamma distributed Table 1 Description of variables and parameters of the model (6) .

Population of susceptible individuals E i (t) Population of exposed individuals in ith exposed stage (i = 1,. . ., m) I j (t)

Population of infected individuals in jth infectious stage (j = 1,. . ., n)

Population of hospitalized individuals in jth hospitalized stage (j = 1,. . ., n) R(t)

Population of recovered individuals Parameter Description P

Effective contact rate

Progression rate from infectious stage j to stage j + 1 (j = 1,. . ., n) c j j Progression rate from hospitalized stage j to stage j + 1 (j = 1,. . ., n) r i

Quarantine rate of exposed individuals in stage i a i a Progression rate from exposed stage i to stage i + 1 (i = 1,. . ., m À 1)

Progression rate of exposed individuals in stage m to first infectious stage Recovery rate of hospitalized individuals in stage n d j (1 6 j 6 n)

Disease-induced death rate of individuals in jth infectious stage d j (n + 1 6 j 6 2n)

Disease-induced death rate of individuals in (n À j)th hospitalized stage l Natural death rate sojourn periods, is given by the following non-linear system of differential equations (a flow diagram of the model is given in Fig. 1 ; and the associated variables and parameters are described in Table 1) :

The model (6) extends the multi-stage model given in [7] by (i) including a term for the loss of infection-acquired immunity (at the rate w). Although the numerical simulations to be carried out in this study are largely based on the 2003 SARS outbreaks (which was a single season epidemic), the model (6) is robust enough to enable the assessment of the transmission dynamics of any arbitrary disease where the infection-acquired immunity is lost either during a single season or in multiple seasons (such as the case of influenza, malaria, and some childhood diseases); (ii) including disease-induced death (at rates d i ; i = 1,2,. . . , 2n). Most diseases, such as HIV, malaria, influenza, TB, etc., have significant disease-induced mortality. Hence, it is crucial that this feature be incorporated in modeling studies; (iii) assuming the average sojourn periods in the exposed, quarantined, infectious and hospitalized classes are distributed (not necessarily equally) among the various stages (these periods are assumed to be equally distributed among each of the aforementioned four infected compartments in [7] ). Although, to our knowledge, there is no definitive epidemiological data to suggest that these periods are equally or unequally distributed, the model (6) is general enough to allow for the assessment of each of the two cases; (iv) assuming varied rates of quarantine and isolation in each quarantine and isolation stage (same rates are used in [7] in all quarantine and isolation stages). This assumption allows for the assessment of progressive refinement of quarantine and isolation measures (this was evident during 2003 SARS outbreaks [8, 17] ).

The model (6) is denoted by GD1 for comparison purposes. It is worth emphasizing that the model (6) reduces to the model in [7] by setting (6) is an extension of the model given in [24] by considering m stages for the exposed (E i ; i = 1,2,. . . , m) and quarantined (Q i ; i = 1,2,. . . , m) individuals and n stages for the infectious (I j ; j = 1,2,. . . , n) and the hospitalized (H j ; j = 1,2,. . . , n) individuals (i.e., the model (6) reduces to the model in [24] by setting n = m = 1, taking into account the assumption that hospitalized individuals do not transmit infection; this assumption is relaxed in [24] ).

In addition to formulating the model in terms of gamma-distributed waiting times for the associated disease stages, this study contributes by way of carrying out a detailed rigorous mathematical analysis of the model (6) . In particular, global asymptotic stability results for the equilibria of the model will be proven (under certain conditions). Furthermore, the model (6) is used to evaluate the impact of the use of quarantine and isolation in combatting the spread of a given communicable disease (such as SARS). This study offers not only important extensions to the model presented in [7] , it also contributes by extending some of the mathematical results presented in [7] (particularly global stability proof of the associated endemic equilibrium of the extended model (6)).

Since the model (6) monitors human populations, all its associated parameters are non-negative. Further, the following basic results can be easily established (see, for instance, [24, 27] ): Theorem 1. The state variables of the model (6) are non-negative for all time. In other words, solutions of the model system (6) with positive initial data will remain positive for all time t > 0.

3. Local stability of disease-free equilibrium (DFE)

The model (6) has a DFE, obtained by setting the right-hand sides of the equations to zero, given by is given by

The next generation operator method [4, 28] will be used to explore the local stability of X 0 . Using the notation in [28] , the non-negative matrix, F, of the new infection terms, and the M-matrix, V, of the transition terms associated with the model (6) , are given, respectively, by

where, A F is 2(m + n) Â m zero matrix, C F , B V are 2(m + n) Â (m + n), (m + n) Â (m + n) zero matrices, respectively. Furthermore, B F is a 2(m + n) Â n matrix, given by

The matrices, A V , C V and D V are (m + n) Â (m + n) are given by

2m þ n þ 1 6 j 6 2ðm þ nÞ:

. . . ; m þ n À 1; B mþn ¼ 1;

. . . ; n; q ¼ n þ 1 À p; C 1;n ¼ 1 and C 2;nÀ1 ¼ jd nÀ1 þ k mþn :

It follows that the control reproduction number [1, 10] , denoted by R c ¼ qðFV À1 Þ, where q is the spectral radius, is given by R c ¼ bD 1 C n;1 k mþn :

Using Theorem 2 in [28] , the following result is established.

The DFE of the model (6), given by (8), is locally-asymptotically stable (LAS) if R c < 1, and unstable if R c > 1.

The epidemiological implication of Lemma 2 is that the disease can be eliminated from the community (when R c < 1) if the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE (X 0 ). For disease elimination to be independent of the initial sizes of sub-populations, the global asymptotic stability of the DFE must be established for R c < 1.

Theorem 2. The DFE of the model (6), given by (8), is globally-asymptotically stable (GAS) in D whenever R c 6 1.

Proof. Consider the following Lyapunov function (with the coefficients B, C and D as defined in (9)):

Lyapunov derivative (where a dot represents differentiation with respect to time) given by

C nÀjþ1;j B mþj k mþj I j À k mþn I n ;

¼ k mþn R c X n j¼1 I j þ X m j¼2 C n;1 ðD jþ1 a j a À D j k j ÞE j þ X nÀ1 j¼1 d j jC nÀj;jþ1 B mþjþ1 À k mþj C nÀjþ1;j B mþj I j À k mþn I n :

It can be shown, after some lengthy algebraic manipulations, that D jþ1 a j a À D j k j ¼ 0;

and, d j jC nÀj;jþ1 B mþjþ1 À k mþj C nÀjþ1;j B mþj ¼ Àk mþn :

Hence, _ F 6 k mþn ðR c À 1Þ X n j¼1 I j 6 0 for R c 6 1:

Since all the parameters of the model (6) and variables are non-negative, it follows that _ F 6 0 for R c 6 1 with _ F ¼ 0 if and only if I 1 = I 2 = Á Á Á = I n = 0. Hence, F is a Lyapunov function on D. Therefore, by the LaSalle's Invariance Principle 

It is clear from (10) that lim sup t?1 E 1 = 0. Thus, for sufficiently small -1 > 0, there exists a constant N 1 > 0 such that lim sup t?1 E 1 6 -1 for all t > N 1 . It follows from the (m + n + 2)th equation of the model (6) that, for t > N 1 , _ Q 1 6 r 1 -1 À k mþnþ1 Q 1 :

Thus, by comparison theorem [26] ,

so that, by letting -1 ? 0,

Similarly (by using lim inf t?1 E 1 = 0), it can be shown that

Thus, it follows from (11) and (12) 

Thus, by combining (10), (13) and (14), it follows that every solution of the equations in the model (6) Theorem 2 shows that if the use of quarantine and isolation can bring (and keep) the threshold quantity, R c , to a value less than unity, then the disease will be eliminated from the community (i.e., the condition R c < 1 is necessary and sufficient for disease elimination). Fig. 2 depicts numerical results obtained by simulating the model (6), with m = 2 and n = 3, using various initial conditions for the case R c < 1. All solutions converged to the DFE, X 0 , (in line with Theorem 2). It should be mentioned that, unless otherwise stated, simulations of the model (6) are carried out using the parameter values in Tables  2 and 4 . These parameter values are consistent with those associated with the 2003 SARS outbreaks [2, 5, 8, 17] . It is worth mentioning that the progressive refinement of quarantine and isolation measures is incorporated in all numerical simulations in this study (unless otherwise stated) by using smaller values of r 1 and r 2 , in comparison to r 3 ; and also smaller values of / 1 and / 2 , in relation to / 3 (see Table 4 ).

In this section, the possible existence and stability of endemic (positive) equilibria of the model (6) (i.e., equilibria where at least one of the infected components of the model is non-zero) will be explored. 

. . . ; m;

The force of infection k, given by (1), can be expressed at endemic steady-state as

As in [24] , the expressions in (15) are re-written in terms of k ⁄⁄ S ⁄⁄ , for computational convenience, as below: Table 3 Values of a i , b i , c i and d i (i = 1, 2, 3) for various number of disease stages (m and n). Table 4 Quarantine and hospitalization rates for various number of disease stages (m and n). 

. . . ; m;

. . . ; n;

where,

and,

Substituting the expressions in (17) into (16) gives

Dividing each term in (18) by k ⁄⁄ S ⁄⁄ (and noting that k ⁄⁄ S ⁄⁄ -0 at the endemic steady-state) gives

The components of the unique endemic equilibrium X 1 can then be obtained by substituting the unique value of k ⁄⁄ , given in (19) , into the expressions in (17) . This result is summarized below.

Lemma 3. The model (6) has a unique endemic (positive) equilibrium, given by X 1 , whenever R c > 1.

Here, the global stability of the endemic equilibrium of the model (6) is given for the special case where the recovered individuals do not lose their infection-acquired immunity (i.e., w = 0) and the associated disease-induced mortality in all classes is negligible (so that, d 1 = d 2 = Á Á Á d 2n = 0). The model (6), with w = d 1 = d 2 = Á Á Á = d 2n = 0, then reduces to: dS dt ¼ P À kS À lS;

Adding the equations of the reduced model (20) gives dN/dt = P À lN. Hence, N ? P/l as t ? 1. Thus, P/l is an upper bound of N(t) provided that N(0) 6 P/l. Further, if N(0) > P/l, then N(t) will decrease to this level. Using N = P/l in (1) gives a limiting (mass action) system given by (20) with

It can be shown that the associated reproduction number of the reduced model, (20) with (21), is given by

It is easy to show, using the technique in Section 4.1, that the reduced model, given by (20) with (21), has a unique EEP whenever R cr > 1.

Lemma 4. The reduced model, given by (20) with (21), has a unique positive endemic equilibrium whenever R cr > 1. Furthermore, we claim the following result (see Appendix B for the proof).

Theorem 3. The unique endemic equilibrium of the reduced model, given by (20) with (21), is GAS in D n D 0 if R cr > 1.

Simulations for the case when R c > 1 are depicted in Fig. 3 , showing convergence of the solutions to the endemic equilibrium (in line with Theorem 3). Fig. 4 depicts the cumulative number of new infections as a function of quarantine rates, from which it is evident that the cumulative number of new infections decreases with increasing quarantine rate. Similar result is obtained by increasing the isolation rate (Fig. 5) . It should be mentioned that the simulation results in Figs. 4 and 5 are consistent with those reported in [7] . Although the global asymptotic stability result given in Appendix B is for a special case (with w = d 1 = d 2 = Á Á Á = d 2n = 0), further extensive numerical simulations suggest that the endemic equilibrium X 1 , of the full model (6), is GAS in D n D 0 whenever R c > 1, suggesting the following conjecture.

Conjecture. The unique endemic equilibrium of the model (6), denoted by X 1 , is GAS in D n D 0 if R c > 1.

The effect of the number of disease stages for the exposed (m) and infectious (n) classes is monitored by simulating the model (6) with various values of m = n. The results obtained, depicted in Fig. 6 , show an increase in the cumulative number of disease-related mortality with increasing values of m = n. Simulations for the cumulative number of probable SARS cases observed during the 2003 outbreaks in the Greater Toronto Area (GTA) of Canada are also carried out. The results obtained, for the case m = n = 3, are compared with those obtained using the exponentially-distributed (ED) equivalent of the model (6) (i.e., model (6) with m = n = 1) and another gamma-distributed version of the model (6) with m = n = 3, denoted by GD2, where the average sojourn time in each of the exposed, quarantined, hospitalized and infectious stages is shared equally among each associated disease stage (this is similar to the model given in [7] ). It should be mentioned that, in such a setting, the standard ED model has the associated reproduction number given by R c ¼ 0:6506. Similarly, the GD2 and GD1 models have R c ¼ 0:6962 and R c ¼ 0:9858, respectively. Furthermore, about 250 probable SARS cases were reported for the GTA (see Fig. 2 in [8] ). The simulation results obtained, depicted in Fig. 7 , show that while the ED and GD2 models underestimated the observed number of probable cases, the GD1 model (6) gave a very good estimate of the observed data. It should be mentioned that the GD2 model is also competitive if the quarantine and isolation rates are distributed (unequally) to incorporate their progressive refinement (as in the case of the model GD1).

Similar comparison are made for the cumulative number of cases recorded for the Hong Kong SARS outbreaks (approximately 1750 cases were recorded in Hong Kong [8] ). Here, too, the GD1 model is more competitive (Fig. 8) . For these simulations, the ED, GD1 and GD2 models have R c given by 0.7345, 0.9710 and 0.7861, respectively. It should be emphasized, however, that the reason why the GD1 model gives different results, compared to the GD2 model (for instance), is that the values of r 1 and r 2 , and also / 1 and / 2 , used in the simulations of the GD1 model are different from the quarantine (r) and isolation (/) rates used in the simulations of the GD2 model. While the values r 1 = 0.0333, r 2 = 0.05, r 3 = 0.1 and / 1 = 0.0666, Table 2 , with b = 0.15, m = 2, n = 3, a 1 = b 1 = 1.5, a 2 = b 2 = 3, c 1 = d 1 = c 2 = d 2 = c 3 = d 3 = 3 and isolation rates as given in Table 4 . The effect of the distribution of sojourn times for the symptomatic period (1/j) is monitored by simulating the GD1 model (6) with the parameters in Table 2 for the case where the periods are either same or varied in each stage (i.e., the case where d j = n = c j versus the case where d j -nc j ). In both cases, the same numerical simulation results were obtained (Fig. 9) . In other words, distributing the average sojourn times equally or unequally between the sub stages of the symptomatic classes (I and H) does not alter the numerical simulation results obtained. The effect of the distribution of sojourn times in the asymptomatic classes (E and Q; given by 1/a) is also monitored by simulating the model with the parameters in Table 2 for three different scenarios. An asymptomatic period 1/a = 6 days is chosen, and distributed as follows:

(I) 2.5 days in E 1 and Q 1 classes (i.e., 1/a 1 a = 1/b 1 a = 2.5 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 a = 1/b 2 a = 2 days) and

1.5 days in E 3 and Q 3 classes (i.e., 1/a 3 a = 1/b 3 a = 1.5 days); Table 2 , with b = 0.15, m = 2, n = 3, a 1 = b 1 = 1.5, a 2 = b 2 = 3, c 1 = d 1 = c 2 = d 2 = c 3 = d 3 = 3 and quarantine rates as given in Table 4 . (II) 2 days in E 1 and Q 1 classes (i.e., 1/a 1 a = 1/b 1 a = 2 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 a = 1/b 2 a = 2 days) and 2 days in E 3 and Q 3 classes (i.e., 1/a 3 a = 1/b 3 a = 2 days); (III) 1.5 days in E 1 and Q 1 classes (i.e., 1/a 1 a = 1/b 1 a = 1.5 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 a = 1/b 2 a = 2 days) and 2.5 days in E 3 and Q 3 classes (i.e., 1/a 3 a = 1/b 3 a = 2.5 days).

The simulation results obtained (Fig. 10 ) clearly show that if the asymptomatic period is distributed such that more time is spent in the early stages of the asymptomatic (latent and quarantine) classes (i.e., more time is spent in E 1 , E 2 , Q 1 , Q 2 classes in comparison to in E 3 and Q 3 classes), the cumulative number of new cases is higher than for the cases where the asymptomatic period is distributed equally among the stages or if more time is spent in the later asymptomatic stages. In other words, unlike for the case of the sojourn time spent in the symptomatic classes (I and H), the way the sojourn time is distributed in the asymptomatic compartments (E and Q) affects the cumulative number of new cases. 

A new deterministic model for disease transmission, subject to the use of quarantine (of asymptomatic cases) and isolation (of individuals with disease symptoms), is presented and rigorously analyzed. The model, which is based on the assumption that the mean waiting periods in all infected classes obey a gamma distribution, adopts a standard incidence formulation for the infection rate. An important feature of this model is that it allows for equal or unequal distribution of the sojourn time in each of the associated infected compartment. Furthermore, it allows for the gradual refinement of quarantine and isolation measures (this was the case during the 2003 SARS outbreaks). The main theoretical findings of the study are given below:

(i) The model (6) has a globally-asymptotically stable disease-free equilibrium whenever the associated reproduction number ðR c Þ is less than unity. (ii) The model has a unique endemic equilibrium whenever the reproduction number exceeds unity. (iii) The unique endemic equilibrium of the model is shown to be globally-asymptotically stable for a special case. Table 2 , with b = 0.2, w = 0, r 1 = r 2 = r 3 = 0.1 and / 1 = / 2 = / 3 = 0.20619.

Numerical simulations of the model (6), using data related to the 2003 SARS outbreaks, show the following:

(a) the cumulative number of new cases of infection decreases with increasing quarantine or isolation rate; (b) the cumulative number of disease-related mortality increases with increasing number of disease stages (m and n); (c) unlike the ED and GD2 models, the model (6) gives numerical results that are consistent with the 2003 SARS outbreaks data for the GTA and Hong Kong; (d) distributing the average sojourn time equally or unequally between the respective symptomatic classes does not alter the numerical simulation result obtained (i.e., the cumulative number of new cases); (e) if the asymptomatic period is distributed such that more time is spent in the early asymptomatic (latent and quarantine) stages, the cumulative number of new cases is higher than for the cases where the period is distributed equally among the asymptomatic stages or if more time is spent in the later asymptomatic stages.

One of the authors (ABG) acknowledges, with thanks, the support in part of the Natural Science and Engineering Research Council (NSERC) and Mathematics of Information Technology and Complex Systems (MITACS) of Canada. MAS gratefully acknowledges the support of the Institute of Industrial Mathematical Sciences and the University of Manitoba Graduate Fellowship. The authors are grateful to Alex Leblanc for helpful discussions on the formulation of the gamma distribution.

Appendix A. Properties of gamma distribution [12] A random variable X that is gamma-distributed with scale h and shape k is denoted by X $ Cðk; hÞ or X $ Gammaðk; hÞ:

Properties:

(i) Summation: if X i has a C(k i , h) distribution for i = 1,2,. . . , N, then

(ii) Scaling: if X $ C(k, h) then for any a > 0; aX $ C k; h a À Á . Proof. Consider the reduced model, given by (20) with (21) . Let R cr > 1, so that the associated unique endemic equilibrium exists. Further, consider the following non-linear Lyapunov function:

y i I j À I ÃÃ j À I ÃÃ j ln

where the coefficients x i (i = 1,. . . , m) and y j (j = 1,. . ., n) are positive constants to be determined. Thus,

Substituting the expressions of the derivatives from the system (20), using (21) , gives (note that the relation P ¼ b 1 S ÃÃ P n j¼1 I ÃÃ j þ lS ÃÃ , at endemic steady-state, has been used) 

The coefficients x i (i = 2,. . ., m) and y j (j = 1,. . . , n) are chosen such that

x 2 a 1 a À k 1 ¼ 0;

x iþ1 a i a À x i f i ¼ 0; for i ¼ 3; . . . ; m À 1; y 1 a m a À x m f m ¼ 0; b 1 S ÃÃ þ y jþ1 d j j À y j f mþj ¼ 0;

b 1 S ÃÃ À y n f mþn ¼ 0;

ð23Þ so that, from (23),

f l a l a ; i ¼ 2; . . . ; m; 

Thus, by combining (31), (34) and (35), it follows that every solution to the equations of the reduced model, with initial condition in D n D 0 , approaches the unique endemic equilibrium of the reduced model (20) with (21) as t ? 1 for R cr > 1 and

Population Biology of Infectious Diseases

Model parameters and outbreak control for SARS

The basic reproductive number of ebola and the effects of public health measures: the cases of Congo and Uganda

On the definition and computation of the basic reproduction ratio R 0 in models for infectious disease in heterogeneous population

Epidemiological determinants of spread of casual agent of severe acute respiratory syndrome in Hong Kong

The influenza pandemic preparedness planning tool InfluSim

Epidemiological models with non-exponentially distributed disease stages and application to disease control

Gumel et al, Modelling strategies for controlling SARS outbreaks

Ordinary Differential Equations

The mathematics of infectious diseases

Effects of quarantine in six endemic models for infectious diseases

Introduction to Mathematical Statistics

Hong Kong in Figs

Quarantine-based disease control in domesticated animal herds

A contribution to the mathematical theory of epidemics

A contribution to the mathematical theory of epidemics, Part II

The epidemiology of severe acute respiratory syndrome in the 2003 Hong Kong epidemic: an analysis of all 1755 patients

Transmission dynamics and control of severe acute respiratory syndrome

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

Sensitivity and uncertainty analyses for a SARS model with time-varying inputs and outputs

Mathematical study of the impact of quarantine, isolation and vaccination in curtailing an epidemic

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

The Prevention of Malaria

Global asymptotic dynamics of a model for quarantine and isolation

Role of incidence function in vaccine-induced backward bifurcation in some HIV models

The Theory of the Chemostat

Mathematics in Population Biology

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

Simulating the SARS outbreak in Beijing with limited data

Appropriate models for the management of infectious diseases

Critical role of nosocomial transmission in the Toronto SARS outbreak

Optimal and sub-optimal quarantine and isolation control in SARS epidemics

Analysis of a model with multiple infectious stages and arbitrary distributed stage durations

R cr > 1. Hence, F is a Lyapunov function for the sub-system of the model (20) consisting of the equations for S, E i (i = 1,. . . , m), I j (j = 1,. . . , n) of the model (20) on D n D 0 . Therefore, by the LaSalle's Invariance Principle [9] , lim t!1 so that, by letting -? 0,Similarly (by using lim inf t!1 E 1 ¼ E ÃÃ 1 ), it can be shown thatThus, it follows from (32) and (33)