Stochastic SIR model predicts the evolution of COVID-19 epidemics from public health and wastewater data in small and medium-sized municipalities: A one year study.

The level of unpredictability of the COVID-19 pandemics poses a challenge to effectively model its dynamic evolution. In this study we incorporate the inherent stochasticity of the SARS-CoV-2 virus spread by reinterpreting the classical compartmental models of infectious diseases (SIR type) as chemical reaction systems modeled via the Chemical Master Equation and solved by Monte Carlo Methods. Our model predicts the evolution of the pandemics at the level of municipalities, incorporating for the first time (i) a variable infection rate to capture the effect of mitigation policies on the dynamic evolution of the pandemics (ii) SIR-with-jumps taking into account the possibility of multiple infections from a single infected person and (iii) data of viral load quantified by RT-qPCR from samples taken from Wastewater Treatment Plants. The model has been successfully employed for the prediction of the COVID-19 pandemics evolution in small and medium size municipalities of Galicia (Northwest of Spain).

Introduction

Since the emergence of the COVID-19 pandemics caused by the SARS-CoV-2 virus, great efforts have been made for the purposes of virus detection and epidemics forecasting. COVID-19 modeling approaches generally fall into one of the following categories, see [1]: statistic models for short-term forecasts; and mechanistic models, whether they are based on differential equations (like compartmental models) or agent-based [2], [3], for analyzing the spread dynamics of SARS-CoV-2, investigating future possible scenarios and/or simulate interventions to control virus spreading. In this work, we introduce a stochastic mechanistic model valid for both testing scenarios and short-term forecasting. The model has been developed to study and predict the evolution of the pandemics at the level of municipalities, and it has been calibrated and tested during a one year study in Galicia (northwest of Spain) using measurements from the health system and viral load from wastewater samples. The model developed for SARS-CoV-2 can be easily adapted to the surveillance of other pathogens and therefore, the methodology presented makes a significant contribution to wastewater based epidemiology [4].

Epidemiological compartmental deterministic models, like the Susceptible–Infected–Recovered (SIR) model firstly described by [5] (and extended versions of it) have been employed to predict COVID-19 spread [6], [7], [8], [9], [10]. However, predictability issues arise and models (whether they are phenomenological, mechanistic, or agent-based) are not efficient to predict the COVID-19 pandemics in the long term [11], [12]. Model predictive control approaches have been proposed to efficiently deal with the uncertainty and predict the effects of mitigation and suppression strategies [13].

The unpredictable nature of the pandemic spread has been tackled, on the one hand from the perspective of deterministic chaos [14], [15], [16] and, on the other hand, using stochastic models [17], [18]. Dynamic stochastic models for COVID-19 spread prediction can be broadly categorized into: (i) stochastic differential equations based in classical SIR models [8], [17], and (ii) compartmental models combined with Mote Carlo methods [6], [19], [20], [21].

Here we use a reinterpretation of the classical compartmental modeling of infectious diseases (SIR type models) as chemical processes, which are inherently stochastic and governed by the Chemical Master Equations (CME). The CME describes the evolution in time of probability distributions [22], [23], and the Stochastic Simulation Algorithm (SSA) by [24] can be used to compute exact realizations (time course trajectories) of the CME. Data series of new infected persons, as provided by the public health systems, are indeed realizations of a stochastic process or random walks on the positive integers. Therefore, approaches based on a CME equation solved by the SSA algorithm are particularly convenient to develop COVID-19 predictive models [25], [26], [27].

Data from the health system have an inherent delay (from the time of infection until the positive case is reported). Moreover, the significant percentage of asymptomatic cases characteristic of the COVID-19 pandemics, hampers the prompt detection by the health systems (efforts including screening or contact tracing have been implemented to overcome this difficulty). In this regard, wastewater have been proven to be a good complementary tool for COVID-19 surveillance. The analysis of SARS-CoV-2 viral load in sewage at wastewater treatment plants can be interpreted as one pooled test for the area where they are located. Pooled or group testing is an excellent tool for surveillance of diseases spread in animals and humans [28], indicated also to expand COVID-19 surveillance [29]. Different studies confirmed wastewater monitoring as a convenient complementary approach to COVID-19 surveillance and testing strategies [30], [31], [32] and, in fact, viral RNA detection in Wastewater Treatment Plants (WWTPs) has shown an anticipative capacity with respect to the cases reported by the health system in several studies [33], [34], [35], [36].

