CUMULATIVE AND MAXIMUM EPIDEMIC SIZES FOR A NONLINEAR SEIR STOCHASTIC MODEL WITH LIMITED RESOURCES

. The paper deals with a stochastic SEIR model with nonlinear incidence rate and limited resources for a treatment. We focus on a long term study of two measures for the severity of an epidemic: the total number of cases of infection and the maximum of individuals simultaneously infected during an outbreak of the communicable disease. Theoretical and computational results are numerically illustrated.


1.
Introduction. The SEIR model is one of the basic compartmental epidemic models, which represents infection dynamics of an epidemic process showing a significant period of time during which the host of the infection is no able to transmit the pathogen.
Commonly, the understanding of the dynamics of the spread of an infectious disease relies on deterministic and stochastic epidemic models. The former are defined by a set of differential equations related to the community proportion of individuals showing each possible health state at any time t. Stochastic models reflect randomness in epidemic events hence they are preferable when infection flows in small communities where individualities do not behave equally likely regarding contacts, infection or incubation times. Therefore, chance fluctuations really affect overall results.
For both type of models the incidence is described by a function measuring the rate at which a susceptible gets the infection. The most frequently used is the bilinear incidence rate or mass action function βSI, proposed by Kermack and McKendrick [28], where S and I denote the number of susceptible and infectious, respectively. Bilinear rate is the natural assumption when dealing with a homogeneous well mixed population. General incidence rates of the form β(S, I) introduces in the model behavioral changes of susceptible individuals, the crowding effect of the infective ones or even the response of the immune system to the concentration of the infectious agent (see [14], [15], [24], [25] or [34]).
After infection, patients recover due to their own immune system and/or to the help of a treatment. Most of classical epidemic models assume a recovery rate proportional to the number of infective individuals, which is in accordance 3138 JULIA AMADOR AND MARIAJESUS LOPEZ-HERRERO to the idea that recovering depends on the individual's immune system. However, being recovered from diseases such us measles, flu or tuberculosis, depends also on medical treatment. One can assume that treatment is available for every infected but it is more realistic to consider the case of limited access to health facilities. Resources for treatment should be large enough, but if a large quantity is fixed, community pays unnecessary cost for it, instead if it is too short, then population is not protected against an eventual outbreak of the disease. Thus, to consider recovery rate depending on a saturated treatment function contributes to make the model closer to reality. Wang [39] introduced a treatment function proportional to the number of infective while there are treatment resources, and otherwise, it takes the maximal saturated level. A saturated treatment of Verhulst type is considered in [38], [42] and [43], that is increasing for small number of infective and tends to the saturated value for large number of infective.
Our research deals with a stochastic SEIR model with general incidence rate and the piecewise treatment, introduced by [39], that is given by a rate proportional to the infective individuals up to a threshold and constant beyond this level. We are interested in the computation of the long term distribution of the final size of the epidemic and in the maximum number of simultaneously infected individuals. Both of them are useful for measuring the severity of an epidemic and for appraising the effect of control measures [20]. The final size, or the cumulative number of cases of infection during the outbreak, is a quantity of interest to statistical epidemiologists who use it to determine the attack rate; i.e., the incidence proportion or frequency of morbidity for an attacking pathogen in a population. The maximum number of infective in the course of the epidemic is the peak size of the epidemic curve and it gives an idea of how large treatment resources should be with the purpose of keeping demand for facilities around the available amount. In the paper we will compute distributions by following the first-step methodology which produces stable algorithmic results when applied to Markovian compartmental models ( [2], [9]). This methodology is a generalization of the one introduced by Neuts and Li in [32] for the numerical implementation of the final size distribution of the stochastic SIR model. A critical review and a systematic comparison of methods for evaluating the probability mass functions of the total and maximum sizes during a stochastic epidemic is made in [26].
The rest of the paper is organized as follows. Section 2 describes the SEIR stochastic model. Section 3 contains theoretical and algorithmic results for the total size of the epidemic. In Section 4, the probability mass function for the peak size is studied. Implementation of theoretical and algorithmic results is illustrated on Section 5. Finally, we give a discussion of results and outline future work.
2. Model description. We consider a homogeneous closed population of N individuals which is affected by a communicable disease conferring immunity once recovered and presenting a latent period. Individuals in the population show no preferences in their relationship and we can assume that contacts occur at random.
Regarding the health state, population can be classified as: susceptible S, exposed or latent E, infective I and removed R (either recovered and immune, or isolated or dead). We assume that, while they are exposed, individuals incubate the disease but they do not transmit it; i.e., they are infected but not infectious yet. Latent periods of simultaneous exposed individuals are independent and we model them as independent exponentially distributed random variables, with expected length given by 1 σ . The pathogenic agent is transmitted by direct contact between an infective individual and any susceptible, at time points of a time homogeneous Poisson process with a general non linear incidence rate depending on infectious and susceptible individuals: β(S, I). Recovering of infective individuals depends upon a treatment that confers immunity but it is not available to all the infective. When individuals receive it, each one recovers independently of the rest in an exponential time of rate γ (i.e., 1/γ is the expected time to an infectious be recovered when the treatment is applied). It follows that recoveries in the population are exponentially distributed with a piecewise recovering function of the form where I 0 is a fixed constant. Let us denote by S(t), E(t), I(t) and R(t) the number of susceptible, exposed, infective and recovered individuals at time t, respectively. For practical purposes we assume that the starting point of any epidemic outbreak is a single infectious individual. As we are dealing with a closed population S(t)+E(t)+I(t)+R(t) = N and the evolution of the epidemic is described by a three dimensional continuous time Markov chain (CTMC ) , where part of states -those corresponding to situations where there are neither exposed nor infective individuals-are absorbing. Let us represent by Then the infinitesimal generator Q, of the CTMC X has a block triangular representation which is suitable for implementing the recursive algorithmic procedures exposed in the rest of the paper.
As X is a CTMC with finite state space, the absorption into S A is certain and it occurs in a finite expected time. Consequently, the epidemic process consists of a single outbreak that ends as soon as X enters into the set S A , and the infection remains active during a period of time having a finite mean length.
3. Total number of infected individuals in an outbreak. The epidemic final size has been the subject matter of many research papers, a revision on the literature can be found in [26]. For Markovian models, Stone et al. [36] studied an SIS model with external source of infection, Artalejo and Lopez-Herrero [7] investigated the number of different individuals that become infected during an outbreak of an infectious process modeled in terms of a Birth-and-Death model, Black and Ross [12] develop a new methodology for computing the epidemic final size distribution for basic SIR models and generalizations of this, Amador in [4] and Amador and Artalejo (see [5], [6]) explore the probabilistic behavior of the number of cases of infection regarding computers security. Recently, an extension of a deterministic SEIR model was studied by Barbarossa et al. [11] for understanding the dynamics of the 2014 Ebola Virus Disease outbreak.
In this section we are interested in characterizing the random nature of Z defined as the total number of individuals that get the infection during an epidemic outbreak; that is the final size of the epidemic process. The objective is to present a recursive algorithm for determining the distribution mass function of the random variable Z.
Let us introduce some auxiliary probabilities associated to the current state of the epidemic process. Let x k sei be the probability of having k additional infections until the end of the outbreak, given that at present the state of the population is (s, e, i) ∈ S.
In order to study Z, we observe that the outbreak starts from the state (N − 1, 0, 1) and consequently Now we focus on the obtainment of the auxiliary probabilities. First we write down a trivial result where δ ij represents Kronecker's function, defined as 1 when i = j and by 0 otherwise.
Let us consider a fixed integer k, 0 ≤ k ≤ N − 1. A first step argument, conditioning on the first transition out of the state (s, e, i) ∈ S T , yields the fundamental equation For each fixed k, equation (2) and the trivial result (1) are combined to compute auxiliary probabilities in the natural order (i.e., the above one in which we partition the state space S). More specifically, we need the sequence {x k N −1,0,1 , 0 ≤ k ≤ N − 1}, associated to the initial state of the outbreak. The recursive scheme can be summarized in the following algorithmic structure. Algorithm 1. Probabilities x k sei , of having k additional infective individuals in the rest of the outbreak given a current state (s, e, i) ∈ S, are computed by using the following recursive scheme: Step1: Set k = 0.
Remark 2. Probabilities P {Z = k} can be numerically determined from Algorithm 1. Moreover, when k = 1, 2 they present the following explicit forms: that comes from exponential distribution properties, by observing that there is only a single case of infection during the outbreak if and only if the initial infectious recovers prior to establish any contact with the rest of the population. And that derives from (2) after some algebra.
Notice that first probability (3) is trivially independent of the latent period and a simple use of algebra will show that probability (4) does not depend on σ for the classical choice on transmission and recovery functions (i.e., β(s, i) = βsi, γ(i) = iγ). This fact is satisfied for any number of infective in the outbreak when the model involves a more general class of transmission functions. We state the result in the corollary of the next theorem.
Theorem 3. In an SEIR model with transmission and recovery functions satisfying β(s, i) = iβ(s, 1) and γ(i) = iγ, probabilities x k sei , for (s, e, i) ∈ S and 0 ≤ k ≤ N −1, of having k additional infective individuals prior the end of the outbreak, given the current state (s, e, i), are given by the following set of explicit formulae, where any empty product is taken to be the multiplicative identity (with the appropriate dimension) for e = k, is a diagonal matrix with main diagonal entries given by the vector v. Finally, M j is a j × (j + 1) matrix, with 0's under the main diagonal and 1's elsewhere. Notation (v) T corresponds to a row vector v and 1 n is an ndimensional all-ones vector.
Proof. The proof lays on a four order induction on the set F = {(s, e, i, k) : (s, e, i) ∈ S, 0 ≤ k ≤ N −1} which can be indexed linearly according to the ordering described in the recursive scheme developed in Algorithm 1.
After some algebra, it is possible to proof the following basic cases dealing with probabilities of having at most three infections in the rest of the outbreak. According to the notation used in the statement of the theorem, these are Now we consider another state (s, e, i, k) ∈ F , with 3 ≤ k ≤ N − 1. Let us assume that the statement is true whenever we consider a state (s , e , i , k ) < (s, e, i, k).
First, we point out that the multiple sums in the expression (6) can be written in the matrix-vector form appearing in the description of the function A(s, e, i, k) introduced in (5). So, we can assume that auxiliary probabilities for any state (s , e , i , k ) indexed before than (s, e, i, k) are expressed as Next we observe the relationship shown in equation (2), where probabilities x k sei depend on three auxiliary probabilities of states ordered before the state (s, e, i, k) regarding the ordered structure of F . So, plugging in (2) the expressions given by the induction hypothesis, (7), we get: .
After some algebra, we write the expression in brackets in equation (8) as q sei A(s, e, i, k)/(β(s, 1) + γ) and rearranging terms we get the inductive result.
As a consequence of the previous result, we get explicit expressions for the distribution of the final size of an epidemic outbreak.
Corollary 4. Under the same hypothesis and notation of Theorem 3, the final size distribution presents the following explicit form: For simplicity, B j denotes the j− dimensional vector b N −1,j (j) introduced in the above theorem. Again any empty product is taken to be the multiplicative identity with the appropriate dimension.
Numerical implementation of this result, when the mass action transmission function β(s, i) = βsi is considered, provides the same final size probabilities of a classic stochastic SIR model (see [1], [12] or [18]). Which is in agreement with the well known result that the final size distribution of the classic stochastic SEIR is insensitive with respect to latency parameter [10].

