Stochastic epidemiological model: Simulations of the SARS-CoV-2 spreading in Mexico

In this paper we model the spreading of the SARS-CoV-2 in Mexico by introducing a new stochastic approximation constructed from first principles, where the number of new infected individuals caused by a single infectious individual per unit time (a day), is a random variable of a time-dependent Poisson distribution. The model, structured on the basis of a Latent-Infectious-(Recovered or Deceased) (LI(RD)) compartmental approximation together with a modulation of the mean number of new infections (the Poisson parameters), provides a good tool to study theoretical and real scenarios.


Introduction
Since the late 2019 to date, the rapid worldwide spread of the SARS-CoV-2 has caused around four and a half million of human deaths [1], placing mankind in one of the most challenging episodes in the recent human history.An extraordinary effort has been made to implement mathematical methods to accurately describe the spreading of the epidemic, looking to forecast and to implement non-pharmaceutical responses to reduce the damage in the society [2].These methods, ranging from standard compartmental models (typically employed to determine the initial epidemiological parameters [3][4][5][6][7][8]), to hybrid methods that incorporates stochastic meta-population network models with local and global mobility patterns [9][10][11][12][13], attempt to overcome the complex behavior of social interaction characterized by the tendency of the population to cluster [14], following quasi-periodic patterns of mobility in large dense urbanized areas [12,13,15].Furthermore, in addition to the complexity for determining the degree of connectivity among individuals (the contact network), regulatory measures such as home lockdown and social distancing were promoted to reduce the transmission of the infection [1], providing an additional degree of complexity in determining the spreading of the disease.
In this regard, how and when to promote regulatory measures became one of the most difficult decisions to follow, because of their effects in public health, the economy and several other social factors.These decisions had to be supported in predictive models possessing a good equilibrium between registered data (reliable readouts about registered confirmed cases), mobility patterns followed by the population, and computational efficiency of the epidemiological models [8,11,13,[16][17][18][19].
Moreover, the efficiency of a given model rely on a good accessibility and characterization of the available data [13], which could be more difficult or impossible to be obtained in the case of less developed countries; in this context, stochastic models, which introduce a randomization about certain unknown elements could provide an alternative guidance.
Recently, some stochastic models have been employed to study the Sars-CoV-2 spreading [20][21][22]; however these models still rely on the law of mass interaction governing the probabilities of infection, an assumption that may not be fulfilled when the dispersion of the disease happens in highly structured social networks following confinement measures; e.g., in [20,21] a compartmental description is used as base model and then additive white Gaussian noise is introduced in the contact parameter β; other approximations are carried out by considering a master equation following transition probabilities which also rely on the standard SIR model dynamics and hence a probability of infection proportional to the infected and susceptible populations.In this paper, we introduce a new stochastic model which attempts to overcome the law of mass interaction underlying in the traditional compartmental models.The model has served us to simulate and follow the spreading of the Sars-CoV-2 in Mexico with a good agreement to the real cases; it is structured on the basis of a LI(RD) compartmental model (Latent-Infectious-(Recovered or Deceased)) where the number of infections caused by a single infected individual, per unit time (a day) is randomized, while the daily mean of infected population is modulated through a weight-like time dependent function.The modulation help us to introduce tendencies in the mean of the daily infections caused by several phenomenological or fundamental behavior such as pharmaceutical or non-pharmaceutical interventions and herd immunity as well.
Finally, through this model we analyze the evolution of the disease in some Mexican states (some of them housing the largest metropolitan areas of Mexico), by deriving an empirical approximation of the weight function which is in turn deeply connected to the effective reproduction number RðtÞ.

The model
The epidemiological model we propose attempts to describe a scenario about how many people can infect one infectious individual per day when the infectious events are considered homogeneously distributed in time and when the probability of infection is affected by pharmaceutical or non-pharmaceutical interventions; it consists of the randomization of the daily number of infections using a time-dependent Poisson processes to generate the new infections caused by each of the infectious individuals, along a given period of time (the time unit).The core of the model is constructed on the basis of a compartmental description: the susceptible population S(t) which serves merely to have a finite resource about the number of new infections, the infected-latent population L(t) which is randomly obtained and the infectious population I(t) which represents the the part of the population capable of infecting.Once the number of infections per infected individual at a single time step (i.e., the daily infected (but not infectious) population per infected individual) has been obtained, they are removed from the susceptible condition an placed into the latent condition L(t) which characterizes the part of the population that is infected but is not capable of transmitting the virus until an incubation period or latency period τ L , has passed.After the latency period, the infected-latent population becomes contagious, passing into the infectious condition I(t) and hence becoming able to transmit the disease to the susceptible population by associating to each member of this group a new random process of infection.A schematic representation of the stochastic model is given in Fig 1 ; in the figure, an infectious population number n will generate a set of random numbers following a Poisson distribution associated with infectious individuals and which are summed to the infected population.The number of the latent and the infectious population at time t j+1 may be written as follows: where τ L and τ I represents the period of latency and the period of infectivity, the former defined as the period of time in which an infectious person is capable of transmitting the disease, also, the theta functions θ(�) are placed in the equations to start counting or removing individuals from the different categories after the latency or infectivity periods have passed.
The set of random variables {χ(t j )} i gives the amount of new infected individuals due to the ith infectious individual at time t j (which with absolute certainty become infected).The random variables are obtained from a Poisson process, i.e.: where the intensities of the Poisson process, (i.e. the parameters λ(t j )), describe the mean number of contagious events at the time t j .In a real scenario, the spreading of a disease depends on the degree of close contact among the individuals and therefore, on the degree of urbanization and mobility of the population [23]; however, we believe that part of these complex aspects could be captured into our model by a proper parametrization of the Poisson processes, i.e., the daily mean number of infections per infectious individual λ i (t j ).In this regard, we introduce a time-dependent function which indirectly serves to modify the mean number of the daily infections by associating at each time t the following mean of the number of infections produced daily: where λ i represents individual rates of contagious which could be represented as additional random variables assigned to the i-th infectious individual following a probability distribution P(% o ), i.e., {λ} i P(% o ) with parameter % o representing an initial estimation about the average of the number of infections that a single infectious individual can cause per unit time, i.e., the ratio between the basic reproduction number (calculated at the beginning of an epidemiological event [5,24,25]) and the infectious period: additionally, W(t) in 4 is a time-dependent function (the weight function) which serves to modulate the mean of the number of infections.In other words, the stochasticity of λ i attempts to simulate the contact patterns followed by the different individuals of the population given a specific mean of contacts in the whole population while the weight function W(t) describes modulations about the probability of getting infected, e.g., herd immunity, social distancing or confinement.Along this paper we will focus on the case where λ i follows a punctual distribution, i.e.P(% o ) ! % o , and we will address the employment of W(t) in the following subsection.Finally, and following within the compartment direction, we consider that the infectious population could pass, either to the recovered R(t) or to the deceased D(t) condition, depending on the development of the disease in the infected individual (see Fig 1).In the former, we define the recovery period τ R , after which the infected population heals with a given probability of recovering 1 − p d , while the deceased population is the part of the infected population which does not heal according to the decease probability p d .For implementing this procedure we make use of an additional random procedure to randomly select from the infected-latent individuals and according a given fatality rate l, the infected population that will pass into the deceased category at the time t j i.e., we count every new set of infected-latent individuals appearing at time t j+1 and for each of the new cases, we use a uniform distribution to generate a random number r 2 unif(0, 1) which is compared to p = 1 − l and if r > p, we then remove in the future time t j+1 + τ L + τ I = t j+1 + τ R this individual from the infectious condition and place him into the deceased condition D(t).In this sense, the number of recovered and deceased population at the time t j+1 is given by: Along this paper, we will simulate the evolution of a disease possessing similar epidemiological parameters to those of the COVID-19.We use a basic reproduction number of R o = 4, an incubation period τ L = 4 days and an infectivity period of τ I = τ R − τ L = 14 [6,9,10,[24][25][26][27].
Modulation of the Poisson parameters.Without the inclusion of the weight function W (t), the model we propose represents a probabilistic model with replacements, i.e., the probability of infecting a certain amount of susceptible per infected individual would be only determined by a stationary given value, independently of the total population being infected or the contact network.Nevertheless, an intuitive behavior is that as the population of susceptible decreases, then also the chances of having large number of susceptible to fall into close contact with the infectious population; in fact, this is exactly the underlying idea in the emergence of the herd immunity effect.On the other hand, the probability of infection is also continuously changing when contingency measures such as social distancing and home lockdown are implemented in the population.In this regard, an appropriate functionality of the weight function could help us to incorporate such effects.
In the context of an exemplification, we make use of a weight function being the product of a function characterizing the herd immunity effect with an additional function characterizing the variability of the probability of infection along the transients of the epidemiological dynamics due to confinement and other related contingency measures, i.e.:

WðtÞ ¼ HðtÞ CðtÞ ð8Þ
where H(t) represents the herd immunity effect which is estimated to emerge when a large proportion of the population (but not all), has gain certain immunity [28,29], while C(t) represents additional changes in the probability of infection along the transients of the dynamics.
For the herd immunity effect, we employ a reversed logistic-like function whose argument depends on the fraction of the population that has become infected along the evolution of the epidemic, i.e.: where c N ðt j Þ ¼ P t j t¼t o P IðtÞ i¼1 w i ðtÞ=N, is the fraction of the cumulative infected population (Latent and Infectious) at time t j ; α is a free parameter serving to adjust the stationary value of the infection in the long time limit, and a is a lower bound at which the probability of infection is reduced sufficiently to reach the stationarity.In our simulations, we set a = 0.1, whereas we have seen that by choosing α = 0.22 the herd immunity is achieved when close to the 80% of the population has been infected [29].
In Fig 2 we present the effect of the weight function to the epidemiological variables when it only characterizes the herd immunity effect (i.e.W(t) = H(t), C(t) = 1).In the panels at the left we present single realizations of the latent, the infected, the recovered the deceased, the cumulative of the incidence (top left panel) and the incidence (bottom left panel), together with the form of the weight function generating the herd immunity effect.At the right panel, the normalized incidence is presented for different population sizes and when averaged over 5000 realizations.
The parameters employed in From Fig 2, one notices that the maximum incidence for a population of one million is obtained from around 3 months after the beginning of the disease spread, reaching at its maximum an amount of roughly 2.5% to 2.7% of the total population.Additionally, if the population is increased in size by one order of magnitude, the maximum is shifted around 20 to 30 days when no contingency measures are implemented in the population.
The effect of non-pharmaceutical strategies to contain the spreading of a disease is another way to modify the probability of infection.In this regard, one could think that some of the most common or intuitive responses of the population under an epidemiological risk: a confinement responding to the daily experience about the development of the disease, e.g., a confinement depending upon the number of active cases.In other words, when a certain fraction of the population has become symptomatic-infected (or deceased), it is more likely that some of the susceptible population has knowledge about infected individuals in their social circles or in the neighboring community, reacting with lockdown due to the fear of becoming infected.Another possibility is that contingency measures are placed over the population (typically by health authorities) along different stages of the evolution of the disease, attempting (primarily) to find an equilibrium between the public health resources and different economical activities that require contact among the population.In this case and as we have experienced with the COVID-19 pandemic, all populations have gone through lockdown and relaxation of the confinements during different stages, which in turn, can be imposed at any time of the epidemiological development by the health authorities.
We use our stochastic model to explore the behavior of the COVID-19 spreading in two different confinement scenarios, a confinement triggered upon the number of infectious population and an idealized confinement regulated by health authorities at different stages of the dispersion of the disease.In the former, we let C(t) to be Gaussian-decaying function of the active cases (the infectious population), triggered once certain part of the population has become infectious-symptomatic or deceased, i.e.: where I(t) are the infectious cases at time t, while I o represents a threshold about the amount of infectious population at which the confinement function is triggered; γ is a decaying-rate parameter describing how strong is the confinement and N is the total population.Fig 3 shows the evolution of the disease for a confinement following a Gaussian decay as described in (10).
The left four frames represent I o = 1% and I o = 10% (from top to bottom) while from the left to the right, γ = 5 and γ = 10.
In the figure one can see that the outcomes of a confinement relying on the number of the infectious population depend on how strong and rigorous is the confinement and at what stage of the dispersion of the disease is implemented.The different outcomes go from a flattening of the epidemic curve, happening when the confinement does not happens abruptly, to revivals in the incidence which become periodic and more pronounced when the confinement is strong and happens at earlier stages of the epidemic.At the right panel we have plotted the latent and the infectious population under an abrupt confinement (γ = 50) at an early stage of infection (I o = 1%), from which several revivals can be seen.These revivals can be explained by looking at the curve of the latent population: in an abrupt confinement, large part of the population remains on the latent condition and when the number of infectious population is reduced, the latent population will tend to break out the confinement, beginning to produce new contagious events.These results exhibit the need to employ correct times and duration of the confinement measures and that abruptly confinements without proper regulatory measures may trigger revivals.
In the context of a confinement based on regulatory measures such as lockdown, social distancing and restrictions on mobility; they could be implemented at any time of the epidemiological development and they will not follow a deterministic behavior (as shown previously).
In Fig 4, we explore the generation of the incidence when the probability of infection is manipulated by a piece-wise time-dependent C(t) function.In this figure, we show single realizations of the behavior of the incidence.In the figure one sees that if confinement is applied at relatively early stages, then a reduction of the C(t) function below the 25% of its initial value produces a deceleration of the incidence, at 25% the incidence is maintained approximately at a constant rate while anything above the 25% will correspond to increments in the incidence with stronger accelerations for larger values of the C(t) function.

Empirical estimation of W(t).
Having a self-consistent mechanism that could provide us with an estimation about the evolution of the weight function based on the real available data, would be desirable.We approximate to this problem by using the information accessible through empirical data, such as the empirical incidence i e (t) and its cumulative c e (t).In our stochastic approximation, the daily synthetic incidence is obtained from a set of random variables following a Poisson distribution; i.e. i s ðt j Þ ¼ P Iðt j Þ j¼1 w i ðt j Þ, hence the statistical mean of the cumulative of the daily incidence may be written as: where M represents the total number of trajectories to which the statistical mean is performed.
If the number of trajectories is large enough, then the statistical mean will approach the expectation value of the random variables, i.e., the Poisson parameters associated to the infectious individuals.Now, if the infection is sustained in the population (i.e. the probability of having none new infections is very low), then fluctuations around the mean number of infectious individuals at a given time t j does not significantly contribute to the averaged incidence over the ensemble; hence by considering the average of the infectious individuals at a given time t, i.e. � IðtÞ ¼ 1=M P M l I l ðtÞ, we can approximate the average of the cumulative to: In order to gain convergence, we fix as stated before, the Poisson parameters λ i to be obtained from a punctual distribution, i.e.P(% o ) ! {λ} i = % o δ ii = R o /t I δ ii , hence we write for the average of the cumulative: Our aim is to give an approximate description about the time dependent reproduction number through empirical quantities; in this regard, we do the replacement of the average synthetic cumulative and the averaged infectious population with their correspondent empirical descriptions; � c s ðt j Þ !c e ðt j Þ and � Iðt j Þ !I e ðt j Þ thus, by expanding the sum to the first steps of propagation, one can recurrently obtain the value of the weight function at the different times,  and c)), the confinement begins at the day 30 and the decreasing along 7 days the function C(t) from its initial value (one) to a 20% of its initial value (panel a)), to a 25% of its initial value (panel b)) and to a 30% from its initial value (panel c)).At the second row (panels d), e) and f)) the confinement begins at the day 60 where the function C(t) is decreased from its initial value to a 20% of its initial value (panel d)), to a 25% of its initial value (panel e)) and to a 30% from its initial value (panel f)).Finally at the third row (panels g), h) and i)), the function C(t) is initially decreased to a 20% of its initial value along a period of 7 days, increased back to a value of a 50% of its initial value along a period of 7 days and finally decreased back again at the day 120 along a period of 7 days to a 20% of its initial value (panel g)), to a 25% (panel h)) and to a 30% (panel i)).The y-axis labeled at the right of panels c), f) and i) corresponds to the values taken by W(t), H(t) and C(t).The figure was done using the COVID-19 parameters in a population of one million.https://doi.org/10.1371/journal.pone.0275216.g004 i.e.: therefore, at any given time, one can write for the weight function; while the number of infectious individuals at time t j can be approximated, according the definitions done earlier, as: I e ðt j Þ ¼ P t I k¼0 i e ðt j À t L À kÞ.

