Qualitative study of a quarantine/isolation model with multiple disease stages.

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.

Introduction

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], [5], [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, Newcastle disease and rabies [11], [14]. It is known, however, that quarantine and isolation measures, especially in the context of a new emerging disease, are initially not administered effectively, but are gradually refined (as more data and knowledge of the disease transmission process becomes available (see, for instance, [8])).

Numerous mathematical modeling work have been carried out to assess the impact of quarantine and isolation in combatting the spread of the diseases (such as some of the aforementioned modeling studies for SARS). However, many of the models used for assessing the impact of the quarantine and isolation measures tend to be built based on the assumption that the disease stages are exponentially distributed. However, some recent studies [7], [30] show that it is more realistic to use gamma distribution assumption for the waiting time in the disease stages (rather than exponential distribution assumption). Furthermore, Feng et al. [7] showed that quarantine and isolation models that assume exponential distribution (for the disease stages) may not be suitable for diseases with relatively long latent and/or infectious periods for the case when isolation is not completely effective (i.e., isolated individuals can transmit infection).

The purpose of the current study is to provide a rigorous qualitative analysis of a new deterministic model for transmission dynamics of a communicable disease, subject to the use of quarantine and isolation, where the waiting time in the associated infected classes are assumed to have gamma distribution. The model to be designed extends the SEIQHR model given in [24] by considering multiple stages of the exposed, infectious, quarantined and hospitalized individuals (unlike in [24], it is assumed here that hospitalized individuals do not transmit the infection). Diseases like HIV [25] and influenza [6] are known to have multiple disease (infection) stages.

The paper is organized as follows. The model is formulated in Section 2. The global asymptotic stability of the disease-free equilibrium (DFE) is established in Section 3. The existence of the endemic equilibrium is analyzed in Section 4. Local and global stability proofs for the endemic equilibrium, for special cases, are also provided using a Krasnoselskii sub-linearity trick and a non-linear Lyapunov function of Goh–Voltera type, respectively.

Model formulation

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

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 Π. Susceptible individuals may acquire infection, following effective contact with infectious individuals (in any of the n infectious stages) at a rate λ, where

It is assumed that infected individuals in the classes E , Q (with i  = 1, 2, … ,  m) and H (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), β 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 μ), and increased when recovered individuals lose their infection-acquired immunity (at a rate ψ). Thus, the rate of change of the susceptible population is given by

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

The population of exposed individuals in stage i (with 2 ⩽  i  ⩽  m) is generated by the progression of individuals in stage E into the stage i (at a rate a α). It is decreased by progression to the next exposed stage (at a rate a α), quarantine (at a rate σ ) and natural death (at the rate μ), 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 α). It is decreased by progression to the next infectious stage (I 2; at a rate d 1 κ), hospitalization (at a rate ϕ 1), natural death (at the rate μ) and disease-induced death (at a rate δ 1). This gives

The population of infectious individuals in stage j (with 2 ⩽  j  ⩽  n) is generated by progression of individuals in stage j  − 1 (at a rate d κ). It is decreased by progression to the next infectious stage (at a rate d κ), hospitalization (at a rate ϕ ), natural death (at the rate μ) and disease-induced death (at a rate δ ). Individuals in the final (n) stage of infectiousness recover (at a rate γ 1  =  d κ). Thus,and,

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

Similarly, the population of quarantined individuals in stage i (with 2 ⩽  i  ⩽  m  − 1) is generated by the quarantine of exposed individuals in stage E (at the rate σ ) and the progression of quarantined individuals in stage Q into the stage Q (at a rate b α). It is decreased by progression to the next quarantined stage (at a rate b α) and natural death (at the rate μ). Thus,It should be mentioned that the parameters σ (i  = 1, 2, … ,  m) can be used to model progressive refinement of quarantine measures in the population, by assuming smaller values of σ at the beginning and higher rates for later stages (e.g., for m  = 3, we can assume smaller values for σ 1 and σ 2, but a higher value for σ 3; i.e., σ 1  <  σ 2  <  σ 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 α) 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 κ), natural death (at the rate μ), and disease-induced death (at a rate δ ). Thus,