Maximum number of simultaneous infective individuals in an outbreak.
In this section we consider again an outbreak of the epidemic process, our interest is to study the peak prevalence or number of individuals that require to be treated from the infection at the same time.
As we are dealing with a model having a limited number of resources the interest in characterizing the maximum of simultaneous infective individuals involved during an outbreak is well justified. In fact control strategies should include the objective of lowering the peak size in the way that demand for facilities stand below the available resources.
In spite of the applicability of the peak prevalence to control the spread of an infectious disease [20], literature dealing the maximum number of infective individuals is very sparse. In [36] the distribution of the maximum number of an SIS model with immigration is studied, Artalejo et al. [8] provide efficient algorithms for computing the distribution of the number of infected in transient and long term regimes for an SIS model, the paper by House et al. [26] also summarizes algorithmic procedures applied to calculate the maximum size distribution.
Let us denote by M the random variable that gives the maximum number of simultaneously infected individuals that can be observed during an outbreak. M is a discrete random variable with a finite number of mass points, due that we are dealing with a closed population of size N . So our goal is to determine the mass distribution function of M P {M = m}, for 1 ≤ m ≤ N.
As we did in the previous section, we are going to introduce a set of auxiliary conditional probabilities. In this particular case, probabilities are subjected to the current state of the epidemic process, updating the general compartmental situation in the population and also recording the maximum number of simultaneous infective individuals historically observed at any intermediate epoch prior to the end of the outbreak. In more detail, for any integer m, 1 ≤ m ≤ N , we define y m seik as the probability of having at least m simultaneously infective individuals in the outbreak, given that the current situation in the population is s susceptible, e exposed, i infective and until now we have recorded a maximum of k simultaneous infective individuals.
All of the epidemic outbreaks start from a single infective individual so the event {M ≥ 1} almost surely occurs and the level crossing probability is The rest of tail probabilities, for 2 ≤ m ≤ N , can be recovered from the auxiliary probabilities by noticing that For any given integer m, 2 ≤ m ≤ N , we denote by The point is to determine the auxiliary probabilities {y m seik : (s, e, i, k) ∈ S m }. For some special states auxiliary probabilities are trivially determined: Result in (10) also holds when the current state of the population do not guarantees enough number of individuals as to reach the level of m simultaneous infected. That is Probabilities conditioned to the remainder group of states in S m are determined recursively from linear equations which are derived from a first step argument. More specifically, given a state (s, e, i, k) where s + e + i ≥ m, by conditioning to the possible events going out of the current state we get The sequence of probabilities {y m N −1,0,1,1 : 2 ≤ m ≤ N } driving to the distribution of M is determined from the above equations (12)- (13) with the help of the trivial results (9)- (11). For each fixed m, with 2 ≤ m ≤ N , auxiliary probabilities y m seik are recursively computed in the natural order, that is 0 ≤ s ≤ N − 1, 0 ≤ e ≤ N − 1 − s, 1 ≤ k ≤ m, 0 ≤ i ≤ k, using the following pseudo-code Algorithm 5. The set of probabilities y m seik , 2 ≤ m ≤ N, conditioned to the current state (s, e, i, k) ∈ S m , are determined as follows: Step 1: Set m = 2.
Step 3a: If IC < m and k < m then set y m seik = 0.
Step 3b: If IC ≥ m, k < m and i < k then compute y m seik via equation (12).
Step 3c: If IC ≥ m, k < m and i = k then compute y m seik via equation (13).
Step 3d: If k = m then set y m seik = 1.
In contrast with the final size's behavior regarding latent periods, maximum size's distribution is really influenced by the latency parameter, even for the bilinear mass action choice, as numerical examples in the following section we will show. 5. Numerical results. In this Section we present several numerical results to illustrate the influence of the incidence function, β(s, i), the level of saturation for resources, I 0 , and the latency rate, σ, on the distribution of the final size, Z, and the maximum simultaneous number of infectious individuals, M , observed in an outbreak.
Numerical scenarios will share the population size, N = 100, and the value γ = 1.0 of the individual recovery rate when treatment is applied; consequently the rate recovering function is γ(i) = min{i, I 0 } for 0 ≤ i ≤ N.
Regarding the mechanism of transmission of the disease, we consider three different incidence functions. The bilinear mass action β 1 (s, i) = βsi/N , which assumes equal contribution of susceptible and infective individuals in the transmission; a rescaled version of the incidence function appearing in [34], β 2 (s, i) = βsi/N (1 + i 2 ), modeling inhibitory effect for susceptible when the number of infective is getting large and the force of infection given by β 3 (s, i) = βsi(1 + i)/N, considered in [16] in the context of an SEIR reaction-diffusion model. These three functions present the same initial infection force βi and satisfy β 2 (s, i) ≤ β 1 (s, i) ≤ β 3 (s, i), for any fixed number of susceptible individuals.
We articulate the exposition of results in two parts: The first one deals with the cumulative size of the epidemic and the second focusses on the maximum number of simultaneously infectious individuals.  Table 1 Notice that, the risk of having an outbreak where all the individuals are affected by the disease is around the 90% assuming a reaction-diffusion transmission, around the 40% under mass action and practically impossible assuming the inhibitory incidence function. This is in agreement to the relationship satisfied by transmission functions. In comparison with the mass action transmission, the inhibitory choice produces a decrease in the number of contacts for a large number of infective while the opposite behavior is produced when reaction-diffusion transmission is considered. Similar comments can be done when comparing distribution's quartile or the expected number of cases of infection, for the three incidence choices. Next, we focus on the effect of the latency rate σ on the distribution of Z. Figure  1 presents the mass distribution function for the final size when σ ∈ {0.00005, 0.05, 1.0}. The incidence function is β 2 (s, i) = βsi/N (1 + i 2 ), with β = 5.0 and the limit for resources is I 0 = 5. For a better visualization, distributions are represented for mass points in [2,100], the probability of having a single infectious individual in the outbreak does not depend on the latency rate and its value is 0.287769 according to equation (3). The distribution shows a single mode at (Z = 1) for σ = 1.0 and shows a bimodal shape for σ = 0.00005 or σ = 0.05. During an outbreak, infections with a large latency phase (i.e., small value for σ) give more chance to a big cumulative number of infective individuals than infective processes showing shorter latent periods.