Development of the COVID-19 in some Mexican states
Along the development of the COVID-19 pandemic, we have used the stochastic model to follow the evolution of the spreading of the COVID-19 in certain Mexican states, some being the largest populated states, (Estado de Me ´xico �17 million, Ciudad de Me ´xico �9 million, Jalisco �8 million, Nuevo Leo ´n �6 million), housing the largest Mexican metropolitan areas (i.e., Me ´xico City, Guadalajara (Jalisco) and Monterrey (Nuevo Leo ´n) having populations around 5 to 8 million), and some middle-size populated states (i.e., Chiapas �5.5 million, Michoacan �5 million and Oaxaca �4 million) and the state of Nayarit which has a relative small population �1 million.
To model the spreading of the COVID-19 in these states, we collect the empirical data of the incidence through the reported cases by the scientific division of the Mexican federal government (CONACyT): https://datos.covid-19.conacyt.mx/and with that, we construct the empirical weight function as described in (Eq (16)), using the fixed epidemiological parameters of τ I = 14 days and R o = 4. Once the empirical weight function is obtained, we generate the synthetic incidence i s as described in our model (following a Poisson distribution modulated through the weight function).For doing this, we fix an initial time of propagation when the empirical data shows a sustained incidence plus the latency period, and by entering the initial latent individuals (the incidence from the past latency period), we initiate the generation of the incidence by fitting the initial number of infectious individuals such that the averaged trajectories of the synthetic incidence is in good agreement to the empirical data (see the supporting information for reference about the used initial conditions and estimation of the infectious of each case included the empirical incidence used in our simulations).
Our simulations are presented from Figs 5 to 7. They show the comparison between the incidence (first rows), the cumulative of the incidence (second rows) and the form of the weight functions for each of the studied cases (third rows), all of them when averaged over 1000 trajectories of the stochastic model.In the general context one can notice same tendencies in the dispersion of the disease and a good agreement to the real cases, with the possible exception of Campeche, (and less noticeable in Chiapas and Michoacan); all the discrepancies starting around the day 350 from the beginning of the dispersion (February 18th, 2020), in which the empirical data lies below our simulations.This discrepancy could be due to i) that the weight function should have taken smaller values than those obtained by the empirical data, which leads us to assume an under-counting on the infectious cases, i.e. there were many infectious imported cases (see Eq (16)), or ii) there was an under-counting in the incidence.The former case seems more likely since the discrepancy occurs at a time when the effect of the second wave at a national level has begun to subside hence possible relaxation of mobility measures may have happened.
At the last row of Figs 5-7 we present the form of the weight function for each of the studied cases.The large values of W(t) shown in the gray shaded areas represent an initial interval of the infection where the occurrence of the incidence, without reported infectious cases require that the majority of the real infectious cases are imported cases while after the shaded area one can expect that the dispersion of the disease relies mostly on the connection network of the local region, the confinement and social distancing regulations.
In this context, one could ask about encoded periodic patterns in the empirical weight function when the local dispersion regime dominates, (i.e., the regime at which confinement is well established and imported infectious cases do not contribute largely to the dispersion and estimated in clear areas of the weight function figures (last rows of Figs 5-7)), which could provide us with information about the development of the pandemic e.g.the local restrictions, the testing, mobility and the contact network, or even emergent periodic behavior encoded in the probability of infection, during the development of the COVID-19.These periodic patterns can be examined by performing the Fourier transform of the weight function for which, the peaks appearing in the spectrum correspond to frequencies associated to these periodic behaviors.
In Fig 8 we present the absolute value of the one-sided Fourier transform of the weight function of each of the studied cases.In the figure we identify different set of frequencies whose corresponding periods (given in days) lie on time scales associated to the weekly agenda, possibly confinements and de-confinements and also to larger patterns such as pandemic waves or even the emergence of seasonality.We identify a first set belonging to periods within a week for which in almost all the cases (with the exception of Nayarit) three peaks are present exactly at the same periods i.e. at 7 days, 3.5 days and 2.3 days which may confirm a global (not local) behavior suggesting these are related to the weekly agenda concerning the testing and readouts.A second time scale lies on the range from one two six months which we believe are connected to pandemic waves or periods of confinement, de-confinement and holiday seasons.At this time scale one observes periods that are shared by several states; the periods of roughly six months shared by Nuevo Leo ´n, Chiapas, Nayarit and Campeche; periods of 3 months shared by Ciudad de Me ´xico, Jalisco, Estado de Me ´xico, Michoacan, Chiapas, Nayarit and Oaxaca, and periods of roughly one and a half months which appearing in the states of Ciudad de Me ´xico, Jalisco, Estado de Me ´xico, Michoacan and Nayarit.The shared periods could be in some cases explained due to the closeness of certain states (such is the case of Ciudad de Mexico and Estado the Mexico or Jalisco and Michoacan) while other shared peaks between not neighboring states may be telling us something about the connection network of those states.The larger time scale corresponds to a period of roughly 9 months and is represented by the first peak appearing in Ciudad de Me ´xico and Estado de Me ´xico.These states share the largest metropolitan area of Mexico (� 21 million), hence this pattern suggest a relation to pandemic waves or even the emergence of a seasonal behavior of the COVID-19.This time scale is only present in these states and it may be a consequence of the amount of the population that has become infected, and the large degree of urbanization of this region.