In this work, we developed a stochastic SIR model with good predictive capacity which incorporates the data from health system and wastewater analysis. The model has been calibrated and tested during a one year study with health system data and wastewater analysis data from five different small-medium size municipalities in Galicia (Northwest of the Iberian peninsula) [36]. We implemented the model in a software package that can be used for predictions in different scenarios and forecasting.

The model is a stochastic SIR: importantly, in stochastic mechanistic models uncertainty increases as we move into the future (like in the case of statistical models and unlike in the case of mechanistic deterministic models).

The model is simple in terms of number of states and parameters, and it shows a good predictive capacity with very few parameters. Simplest SIR models were previously reported to perform better than other models with greater complexities [9].

The model integrates, for the first time to the best of authors knowledge, SARS-CoV-2 viral load data from WWTPs within a mechanistic model with predictive capacity.

The model is very robust (we obtain similar values of the parameters for the WWTPs tested).

The model is capable to predict superspread events.

The model can be used to test different scenarios and for forecasting in small and medium size municipalities with WWTPs in time horizons of 7–10 days.

In the next section we introduce the model with its different extensions. First, we develop a simple stochastic SIR (which proves to correctly predict the evolution of infected individuals from health system in small and medium-size municipalities). Then, we extend the model with the capability of incorporating variable degradation rates (SIRv). In order to be able to predict superspread events, we also introduce a SIR model with jumps (SIRj) allowing for the infection of various persons at the same time from one infected individual. Finally, we incorporate the WWTP data together with the public health system data into an integrated model (SIRO).

The SIR model and its extensions (SIRv, SIRj and SIRO) have been developed and tested in the context of COVID-19 surveillance of small and medium size municipalities in the Northwest of the Iberian peninsula [36] from June 2020 to December 2021.

Methods

In this section we introduce the models developed to predict the evolution of the COVID-19 pandemics. Our aim was to obtain a good predictive capability with the simplest possible model. Importantly, our model is intended to predict the evolution at the local level (municipalities), and for short term time horizons (7–10 days maximum). First, we start with a stochastic SIR model with standard rate expressions.

SIR model

We consider a total population of individuals, distributed in the following types (or compartments):

Susceptible, , susceptible to acquire the disease.

Infected, , capable of spreading the disease to a susceptible person.

Recovered, , recuperated/dead from disease.

Considering the three categories , and , we can represent their interactions as follows: with and being parameters (infection and recovery rate constants, respectively). The equation governing the dynamics of this process taking into account its inherent stochasticity is the Chemical Master Equation (CME) [22]. Let be the states of the system and be a probability density function of the state at time . The CME governing the process (1) reads: We use the Stochastic Simulation Algorithm (SSA) [24] to simulate realizations of the SIR process (1). A trajectory of a single SSA simulation (or realization) is an exact sample from the probability function that is the solution of the CME (2), therefore, the solution of the CME can be approximated by a reasonable number of realizations.

In this way, the number of infected persons reported daily by the health system can be considered as single realization of the SSA for the SIR process (1). To calibrate with data, make predictions or forecasting for a given time horizon, we use the mean and standard deviation of SSA realizations.

SIR with variable infection rate

The infection rate for the SIR model is constant. The SIR model with constant parameters cannot capture a turning point in the evolution of the pandemics unless the number of susceptible remains a limiting factor (which was never the case at least in the first two years of the COVID-19 pandemics). However, in many territories (including the municipalities under study) a series of policies (including lockdowns, travel bans, capacity limitations for social gatherings and other restrictions) have been applied by the authorities, at different stages of the pandemics, to decrease the infection rate. Therefore, time course data of infected people used for calibration might show changes of tendencies (turning points) that cannot be captured correctly by the standard SIR. Here we define a turning point () as a point in time in which a (sustained) change in the sign of the slope of the total number of infected individuals is detected.

