# Stochastic Models for the Infectivity Function in an Infinite Population of Susceptible Individuals

Authors: Viswanathan Arunachalam; Liliana Blanco
Journal: Journal of Probability and Statistics
 (2017)

---

## Abstract

Two stochastic models to study the course of the transient behaviour of the total infectivity present in an infinite population of susceptible individuals are developed. The conditional intensity function of the contagion comprises two components: one is due to the external sources only and the other is the contribution of each of the infected persons which is nonstationary in nature. The statistical characteristics of the number of infected individuals at any time are explicitly obtained. Estimation of the model parameters is also indicated.


---

## Body

## 1. Introduction

Infectious diseases impact public health and consequently a nation’s progress in several ways. Diseases such as influenza, gastroenteritis among humans and bird flu, and foot and mouth diseases among animals are a matter of great public concern. In the face of such emerging diseases, time is of the essence in developing and implementing interventions. The ability to predict the transient infectivity function (infected population) imposed by such epidemics is essential for intervention design. Stochastic modeling of such infectious diseases is the only effective tool available for the health planners for optimal intervention design and their implementation.From the time Kermack and McKendrick [1] published their classic in the proceedings of the Royal Society, there have been various approaches to describe the dynamics of the spread of the infectious diseases. They include stochastic compartment models, conditional intensity models [2, 3], simulation models, inference models [4], Monte Carlo methods [5], and SIR model [6–9]. Recently, Lie algebra approach was used to study susceptible-infected-susceptible (SIS) epidemic models [10, 11].Kermack and McKendrick [1] model is characterised by the deterministic equation:(1)S′t=St∫0tS′t-τgτdτ-ft,where Strepresents the number of susceptible individuals at time t. This model is completely characterized by S0and the functions ftand gt.The function ftmodels a possible time dependent outside source of infection or represents the history of the epidemic up to time zero. The function gt describes the infectivity of an individual as a function of his age of the infection t, that is, the time elapsed since the moment of his own infection.Gielen [12] proposed a natural stochastic version of the above model. He observed that his model fills the gap between Bailey’s simple stochastic epidemic model [13] (which is obtained as a special case of Kermack and McKendrik’s model when the infectivity parameters ft and gt are constants) and the Reed-Frost model [12] (when the functions ftand gt are Dirac delta functions). Gielen [14] further proposed an interesting variant of the above model for a small epidemic in a large population of susceptible individuals. This model can be seen to be the direct stochastic version of the Keyfitz [15] deterministic model, governed by the renewal equation:(2)I′t=∫0tI′t-τgτdτ+ft,where It indicates the number of infective individuals at time tin an infinite population of susceptible individuals. Gielen [14] explicitly derived the probability distribution of the number of infective individuals at any time t, using which he arrives at several interesting results. However, the expression for the distribution of the number of infective individuals which involves n-fold integrals can at best be thought of as a formal solution and is not useful for explicit evaluations. Consequently, the statistical characteristics of the number distribution which are very important in the optimal decision making problems cannot be arrived at that in the general case. One of the purposes of the present work is to analyse the stochastic version of the Keyfitz model [15] and obtain explicitly the statistical characteristics of the number of infective individuals in an infinite population.The traditional epidemic models assumed the contagion rate to depend on the number of infected persons in the population, thereby tacitly assuming that all individuals exhibit the same pattern of infectiousness. However, this being not so in the reality led to the notion of infectivity function. The infectivity functionft is the force of infection at any time twhich is the probability/unit of time that a susceptible becomes infected. In the dynamic contagion models, the infectivity functions impact endogenous and exogenous factors of the underlying infection. The proposed stochastic infection models in the present paper split the conditional intensity function of infection into endemic and nonstationary epidemic components. The latter, in our models, is the spread of the disease by person to person contact and provides a description preserve at time t caused by each of the infected individuals. The endemic component models imported cases and is contrary to the epidemic component. The number of infected persons is modeled as a point process with its intensity function dependent on the point process itself, thus providing a more effective model to capture the contagion phenomena. The mean number of infected persons obtained from our stochastic epidemic model, for the special case when both the endemic and epidemic components of the intensity function are constants, satisfies Keyfitz’s renewal equation (2). Gielen’s stochastic version of Keyfitz’s model predicts a negative binomial distribution for the number of infected persons (see [14]) with the same mean thus validating our approach. The proposed formulation is broad enough to handle a wide variety of disease progression models for which infectivity varies both within and external to factors, including Ebola virus disease (EVD) which has generated a large epidemic in West Africa and other parts of the world [16, 17].The present paper is organised as follows. In the next section we propose two stochastic models: the first one is an endemic model in which each person infected from outside source contributes a random amount of infectivity to the force of infectivity (or infectivity function) which decreases at a constant rate over the time. The total infectivity present at any timet is the quantity of interest whose statistical characteristics are obtained. The second model which captures an epidemic has both endemic and epidemic components in the intensity function. The statistical characteristics of the number of infected persons are explicitly obtained. Section 3 presents a discussion on the developed models along with estimation of maximum likelihood parameters of the model.


