1.Introduction

Epidemic spread in human populations is a complex phenomenon, throughout human history infectious diseases have caused debilitation and premature death to large portions of the human population, leading to serious social-economic concerns. The bubonic plague, smallpox, dengue, ebola, yellow fever, influenza, malaria, polio, avian flu, leprosy, tuberculosis, whooping cough, severe acute respiratory syndrome (SARS), among others, are examples of the most devastating ones (see ¹). Although, some viruses and bacterias causing these diseases have been completely eradicated others have only been possible to be controlled and are a health problem yet. A recent case took place in the late March and early April 2009, when reports of respiratory hospitalizations and deaths among young adults in Mexico alerted local health officials to the occurrence of atypical rates of respiratory illness at a time when influenza was not expected to reach epidemic levels. Infections with novel swine-origin influenza AH1N1 virus were confirmed in California, (United States), on April 21 (see ²) and in Mexico on April 23 (see ³).

The formulation of theories on the nature of infections diseases dates back to ancient times; they have evolved and currently modern theories of differential equations have been widely adopted in the analyses, and simulation of biological populations including the evolution of epidemics. The first mathematical model to describe an epidemic was developed by D. Bernoulli which was applied to smallpox. Later, W. Heaton formulated a model to analyze the propagation of measles and R. Ross developed a work on malaria transmission. On the basis of a scheme with three elementary compartments: susceptible, infected, and removed, a mathematical model more general commonly known as SIR, was developed by ⁴ for describing any disease of viral kind.

Most epidemic models are based on dividing the host population into a small number of compartments, each containing individuals that are identical in terms of their status with respect to the disease in question. For instance, the standard approaches are the so-called susceptible-infected-removed (SIR) and susceptible-exposed-infected-removed (SEIR) models have contributed much in the present for understanding the nature of recurrent epidemics. Traditionally, deterministic models have been the base of mathematical epidemiology, however, stochasticity has been recognized as an important tool within epidemic modelling, but there is still much concerning its precise dynamical role and how it can be understood from a theoretical point of view. Nevertheless, stochastic models have recently been considered as a more realistic descriptions of epidemics. This is evidenced by the large amount of systematic treatments of the subject, for instance, the first stochastic simulation on epidemiological system, was implemented by Bartlett, who was interested in extinction dynamics (see ⁵). In more recent studies,⁶ , have applied the well known Van Kampen method to the SEIR model with distributed exposed and infectious periods for studying whooping cough. Moreover, McKane et al. developed a method for fluctuations around cycles of forced or unforced systems, a method which was extended by ⁷ to a seasonally forced SEIR model in a realistic parameter region for childhood infectious diseases where different attractors exist or coexist. According to studies on avian flu ⁸, its is argued that a fully stochastic theory based on environmental transmission provides a simple, plausible explanation for the phenomenon of multi-year periodic outbreaks. On the other hand, a stochastic theory for the major dynamical transitions in epidemics from regular to irregular cycles, which relies on the discrete nature of disease transmission and low spatial coupling was proposed by ⁹, while a study on the estimation the parameters in dengue fever case studies was performed by ¹⁰. Among the studies on epidemics in Mexico, we can mention, for instance, the study performed by ¹¹, who used a time dependent modification of the ⁴ model to study the evolution of the influenza AH1N1 epidemic reported in the Mexico City area under the control measures used during April and May 2009. Meanwhile, ¹² analyzed the epidemiological patterns of the pandemic during April-December 2009 in Mexico and evaluated the impact of nonmedical interventions, school cycles, and demographic factors on influenza transmission.

This work is organized as follows. In order to establish the terminology and nomenclature used in our analysis, in Sec. 2 we summarize a deterministic mathematical model SIR and we give a brief description of the general master equation. An application of the stochastic formalism that we have developed on a realistic case is presented in Sec. 3. In Sec. 4, we present an analysis where the important parameters and numerical comparison are discussed. Finally in Sec. 5 we address our concluding remarks and conclusions.

2.The SIR model

In Mathematical Epidemiology, one of the most widely used models still is the simple SIR model developed in Ref. ⁴. In Subsec. 2.3, we will present a stochastic extension of such model. For context, we will briefly review the deterministic SIR model. For more details on this mathematical formulation, Ref. ¹ can be consulted.

2.4.One step master equation

In what follows, the description done so far on the master equation, will be limited only to the three population groups: susceptible, infected, and removed, as was discussed above are denoted by S, I, and R, respectively. We also will consider that an individual only has interaction with their nearest neighbour.