For calibration purposes, when the data show a turning point, we propose an extended SIR model starting from (1) with a variable infection rate defined as: where the is the simulation time, and is a new parameter of the SIR model representing the turning point date. Note that, if , SIRv is equivalent to the SIR model.

SIR with jumps

Standard SIR models consider the infection of one person at a time (the stoichiometric S/I ratio is ). This is the reason why SIR models cannot capture the jumps observed in time course data from SARS-CoV-2 infected individuals. One infected individual can infect more than one susceptible person as a result of the same interaction event (in this case the stoichiometric S/I ratio is with , being a positive integer). In order to consider these bursting (superspreading) events, we generalize the reactions associated with the SIR model (1) as follows: In this work we consider which is compatible with the mass gathering restrictions established by the health authorities in the period and geographic area under study. Note that, for , SIRj is equivalent to the simplest SIR model proposed.

Let be the system states and be a probability density function of being in state at time . The associated CME to the reaction set (4) reads: Note that this CME is an extension of the one obtained for the simplest SIR model, Eq. (2), by adding the new terms related to the new reactions (more than one infected at the same time) in the second line of expression (5).

SIRO model

The SIRO model is an adaptation of the SIR formulation described in Section 2.1 to incorporate the number of infected people (I) observed through the viral load monitoring in WWTPs, together with the subset of the total infected people (O) observed by the public health system, through the following set of reactions: where , and are the susceptible, total infected and total recovered individuals in a given municipality. The infected individuals detected and reported by the health system are denoted by , whereas are recovered individuals that had been previously reported as infected. Let be the states of the system and be a probability density function of the state at time . The CME governing the dynamics of the process described by the set of reactions in (6) reads: We consider two observables of the model, denoted by ; on the one hand, the subset of the total infected people being reported by the health system is , and on the other hand the viral number of gene copies detected in wastewater is: where , , and are defined in Table 1. This formula is considered to properly normalize the wastewater data. In weeks where the number of total infected detected by viral load in WWTP at first day is lower than the infected number reported by the health system at first day , we re-normalize as follows:

Table 1

Notation and units.

Notation	Units	Description
C_w	copies/L	Copies per Liter of wastewater [37], [38]
Q	m^3/day	Inlet flow by GESECO Aguas https://www.gesecoaguas.es/
ρ	g faeces/mL faeces	Faeces density [39]
C_f	copies/mL faeces	Average number of copies per mL of faeces in infected person [37], [38]
ξ	g faeces/person/day	Average weight of faeces per person per day [39]

Results and discussion

In this section, we first describe the scripts implemented for the proposed SIR and SIRO models. Then, we illustrate how the simplest SIR model has the capacity to represent the data provided by the public health system. After that, we show the capacity of the SIRO model to predict in a week time horizon the evolution of the number of infected people, starting from the WWTP viral load levels detected during the previous week.

Data and scripts

The models proposed in the previous section have been implemented in MATLAB and are available at https://github.com/manuelpajaro/stochasticSIR_O inside the folder stochasticSIR_O. The main programs are SIR_ppal.m for the simplest SIR model, SIRv_ppal.m and SIRj_ppal.m for the SIRv and SIRj model extensions respectively, and SIRO_ppal.m for the SIRO model. The data used for this study are saved in the DATA.mat archive which is available within the folder sotchasticSIR_O. The data of new infected persons per day and municipality were provided by the Galician Health System (Servizo Galego de Saúde SERGAS). The data from SARS-CoV-2 viral load in the WWTPs under study were obtained within the DIMCoVAR project consortium [36]. Specifically, we provide the number of infected persons detected by the health system per day (I_Locality variables) and the corresponding cumulative infected cases for fourteen days (Icum14_Locality variables). The variable in (8) computed from the measurements of viral load in sewage is stored in WTP_Locality.

