An Application of Queuing Theory to Sis and Seis Epidemic Models

In this work we consider every individual of a population to be a server whose state can be either busy (infected) or idle (susceptible). This server approach allows to consider a general distribution for the duration of the infectious state, instead of being restricted to exponential distributions. In order to achieve this we first derive new approximations to quasistation-ary distribution (QSD) of SIS (Susceptible-Infected-Susceptible) and SEIS (Susceptible-Latent-Infected-Susceptible) stochastic epidemic models. We give an expression that relates the basic reproductive number, R 0 and the server utilization, ρ. 809 810 HERN´ANDEZ, CASTILLO-CHAVEZ, MONTESINOS AND HERN´ANDEZ 1. Introduction. Queuing theory deals with the analysis of serving customers arriving to a facility with a fixed number of servers. Basic queuing models include the inter-arrival distribution, the number of servers and the distribution of the service time, once an individual reaches a server. The notation M/G/n is used to indicate that customers arrive according to a Poisson process, that the service time once in the server follows a general distribution and that there are n servers. Some important questions in queuing theory relate to the average time of a customer in the system, the average time to service, the distribution of the number of customers in the system and so forth. One important parameter in queuing systems is the server utilization, which is the long run proportion of time a server is busy. For an introduction to queuing models see [1]. Epidemic models deal with the dynamics of a population where individuals can be classified in one of several possible states, each describing some type or level of illness. The main use of epidemic models is to find out how the different parameters relate to the progression of the epidemics, and thus, they become a valuable tool to control disease. One of the most important parameters in epidemic models is the basic reproduction number or R 0 , which is basically a measure of the infection potential of the disease. If R 0 ≤ 1 then the disease will fade out eventually with probability one, so, R 0 becomes an important threhsold parameter. There is relatively little work relating queuing theory with epidemiology. Kendall [2] and [3] discusses the relationship between M/G/1 queues and birth-death processes. Kitaev [4] works on the relation between birth and death processes and the M/G/1 queues with processor sharing, and the SIR epidemic model constructions used in [5] and …