## 2. The Models

In this section we propose an endemic and an epidemic stochastic model to describe the dynamics of a disease in a large population of susceptible individuals.


### 2.1. An Endemic Model

We assume that at the time of origin, there are infinitely many susceptible individuals. A susceptible may be infected at timet>0, after which it is classified as infected. It is assumed that the random variable denoting the time between two successive infections is exponentially distributed and that infections are due to a source of infection from outside of the population. Our interest is analysing the statistical characteristics of the total force of infectivity ft present in the population at any time t due to all the infections until time t.In this regard, the contribution to ftdue to an infection at time ti<t is composed of two components, a discrete jump ati and a nonstationary component ht-ti. Thus the force of infectivity at any time t can be written as follows:(3)ft=∑i=1Ntatiht-ti,where Ntis the number of infections until time t.From our assumptions,Ntis easily seen to be Poisson process. In order to derive the probability frequency functions of the Markov process ft we define(4)πy,tdy=Py<ft<y+dy∣f0=0.The Chapman-Kolmogorov backward equation yields(5)πy,t=1-λΔπy,t-Δ+λΔπy-aht,t-Δ+oΔ,where we have assumed ati=a.The above equation after simplification yields(6)∂πy,t∂t=-λπy,t+λπy-aht,twith initial condition(7)π0,t=1.Applying the Fourier transform(8)π~s,t=12π∫-∞∞πy,teisydyto both sides of (6), we obtain(9)∂π~s,t∂t=-λπ~s,t+λeisahtπ~s,t=-λ1-eisahtπ~s,t.The above partial differential equation with (7) can be solved to give the solution:(10)π~s,t=exp⁡-λt+λ∫0teisahudu.Now πy,t could easily be obtained by inverting (6)for specific forms of the function h·. In the foregoing model, we have assumed that the contribution of each infected individual to the force of infectivity ft is the same constant a. However in practice, this may not be so. Thus if we assume the jumps “a” to be a random variable with probability density pa, then πy,t in this case can be obtained with minor modifications to the above model. The backward equation (3) now reduces to(11)∂π∂t=-λπ+λ∫aπy-aht,tpadaand the corresponding equation after applying Fourier transform on both sides turns out to be(12)∂π~∂t=-λπ~+λ∫aeisahtπ~s,tpada.Now letting(13)ϕu=∫aeiaupadawhich is nothing but the characteristic function of the random variable a and using (13) in (12) we obtain(14)∂π~∂t=-λπ~+λπ~ϕsht.Thus we obtain(15)π~=exp⁡-λt+λ∫0tϕshτdτ.πy,t can be obtained by inverting (15) for specific forms of the function ht.


### 2.2. An Epidemic Model