Discussion
The good agreement between the stochastic model here proposed and the studied cases rely on considering the daily infectious events to follow Poisson distribution, together with the approximations when deriving the empirical weight function; in fact, since the Poisson distribution possess the property that if a set of random variables χ 1 , . .., χ n , each following independent Poisson distributions with parameters λ 1 , . .., λ n , then the sum of the random variables w ¼ P n i¼1 w i also follows a Poisson distribution with parameter, l ¼ P n i¼1 l i .In our context, the total number of the daily new infections is a random variable of a Poisson distribution, i.e.
which, under the assumption that all the individuals have the same probabilities of transmission, (this point refers to the employment of a punctual distribution about the λ i 's, i.e. λ i = % o δ ii ), then one can state that the probability of having n infections at the day t j should be given by p n (t j ) = (% o I(t j )W(t j )) n exp(−% o I(t j )W(t j ))/n!, or by considering the empirical estimation of the weight function (Eq (16)): i.e., the good agreement relies on the fact that we generate the daily incidence from a Poisson distribution whose parameter is the empirical daily incidence.However, in the most general case the assumption of considering same probabilities of transmission for any individual might be difficult to meet and instead one could attempt to connect the weight function to additional empirical quantities.In this regard, one could ask about the relation between W(t) and an effective reproduction number RðtÞ, the later representing the statistical mean of the infections caused by single individuals once the disease has begun to disperse.To answer this question, lets consider the statistical mean of the number of infected at the time t j due to the i- th infectious individual: where the [χ i ] k represents the possible outcomes of the random variable χ i of the i-th infectious at time t j , i.e., the possible number of infections that the i − th infectious could produce with probability p k [λ i W(t j )] of the k-th event and λ i W(t j ) is the statistical mean of all possible outcomes of the i-th contagious individual at time t j .In the case where the Poisson parameters λ i are also distributed according to a probability distribution P(% o ) (i.e., the rate at which each infectious individuals infect is also distributed around a mean % o ), then the average of the total number of infected individuals at a fixed time t j may therefore be given by: Let us consider now that at certain given time t j � t 0 , the number of infectious has become large and representative about the dispersion of the disease in the population, i.e., the number of infectious can be found homogeneously distributed in the population and the quantity P Iðt j Þ i¼1 l i =Iðt j Þ becomes representative about the mean % o .In other words, for times t < t 0 fluctuations are expected to dominate and as the number of infectious increases, the fluctuations reduce yielding a more localized value of the probability of infection.Therefore, at time t j > t 0 one could find a close relation between the weight function defined earlier and the timedependent effective reproduction number: Finally, the stochastic model we are presenting has the advantage that it does not require large computational resources to simulate the dispersion in high populated areas and provide us with a tool to simulate and study idealizations about the changes in the structure of the contact network among the populations.In the context of simulating real scenarios, the success of the model relies on certain information about the evolution incidence, something that cannot be known a priori, although the incorporation of agent-based or complex network models could bring insights about the tendency of the incidence and forecast the spreading of the disease in a given interval of time.Moreover, the shared frequencies shown in Fig 8 suggest a possible synchronization of infectious events between populated areas, revealing deeper complex connections among those regions which could be explored through the implementation of hybrid stochastic-complex network models.

Conclusion
In this paper we have derived an stochastic compartmental epidemiological model constructed from first principles consisting on a randomization about the number of the new infected population caused daily and following a Poisson process.We have shown that under this assumption, one can reconstruct and simulate the evolution of the development of the COVID-19 pandemic in Mexico, by introducing an additional time-dependent function (the weight function) which is in turn connected to a normalized effective reproduction number.Along this paper, we have focused on the epidemiological parameters corresponding to the COVID-19 disease and through the employment of the weight function, the model is capable of introducing and studying some conceptual behaviors such as herd immunity or certain idealized confinement scenarios.In the former we have employed an inverse-like logistic function of the fraction of the total infected population, which for the COVID-19 epidemiological parameters and without any confinement measures, we have found that the peak of the incidence scales by 20 to 30 days when the total population is increased by one order of magnitude while independently of the population sizes, the maximum incidence reaches from 2.5% to 2.7% of the total population; in the latter, we have explored the reaction of the dispersion of the disease when the population reacts to the infectious population (representing an intuitive reaction of the population under an epidemiological emergence), finding revivals in the incidence (infective waves) if confinement is abrupt and happens at earlier stages in the dispersion of the disease, and a flattening of the epidemic curve on the contrary situation.In this context, we have shown an acceleration of the generation of the incidence when the weight function takes values above 25% of its initial value, a steady behavior in the generation of the incidence for a 25% of its initial value and a deceleration in the incidence when the weight function take values below the 25% of its initial value.
In addition, we have employed our stochastic model together with the definition of the empirical weight function, to simulate the dispersion of the COVID-19 in some Mexican states, some of them housing the major metropolitan areas in Mexico and finding a very good agreement to the real scenarios implying that the infectious events in Mexico could be interpreted as homogeneously distributed events, providing us with an indirect mechanism to estimate dates of super-infectious or anomalous events.
Finally, we have applied the one-sided Fourier transform to the empirical description of the weight function with the intention to look at periodic patterns emerging in the mean of the daily infection which may give us insights about the evolution of the pandemic in Mexico.In this regard, we have found three different set of frequencies corresponding to different timescales which we identify to a weekly agenda about the capture of the readouts of the testings, confinement and also larger patterns which may be related to pandemic waves or even seasonality.

Fig 1 .
Fig 1.Schematic representation of the model.The figure shows an schematic representation of the stochastic model, where according to the number of infectious population a random number associated to each infectious individual provides the new number of infections.After a latency period, the amount of infected population at time t j passes into the infectious population generating a new set of random numbers associated to the new infections.Once the recovery period has passed since the beginning of the infections, the infectious population is removed from this condition and moved either to the recovered population with probability 1 − p d , or to the deceased population with probability p d , where p d represents a probability of decease.https://doi.org/10.1371/journal.pone.0275216.g001 Fig 2 are fixed to the estimated values of the COVID-19 disease described above, although this choice is done only for demonstrative purposes to exemplify the manipulation of the incidence through the weight function and the Fig 2 does not reflects a herd immunity effect emerging in the COVID-19 pandemic.

Fig 2 .
Fig 2. Effects of the weight function representing the herd immunity over the spreading of a disease based on the estimated COVID-19 parameters.At the first column (from left to right), the development of the latent, the infectious, the recovered, the deceased population and the cumulative of the incidence are plotted for a single trajectory.At the left bottom panel is plotted the incidence with the form of the weight function superimposed on the incidence, and representing the herd immunity effect.The y-axis labeled at the right of this panel corresponds to the values taken by H(t).In the panel at the right, the averages over 5000 trajectories of the normalized incidence for different sizes of the population are shown.The epidemiological parameters employed in the figure corresponds to the estimated values of the COVID-19.https://doi.org/10.1371/journal.pone.0275216.g002

Fig 3 .
Fig 3. Effect of a confinement based on a Gaussian decay as given in (10) for the COVID-19 parameters.At the four left panels, the incidence is shown with the weight function (containing the effects of herd immunity as presented before, and the effects of the confinement), over-imposed on the incidence.The rows from top to bottom show increasing values of I o : I o = 5% (top), I o = 10% (bottom), while the columns from left to right show an increasing decaying rate parameter: γ = 5 (left), γ = 10 (right).At the right, the figure shows the effect of a strong and early confinement, γ = 50 and I o = 1%, to the latent and the infectious population.In all the cases, the heard immunity is incorporated with the same parameters as used in Fig 2 while we use a total population of 10 million.https://doi.org/10.1371/journal.pone.0275216.g003

Fig 4 .
Fig 4. Effects of the piece-wise confinement in the incidence.At the first row (panels a), b)and c)), the confinement begins at the day 30 and the decreasing along 7 days the function C(t) from its initial value (one) to a 20% of its initial value (panel a)), to a 25% of its initial value (panel b)) and to a 30% from its initial value (panel c)).At the second row (panels d), e) and f)) the confinement begins at the day 60 where the function C(t) is decreased from its initial value to a 20% of its initial value (panel d)), to a 25% of its initial value (panel e)) and to a 30% from its initial value (panel f)).Finally at the third row (panels g), h) and i)), the function C(t) is initially decreased to a 20% of its initial value along a period of 7 days, increased back to a value of a 50% of its initial value along a period of 7 days and finally decreased back again at the day 120 along a period of 7 days to a 20% of its initial value (panel g)), to a 25% (panel h)) and to a 30% (panel i)).The y-axis labeled at the right of panels c), f) and i) corresponds to the values taken by W(t), H(t) and C(t).The figure was done using the COVID-19 parameters in a population of one million.