We suppose that the transition rule of the population members between the groups S, I, R is as follows: The susceptible individuals, initially healthy, may become infected with a rate ΒSI/N, whereas the infected individuals, have the disease and are able to transmit it to other individuals at a rate γI, the removed individuals are immune to the disease after recovery. For simplicity in this model, it is assumed that during the epidemic, the virus is transmitted only between close neighbours with the possibility that an infected individual may be removed. Figure 1 illustrates schematically the transition flows from one group to another.

Figure 1 Flow chart for the SIR model. 

The elementary processes occurring into the system can be divided in two elementary groups; infection and recovery, which are defined as follows: An infection process implies that a susceptible individual enters in contact with an infected individual, generating therefore a newly infected individual. The recovery process signifies that an infected individual who is able to recover from the disease, cannot be susceptible again, because he has acquired immunity. In summary, the transition rates expressed in terms of the SIR model variables are defined by the following events,

{S,I}⟶{S-1,I+1}{S,I}⟶{S,I-1}, (3)

which are represented schematically in Fig. 2 together with their respective transition rates. We take just two variables n={S,I}, because the variable 𝑅 can be completely determined by the restriction N = S + I + R.. In this way Eq.(2), can be written explicitly as

dP(S,I,t)dt=TS,IS+1,IPS+1,I-1,t+TS,IS,I+1PS,I+1,t-TS-1,I+1S,IPS,I,t-T()(S,I-1|S,I)P(S,I,t) (4)

Figure 2 Processes of interaction between individuals. 

Equation (4) is a master equation describing a one-step process, whose corresponding associated transition rates we shall assume by proposing the following expressions,

TS-1,I+1∣S,I=βSIN, (5)

TS,I-1∣S,I=γI. (6)

Here, T (S -1,I +1 | S,I) measures the probability rate for an infected individual to become susceptible I, while the number of removed R is constant. We take this quantity to be proportional to the temporal rate of change of S provided by the Kermack McKendrick model, as given by the first equation of Eq. (1), which is conceptually resembling. Similarly, T(S,I -1 | S,I) gives the transition probability at which the number of infected individuals increase by one while the susceptible ones is kept constant. Because the total number is constant, this is related with the reduction of the recovery individuals R whose temporal ratio is provided by the last equation of Eq. (1). We identify both quantities as proportional to construct our master equation. For convenience, the one-step operator 𝐸, is defined through its effect on an arbitrary function f(m) as E±1f(m)=f(m±1), thus in terms of the step operator, the probability transitions given by Eqs. (5) and (6), can be rewritten as follows

E±1fm=TS,IS+1,I-1PS+1,I-1=ES+1EI-1βSINP(S,I,t) (7)

and

E±1fm=TS,IS,I+1PS,I+1,t=EI+1γIPS,I,t. (8)

From transition rates represented by expressions (7) and (8), Eq. (4) can be written in a compact form as

dPdt=ES+1EI-1-1,βSIN+EI+1-1γIP. (9)

In order to solve Eq. (9), we use the method that has been proposed by ¹³. This method suppose that the probability distribution P behaves like a sharp Gaussian statistical distribution, where the random variables are rewritten as the superposition of a macroscopic deterministic variable and a mesoscopic random variable. According to this formulation, the following ansatz for the independent variables are suggested, namely,

S(t)=N(ϕ(t)+qx) (10)

I(t)=N(ψ(t)+qy), (11)

where ϕ and ψ are the macroscopic variables which represent both to susceptible individuals and infected individuals, respectively, and x and y are small stochastic corrections to these variables, q=1/N is a smallness parameter. The Van Kampen approximation allows us to find the deterministic evolution of Φ(t) and ψ(t) as well as the two first momenta of the distribution P in terms of the corrections x and y. If we apply the step operators ES±1 and EI±1 to expressions (10) and (11), we have to replace x → x + q and y → y + q, so the resulting expressions ES±1(x+q) and EI±1(y+q) are expanded in Taylor’s series,

ES±1=1±q∂∂x+q22∂2∂x2+O(q3), (12)

EI±1=1±q∂∂y+q22∂2∂y2+O(q3). (13)

Inserting expressions (10-13) into Eq.9, and multiplying the hand right side by NP, it becomes

β1+q∂∂x+q22∂2∂x21-q∂∂y+q22∂2∂y2-1×(ϕ+qx)(ψ+qy)+γq∂∂y+q2∂2∂y2(ψ+qy). (14)