In the earlier endemic model, it was assumed that the spread of the infection is due to external source only while each infection contributed a nonstationary input to the force of the infection. The present section develops an epidemic model which spreads due to an external source as well as through contagion of the infected population. Instead of the number of infected persons, we use the more realistic concept of the force of infection. Thus the number of infected at any timet is a point process whose conditional intensity function comprises both endemic and epidemic contributions. More specifically, we assume a conditional intensity function of the form(16)ft=λ+∑ibexp⁡-at-ti.The above form is motivated by the following interpretation. The terms on the RHS could be interpreted as follows: the intensity function ft which consists of two parts:(1)
The external rate of infection which occurs at a constant rateλ.(2)
The contagion rate due to the already infected persons consists of a constant jumpb at the time of infection which decays at a constant rate a so that the contribution to the effective contagion rate at later time t due to an infection at time ti is bexp⁡-at-ti. Thus the net contribution of the infection rate at time t due to all previous infections at t1,t2,…,tnis(17)∑ibexp⁡-at-ti.The rate could as well be interpreted as a rate at which the instantaneous jump b decreases due to intervention action.Letπy,tdy be the probability that ft has a value between y and y+dyat time t.Writing the forward equation,(18)πy,t+Δdy=πy-b;ty-bΔdy+1-yΔπy-aλ-yΔ;t+oΔ.The above equation is obtained by noting that the rate of decrease of the infection rate is specified by(19)dydt=aλ-y.Thus (18) reduces to(20)∂π∂t=y-bπy-b;t+a-yπy,t-aλ-y∂π∂y.Letting(21)Pn,t=∫0∞ynπy,tdy(20) reduces to(22)∂Pn,t∂t=anλPn-1,t-anPn,t+∑inibiPn-i+1,twith initial condition(23)P0,1=1,Pn,0=λn.Letting n=1 and solving (22), we obtain(24)P1,t=aλa-b-bλa-bexp⁡-a-bt.When b=0, so that the contribution of the second term of the contagion rate is not there and only external factor is present, we obtain P1,t=λand ENt=λt which is consistent with the assumptions of the model.From the solutionP1,t, we are in a position to obtain the expected number of infected persons up to time t. To see this we note that γ1t, the product density of order one of infections (i.e., probability an infection occurs between tand t+Δ), is straight away given by(25)γ1t=P1,t=∫0∞yπy,tdy.Thus the expected number of infections in 0,t is given by(26)ENt=∫0tγ1udu=1a-baλt-bλa-b1-e-a-bt.In a similar form, writing the equation satisfied by the conditional density function Πy2,t2∣y1,t1,λ representing the probability that ft has a value between y2 and y2+dy2 at time t2 given that it has a value between y1 and y1+dy1 at time t1 and initially λ, we obtain the second order product density γ2t1,t2.Then using the relation (27)ENt2=∫0tγ1t1dt1+2∫0t∫0tγ2t1,t2dt1dt2we obtain explicit expression for the above equation after some calculations, which we do not give as they are lengthy and unwieldy.


## 2.1. An Endemic Model

We assume that at the time of origin, there are infinitely many susceptible individuals. A susceptible may be infected at timet>0, after which it is classified as infected. It is assumed that the random variable denoting the time between two successive infections is exponentially distributed and that infections are due to a source of infection from outside of the population. Our interest is analysing the statistical characteristics of the total force of infectivity ft present in the population at any time t due to all the infections until time t.In this regard, the contribution to ftdue to an infection at time ti<t is composed of two components, a discrete jump ati and a nonstationary component ht-ti. Thus the force of infectivity at any time t can be written as follows:(3)ft=∑i=1Ntatiht-ti,where Ntis the number of infections until time t.From our assumptions,Ntis easily seen to be Poisson process. In order to derive the probability frequency functions of the Markov process ft we define(4)πy,tdy=Py<ft<y+dy∣f0=0.The Chapman-Kolmogorov backward equation yields(5)πy,t=1-λΔπy,t-Δ+λΔπy-aht,t-Δ+oΔ,where we have assumed ati=a.The above equation after simplification yields(6)∂πy,t∂t=-λπy,t+λπy-aht,twith initial condition(7)π0,t=1.Applying the Fourier transform(8)π~s,t=12π∫-∞∞πy,teisydyto both sides of (6), we obtain(9)∂π~s,t∂t=-λπ~s,t+λeisahtπ~s,t=-λ1-eisahtπ~s,t.The above partial differential equation with (7) can be solved to give the solution:(10)π~s,t=exp⁡-λt+λ∫0teisahudu.Now πy,t could easily be obtained by inverting (6)for specific forms of the function h·. In the foregoing model, we have assumed that the contribution of each infected individual to the force of infectivity ft is the same constant a. However in practice, this may not be so. Thus if we assume the jumps “a” to be a random variable with probability density pa, then πy,t in this case can be obtained with minor modifications to the above model. The backward equation (3) now reduces to(11)∂π∂t=-λπ+λ∫aπy-aht,tpadaand the corresponding equation after applying Fourier transform on both sides turns out to be(12)∂π~∂t=-λπ~+λ∫aeisahtπ~s,tpada.Now letting(13)ϕu=∫aeiaupadawhich is nothing but the characteristic function of the random variable a and using (13) in (12) we obtain(14)∂π~∂t=-λπ~+λπ~ϕsht.Thus we obtain(15)π~=exp⁡-λt+λ∫0tϕshτdτ.πy,t can be obtained by inverting (15) for specific forms of the function ht.