Fig 5 .
Fig 5. Comparison between the synthetic data generated by the stochastic model and real scenarios happening in some Mexican states and the form of their corresponding weight functions.The figure shows a comparison of the incidence (first row) and its cumulative (second row), between the synthetic data ( � i s and � c s ) generated by the stochastic model when averaged over 1000 trajectories to the real scenarios (i e and c e ), happening in Mexico City(CDMX), Jalisco(Jal) and Nuevo Leo ´n(NL).The shaded region for the cumulative of the incidence represents the 1st and the 3th quartil of the cumulative of the synthetic incidence generated randomly.The synthetic data was generated by employing the empirical estimation of the weight function (last row) obtained form Eq (16).The figure shows a period of roughly a year and a half of the spreading of the SARS-CoV-2 (from February 18th, 2020 to August 20th, 2021).https://doi.org/10.1371/journal.pone.0275216.g005

Fig 6 .
Fig 6.Comparison between the synthetic data generated by the stochastic model and real scenarios happening in some Mexican states and the form of their corresponding weight functions.The figure shows a comparison of the incidence (first row) and its cumulative (second row), between the synthetic data ( � i s and � c s ) generated by the stochastic model when averaged over 1000 trajectories to the real scenarios (i e and c e ), happening in Estado de Mexico(Edmx), Michoacan(Mich) and Chiapas(Chis).The shaded region for the cumulative of the incidence represents the 1st and the 3th quartil of the cumulative of the synthetic incidence generated randomly.The synthetic data was generated by employing the empirical estimation of the weight function (last row) obtained form Eq (16).The figure shows a period of roughly a year and a half of the spreading of the SARS-CoV-2 (from February 18th, 2020 to August 20th, 2021).https://doi.org/10.1371/journal.pone.0275216.g006