(Communicated by Zhilan Feng)
Abstract.In this work we consider every individual of a population to be a server whose state can be either busy (infected) or idle (susceptible).This server approach allows to consider a general distribution for the duration of the infectious state, instead of being restricted to exponential distributions.In order to achieve this we first derive new approximations to quasistationary distribution (QSD) of SIS (Susceptible-Infected-Susceptible) and SEIS (Susceptible-Latent-Infected-Susceptible) stochastic epidemic models.We give an expression that relates the basic reproductive number, R 0 and the server utilization, ρ.
1. Introduction.Queuing theory deals with the analysis of serving customers arriving to a facility with a fixed number of servers.Basic queuing models include the inter-arrival distribution, the number of servers and the distribution of the service time, once an individual reaches a server.The notation M/G/n is used to indicate that customers arrive according to a Poisson process, that the service time once in the server follows a general distribution and that there are n servers.Some important questions in queuing theory relate to the average time of a customer in the system, the average time to service, the distribution of the number of customers in the system and so forth.One important parameter in queuing systems is the server utilization, which is the long run proportion of time a server is busy.For an introduction to queuing models see [1].
Epidemic models deal with the dynamics of a population where individuals can be classified in one of several possible states, each describing some type or level of illness.The main use of epidemic models is to find out how the different parameters relate to the progression of the epidemics, and thus, they become a valuable tool to control disease.One of the most important parameters in epidemic models is the basic reproduction number or R 0 , which is basically a measure of the infection potential of the disease.If R 0 ≤ 1 then the disease will fade out eventually with probability one, so, R 0 becomes an important threhsold parameter.
There is relatively little work relating queuing theory with epidemiology.Kendall [2] and [3] discusses the relationship between M/G/1 queues and birth-death processes.Kitaev [4] works on the relation between birth and death processes and the M/G/1 queues with processor sharing, and the SIR epidemic model constructions used in [5] and [6] can be seen as M/G/1 queues with processor sharing.Ball and Donnelly [7] use M/G/1 theory to find the total cost of the epidemic and most recently Trapman et al [8] use M/G/1 queues with processor sharing to model an SIR epidemic with detections.We are not aware of work in which the epidemic process is considered an M/G/N queuing process implying that every individual is a server that can be busy (infected) or idle (susceptible).In this work we use this approach for the SIS model, an we extend the results to an SEIS model.
In closed population stochastic epidemic models, the process reaches an absorbing state once the population is free of infected individuals.Absorption into this state will occur eventually, with the time to reach this state depending strongly on the infection potential ρ = λ/µ, where λ and µ are the individual contact and recovery rates respectively, therefore, the process is said to have a degenerate limiting distribution.
The expected time to extinction generally increases with ρ (see [9], [10]); thus, an interesting property of these process is its behavior before going to absorption.Quasi-stationary distributions provide information on the limiting distribution of the process conditioning on non-absorption (see [11] and [12]) that is, they allow to analyze the behavior of the disease while it is in the endemic state.These distributions are difficult to find, and instead, some approximations were suggested like the "reflecting state 0" or the "one permanently infected" (see [13]).These approximations will be explained later.
Although the current approximations to the quasi-stationary distribution of the number of infectives seem not to have a closed "known" distribution, we show that a good approximation to the quasi-stationary distribution for the susceptible is Poisson distributed for a general distribution of the duration of the infectious state.Since the next natural step is to add a latent period that allows for the latency of the infection, we have applied this approach to approximate the quasi-stationary distribution (bivariate) of an SEIS epidemic model.Quasi-stationary distributions assume in general that the distribution of the infectious time is exponential.In this paper all approximations assume a general distribution of the duration of the infectious period.
Here we derive the distribution of the approximation "one permanently infected" for SIS and SEIS models [13].We use the term "conditional endemic distribution" for this type of distributions in stochastic epidemic models.The paper is organized as follows: sections 2 and 3 deal with SIS and SEIS models respectively.Section 2.1 introduces the SIS model; Section 2.2 introduce the quasi-stationary distribution while 2.3 deals with the quasi-stationary distribution for the SIS epidemic model.In Section 2.4 the case of a general distribution of the duration of the infectious state is analyzed using standard results from queuing theory.Numerical comparisons of these results via simulations are presented in Section 2.5 Section 3.1 introduces the SEIS model.In Section 3.2 it is shown that the quasi-stationary distribution of the number of infected individuals (latent + infective) can be approximated with that of an SIS model with appropriate parameters.In this section we also derive an approximation to the joint quasi-stationary distribution of latent and infective, and numerical comparisons for the SEIS are presented in Section 3.3.

