Table of contents:
2.1. The stochastic SIR epidemic model
3. The SHAR modelling framework
4. The SHARUCD model framework
4.1. Modeling the effects of the control measures (first wave)
4.2. Model simulations with control and data
4.3. Model Validation and short-term prediction exercise with control measure
4.4. Growth rates and reproduction ratio
5. The refined SHARUCD MODEL with import and seasonality (second phase)
5.1. Model simulations and short term predictions
5.2. Short term predictions of COVID-19 in the Basque Country
6. Acknowledgments and contributions
In December 2019, a new respiratory syndrome (COVID-19) caused by a new coronavirus (SARS-CoV-2), was identified in China  and spread rapidly around the globe. COVID-19 was declared a pandemic by the World Health Organization (WHO) in March, 2020 . Up to date, more than 230 million cases were confirmed with more than 4.7 million deaths, with a global case fatality ratio (CFR) of approximate 2% . One year after the pandemic started, four vaccines are now licensed for emergence use in Europe:
- Pfizer-BioNTech COVID-19 Vaccine and Moderna COVID-19 Vaccine (mRNA ~ 95% VE)
- Oxford-AstraZeneca and Johnson & Johnson’s Janssen COVID-19 Vaccine (Viral vector ~ 70% VE) .
Figure 1: COVID-19 pandemic figures. In a) bar chart race for the confirmed cases around the globe, in b) bar chart race for the deceased cases around the globe and in c) bar chart race for the vaccine doses administrated per 1 million people around the globe. As of 7 April 2021, a total of 650.382.819 vaccine doses have been administered globally .
In March 2020, a Multidisciplinary Task Force (so-called Basque Modelling Task Force, BMTF) was created to assist the Basque Health managers and the Basque Government during the COVID-19 responses. BMTF is a modeling team, working on different approaches, including stochastic processes, statistical methods and artificial intelligence. Members were collaborating taking into consideration all information provided by the public health frontline and using different available datasets in respect to the COVID-19 outbreak in the Basque Country. The objectives were, besides projections on the national health system’s necessities during the increased population demand on hospital admissions, the description of the epidemic in terms of disease spreading and control, as well as monitoring the disease transmission when the country lockdown was gradually lifted. All modeling approaches are complementary and are able to provide coherent results, assuring that the decisions made using the modeling results were sound and, in fact, adjusted to the current epidemiological data.
In this dashboard we describe and present the results obtained by one of the modeling approaches developed within the BMTF, specifically using extended versions of the basic epidemiological SIR-type models, able to describe the dynamics observed for tested positive cases, hospitalizations, intensive care units (ICUs) admissions, deceased and finally the hospital discharges used as a proxy for a proportion of recovered individuals. We start presenting the properties of the basic SIR epidemic model for infectious diseases , with the goal to introduce notation, terminology, and results that will be generalized in later sections on more advanced models describing the COVID-19 epidemiology. Graphics will be updated every month.
2. The SIR epidemic model:
The SIR epidemic model divides the population into three classes: susceptible (S), Infected (I) and Recovered (R). It can be applied to infectious diseases where waning immunity can happen, and assuming that the transmission of the disease is contagious from person to person, the susceptibles become infected and infectious, are cured and become recovered. After a waning immunity period, the recovered individual can become susceptible again .
In the simple SIR epidemics without strain structure of the pathogens we have the following reaction scheme, for the possible transitions from one to another disease related state, susceptibles S, infected I and recovered R is shown in Fig. 2 a). For a host population of N individuals, with contact and infection rate β, recovery rate γ and waning immunity rate α, the dynamic model in terms of ordinary differential equations are shown in 2b), and the dynamical behavior for each variable is shown in Figure 2c).
Figure 2: For the basic SIR type model, in a) reaction scheme, in b) ODE system and in c) time dependent solution simulation, with a population N=100 (for the interpretation in percentages), and starting values I=40, S=60 and R=0 we fixed β=2.5, α=0.1, and γ=1. In green the dynamics for the susceptible S(t), in pink the dynamics for the infected I(t) and in blue the dynamics of the recovered R(t).
2.1. The stochastic SIR epidemic model:
The stochastic SIR epidemic is modeled as a time-continuous Markov process to capture population noise. The dynamics of the probability of integer infected and integer susceptibles, while the recovered follow from this due to constant population size, can be give as a master equation  in the following form:
This process can be simulated by the Gillespie algorithm giving stochastic realizations of infected and susceptibles in time [11,12]. The deterministic approach is obtained via the mean field approximation. For more details on the calculations see e.g. [7,10].
How does the program work: a short guide
Choose the parameters values for your simulation. Each susceptible individual is surrounded by 8 neighbors. “P” defines the probability of infection for each iteration. “Population Immunity” defines how many susceptible people exist in the population, and “Mortality” refers to the disease induced death probability.
- Color code: green cells represent susceptible individuals, red cells represent infected individuals, blue cells represent recovered and immune individuals, and black cells represent deceased individuals.
When the simulation starts, susceptible individuals (green cells) become infected randomly and cells turn red. For each iteration, new infection might occur within the 8 neighbors, with probability “P” of infecting someone. For example, if P = 50%, on average each cell will infect 4 neighbors per iteration.
Individuals are infected during the “Infected Period” and become recovered and immune (blue cells). For disease induced death (black cells), infected cells have a probability of “dying” at each iteration given by 1 – (1 – μ)^ (1 / p), where μ is the mortality rate and p is the infected period.
- Output: On the right hand side, model output in terms of time series, for each variable, susceptible, infected, recovered/immune and deceased.
3. The SHAR Modeling framework:
We develop a basic SHAR-model including now a severity ratio η for susceptible individuals (S) developing severe disease and possibly being hospitalized H or (1- η ) for milder disease, including sub-clinical and eventually asymptomatic (A) infections, where mild infected A have different infectivity from severe hospitalized disease H, parametrized by a ratio Φ to be smaller or larger than 1 comparing to baseline infectivity rate β of the “hospitalized” H class and altered infectivity rate Φβ for mild/“symptomatic” A class. Hence we obtain a SHAR-type model
that needs to be further refined to describe COVID-19 dynamics, adjusting the modelling framework to the available empirical data.
4. The SHARUCD model framework:
To model COVID19 dynamics, we use SHARUCD-type models, an extension of the well known simple SIR (susceptible-infected-recovered) model, with infected class (I) partitioned into severe infections prone to hospitalization (H) and mild, sub-clinical or asymptomatic infections (A). For severe infections prone to hospitalization, we assume the following dynamics: severe hospitalized individuals H could either recover, with a recovery rate γ , be admitted to the ICU facilities U, with a rate ν, or eventually decease into class D before being admitted to the ICU facilities, with a disease induced death rate μ. The ICU admitted patients could recover or die. For completeness of the system and to be able to describe the initial introductory phase of the epidemic, an import term ρ should be also included into the force of infection.
As we investigate cumulative data on the infection classes and not prevalence, we also include classes C to count cumulatively the new cases for “hospitalized” CH, “asymptomatic” CA, recovered CR and ICU patients CU. In this way we can easily include a ratio ξ of under-notification of mild/asymptomatic cases. The deceased cases are automatically collecting cumulative cases, since there is no exit transition form the death class D. We consider primarily SHARUCD model versions as stochastic processes in order to compare with the available data which are often noisy and to include population fluctuations, since at times we have relatively low numbers of infected in the various classes. The stochastic version can be formulated through the master equation in application to epidemiology in a generic form using densities of all variables x1:=S/N , x2:=H/N, x3:=A/N , x4:=R/N , x5:=U/N , x6:=CH/N, x7:=CA/N , x8:=CU/N and x9:=D/N and x10:=CR/N, hence state vector x := (x1,…,x10)tr, giving the dynamics for the probabilities p(x,t) as
For a more detailed description of the SHARUCD model framework, see [8,9]. The SHARUCD model is able to describe the dynamics observed for PCR tested positive cases, hospitalizations, intensive care units (ICUs) admissions, deceased and finally the recovered. Keeping the biological parameters for COVID-19 in the range of the recent research findings, but adjusting to the phenomenological data description, the models were able to explain well the exponential phase of the epidemic and fixed to evaluate the effect of the imposed control measures. With good predictability so far, consistent with the updated data, we continue this work while the imposed restrictions are relaxed and closely monitored. The deterministic version of the refined model is given by a differential equation system for all classes, including the recording classes of cumulative cases CH, CA, CR and CU by
Note, that in original versions of the SHARUCD model we had included ICU admissions similar to other transitions like progression to deceased, but then due to the observed synchronization of transitions to hospitalization and asymptomatic cases (see section 4.4), we refined the model, changing the transition into ICU admissions from the previous assumption of hospitalized patients recovering with recovery rate γ, being admitted to ICU facilities with rate ν or dying with disease induced death rate μ by the assumption that ICU admissions are consequence of disease severity prone to hospitalization analogously to rate η.
4.1. Modeling the effects of the control measures (First phase):
With the initial parameters estimated and fixed on the exponential phase of the epidemic, we model the effect of the disease control measures introduced using a standard sigmoid function σ(x)=1/(1+ex) . This function was able to describe well the gradual slowing down of the epidemics, as it turned out later in the response of the disease curves to the control measures. The current SHARUCD framework evaluates the seasonal effect on COVID-19 dynamics in the Basque Country, after social distancing measures start to be lifted (from May 4 – phase 0 and from May 11 -phase 1 towards the “new normality”). Model simulations with a single parameter set include now 2 control functions, shown in Fig. 3: one for the lockdown implementation and a second one for gradual lockdown lifting. Low seasonality is assumed to play a role, helping to keep transmission at low levels (see Fig. 5).
The seasonal forcing is implemented as a simple cosine function
β(t)=β0 (1+ θ cos(ω(t+φ))) on top of the already implemented sigmoidal control function σ(x)=1/(1+ex) from what we obtain the final format β(t)=β0 σ–(x) + β1 σ+(x) σ–(x1) + β2 σ+(x1 ) with σ–=1/(1+ex); σ+=1/(1+e-x); x=a(t-tc) and x1= a1(t-tc1).
Figure 3: Sigmoid function as effect of control measures in infection rate β = β (t)
4.2. Model simulations
with control and data:
The results presented here were obtained using two control functions described above and shown in see Fig. 3. Although still difficult to disentangle the possible seasonal effect from the draconian social distancing measures implemented during the lockdown, the SHARUCD model is unable to describe the most recent data without assuming low seasonal effect, showing a considerable overestimation of hospitalizations and deceased cases by assuming differences in transmission rate restricted to the social behaviour only.
In good agreement, the refined model can describe well the hospitalizations, the ICU admissions and the deceased cases, well matched within the median of the 200 stochastic realizations (incidences are shown in Fig. 4 and cumulative cases are shown in Fig. 5). The cumulative incidence for all positive cases follows the higher stochastic realizations range and the deviation observed between model simulations and data can be explained by the increasing testing capacities since March 22, 2020, followed by the introduction of rapid tests, mostly used as screening tools in nursing homes and public health workers. The current model (without testing feedback) can not describe this variable quantitatively, and would need further refinements and parametrization. The cumulative incidences for alive hospital discharges (black dots), used as a proxy for notified COVID-19 recovered individuals who were hospitalized, keeps following the lower realizations range when model simulations are obtained for the cumulative notified recovered (H+U+ξA), shown in Fig. 5 e) Figure 5 f) shows model simulations obtained for the cumulative notified recovered from hospital (H+U) only. A new transition is needed to record a proportion of recovered individuals tested positive without being admitted to the hospital (χA) in order to describe quantitatively describe the available data.
Figure 4: SHARUCD Model validation from March 4 to June 30 (incidences) Data are plotted as black line. The mean of 200 stochastic realizations are plotted as a blue line. The 95% confidence intervals are obtained empirically from 200 stochastic realizations and are plotted as purple shadow.
Figure 5: Ensemble of stochastic realizations of the refined SHARUCD-model. Model validation via data matching from March 4, 2020 to June 30, 2020. In a) cumulative tested positive cases Icum(t) (PCR+ only), b) cumulative hospitalized cases CH (t), c) cumulative ICU admissions CU(t), and d) cumulative deceased cases D (t). Using data for “hospital discharges alive”, we match in e) the simulations for the cumulative notified recovered CR(t) for H+U+ξA. In f) we plot the simulations for the cumulative notified recovered CR(t) for H+U only.
4.3. Model Validation and short-term prediction exercise with control measure:
Model validation and short-term predictions considering the effective control measures described above are shown in Fig. 6. For this exercise, empirical data available up to June 30, 2020, and model simulations are obtained for a four weeks longer run than the available data. New data will be included to check the quality of the short prediction exercise.
Figure 6: Mean and the 95% CI for the ensembles of stochastic realizations of the SHARUCD-model with control and data matching, starting from March 4, 2020. In a) cumulative tested positive cases Icum(t) (PCR+ only), b) cumulative hospitalized cases CH(t), c) cumulative ICU admissions CU(t), and d) cumulative deceased cases D (t). Using data for “hospital discharges alive”, we match in e) the simulations for the cumulative recovered CR(t) for H+U+ξA and in f) we plot the simulations for the cumulative recovered CR(t) for H+U only.
4.4. Growth rate and reproduction ratio:
The initial exponential growth rate λ of an epidemic is an important measure that follows directly from data at hand, commonly used to infer the basic reproduction number r, often called R0. While the momentary growth rate follows directly from the time continuous data at hand (see Fig. 7), the momentary reproduction ratio depends on the notion of a generation time γ-1.
Figure 7: Growth rate estimation from the data on all PCR positive tested cases for various variables: tested positive cases in yellow, hospitalizations in red, ICU admissions in purple, hospital discharges alive (proxy for recovered) in green and deceased in black.
The momentary growth rates and momentary reproduction numbers for the tested positive cases are calculated using data at hand available, from March 4 to June 39, 2020, shown in Fig. 8 a-b). Fig. 8 c) shows the behavior of the three synchronized variables, Icum, CH, and CU, crossing the threshold to a negative growth rate on April 1st, 2020, confirming the observed r(t) trend obtained by looking at data on Icum alone. Fig. 8 d) shows the synchronization of recovered and deceased cases, following toward the negative growth rate 1-2 weeks later, due to the delay between onset of symptoms, hospitalization and eventually death, reaching negative growth rate on April 7 and April 11, 2020, respectively. These results led to information about how to refine the model in order to capture the dynamics of the ICU admissions, better than in the original model, described in section 4.
Figure 8: a) Growth rate estimation from the data on positive tested infected cases, and b) reproduction ratio from the same data, with different recovery periods. In c) combined growth rate for tested positive cases (yellow), hospitalizations (red) and ICU (purple) and in d) combined growth rate for recovered (green) and deceased (black). Clearly, the two groups can be distinguished.
Although the absolute value of r(t) can vary given the modeling assumptions considered during its estimation, it is advised to rather look at the threshold behavior, using the growth rates primarily as it is independent of modeling uncertainties and clearly indicates if the outbreak is under control or not when estimations are below or above threshold.
The growth rates for the tested positive cases have oscillated, crossing the threshold in middle June when some isolated outbreaks were identified and controlled. At the moment, the growth rate is negative, indicating a decrease of transmission. The reproduction ratio r(t) is also estimated to be below the threshold behavior of r=1, and careful monitoring of the development of the outbreak is needed. Important to note that the increased testing capacity allows the identification of mild/asymptomatic cases, is influencing the reproduction ratio behaviour.
Complementary measures of growth rates, shown in Fig. 9, for different variables such as hospitalization, intensive care units (ICU) admissions and deceased, can be also evaluated as data is consistently collected and defined. The evaluation of those measures together balance the interpretation of the behaviour of disease transmission in the Basque Country (see Fig. 9).
Figure 9: Growth rate estimation for various variables. In a) hospitalizations (red) and in c) ICU admission cases (purple). In c) growth rate for recovered (green) and in d) deceased (black).
5. The refined SHARUCD MODEL with import and seasonality (Second phase):
The current SHARUCD framework evaluates the seasonal effect on COVID-19 dynamics in the Basque Country, after social distancing measures start to be lifted (from May 4 – phase 0 and from May 11 -phase 1 towards the “new normality”). An import factor is also included in the model dynamics, after the full lockdown lifting in July. Although this factor was not important during the exponential growth phase, imported infections play a major role when small numbers are detected (stochastic phase pre/post exponential phase). It refers to infected individuals (most likely asymptomatic) coming from outside the studied population (either an infected foreigner visiting the region or an infected Basque returning to the country and that are not detected by the current testing strategy.)
Note that when the community transmission is under control (social distancing, masks and hygienic measures), the import factor does not contribute significantly to the epidemic, only starting isolated outbreaks (variable sizes depending on the momentary infection rate), but not driving the current epidemic into a new exponential growth phase.
Figure 10: Figure 10: In a) COVID-19 in the Basque Country figure From March 4 to April 7, 2021, positive PCR cases (in yellow), Hospitalizations (in red), ICU admission (in purple) and deceased cases (in black) are shown. In b) COVID-19 tests performed in the Basque Country (in light blue) and positivity rates (red line).
5.1 Model simulations and short term predictions:
The model shown in section 4 analyses now isolated outbreaks (assuming now import to asymptomatic infection) after lifting of lockdowns and increased detection of asymptomatic due to now intensified contact tracing with time dependent infection rate β, and now also increased import ρ and increased detection of asymptomatic cases ξ over time.
Model parametrization is now also considering the detected positive cases data, in order to disentangle the role of import ρ from community spreading β. From September 15 onwards the daily testing capacity is varying significantly, however the cumulative positivity rates are stabilized (≈ 6.3%). The model predicts an endemic state for positive cases and it is, at the moment, overestimating the mean of positive cases detected using PCR tests. That is because the current assumed detection rate parameter has decreased (and keeps varying) and must be adjusted. More tests would eventually identify those asymptomatic cases predicted by the model, as observed in the last few days with the number of tests increasing significantly. Antigen tests were introduced on Oct. 17, 2020, with an important increase in testing capacity (See Fig. 10b)).
As the trend of stationarity continues with import cases (i.e. non detected asymptomatic infections) possibly causing isolated outbreaks, the community transmission is below the threshold, i.e controlled and not increasing to a new exponential phase in the next few days. Note that on October 27, a new lockdown was implemented and transmission rate has decreased even more. Although the trend of stationarity does not appear to be changed, a new control function was included in order to describe the current epidemiological scenario of decreasing detected positive cases.
5.2 Short term predictions of COVID-19 in the Basque Country:
Fig. 11: Ensemble of stochastic realizations of the refined SHARUCD-model with import. In a) cumulative positive cases I cum (t), in b) cumulative hospitalized cases C H (t), in c) cumulative ICU admissions C U (t) and in d) cumulative deceased cases D(t). The mean of 500 stochastic realizations are plotted as a blue line._
Fig. 12: 20 days predictions. In a) hospitalizations and in b) ICU admission. In c) deceased cases and in d) detected PCR positive cases.
Predictions are made with the refined SHARUCD model parameterized with empirical data described above. Deceased cases are estimated from the current data referring to detected positive cases and severe cases prone to hospitalization and are possibly been overestimated by the model (as PCR method can identify nucleic acid from the virus but can not differentiate between active or residual viral particles). Model predictions are shown as light blue curves. The 95% confidence intervals (CI) are plotted in light purple (shadow). Empirical data are plotted on top of the prediction curve in black.
6. Acknowledgments and contribution:
Maíra Aguiar, BCAM – Ikerbasque, MTB research line leader & University of Trento Marie Curie Fellow (COMPLEXDYNAMICS-PHIM Marie Skłodowska-Curie grant agreement No 792494): conceived of the study, developed and analyzed the models.
Nico Stollenwerk, BCAM, MTB Group & University of Trento: conceived of the study, developed and analyzed the models.
Nicole Cusimano, BCAM – MTB Postdoctoral researcher: data collection, analysis and graphical illustration.
Marina Echeverría Ferrero, BCAM – MTB PhD Student (Predoc Severo Ochoa 2018): data collection and graphical illustration.
Damián Knopoff, BCAM – MTB Postdoctoral researcher: data collection, analysis and graphical illustration.
Francisco P. Rodrigues, student at Instituto Superior Técnico, University of Lisbon, Portugal, has developed the SIR simulation program.
Eduardo Millán, Osakidetza Basque Health Service, has collected for the Basque Country and prepared the data sets used for the SHARUCD modeling exercise.
We thank the huge efforts of the whole COVID-19 BMTF, specially to Joseba Bidaurrazaga Van-Dierdonck, Public Health, Basque Health Department, Javier Mar, Biodonostia Health Research Institute and Adolfo Morais Ezquerro, Vice Minister of Universities and Research of the Basque Government, for the fruitful discussions.
 World Health Organization. Emergencies preparedness, response. Novel Coronavirus – China.
 Maíra Aguiar, Nico Stollenwerk and BobW. Kooi. (2012). Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study. Chapter In book: Epidemiology Insights.
 FredBrauer. (2005). The Kermack–McKendrick epidemic model revisited. Mathematical Biosciences, 198(2), 119-131.
 Stollenwerk, N. & Jansen, V. (2010). Population biology and criticality, ISBN 978-1-84816-401-7. (Imperial College Press, London).
 Aguiar, M., Ortuondo, E.M., Bidaurrazaga Van-Dierdonck, J. et al. Modelling COVID 19 in the Basque Country from introduction to control measure response. Nature Scientific Reports 10, 17306 (2020).
 Aguiar, M. & Stollenwerk, N. SHAR and effective SIR models: from dengue fever toy models to a COVID-19 fully parametrized SHARUCD framework. Commun. Biomath. Sci. 3(1), 60–89 (2020).
 van Kampen, N. G. (1992). Stochastic Processes in Physics and Chemistry, ISBN 978-0-44452-965-7. (North-Holland, Amsterdam).
 Gillespie, D.T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics, 22, 403-434, ISSN 0021-9991.
 Gillespie, D.T. (1978). Monte Carlo simulation of random walks with residence time dependent transition probability rates. Journal of Computational Physics, 28, 395-407, ISSN 0021-9991.
 Aguiar, M., Van-Dierdonck, J. B. & Stollenwerk, N. (2020). Reproduction ratio and growth rates: measures for an unfolding pandemic. PLoS ONE 15(7), e0236620.
 Aguiar, M. & Stollenwerk, N. (2020). Condition-specific mortality risk can explain differences in COVID-19 case fatality ratios around the globe. J. Public Health, 188, 18-20.
 Maíra Aguiar, Joseba Bidaurrazaga Van-Dierdonck, Javier Mar, Nicole Cusimano, Damián Knopoff, Vizda Anam & Nico Stollenwerk. (2021).The interplay between subcritical fluctuations and import: understanding COVID-19 epidemiology dynamics. Pre-print available Repository: medRxiv