Notice that the expression between square brackets up to second order in q² can be approximated as

1+q∂∂x+q22∂2∂x21-q∂∂y+q22∂2∂y2-1=q∂∂x-∂∂y+q22∂2∂x2+∂2∂y2-2∂2∂x∂y+O(q3). (15)

Thus, if we expand the other terms and group the resulting ones in powers of q we can obtain

dPdt=1qβϕψ∂P∂x+γψ-βψϕ∂P∂y+∂∂x-∂∂y×xψ+ϕyP+ψϕ2∂2∂x2+∂2∂y2-2∂2∂x∂yP+γ2ψ∂2P∂y2+γ∂∂yyP, (16)

where the terms of order q have been neglected. It is important to remark that Eq. (16) is the Fokker-Planck’s equation (see ¹³). Moreover, it is useful to express the probability distribution in terms of S and I, which in turn can be written as a function of the continuous variables x and y, namely, P(S,I,t)=Π(x,y,t), it yields

dPdt=∂Π∂t-1qdϕdt∂Π∂x-1qdψdt∂Π∂y. (17)

By comparing the right hand side of the Eq. (16) and that of the Eq. (17), it can be infered by equating the lowest power in q that

dϕdt=-βϕψ, (18)

dψdt=βϕψ-γψ. (19)

Notice that the system of Eqs. (18-19), has the same functional form that SIR model described by Eq. (1), whose solutions ϕ=E[S] and ψ=E[I], are the expected value of the Eqs. (10-11). On the other hand, up to zero order q the following equation is fulfiled.

∂Π∂t=β∂∂x-∂∂yxψ+ϕyΠ+βψϕ2∂2∂x2+∂2∂y2-2∂2∂x∂yΠ+γ2ψ∂2Π∂y2+∂∂yyΠ (20)

which can be expressed as a continuity balance equation for the probability ∏ given by

∂Π∂t+∇⋅J=0, (21)

where the two dimensional probability current vector J is defined as

J=β1-1(xψ+ϕy)Π+βψϕ21-1-11∂Π∂x∂Π∂y+01γ2+∂Π∂y+γyΠ. (22)

We remark that Eq.(21) guaratees the conservation of probability.

To analyze the dynamics of the expected values of the ramdom variables x and y we multiply by x or by y Eq.(20) and integrate over the whole domain. It yields

d<x>dt=-β<xψ+ϕy>, (23)

d<y>dt=β<xψ+ϕy-γy>, (24)

where as usual the expected value is given by

<u>=∫-∞∞∫∞∞uΠ(x,y)dxdy. (25)

Here u is a variable which can be a function of x and y. By adding Eqs. (23) and (24) we get the expression

ddt<x>+<y>-γy=0, (26)

which allows us to find a time invariant of the system. To find the dynamics of the second momenta, we proceed similarly and integrate by parts the resulting expressions by assuming that the density of probability ∏ and its partial derivatives with respect to x and y vanish at x=±∞. It leads to

d<x2>dt=-2β<xxψ+ϕy>βψϕ, (27)

d<y2>dt=2β<yxψ+ϕy>+βψϕ+yψ-2γ<y2> (28)

d<xy>dt=β<y-xxψ+ϕy>-γ<xy>-βψϕ (29)

After adding the first and second expressions and substracting the third one we obtain following relation

ddt<(x-y)2>=γψ-2γ<y(x+y)>, (30)

It is convenient to stress that the model described by the set of Eqs.(23), (24), (27), (28), and (29) is general and therefore, it can be applied to any epidemic of viral transmission fulfilling the restriction of having a narrow statistical distribution; it means that he solution P(S,I,t) of Eq. (9) obeys a Gaussian distribution with averages Φ(t) and ψ(t), whose fluctuations bands ϕ(t)±σS and ψ(t)±σI, are given by the solutions of the Eqs.(23), (24), (27), (28), and (29). Here the respective standard deviations of each variable are defined as usual: σS2=E[x2]-E[x]2 and σI2=E[y2]-E[y]2.

3.Application of this model to study a realistic case

Before applying our stochastic formulation, it is good to mention that although the behavior of the groups S and R can be predicted by using a deterministic model, in contrast with infected group, its quantification for the whole population is very difficult, and, in consequence, there is not enough field data. The infected group is considered as the variable with more impact and as such, it is therefore the most commonly used in epidemiological predictions. For these reasons, we focus our discussion mainly in infected cases I.