## 2.2. An Epidemic Model

In the earlier endemic model, it was assumed that the spread of the infection is due to external source only while each infection contributed a nonstationary input to the force of the infection. The present section develops an epidemic model which spreads due to an external source as well as through contagion of the infected population. Instead of the number of infected persons, we use the more realistic concept of the force of infection. Thus the number of infected at any timet is a point process whose conditional intensity function comprises both endemic and epidemic contributions. More specifically, we assume a conditional intensity function of the form(16)ft=λ+∑ibexp⁡-at-ti.The above form is motivated by the following interpretation. The terms on the RHS could be interpreted as follows: the intensity function ft which consists of two parts:(1)
The external rate of infection which occurs at a constant rateλ.(2)
The contagion rate due to the already infected persons consists of a constant jumpb at the time of infection which decays at a constant rate a so that the contribution to the effective contagion rate at later time t due to an infection at time ti is bexp⁡-at-ti. Thus the net contribution of the infection rate at time t due to all previous infections at t1,t2,…,tnis(17)∑ibexp⁡-at-ti.The rate could as well be interpreted as a rate at which the instantaneous jump b decreases due to intervention action.Letπy,tdy be the probability that ft has a value between y and y+dyat time t.Writing the forward equation,(18)πy,t+Δdy=πy-b;ty-bΔdy+1-yΔπy-aλ-yΔ;t+oΔ.The above equation is obtained by noting that the rate of decrease of the infection rate is specified by(19)dydt=aλ-y.Thus (18) reduces to(20)∂π∂t=y-bπy-b;t+a-yπy,t-aλ-y∂π∂y.Letting(21)Pn,t=∫0∞ynπy,tdy(20) reduces to(22)∂Pn,t∂t=anλPn-1,t-anPn,t+∑inibiPn-i+1,twith initial condition(23)P0,1=1,Pn,0=λn.Letting n=1 and solving (22), we obtain(24)P1,t=aλa-b-bλa-bexp⁡-a-bt.When b=0, so that the contribution of the second term of the contagion rate is not there and only external factor is present, we obtain P1,t=λand ENt=λt which is consistent with the assumptions of the model.From the solutionP1,t, we are in a position to obtain the expected number of infected persons up to time t. To see this we note that γ1t, the product density of order one of infections (i.e., probability an infection occurs between tand t+Δ), is straight away given by(25)γ1t=P1,t=∫0∞yπy,tdy.Thus the expected number of infections in 0,t is given by(26)ENt=∫0tγ1udu=1a-baλt-bλa-b1-e-a-bt.In a similar form, writing the equation satisfied by the conditional density function Πy2,t2∣y1,t1,λ representing the probability that ft has a value between y2 and y2+dy2 at time t2 given that it has a value between y1 and y1+dy1 at time t1 and initially λ, we obtain the second order product density γ2t1,t2.Then using the relation (27)ENt2=∫0tγ1t1dt1+2∫0t∫0tγ2t1,t2dt1dt2we obtain explicit expression for the above equation after some calculations, which we do not give as they are lengthy and unwieldy.


## 3. Discussion