The main programs (SIR_ppal.m, SIRv_ppal.m, SIRj_ppal.m, and SIRO_ppal.m) share the same structure. The user can choose the municipality by modifying the value of the variable in the following code: where currently to select Ares, Baiona, Gondomar, Melide or Nigrán, respectively (of course the list can be extended to other municipalities of interest if access to data is provided). The starting date of the simulation can be chosen by assigning to variable the date in the format year month day, , as follows: Finally, the free parameters for each model are indicated next. In the SIR model (SIR_ppal.m) the parameters (recovery rate constant) and (infection rate constant) are free a priori. To avoid autocorrelation problems, we fix the parameter (cumulated incidence is calculated for a time interval of 14 days). Importantly, the mitigation policies do not affect the value of the recovery rate constants of individuals (which is coherent with fixing the parameter ), but the infection rate constant (which should be therefore calibrated from data). For the SIRv model (SIRv_ppal.m), the parameters , and can be calibrated. is the point in time where a change in the sign of the slope of the dynamics occurs (i.e. a turning point as defined in the previous section). For the SIRj model (SIRj_ppal.m) the parameter can be used for calibration purposes. In the SIRO model (SIRO_ppal.m), the parameters and can be used for calibration. Generally, the recovery rate constant should be equivalent, i.e. .

SIR predictions

The SIR model is calibrated from the data of infected persons provided by the public health system. The parameter for the SIR model described in Section 2.1 is estimated from the data of infected persons provided by the public health system in order to fit the trajectory of fourteen days cumulative infected cases. Calibrations are done per week. As justified in the previous section, the recovering rate constant is fixed at .

Remarkably, the evolution of the number of infected persons detected by the health system could be accurately predicted for most weeks (more than 70% of the total of 265) during the frame of the study with the simplest stochastic SIR model with a fixed and only four different values of the infection rate constant . Specifically, , see Fig. 1 and Table 2. This has important implications for assessment and forecasting. First, the model can be used for quantitative assessment of COVID-19 mitigation policies, i.e. to quantify the impact of restrictions on the infection rate constant. Second, this facilitates the use of the model for forecasting purposes, to the point of needing only the one parameter (which can be fixed to the value obtained for the previous week) and the initial condition (number of infected reported by the health system) at the starting date, in order to obtain predictions at one week time horizon.

Fig. 1

From SSA realizations of the SIR model we compute the infected mean (blue line) and standard deviation (dotted blue lines). Black squares represent the accumulated infected cases reported by the health system. for all cases. The infection rate constants are: for Gondomar, for Baiona, for Ares and for Nigrán.

Table 2

Total number and percentage of weeks for which the specified model (left column) provided an accurate prediction.

Model	Ares	Baiona	Gondomar	Melide	Nigrán	Total	%
SIR (α=114)	40	39	36	33	39	187	70.6
SIR (α<114)	11	11	9	13	10	54	20.4
SIRv	1	1	2	2	5	11	4.1
SIRj	1	1	4	4	3	13	4.9
Total	53	52	51	52	57	265	100

As it is indicated Table 2, the SIR model accurately predicts the evolution in time of the infected cases detected by the health system (90% of the weeks by the simplest stochastic SIR with fixed and fixed , 4.1% of the weeks showed a variable , and 4.9% of the weeks showed jumps).

In Fig. 2 we show two exceptional cases in which the best fits were obtained with (this might happen for example for those weeks in which many infected persons are recovered in the first days) as it happened in Melide for the week of 18 April 2021.

Fig. 2

From SSA realizations of the SIR model we compute the infected mean (blue line) and standard deviation (dotted blue lines). Black squares represent the accumulated infected cases reported by health system. The infection rate is for both cases. The recovered parameters are and for the 40th and 60th weeks of Melide, respectively.

SIRv predictions