Introduction.
The susceptible-infected-susceptible (SIS) stochastic epidemic model attempts to reproduce the behavior of epidemics running through a population with no vital dynamics; that is, no births and deaths occur.In this model it is assumed that no individuals are removed from circulation either by immunization or isolation.The deterministic version of the SIS model was introduced in [14] and has been fully analyzed since then.Its stochastic counterpart, also called the stochastic logistic epidemic model (see [15] and [16]) was introduced early in [17], and has been applied similarly to study the transmission of rumors (see [18]).However, most of the relevant results concerning this model have been in the epidemics context.
In the SIS model, susceptible individuals may become infected by contact with infective individuals, and hence, it is assumed that there is no incubation period for the disease; that is, infected individuals become infectious immediately.After some time, they become healthy and susceptible again.The process of contagion is assumed to be driven by the homogeneous mixing of individuals in the population.
We let I(t) and S(t) be the number of infective and susceptible respectively at time t, and note that since the population is closed, the state of the process at time t can be fully described by either I(t) or S(t).It is customary to follow I(t).I(t) takes values on Ω = {0, 1, 2, . . ., N }, N being the population size.The SIS stochastic epidemic model is a discrete space, continuous time Markov Chain.In particular, it is a unidimensional continuous time birth-death process.One way to dissect the SIS epidemic process by specifying the following two rules: i ) (Homogeneous mixing) Every individual comes into contact with another at random intervals which are independent and identically distributed random variables.The distribution of these intervals is exponential with parameter β.If a contact involves a susceptible individual and an infectious one, the probability of infection is θ.That is, β is the contact rate and λ = βθ is the effective contact rate.
ii ) The duration of the infectious state is an exponential random variable with parameter µ, that is, µ is the recovery rate.
From i), all k infective individuals come into close contact with others according to a Poisson process with parameter β k.Since every one of these contacts would be with a susceptible individual with probability (N − k)/N , and will result in an infection with probability θ, by thinning the Poisson process, we conclude that for fixed k, contacts between infective and susceptible individuals that will end up in a infection of the susceptible occur according to a Poisson process with parameter λk(N − k)/N .Upon defining the transition probabilities are Here, o(δ) is a function such that lim δ→0 o(δ)/δ = 0 2.2.The quasi-stationary distribution.We are interested in the future of the SIS epidemic, which depends strongly on the ratio ρ = λ/µ, called the transmission factor, basic reproduction ratio, basic reproduction number or infection potential.
The deterministic model has a threshold at ρ = 1, and results in an endemic infection of size N (1 − 1/ρ) if this value is greater than 1.In the stochastic model, since {0}, the disease free state is an absorbing state that can be reached with positive probability from every state, the process will end up in this state for any value of ρ given a sufficient amount of time, for any initial number of infected as long as N is finite.The certainty of extinction imposes a problem if our interest is in the long-time behavior of the epidemic, conditioning in the process not being in the absorbing state.Let X j (t) being the time spent in state j up to time t, j ∈ Ω.Thus, π j denotes the proportion of time spent in state j as t → ∞.Π is called the stationary distribution of the process.If p j (t) denotes the probability that the process is in state j at time t then it is possible to give an alternative representation for π j π j = lim t→∞ p j (t).
The stationary distribution for the SIS epidemic model is also degenerate, that is, will eventually be absorbed at state 0, which is true for any finite value of ρ.However, when ρ > 1 it is reasonable to assume that the disease will be in an endemic state for a while, and any information on the behavior of the process previous to extinction would be useful in the understanding of the epidemic.This interest led to the development of the concept of quasi-stationary distributions.Quasi-stationary distributions (QSD) are limiting distributions conditioning on the process not being in an absorbing state.If we let Q = {q 1 , q 2 , ..., q N } denote the quasi-stationary distribution of the SIS epidemics, then .
Observe that when ρ > 1, Q is a conditional endemic state distribution.No simple expression exists to calculate Q, although Nåsell (see [19]) proposed a numerical algorithm for its calculation.Most of the relevant work regarding the calculation of analytical expressions for the QSD of the SIS model is based on approximation methods.Two of these approximations were suggested by Kryscio and Lefèvre (see [13]) and analyzed in detail in [19], see also [16].The common characteristic of these approximations is that the process is modified in such a way that it lacks the absorbing state {0} and thus the possibility of degenerate distributions is avoided.
In one approximation the number of infected in the population is at least one, and it is called the SIS model with one permanently infected individual.In this process every recovery rate µ j = µ j is replaced by (j − 1)µ, while the infection rates are unchanged.In the second approximation, the only rate that is changed is µ 1 , which is replaced by zero.This latter approximation is referred to as the SIS model with the origin removed or reflecting state approximation model (see [20]).Andersson et al [10] studied the SIR and SEIR epidemic models with demography and obtained expressions for the time to extinction starting from the QSD.In their work, the duration of the infectious state was not restricted to be an exponential distribution but instead a gamma distribution.Ovaskainen [19] improved previous approximations for the QSD and the time to extinction for two cases: when N → ∞ and when the basic reproduction ratio R 0 → ∞.Nåsell [21] gave approximations for the QSD and the time to extinction for the stochastic version of the Verhulst epidemic model, depending on R 0 being greater, smaller or in the transition region.For the first two regions, the approximations yielded normal and geometric distributions respectively.Nåsell [22] derived approximations to the QSD for SI, SIS SIR and SEIR epidemic models with demography for N large.For the region R 0 > 1, Nåsell concluded that the normal distribution yielded reasonable approximations to the QSD.
A different approach to obtain the reflecting state approximation model was suggested in [23] to obtain the time to extinction starting from a general state k.Stefanov and Wang [24] and Ball and Stefanov [25] derived higher moments for the time to extinction for the SIS starting from a general state k.
Here we use q (1) j to denote the approximation to q j when using one permanently infected individual and q (0) j for the reflecting state approximation model.For the SIS, it has been shown in [13] that when ρ > 1 and N → ∞ that is, both approximations converge in distribution (see also [20]).Also, in [19] it is proved that for ρ > 1, the distribution of the number of infectives under both the reflecting state 0 approximation and the one with one permanently infected individual are approximately normal with mean N (1 − 1/ρ) and variance N/ρ when N → ∞ and ρ is fixed.