We have proposed two stochastic models to analyse the dynamics of the spread of a disease. The core of the models is the force of the infection. Each of the infected individual contributes a nonstationary term to the total force of the infection. This comprises a discrete jump at the instant of infection which could be deterministic or random and a response function. The statistical characteristics of the number of infected persons are obtained from the conditional intensity function. The developed models subsume some interesting models available in the literature.In the epidemic model, lettingb=0 in (16) so that the contribution to the intensity function is only from external source at a rate λ, we note that the mean and variance of the number of infected persons are equal to λt showing that the infection grows according to Poisson process with rate λ. On the other hand letting a=0 so that the contribution to the intensity function from the external source and infected individuals are constants λand b, respectively, we are in the domain of a simple stochastic epidemic [13, 14]. In such a scenario, our model predicts (28)ENt=1bλebt-1and Gielen [14] shows the number of infected persons has a negative binomial distribution with the same mean. In this case with probability 1, there are ultimately infinitely many infective individuals. Also the mean number of infective individuals obtained using our model satisfies the renewal equation (2) of Keyfitz’s deterministic model [15].With the usual notations, the equations of the susceptible-infected-recovered (SIR) model are given by(29)dSdt=-βISN,dIdt=βISN-ϑI,dRdt=ϑI.These equations assume that infected people cause infections at rate β and recover at rate ϑ, giving an exponentially distributed infection duration. The process is obviously memoryless. In contrast, in real life situations, for diseases the “age” of an individual’s infection affects his/her infectiousness and consequently probability of recovery. Thus our modeling could be visualised as transcending the memoryless regime.Gielen [14], in his model, considered the intensity function of the form (30)ft=at+∑gt-ti.In such case, the total contribution to the infectivity function by the external source and by an infective during the entire infective period, respectively, is given by(31)∫0∞ftdt=F∞,∫0∞gtdt=G∞.These two quantities are dimensionless. If F∞=∞, then, with probability 1, there are finally infinitely many defectives which is the case in our model (see page 8 of the paper, [14]). However in Gielen’s model, it is possible that 0<F∞≕A<∞ and G∞=∞.In such a scenario, there is a nonzero probability that infection may not start at all so that the probability of infinite infective individuals ultimately is defective with probability 1-e-A. In our model we start with f0=λ which preclude such an eventuality.In our endemic model, lettinga=1and ht=δt-1 where δ is the Dirac delta function, we note that each infection increases the force of the infectivity by one unit. This is reflected in the final solution as we get(32)Πn,t=Pft=n=e-λtλtnn!,n=0,1,….Thus we arrive at a Poisson process of infections. In the endemic model, the contribution to the force of infectivity arises only from the external source. If we further assume that infections also occur due to the already infected persons but such infected individuals are removed from the pool of infected, then the probability of an infection due to the infected population in t,t+dt is given by(33)∫0∞yΠy,tdydt.This can be identified with the product density of order one, namely, γ1t. Hence the expected number of infected individuals at time t is (34)∫0tγ1uduwhile the expected number of infected persons due to external source is λt.Finally, we wish to mention the useful tools that a planner needs in designing interventions is the estimated values of the model parameters. In order to estimate the parameters of our model,a and b, we need the function (35)γnt1,t2,…,tn;a,bwhich is the product density of order n.This function which can serve as the likelihood function specifies the probability that infections occur at t1,t2,…,tn under the assumed values of a and b. Although such a function could be obtained, it may not be of much value in actual use due to its unwieldy form. Thus we make the reasonable approximation (36)γnt1,t2,…,tn;a,b=γ1t1;a,bγ1t2-t1;a,b⋯γ1tn-tn-1;a,b.The estimate of the parameters a and b yields the maximum likelihood parameters a^ and b^, by using methods proposed by Rubin [18] and Ozaki [19].This paper generalizes the stochastic version of Kermack and McKendrick model of Gielen [14] by considering the transient behaviour of the total infectivity present in an infinite population of susceptible individuals. A noteworthy aspect is the explicit analytical expressions for the probability distribution and statistical characteristics of the number of infected individuals. The proposed models can also be solved using Lie algebra approach for finding analytical solutions. The proposed model could be useful in order to clarify theoretical concepts on infectiousness, the force of the infection, in order to compare the infectiousness of different communicable diseases such as Ebola virus disease. The future work is proposed to carry out specific statistical analysis and estimation of parameters using real-life data for epidemic disease.



---

*Source: 9064821-2017-01-11.xml*