Those exceptional weeks for the overall time of the study showing turning points are specified in Table 3. Two selected examples are depicted in Fig. 3 where the predictions obtained by SIR and SIRv are compared. We show the 36th week of Ares and the 43rd week of Nigrán. As it can be seen in Fig. 3, the SIR model with two infection rates, and , generates much more precise predictions of the cumulative level of infection. Note that SIRv can be used for calibration, and to analyze the conditions under which the turning point is produced, but, without additional information on the time of the turning, it cannot be used for forecasting (in the conclusion we propose as future work the implementation of real time calibration and machine learning methods to overcome this limitation).

Table 3

Parameters for the SIR model with 2 different beta, , .

Week	Locality	Date	β_0	β_1	α	T_change
36	Ares	01/11/20	0.14	0.03	1/14	3
36	Baiona	01/11/20	0.14	0.03	1/14	4
35	Gondomar	25/10/20	0.14	0.03	1/14	3
46	Gondomar	10/01/21	0.22	0.03	1/14	5
46	Melide	10/01/21	0.22	0.14	1/14	3
48	Melide	24/01/21	0.03	0.07	1/14	3
31	Nigrán	27/09/20	0.14	0.03	1/14	2
35	Nigrán	25/10/20	0.14	0.03	1/14	3
36	Nigrán	01/11/20	0.14	0.03	1/14	3
43	Nigrán	20/12/20	0.14	0.03	1/14	4
60	Nigrán	18/04/21	0.22	0.03	1/14	5

Fig. 3

From SSA realizations of the SIR model we compute the infected mean (blue line) and standard deviation (dotted blue lines). Black squares represent the accumulated infected people reported by the health system. The results of SIR and SIRv models are shown in the first and second columns, respectively. for all cases. The infection parameters for SIR is , for both localities and for the SIRv model. is respectively 3 and 4 for the weeks under study.

SIRj predictions

The SIR model fails to predict superspreading events. Using the SIRj model developed in the Methods section, we not only reproduce the dynamics of cumulative cases, but also obtain better realizations than those obtained using the simplest SIR model. In order to measure the accuracy of the realizations for new infected cases we define the following metrics: where is a vector of real data provided by the health system whereas is one realization obtained from SSA simulation. We define as the number of realizations with , where is the number of exact realizations (see for example the last plot in Fig. 4) and is the number of realizations which have exact number of infected for all days except one for which there is only a difference of one (see for example the first plot of the second row in Fig. 4).

Fig. 4

From SSA realizations we compute the infected mean (blue line) and standard deviation (dotted blue lines). Black squares represent the accumulated infected (first row) or the new infected cases (second row) reported by the health system. Red squares are the new infected cases obtained from the best SSA realization. The results obtained with SIR and SIRj models are shown in the first and second columns, respectively. The parameters are, and for SIR and and for SIRj.

In Fig. 4 we compare the results obtained after SSA realizations for the SIR and the SIRj models. We selected the 35th week of Melide which stars at 25 November 2020 for which the SIR model does not capture the infection spread using the four values of proposed as reference. So, we use and for the simplest SIR, and for the model with jumps the parameters are given in Table 4, and . For the parameters chosen, both models reproduce the real cumulative cases accurately, as it can be seen in the first row of Fig. 4. However, when we observe the number of new infected persons per day, the SIRj model generates the best realizations with three exact realizations, (last plot in Fig. 4) and several realizations with from the computed, . For the SIR model there are not exact realizations () and only 7 of the realizations are obtained with (), one of them is shown in the first plot in the second row in Fig. 4. Again, SIR models cannot predict superspreader events in absence of a priori information, but we can use SIRj model to detect and quantify a posteriori the occurrence of a superspreader event.

Table 4

Parameters for the SIR with jumps model.