In what follows, we apply the mathematical model developed in Sec. 2.1 and summarized by Eqs. (23)-(29) to analize the influenza AH1N1 spread in Mexico in 2009. We use real data reported by Mexican Institute for Social Security (IMSS) in that year. In Table I it is presented the total number of infected people, this information includes 35 geographical regions stipulated by the . We use the following methodology: even though it is well known that one epidemiological year encompasses a whole period of 52 weeks, of which we only consider those weeks in which there was a significant number of infected cases N₀. In regard to the proportion of the susceptible members S(t), in the relation N(t) = S(t) + I(t) + R(t), it is taken those individuals who finally became infected, which implies the initial condition I(t) = S(0) . In addition, to estimate the order of magnitude of the parameters β and γ we apply the Levenberg-Marquardt algorithm, whereas the value R0 it is computed directly from the ratio β/ γ. In order to complete the problem formulation, we have also imposed the initial condition S(0) = 1, which implies that final epidemic size for the deterministic SIR epidemic model can be computed from the implicit relation

Hence, the total number infected during the outbreak is I=N-S(∞) and the above expression can be expressed as

this implicit relation is used to calculate the number of susceptible people, by using the Generalized Reduced Gradient (GRG) nonlinear algorithm. Also, for the other expected variables we assume the following initial conditions: E[x=0]=E[y=0]=E[x2]=E[y2]=E[xy]=0. Under these conditions, we have solved numerically the set of equations given by Eqs.(23), (24), (27), (28), and (29) by using a standard Runge-Kutta algorithm (see ¹⁵).

4.Results and discussion

In Table I, we present the corresponding values of the typical parameters β, γ, and R0 that we have estimated. The main results for parameter R0 are as follow: we have estimated the parameter magnitude R0=1.05 and R0=1.08 for the geographical region Distrito Federal Norte and Distrito Federal Sur respectively, whereas a value R0=1.72 has been estimated by for the same geographical region so called Distrito Federal (include both regions Distrito Federal Norte and Distrito Federal Sur). On the other hand, for the whole country, we have estimated the value R0=1.02, while Fraser et al. estimated this value within the interval (1.4,1.6); thus, in a similar analysis done for members belonging IMSS, estimated season values to R0 between (1.8 - 2.1) during Spring, of (1.6 -1,9) for Summer and (1.2 - 1.3) for Fall. We must notice that the order of magnitude of the parameter R0 for all the geographical regions is very close to the value that has been estimated by those authors, who have used different techniques than those we have discussed here.

Figures 3-8 exhibit all the events that we have analyzed. The corresponding comparison between real infected cases I and its respective theoretical average E[I] has been presented for each geographical region. In the same Figs. 3-8, also is illustrated the corresponding associated curves with the extinction probability of the epidemic as a function of the time.

Figure 3 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 4 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 5 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 6 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 7 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 8 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 9 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 10 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

The results show that the best fit correspond to geographical regions Distrito Federal Sur, Guanajuato, Hidalgo, Michoacan, Estado de Mexico Poniente, Tlaxcala and Sonora, because real data are accumulated around the expectation value E[I]. Moreover, the geographical regions Aguascalientes, Baja California Sur, Campeche, Coahuila, Chihuahua, Durango, Guerrero, Morelos, Puebla, Sinaloa, Veracruz Norte, Zacatecas, Nayarit, and San Luis Potosi, also exhibit a good fit with real data more scattered around its theoretical average. However all of them are still within the calculated probabilistic bands given by E[I]±σI. Similar behavior around the expected value E[I] is observed in the regions Jalisco, Quintana Roo, Tabasco, and Tamaulipas, but some real cases are slightly out of the bands E[I]±σI.

Figure 11 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 12 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 13 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value. 

Figure 14 a) Probability distribution theoretical and observed. b) Ep=Extinction probability. Here ψ = E[I] is the expected value of the theoretical distribution and σ _I is the standard deviation around the expected value 

In contrast, the results show that the largest spread of real cases around the theoretical average E[I] correspond to geographical regions Chiapas, Colima, Estado de Mexico Oriente, Yucatan, Baja California, Distrito Federal Norte, Nuevo Leon, Oaxaca, Queretaro, and Veracruz Sur as well as for the whole country. It can seen that in some of these regions and in the whole country, a considerable number of real cases are outside of the theoretical standard deviation bands E[I]±σI.