5.2.
Maximum of simultaneous infectious individuals. Now we deal with the random variable M , the peak prevalence of the infectious disease. First we consider an infection transmitted according to the mass action function β 1 (s, i) = βsi/N , with contact rate β = 10.0, the threshold for treatment resources is of I 0 = 5 units. Figure 2 shows the mass distribution function of the maximum number of simultaneous infected in a population of 100 individuals for latency rate σ ∈ {0.05, 0.5, 1.0, 10.0}.
We notice that the 10% of the epidemics outbreaks involves a single infective individual. The remainder 90% is located, when latency parameter takes small values (i.e., for epidemics showing large latency periods), in intervals that are lower than these observed for large latency parameters (i.e., short latency periods).
Under mass action transmission, latency does not affect the cumulative number of cases but it does influence on the peak prevalence of an outbreak.
Next, we describe the distribution of M using a Box-Whiskers plot diagram. The objective is to compare the patterns of the epidemic peak prevalence when we vary the available level of treatment resources. The box encloses the middle central part of the distribution, lower and upper edges of the box correspond to the lower and the upper quartile, respectively, and the line drawn across the box shows the median of the distribution. Finally, whiskers above and below the box cover the 95% of the distribution.
In Figure 3, we display the distribution of the peak prevalence for I 0 ∈ {5, 25, 50, 75, 100} when considering a mass action transmission with contact rate β = 10.0 and σ = 1.0. Additional numerical results, not included here, for another choices for the latency rate and/or saturated transmission functions present a similar behavior. We observe that a short level of resources conducts to a large size groups of simultaneously infective individuals; the peak prevalence decreases, in general, when I 0 increases. However, to guarantee treatment resources for more of the 50% of the population produces no significative changes on the maximum number of infective. This fact shows the importance of the random variable M to design an appropriate control strategy on the treatment facilities. Figure 4 is a contour graph for the probability P (M ≤ I 0 ) of having enough treatment resources for the peak prevalence, when we vary the initial force of infection, β, and the latency rate. The mechanism of transmission is modeled according to  Main results driven from theoretical and numerical results are summarized as follows: • For transmission functions of the form β(s, i) = iβ(s, 1) and linear recovery functions γ(i) = iγ(1), the mass distribution of the final size of infective individuals presents an explicit form expression which is not influenced by the latency rate σ. • Algorithmic procedures have been checked for moderate population (up to 1000 individuals) and produce numerical results at a low computational cost in a reasonable lapse of time. • Numerical results evidence that marginal distribution of Z (except for the above mentioned scenario) and of M depend on the latency phase of the infective process. • The probability of having enough treatment resources for the peak prevalence, depends on the initial force of infection and also on the extension of the latency period. The study can be continued in several directions. It would be of interest to generalize results to more complex models regarding distribution assumptions or the structure of the population. Relaxing the hypothesis of exponential distribution could be done by considering general infectious period distribution as in [22] or by employing the block-state-dependent-event approach as in [6] or [21]. As examples of future extensions of our work we mention the possibility of consider epidemic models at household levels [13], with population migration [17] or with varying population [23], that relax the homogeneity and constant size hypothesis involved in the present research.