2.3.
The approximation to the quasi-stationary distribution of the number of infectives.In the following calculations we derive an approximation for the QSD of the susceptible individuals.We establish a new result: that the distribution of the QSD of the number of susceptible can be approximated with a Poisson random variable when N is large and ρ is constant.
N } be the stationary distribution of the approximation to the QSD when there is one permanently infected individual.Here µ j = µ(j − 1) and λ j = λj(N − j)/N , j = 1, 2, . . ., N .It can be shown (see [16]) that the q (1) n follow the relationship: N −k as the QSD approximation to the number of susceptible.Thus Solving for q (0) For ρ > 1, which is the epidemiological relevant case, when N → ∞ we have a Poisson random variable with parameter N ρ −1 .Nåsell [19] had already established that the approximations to the QSD for the infective individuals yielded a normal distribution with mean N (1 − ρ −1 ) and variance N/ρ for N → ∞ and ρ constant.
From (2) we can see that when N → ∞ using π N −k = p k , the approximation to the QSD for the number of susceptible is a Poisson random variable, which in turn can be approximated with a normal distribution with mean and variance N/ρ; hence, that of the number of infected is normal with mean N (1 − ρ −1 ) and variance N/ρ.

2.4.
The approximation to the quasi-stationary distribution of the number of susceptible for non-exponential duration of the illness state: An application of queuing theory.An SIS epidemic process can also be seen as a queuing process with state dependent arrival rate.We can think of N servers, where a busy server corresponds to an infected individual.Individuals are served at a rate µ, and the arrival rate is λj(N − j)/N when the number of busy servers is j.We will use this analogy to derive an approximation to the QSD of the number of susceptible individuals .The constant hazard rate characteristic of the exponential distribution makes it difficult to adapt for most diseases.The idea of assuming that the probability that a person will be cured in the next s units of time given that she/he has been for t units of time is independent of t is not very realistic.In this section we derive an approximation to the QSD of SIS models when the illness state is non-exponential.
From queuing theory, (see [26]) the limiting distribution of the number of busy servers in a system M/G/N , with state dependent arrival rate is given by where E[S] is the expected value of the service time.Hence, the approximation to the QSD for the number of busy servers (infectives) depends on the duration of the service (illness) state only through its first moment.The general expression for the approximation to the QSD of the number of susceptible becomes Observe that, for N (λE[S]) −1 moderately large, p k can be approximated by a Poisson distribution with parameter N (λE[S]) −1 .
In the equilibria, the rate of new infections equates the recovery rate, thus, in the equilibria the proportion of the number of infected is (1 − µ/λ).This can also be interpreted as the long-run proportion of time an individual is infected.In our analogy with queuing models, this is the long-run proportion of time a server is busy or server utilization.From here, we arrive to an expression that relates the basic reproductive number in epidemiology, R 0 , with ρ, the server utilization of the system: Simulations.In this section, we simulate the approximation 'one permanently infected' to the SIS epidemic model with N = 50 and ρ > 1.We assume a gamma distribution for the duration of the infectious state with three sets of  parameters.Only the distribution of susceptibles is plotted against the approximation given in (4). Figure 1 shows the gamma distributions used for the duration of the illness state for a shape parameter α and scale parameter β.The goal was to analyze the effect of different distributions for the duration of the infectious state keeping the average time of the infectious state (and ρ) constant.As expected from (4), the three gamma distributions yielded the same distribution and only one is plotted.Figure 2 shows the time spent in every state in the SIS epidemic model with a permanent infected individual (histogram) and the Poisson approximation (solid line).
The simulations were performed in MATLAB.These were implemented as follows: for each parameter set, a stochastic SIS epidemic was simulated during 2000 units of time, and the time spent in each state during this interval was recorded.If the epidemic went to extinction before reaching 2000 units of time, it was discarded, although due to the ρ used this would not happen very often.The distribution of the time spent in each state is plotted against the Poisson approximation given by (4).It can be seen in Figure 2 that the approximation is good in every case except in the right tails, which is due to the time at which the simulations stopped, and will tend to be more precise when the simulation time is increased.As expected from (4), the approximation to the QSD is insensitive to higher moments of the duration of the infectious state.