On the other hand, it is also important to point out that in the geographical regions Puebla, Sinaloa, Colima, Jalisco, Estado de Mexico Oriente, Nayarit, San Luis Potosi, Yucatan, Distrito Federal Norte, Nuevo Leon, Oaxaca, Queretaro, and Veracruz Sur an oscillatory trend is observed particularly notorious at the level of the entire country. The origin of these discrepancies probably, is owing to two main reasons: 1) the simplicity of the model that we have utilized and / or 2) inaccuracies in the acquisition of field data. A possible solution in order to have a better fit, could be the implementation of a more realistic model that takes into account infected death, vaccination, births, population migration, incubation period of the virus, among other factors which we have neglected. Moreover, the oscillatory behavior observed, characteristic in large populations, suggests the necessity of introducing seasonal forcing into the model. In parallel, suitably trained staff would represent a direct improvement on the quality and reliability of the data.

In regard to extinction probability as a function of time, we have identified the week for which there is a 90% of extinction probability of the epidemic -a threshold that has been chosen arbitrarily and could easily be changed-. Under this criteria, the average extinction time for each geographical region as well as for whole country was around 13 and 30 weeks respectively. Moreover, geographical regions in which it is observed a longer period of extinction corresponding to 90% are, Colima, Nayarit, and Queretaro with 19 weeks, Jalisco and Baja California with 20 weeks and Oaxaca with 21 weeks. Complementary, shorter periods were found for geographical regions Quintana Roo and Tabasco with 4 weeks.

As already mentioned, these estimates values could be improved with the implementation of a more realistic model and a more reliable registration of infected cases. However, a lot of the neglected processes and unrevealed variables can be incorporated to some extent, within the standard deviation band that we have computed. This is indeed the power of our stochastic approach because even though there exist various additional variables, not considered here, which modify our epidemic evolution, we can embody all of them within the statistical distribution of the parameters of a quite simple dynamic model which describes the dominant behavior of the phenomenon.

5.Concluding remarks

We have established a stochastic model based on a simple compartmental SIR formalism. Using this approach we have studied the propagation of the AH1N1 influenza that took place in Mexico in 2009. We have restricted this study to

fected cases reported by Mexican Institute for Social Security (IMSS) in that year.

Our main conclusions are the following: the order of magnitude estimated for the basic reproductive number of infection R0 is within the interval (1.03 - 1.23) which is shorther than those previously has been reported in the literature. The model that we have developed is quite general and therefore is not only applicable to influenza viruses, but it can also be applied to any epidemic in which the following conditions are fulfilled: the infectious agent is a virus, there are not demographic dynamics, the population is homogeneous, the seasonal period is negligible, removed individuals may not be susceptible again, the incubation time of the virus is not taken into account, any two individuals have the same probability of contacting each other, and the duration of the epidemic is short compared with the life expectancy of the individuals. This type of study, applied to real cases, the first of its type in Mexico, as far as we know. In spite of the simplicity of the model that we have implemented, the results show clearly that the probability distribution predicted by the stochastic method for the infected individuals is in good agreement with the distribution obtained from the dataset provided by IMSS. This is very important because provides a solid confidence in our formalism for working in subsequent developments. Thus, the results of this research may have immediate practical application because they constitute a tool for taking decisions and for the optimization of resources against possible epidemic scenarios. The small discrepancies between the forecast and observed data in some geographical regions, might have their origin in: 1) the simplicity of the model implemented; in this case, a more complete formulation is needed in order to provide a better description with respect to the real data. For instance, the implementation of a model with seasonal forcing would be suitable for those regionswhich do not evidently exhibit a Gaussian distribution around the average and 2) data reliability. It is advisable to implement strategies that contribute to train staff effectively; hence, achieving a better quality on the information gathered.

As future work, it will be interesting to formally extend the present study by including the main health agencies that have a presence throughout the country, such as the Mexican Secretariat of Health (SSA), Institute for Social Security and Services for State Workers (ISSSTE) and Navy Secretariat (SEMAR). The data set coming from these institutions would complement the information provided by the Mexican Institute for Social Security (IMSS) for a nationwide study. We can even consider other kind of epidemics, as well as an extension to our model by means of incorporating important variables such as: vaccination, demographic dynamics (emigration, immigration, births, and deaths), information on passive immunity, age and sex of the population, quarantine period, incubation period of the infectious agent, loss of immunity, forcing seasonal, among other factors. It would also be interesting to consider an interaction between individuals beyond their immediate neighbors, that is, second and third neighbors. Nevertheless, we should be aware of including only a few additional important variables in order to keep a very simple and easy-to-use description.