The population of hospitalized individuals in stage j (with 2 ⩽  j  ⩽  n) is generated by the hospitalization of infectious individuals in stage j (I ) (at the rate ϕ ) and the progression of hospitalized individuals in stage j  − 1 (H ) into the H class (at a rate c κ). It is decreased by the progression to the next hospitalized stage (at a rate c κ), natural death (at the rate μ) and disease-induced death (at a rate δ ). Individuals in the final n stage of hospitalized recover (at a rate γ 2  =  c κ). Thus,and,As in the case for the of quarantine measures discussed above, the parameters ϕ (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 γ 1 and γ 2, respectively). It is decreased by the loss of natural immunity (at the rate ψ) and natural death (at the rate μ), so thatIt should be stated that, in the above formulation, a , b , c , d (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 byIn (2), , , and 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/κ for the hospitalized compartment, H) is shared among the various stages in that compartment. In other words, the time period 1/κ is distributed equally (if c 1  =  c 2  = ⋯ =  c  =  n) or unequally (if c 1  ≠  c 2  ≠ ⋯ ≠  c  ≠  n) between all the H (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/α) and symptomatic (1/κ) compartments. In line with [7], it is assumed that the mean exposed and quarantined periods are the same (1/α) and the mean infectious and hospitalized periods are the same (1/κ).

Let,It follows from (2), (4), using the properties of gamma distribution ([12]; see also Appendix A for a brief description), that the compartments E, I, Q and H indeed have gamma distributions, given, respectively, bywhere the associated exposed, infectious, quarantined and hospitalized periods are given, respectively, by (see also [7], [33])The above formulation ((3), (4)) reduces to that in [7] by setting a  =  b  =  m (for i  = 1, … ,  m) and c  =  d  =  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/α, 1/κ, 1/α and 1/κ, 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, byare 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 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

Fig. 1

Flow diagram of the model (6).

Table 1

Description of variables and parameters of the model (6).

Variable	Description
S(t)	Population of susceptible individuals
Ei(t)	Population of exposed individuals in ith exposed stage (i = 1, … , m)
Ij(t)	Population of infected individuals in jth infectious stage (j = 1, … , n)
Qi(t)	Population of quarantined individuals in ith quarantined stage (i = 1, … , m)
Hj(t)	Population of hospitalized individuals in jth hospitalized stage (j = 1, … , n)
R(t)	Population of recovered individuals
Parameter	Description
Π	Recruitment rate
β	Effective contact rate
djκ	Progression rate from infectious stage j to stage j + 1 (j = 1, … , n)
cjκ	Progression rate from hospitalized stage j to stage j + 1 (j = 1, … , n)
σi	Quarantine rate of exposed individuals in stage i
aiα	Progression rate from exposed stage i to stage i + 1 (i = 1, … , m − 1)
amα	Progression rate of exposed individuals in stage m to first infectious stage
biα	Progression rate from quarantined stage i to stage i + 1 (i = 1, … , m − 1)
bmα	Hospitalization rate of quarantined individuals in stage m
ϕj	Hospitalization rate of infectious individuals in jth infectious stage (j = 1, … , n)
ψ	Rate of loss of infection-acquired immunity
γ1	Recovery rate of infectious individuals in stage n
γ2	Recovery rate of hospitalized individuals in stage n
δj (1 ⩽ j ⩽ n)	Disease-induced death rate of individuals in jth infectious stage
δj (n + 1 ⩽ j ⩽ 2n)	Disease-induced death rate of individuals in (n − j)th hospitalized stage
μ	Natural death rate

including a term for the loss of infection-acquired immunity (at the rate ψ). 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);

including disease-induced death (at rates δ ; 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;

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;

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 ψ  =  δ 1  =  δ 2  = ⋯ =  δ 2  = 0, a 1  =  a 2  = ⋯ =  a  =  b 1  =  b 2  = ⋯ =  b  =  m, c 1  =  c 2  = ⋯ =  c  =  d 1  =  d 2  = ⋯ =  d  =  n, ϕ 1  = ⋯ =  ϕ  =  ϕ and σ 1  = ⋯ =  σ. Also, the model (6) is an extension of the model given in [24] by considering m stages for the exposed (E ; i  = 1, 2, … ,  m) and quarantined (Q ; i  = 1, 2, … ,  m) individuals and n stages for the infectious (I ; j  = 1, 2, … ,  n) and the hospitalized (H ; 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)).