Fig 7 .
Fig 7. Comparison between the synthetic data generated by the stochastic model and real scenarios happening in some Mexican states and the form of their corresponding weight functions.The figure shows a comparison of the incidence (first row) and its cumulative (second row), between the synthetic data ( � i s and � c s ) generated by the stochastic model when averaged over 1000 trajectories to the real scenarios (i e and c e ), happening in Campeche(Camp), Nayarit(Nay) and Oaxaca(Oax).The shaded region for the cumulative of the incidence represents the 1st and the 3th quartil of the cumulative of the synthetic incidence generated randomly.The synthetic data was generated by employing the empirical estimation of the weight function (last row) obtained form Eq (16).The figure shows a period of roughly a year and a half of the spreading of the SARS-CoV-2 (from February 18th, 2020 to August 20th, 2021).https://doi.org/10.1371/journal.pone.0275216.g007

Fig 8 .
Fig 8. One-sided Fourier transform of the empirical weight function of some mexican federal entities.The figure shows the absolute value of the one-sided Fourier transform of the weight function derived from the empirical weight function of: a) Ciudad de Me ´xico, b) Jalisco, c) Nuevo Leo ´n, d) Estado de Me ´xico, e) Michoacan, f) Chiapas, h) Campeche, i) Nayarit and j) Oaxaca.The colored dots mark the most relevant frequencies associated to periodic patterns in the empirical weight function.The figure shows a period of roughly a year and a half of the spreading of the SARS-CoV-2 (from February 18th, 2020 to August 20th, 2021).https://doi.org/10.1371/journal.pone.0275216.g008