Introduction.
A natural generalization of the SIS model is to consider additional epidemiological states.Here we consider the possibility that an infective individual undergoes a latent state before becoming infectious.Thus, individuals who become infected are not able to transmit the disease for a period of time.Since infected individuals in the latent state can not transmit the disease or acquire additional infection, they play no role in the transmission of the epidemic but rather serve as buffers or reservoirs of infection.This model is referred as the SEIS model or SLIS model.
Quasi-stationary distributions for SEIS models are two-dimensional arrays, π m,n , where π m,n denotes the limiting proportion of time in which there were m latent and n infective individuals conditioning in the process not being in the absorbing state.The first effort to analyze the joint QSD of an epidemic model is reported in [20] with a stochastic version of the Ross malaria model.
It is possible to define a QSD for the total number of infective individuals π k where k = m + n.In this section it is shown how the QSD for the total number of infected in a population in an SEIS model can be derived from that of an SIS model.From this last result the two-dimensional QSD follows directly.
The state space of the process can be stated in terms of the number of latent and infective individuals, (E,I) as: The following diagram shows the transitions between the different epidemiological states: The quasi-stationary distribution of the SEIS model.Define: where E(t) and I(t) are random variables that denote the number of latent and infective individuals at time t, respectively.Clearly {k, m}, {i, j} ∈ Ω.The instantaneous transition probabilities are while all other events are assumed to occur with probability o(δ).In section 3.2.1 an approximation to the QSD of the total number of infected (latent + infective) is derived.In section 3.2.2 the joint QSD distribution of the number of latent and infective is analyzed.

3.2.1.
The quasi-stationary distribution of the total number of infected in the SEIS model.We introduce the random variable I * (t) = E(t) + I(t) which denotes the total number of infected individuals at time t.Define also: and call p * k the k-th element of the QSD of this process.To derive the approximation to the quasi-stationary distribution, it suffices to start from the instantaneous transition probabilities when t → ∞ in (5), and considering that j in this expression is k, k − 1 or k + 1.
Now choose an individual at random.From renewal theory (see Ref. [27], pp.80-86) as t → ∞, the probability that this particular individual is in the infective state given that he is infected (latent or infective) is Observe that, if at some fixed time t two individuals are not in state S, the probability that one of them is in state E (or I) is independent of the other being in state E (or I) then, conditioning in exactly k infected, the number of infective individuals follows a Binomial distribution with parameters k and θ I , that is Using these results,the instantaneous transition probabilities are Equations ( 8) and ( 9) are formally equivalent to those given in (1) with λ θ I and µ 2 θ I replacing λ and µ respectively.Therefore, an approximation to the QSD of the total number of susceptible individuals in an SEIS model can be obtained directly from that of an SIS model (Sec.2.3), which is given by a Poisson distribution with parameter N µ 2 /λ: For SEIR (Susceptible-Latent-Infective-Recovered) models, Trapman et al [8] had found that the sum of the number of infective and latent individuals was the same as the number of infected individuals in an SIR (Susceptible-Infective-Recovered) model up to a random time change.In Section 2.4 we found that the QSD for the number of infective individuals (busy servers) in an SIS epidemic model depends on the duration of the infectious period (service time) only through the first moment.Since an SEIS epidemic model can be seen as a queuing system with two phases on each server, the average service time is µ −1 1 + µ −1 2 while the arrival rate is λ θ I .Thus, using (4), the parameter of the Poisson distribution for the number of susceptible becomes ).Using (6), the mean number of infected becomes N µ 2 /λ, as expected.Thus, the joint QSD has its mean at It is important to observe the resulting lack of dependence of the above result on higher moments of the duration of the latent and infectious period.