Week	Locality	Date	β	α
34	Ares	18/10/20	0.03	1/14
55	Baiona	14/03/21	0.02	1/14
23	Gondomar	07/08/20	0.02	1/14
33	Gondomar	11/10/20	0.02	1/14
39	Gondomar	22/11/20	0.02	1/14
55	Gondomar	14/03/21	0.02	1/14
3	Melide	20/03/20	0.03	1/14
25	Melide	16/08/20	0.02	1/14
35	Melide	25/10/20	0.02	1/14
45	Melide	03/01/21	0.02	1/14
2	Nigrán	12/03/20	0.02	1/14
23	Nigrán	02/08/20	0.03	1/14
55	Nigrán	14/03/21	0.02	1/14

Probability density distributions

The realizations of the SSA algorithm are used to compute the mean number of infected cases and the standard deviation to assess the accuracy of the predictions of the stochastic model. Besides this, for all SIR models we also obtain the probability density distribution of (new and accumulated) infected cases from the SSA realizations which is an approximation of the solution of the associated CME (2). In Fig. 5 we depict the probabilities of a given number of cumulative cases for a specific week (we chose the 36th week of Ares, already discussed, see first row of Eq. (3)). As it can be observed in Fig. 5 the cases reported by the health system are close to the mode of the obtained probability distributions, which means that the proposed model captures well the evolution of the epidemics. Moreover, the SSA algorithm can be used to estimate the new infected cases per day. The probability distribution of these new cases is shown in Fig. 6 where the data provided by the health system (vertical black dotted lines) fall within the most probable cases predicted.

Fig. 5

Probability density distribution for accumulated cases obtained from SSA realizations of the SIR with variable infection rates for Ares starting the 1 November 2020. The vertical black dotted line represent the real cases reported by the health system. The parameters considered are, , , and .

Fig. 6

Probability density distribution for new cases obtained from SSA realizations of the SIR with variable infection rates for Ares starting the 1 November 2020. The vertical black dotted line represent the real cases reported by the health system. The parameters considered are, , , and .

SIRO predictions

The SIRO model can be used to predict, via SSA simulations, the cumulative cases reported by the public health system (observed cases ) and those reported by wastewater treatment plant samples (we assume they are a proxy of the total infected cases ). The viral load in the WWTP samples was quantified using by RT-qPCR [36]. We consider periods of ten days as depicted in Table 5 together with the parameters obtained for the SIRO model. The mathematical model shows a good predictive capacity allowing us to forecast the evolution of infected persons, both total () and observed () by the health system in the municipalities within a horizon of 10 days. In Fig. 7 we show, as a representative illustration of the model outcome, the predictions of the SIRO model using SSA realizations for the five municipalities: Melide, Nigrán, Baiona, Gondomar and Ares. The model predictions, mean and standard deviations (continuous and dashed lines, respectively) for the total number of infected (blue lines) and the observed number of infected (black lines) are depicted in Fig. 7 together with the real data obtained from viral load in sewage (blue circles) and health system (black squares).

Table 5

Parameters obtained for the SIRO model.

	Locality	Date	β	γ	α_I	α_O