Basic properties

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]):

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.

The closed set

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 byThe next generation operator method [4], [28] will be used to explore the local stability of Ω 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, byand,where, A is 2(m  +  n) ×  m zero matrix, C , B are 2(m  +  n) × (m  +  n), (m  +  n) × (m  +  n) zero matrices, respectively. Furthermore, B is a 2(m  +  n) ×  n matrix, given byThe matrices, A , C and D are (m  +  n) × (m  +  n) are given by and,with,Let, It follows that the control reproduction number [1], [10], denoted by , where ρ is the spectral radius, is given byUsing Theorem 2 in [28], the following result is established.

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

The epidemiological implication of Lemma 2 is that the disease can be eliminated from the community (when ) if the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE (Ω 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 .

Global stability of DFE

The DFE of the model (6), given by (8), is globally-asymptotically stable (GAS) in whenever .

Consider the following Lyapunov function (with the coefficients B, C and D as defined in (9)):with Lyapunov derivative (where a dot represents differentiation with respect to time) given byIt can be shown, after some lengthy algebraic manipulations, thatand,Hence,Since all the parameters of the model (6) and variables are non-negative, it follows that for with if and only if I 1  =  I 2  = ⋯ =  I  = 0. Hence, is a Lyapunov function on . Therefore, by the LaSalle’s Invariance Principle [9],It is clear from (10) that lim sup E 1  = 0. Thus, for sufficiently small ϖ 1  > 0, there exists a constant N 1  > 0 such that lim sup E 1  ⩽  ϖ 1 for all t  >  N 1. It follows from the (m  +  n  + 2)th equation of the model (6) that, for t  >  N 1,Thus, by comparison theorem [26],so that, by letting ϖ 1  → 0,Similarly (by using lim inf E 1  = 0), it can be shown thatThus, it follows from (11), (12) thatHence,Similarly, it can be shown thatThus, by combining (10), (13), (14), it follows that every solution of the equations in the model (6), with initial conditions in , approaches the DFE, Ω 0, as t  → ∞ when . □

Theorem 2 shows that if the use of quarantine and isolation can bring (and keep) the threshold quantity, , to a value less than unity, then the disease will be eliminated from the community (i.e., the condition 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 . All solutions converged to the DFE, Ω 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 Table 2, Table 4 . These parameter values are consistent with those associated with the 2003 SARS outbreaks [2], [5], [22], [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 σ 1 and σ 2, in comparison to σ 3; and also smaller values of ϕ 1 and ϕ 2, in relation to ϕ 3 (see Table 4 ).

Fig. 2

Simulation of the model (6) showing the total number of infected individuals as a function of time for . Parameter values used are as in Table 2, Table 4 with β = 0.2, m = 2, n = 3, a1 = b1 = 1.5, a2 = b2 = 3, c1 = d1 = c2 = d2 = c3 = d3 = 3 (so that, ).

Table 2

Estimated values of the parameters of the model (6).

Parameters	Values (per day)	Sources
β	[0.1, 0.5]	[8], [21]
μ	0.0000351	[13]
κ	0.042553	[2]
δi; i = 1, … , n	0.04227	[17]
δi; i = n + 1, … , 2n	0.027855	[2]
α	0.156986	[5], [22]
ϕn	0.20619	[2]
Π	136	[8]
σm	0.1	[8]
ψ	0.005	Assumed

Table 4

Quarantine and hospitalization rates for various number of disease stages (m and n).

Number of stages	Quarantine rates	Hospitalization rates
m = n = 1	σ1 = 0.1	ϕ1 = 0.20619
m = n = 2	σ1 = 0.05, σ2 = 0.1	ϕ1 = 0.1, ϕ2 = 0.20619
m = n = 3	σ1 = 0.03333, σ2 = 0.05, σ3 = 0.1	ϕ1 = 0.0666, ϕ2 = 0.1, ϕ3 = 0.20619

Existence and stability of endemic equilibria

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.

Existence of endemic equilibrium point (EEP)

Defineto be any arbitrary endemic equilibrium of the model (6). Solving the equations of the model at endemic steady-state givesThe force of infection λ, given by (1), can be expressed at endemic steady-state asAs in [24], the expressions in (15) are re-written in terms of λ ∗∗ S ∗∗, for computational convenience, as below:where,and,Substituting the expressions in (17) into (16) givesDividing each term in (18) by λ ∗∗ S ∗∗ (and noting that λ ∗∗ S ∗∗  ≠ 0 at the endemic steady-state) giveswhere,Hence,The components of the unique endemic equilibrium Ω 1 can then be obtained by substituting the unique value of λ ∗∗, given in (19), into the expressions in (17). This result is summarized below.

The model (6) has a unique endemic (positive) equilibrium, given by Ω 1, whenever .

Global stability of endemic equilibrium for special case

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., ψ  = 0) and the associated disease-induced mortality in all classes is negligible (so that, δ 1  =  δ 2  = ⋯  δ 2  = 0). The model (6), with ψ  =  δ 1  =  δ 2  = ⋯ =  δ 2  = 0, then reduces to:

Adding the equations of the reduced model (20) gives dN/dt  =  Π  −  μN. Hence, N  →  Π/μ as t  → ∞. Thus, Π/μ is an upper bound of N(t) provided that N(0) ⩽  Π/μ. Further, if N(0) >  Π/μ, then N(t) will decrease to this level. Using N  =  Π/μ in (1) gives a limiting (mass action) system given by (20) withIt can be shown that the associated reproduction number of the reduced model, (20) with (21), is given bywhere,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 .

The reduced model, given by (20) with (21), has a unique positive endemic equilibrium whenever .

Furthermore, we claim the following result (see Appendix B for the proof).

The unique endemic equilibrium of the reduced model, given by (20) with (21), is GAS in if .

Simulations for the case when 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 Fig. 4, Fig. 5 are consistent with those reported in [7]. Although the global asymptotic stability result given in Appendix B is for a special case (with ψ  =  δ 1  =  δ 2  = ⋯ =  δ 2  = 0), further extensive numerical simulations suggest that the endemic equilibrium Ω 1, of the full model (6), is GAS in whenever , suggesting the following conjecture.

Fig. 3

Simulation of the model (6) showing the total number of infected individuals as a function of time for . Parameter values used are as in Table 2, Table 4 with β = 0.5, m = 2, n = 3, a1 = b1 = 1.5, a2 = b2 = 3, c1 = d1 = c2 = d2 = c3 = d3 = 3 (so that, ).

Fig. 4

Numerical simulations of the model (6) showing the cumulative number of new infections for various values of the quarantine parameters (σ1 and σ2). Parameter values used are as in Table 2, with β = 0.15, m = 2, n = 3, a1 = b1 = 1.5, a2 = b2 = 3, c1 = d1 = c2 = d2 = c3 = d3 = 3 and isolation rates as given in Table 4.

Fig. 5

Numerical simulations of the model (6) showing the cumulative number of new infections for various values of the isolation parameters (ϕ1, ϕ2 and ϕ3). Parameter values used are as in Table 2, with β = 0.15, m = 2, n = 3, a1 = b1 = 1.5, a2 = b2 = 3, c1 = d1 = c2 = d2 = c3 = d3 = 3 and quarantine rates as given in Table 4.

The unique endemic equilibrium of the model (6), denoted by Ω 1, is GAS in if .

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 . Similarly, the GD2 and GD1 models have and , 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 under-estimated 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).

Fig. 6

Numerical simulations of the model (6) showing the cumulative number of disease-induced mortality for various disease stages (m = n). Parameter values used are as in Table 2, Table 3, Table 4, with β = 0.15.

Fig. 7

Numerical simulations of the model (6) showing the cumulative number of probable SARS for the GTA generated using the GD1, GD2 and ED models. Parameter values used are as in Table 2, Table 3, Table 4, with β = 0.2, ψ = 0. GD1 model: m = n = 3, GD2 model: m = n = 3; σ1 = σ2 = σ3 = 0.1 and ϕ1 = ϕ2 = ϕ3 = 0.20619; ED model: m = n = 1.

Table 3

Values of a, b, c and d (i = 1, 2, 3) for various number of disease stages (m and n).

Number of stages	Values of ai, bi, ci, di
m = n = 1	a1 = b1 = 1, c1 = d1 = 1
m = n = 2	a1 = b1 = 1.5, a2 = b2 = 3, c1 = d1 = 1.5, c2 = d2 = 3
m = n = 3	a1 = b1 = 2, a2 = b2 = 3, a3 = b3 = 6, c1 = d1 = 2, c2 = d2 = 3, c3 = d3 = 6

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 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 σ 1 and σ 2, and also ϕ 1 and ϕ 2, used in the simulations of the GD1 model are different from the quarantine (σ) and isolation (ϕ) rates used in the simulations of the GD2 model. While the values σ 1  = 0.0333, σ 2  = 0.05, σ 3  = 0.1 and ϕ 1  = 0.0666, ϕ 2  = 0.1, ϕ 3  = 0.20619 were used in the simulations of the GD1 model (to account for the gradual refinement of quarantine and isolation), the values σ 1  =  σ 2  =  σ 3  = 0.1 and ϕ 1  =  ϕ 2  =  ϕ 3  = 0.20619 were used in the simulations of the GD2 model (that is why the value for the GD1 model is 0.9710, while that of the GD2 model is 0.7861 for this setting).

Fig. 8

Numerical simulations of the model (6) showing the cumulative number of probable SARS for the Hong Kong generated using the GD1, GD2 and ED models. Parameter values used are as in Table 2, Table 3, Table 4, with β = 0.2, ψ = 0 and Π = 122. GD1 model: m = n = 3, GD2 model: m = n = 3; σ1 = σ2 = σ3 = 0.1 and ϕ1 = ϕ2 = ϕ3 = 0.20619; ED model: m = n = 1.

The effect of the distribution of sojourn times for the symptomatic period (1/κ) 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  =  n  =  c versus the case where d  ≠  n  ≠  c ). 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/α) is also monitored by simulating the model with the parameters in Table 2 for three different scenarios. An asymptomatic period 1/α  = 6 days is chosen, and distributed as follows:

Fig. 9

Numerical simulations of the model (6) showing the cumulative number of new cases for various distributions of the symptomatic period (1/κ) using different values of c1 = d1, c2 = d2, and c3 = d3. Parameter values used are as in Table 2, with β = 0.2, ψ = 0, σ1 = σ2 = σ3 = 0.1 and ϕ1 = ϕ2 = ϕ3 = 0.20619.

2.5 days in E 1 and Q 1 classes (i.e., 1/a 1 α  = 1/b 1 α  = 2.5 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 α  = 1/b 2 α  = 2 days) and 1.5 days in E 3 and Q 3 classes (i.e., 1/a 3 α  = 1/b 3 α  = 1.5 days);

2 days in E 1 and Q 1 classes (i.e., 1/a 1 α  = 1/b 1 α  = 2 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 α  = 1/b 2 α  = 2 days) and 2 days in E 3 and Q 3 classes (i.e., 1/a 3 α  = 1/b 3 α  = 2 days);

1.5 days in E 1 and Q 1 classes (i.e., 1/a 1 α  = 1/b 1 α  = 1.5 days), 2 days in E 2 and Q 2 classes (i.e., 1/a 2 α  = 1/b 2 α  = 2 days) and 2.5 days in E 3 and Q 3 classes (i.e., 1/a 3 α  = 1/b 3 α  = 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.

Fig. 10

Numerical simulations of the model (6) showing the cumulative number of new cases for various distributions of the asymptomatic period (1/α) using different values of a1 = b1, a2 = b2, and a3 = b3. Parameter values used are as in Table 2, with β = 0.2, ψ = 0, σ1 = σ2 = σ3 = 0.1 and ϕ1 = ϕ2 = ϕ3 = 0.20619.

Conclusions

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:

The model (6) has a globally-asymptotically stable disease-free equilibrium whenever the associated reproduction number is less than unity.

The model has a unique endemic equilibrium whenever the reproduction number exceeds unity.

The unique endemic equilibrium of the model is shown to be globally-asymptotically stable for a special case.

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

the cumulative number of new cases of infection decreases with increasing quarantine or isolation rate;

the cumulative number of disease-related mortality increases with increasing number of disease stages (m and n);

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;

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);

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.