3.2.2.
The marginal joint quasi-stationary distribution of the number of latent and infective individuals in the SEIS model.The QSD joint distribution for the number of latent and infective individuals can be derived from expression (10).When N µ 2 /λ is large the number of susceptible is given approximately by a Poisson distribution with parameter N µ 2 /λ.Under the assumption that N µ 2 /λ is large we derive an approximation to p m,n , the joint QSD for the number of latent and infective individuals respectively.
The joint distribution is defined as which can be rewritten as (10).Using (7) we conclude that The computation of the marginal QSD of the number of latent and infective is straightforward.Define the marginal QSD's for the number of latent and infective by α m = lim These results are evaluated through simulations in the next section.
3.3.Simulations.In this section, we simulate several SEIS epidemic models with different values of λ, µ 1 and µ 2 , with N = 100, allowing for one permanently infected individual.We use an exponential distribution for the duration of the latent and infectious states.Both observed and theoretical distributions are plotted for a) the number of susceptible and b) the joint number of infected and latent.A contour diagram proves to be very useful for comparing the theoretical and observed joint distributions.For every set of parameters, a simulation of an SEIS epidemic model was run for 5000 units of time, and the proportion of time spent in each state was recorded.Again, if the epidemic went to extinction before reaching 5000 units of time, it was discarded.The approximation to the QSD for the total number  of susceptible is shown in a histogram with its Poisson approximation.An insert shows the joint distribution of latent and infected using a contour diagram together with the approximation (11).
The approximation to the QSD for the total number of susceptible and for the joint distribution in Figures 3-5 seems very good.In Figure 3 the population size is N = 100,with λ = 4, µ 1 = 1/2 and µ 2 = 1. Figure 4 was constructed using the same population size but with λ = 6, µ 1 = 1 and µ 2 = 1/2, that is, the relationship between the duration of the latent and infective states is reversed compared to the previous case.Since the duration of the infected state has been increased as well as the infection rate λ, the joint distribution has moved towards a higher proportion of infected individuals than in the previous case.Figure 5 uses N = 100, λ = 6, µ 1 = 3 and µ 2 = 1/2, that is, the duration of the latent state has been reduced compared to the previous case, which increases the number of infective.

4.
Discussion.Here we derived the exact distribution of the approximation "one permanently infected individual" to the QSD of the number of susceptible in an SIS model for N → ∞ keeping N/ρ constant.We showed the QSD is a Poisson random variable with parameter N/ρ.The QSD had a similar distribution than that of the number of busy servers in an M/G/N queuing process, and thus we showed that the QSD holds for a general distribution for the duration of the latent state with finite mean.This is important since while the assumption of Markovian contacts  For this N = 100, λ = 6, µ 1 = 1, µ 2 = 1/2, ρ = 12.Time = 5000 units.Insert shows the contour diagram of the approximation to the joint QSD (solid line) and its approximation (dotted line).may be feasible, assuming that the duration of the infectious or latent states is exponential may be questionable, due to the memoryless assumption implicit.The results borrowed from queuing theory allowed for an easy generalization to the SIS and SEIS models.
It is interesting to notice that while the asymptotic distribution of the QSD for the number of infected had been found to be Normal distribution, we can see now that this results as an approximation to the distribution of susceptible, which is a Poisson distribution.Since Poisson distributions with a large mean can be approximated with a normal distribution, we can approximate the QSD of the number of susceptible with a normal distribution, and from here the approximation for the QSD of the number of infective is also normal.
In approximations to the SIS model that lack of absorbing states, like the one used here, in which there is always one permanently infected, every individual goes through susceptible and infected states forever.This implies that even if λ < µ the disease will not vanish, thus, (1 − R −1 0 ) is the proportion of time an individual spends infected.
Regarding the SEIS epidemic model, suitable good approximations were derived for both the total number of infected and the joint distribution of the latent and infected.It was shown that the approximation to the QSD of the SEIS epidemic model with infection rate λ and mean time µ −1 1 and µ −1 2 respectively for the number of individuals in the latent and infected states is the same than that of an SIS  epidemic model with infection rate λ and mean recovery rate µ −1 2 .The simulations show that the approximation is accurate in both SIS and SEIS epidemic models, as expected from (10) for the number of individuals in the latent or infected state and from (11) for the joint distribution.Since these are limit distributions, it is important to run the simulations for a long time in order to obtain an appropriate sample from the targeted distributions.This can be easily seen in Figures 4 and 5, that differ only in the parameter µ 1 , but are almost identical.It is a matter of further study to analyze if epidemic models of the form S − E 1 − E 2 − E 3 − ... − E k − I − S behave similarly, that is, if they only depend on the infection and recovery rate.