1	Ares	29/08/20	0.005	0.5	1/5	1/14
2	Ares	10/09/20	0.005	0.05	1/7	1/14
3	Ares	17/10/20	0.35	0.035	1/14	1/14
4	Ares	24/10/20	0.005	0.08	1/5	1/5
5	Ares	31/10/20	0.03	0.005	1/14	1/14
6	Ares	10/01/21	0.005	0.5	1/5	1/14
7	Ares	16/01/21	0.2	0.08	1/14	1/14
8	Ares	24/01/21	0.02	0.12	1/10	1/10
9	Ares	31/01/21	0.005	0.17	1/5	1/5
10	Ares	07/02/21	0.02	0.005	1/7	1/7
11	Ares	28/02/21	0.005	0.0005	1/5	1/5
12	Ares	07/03/21	0.005	0.025	1/7	1/7
13	Ares	14/03/21	0.25	0.015	1/14	1/14
14	Ares	21/03/21	0.02	0.0005	1/7	1/7
1	Baiona	16/08/20	0.3	0.02	1/14	1/14
2	Baiona	11/10/20	0.3	0.1	1/14	1/14
3	Baiona	18/10/20	0.12	0.1	1/14	1/14
4	Baiona	25/10/20	0.12	0.01	1/14	1/14
5	Baiona	02/11/20	0.005	0.025	1/7	1/7
6	Baiona	08/11/20	0.08	0.1	1/14	1/14
7	Baiona	17/11/20	0.05	0.17	1/14	1/14
8	Baiona	25/11/20	0.07	0.12	1/14	1/14
9	Baiona	01/12/20	0.03	0.17	1/14	1/14
10	Baiona	05/12/20	0.1	0.07	1/14	1/14
11	Baiona	14/12/20	0.05	0.05	1/14	1/14
12	Baiona	21/12/20	0.22	0.025	1/14	1/14
13	Baiona	29/12/20	0.12	0.05	1/14	1/14
14	Baiona	04/01/21	0.05	0.18	1/5	1/14
15	Baiona	16/01/21	0.005	0.02	1/14	1/14
16	Baiona	25/01/21	0.005	0.02	1/14	1/14
17	Baiona	30/01/21	0.005	0.25	1/5	1/14
18	Baiona	06/02/21	0.1	0.05	1/14	1/14
19	Baiona	20/02/21	0.1	0.005	1/14	1/14
20	Baiona	01/03/21	0.005	0.12	1/10	1/10
21	Baiona	07/03/21	0.005	3.2	1/7	1/14
22	Baiona	13/03/21	0.25	0.07	1/14	1/14
23	Baiona	20/03/21	0.08	0.005	1/10	1/10
1	Gondomar	11/01/21	0.3	0.02	1/14	1/14
2	Gondomar	18/01/21	0.005	0.03	1/7	1/7
3	Gondomar	25/01/21	0.3	0.0005	1/14	1/14
4	Gondomar	01/02/21	0.005	0.12	1/7	1/7
1	Melide	23/08/20	0.005	0.07	1/5	1/5
2	Melide	30/08/20	0.02	0.01	1/14	1/14
3	Melide	07/09/20	0.05	0.01	1/14	1/14
4	Melide	12/09/20	0.003	0.001	1/14	1/14
5	Melide	17/10/20	0.4	0.05	1/14	1/14
6	Melide	24/10/20	0.02	0.03	1/10	1/10
7	Melide	31/10/20	0.02	0.03	1/14	1/14
8	Melide	07/11/20	0.4	0.01	1/14	1/14
9	Melide	21/12/20	0.1	0.025	1/14	1/14
10	Melide	27/12/20	0.05	0.05	1/10	1/14
11	Melide	03/01/21	0.05	0.1	1/7	1/14
12	Melide	10/01/21	0.25	0.025	1/14	1/14
13	Melide	17/01/21	0.02	0.06	1/7	1/7
14	Melide	25/01/21	0.005	0.07	1/5	1/7
15	Melide	31/01/21	0.005	0.07	1/5	1/5
1	Nigrán	15/08/20	0.06	0.25	1/14	1/14
2	Nigrán	19/09/20	0.16	0.005	1/7	1/7
3	Nigrán	12/10/20	0.17	0.1	1/14	1/14
4	Nigrán	18/10/20	0.4	0.02	1/14	1/14
5	Nigrán	25/10/20	0.0005	0.25	1/5	1/14
6	Nigrán	01/11/20	0.2	0.03	1/14	1/14
7	Nigrán	21/11/20	0.1	0.03	1/14	1/14
8	Nigrán	28/11/20	0.01	0.22	1/14	1/14
9	Nigrán	05/12/20	0.05	0.13	1/14	1/14
10	Nigrán	14/12/20	0.35	0.005	1/14	1/14
11	Nigrán	21/12/20	0.13	0.04	1/14	1/14
12	Nigrán	09/01/21	0.12	0.07	1/14	1/14
13	Nigrán	18/01/21	0.005	0.15	1/7	1/7
14	Nigrán	25/01/21	0.17	0.005	1/14	1/14
15	Nigrán	01/02/21	0.02	0.02	1/14	1/14
16	Nigrán	08/02/21	0.13	0.005	1/14	1/14
17	Nigrán	15/02/21	0.06	0.005	1/7	1/7
18	Nigrán	22/02/21	0.1	0.01	1/14	1/14
19	Nigrán	01/03/21	0.1	0.06	1/14	1/14
20	Nigrán	22/03/21	0.17	0.06	1/14	1/14

Fig. 7

Predictions of the SIRO model using SSA realizations for each locality, Melide (first row), Nigrán (second row), Baiona (third row), Gondomar (fourth row) and Ares (fifth row). The estimated infected from viral load in sewage (blue points) and mean number of infected individuals, , obtained using the SIRO model (blue lines) are shown at the left column. The number of infected persons reported by the health system (black points) and the SIRO model prediction, , (black lines) are shown at the right column. The standard deviations obtained for and from the SIRO model simulations are represented with blue and black dashed lines, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Conclusions

In this work we present a stochastic model based in the classical compartmental models (SIR type) for which we have incorporated the stochastic character of viral spread. We consider the transitions between each group of persons in which the total population is subdivided (for example, Susceptible, Infected and Recovered) as reactions of chemical species. The Chemical Master Equation (CME) is the model that incorporates the inherent stochasticity of chemical reaction systems and therefore, by using this new formulation and proposing the corresponding CME, we are able to incorporate the noise of the viral infection propagation in a natural form. Moreover, in this article we solve the different models proposed using the Stochastic Simulation Algorithm (SSA) of Gillespie, which allow us to obtain the solution of the CME (the time evolution of the probability density function of infected persons) together with realizations (possible trajectories of the time evolution of the number of infected persons). Whereas deterministic SIR models only capture a unique trajectory for a set of parameters, their stochastic versions produce a high number of realizations providing us automatically with a measure of the noise in the epidemics spread. We can observe how the uncertainty grows with time as the standard deviations reported or the tails of the distributions obtained are higher as we move forward in time.

The stochastic version of the classical SIR model presented in this work has been developed with the aim to analyze and predict the evolution of the COVID-19 epidemics at the level of municipalities in Galicia (Northwest of Spain) integrating data from wastewater and the public health system in Galicia. The model has been calibrated and successfully applied for SARS-CoV-2 surveillance in different municipalities [36], using data from the galician public health system (SERGAS), and RT-qPCR measurements from wastewater samples. The SIRO model presented is fed by the total number of infected persons estimated using the viral load in sewage to predict the number of positives observed by the health system during the following week.

The software with the model implementation is provided, such that it can be calibrated and used to make analysis and predictions in other municipalities (with WWTP data available), and for different purposes. The model is easy to calibrate (only 2 parameters for calibration from time course data) and, as it has been shown in Section 3.2, it is very robust, as it allows estimating the number of infected persons in all the municipalities studied during more than one year by using only four different infection rate parameters. Remarkably, we have obtained a set of 4 parameters for the stochastic SIR model which are enough to describe the viral infection in the hole period of time for the five municipalities considered. These parameters can be explained by the restrictions imposed (lockdowns) during the pandemic. Therefore, we provide evidences of the usefulness of this kind of stochastic models, based in CME, to predict the viral spread of infectious diseases, in particular for the case of COVID-19 spread.

Due to the robustness of the model, its easy calibration, its flexibility to be adapted to the surveillance of other pathogens, we believe it makes a significant contribution to wastewater based epidemiology. As a future work, we propose to incorporate artificial intelligence techniques for automated real time calibration of the parameters, which can significantly facilitate the forecasting of turning points.

Software availability

The scripts for the models used are available under GPLv3 license at https://github.com/manuelpajaro/stochasticSIR_O.

CRediT authorship contribution statement

Manuel Pájaro: Conceptualization, Methodology, Software, Data curation, Writing – original draft, Writing – review & editing, Supervision. Noelia M. Fajar: Data curation, Validation, Writing – original draft. Antonio A. Alonso: Conceptualization, Methodology, Funding acquisition. Irene Otero-Muras: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.