Table of contents:
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 11 million cases were confirmed with about 530 thousand deaths, with a global case fatality ratio (CFR) of approximate 5% .
Figure 1: COVID-19 pandemic figures. In a) bar chart race for the confirmed cases around the globe. The global COVID-19 deceased cases are shown in b) daily deceased cases and in c) total deceased cases.
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 webpage 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 usd 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 COVID19 epidemiology. Graphics will be updated every two weeks.
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).
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 [10,11]. The deterministic approach is obtained via the mean field approximation. For more details on the calculations see e.g. .
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.
- Interpretation: 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.
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.
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 famework, see . 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 progession to deceased, but then due to the observed synchronization of transitions to hospitalization and asymptomatic cases (see section 4.4), we refined the model, changimg 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 η. The improved transitions are given in blue in the ODE system above.
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 σ=1/(1+ex), shown in Fig. 3. 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 includes 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 the 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).
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 tool 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 the available data.
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 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 deceases 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.
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 he 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 let 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.
While we do not take responsibility for the absolute value of r(t) as it is bound to many internal assumptions, recovery period, smoothing and approximations, which are not valid for time dependent parameters, the threshold behavior is independent of those uncertainties and clearly indicates that the outbreak is under control or not, when estimations are respectively below or above 1. The growth rates for the tested positive cases have oscillated, crossing the threshold in midle 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 unites (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).
Maira Aguiar, BCAM – Ikerbasque, MTB research line leader/ Marie Curie Fellow at University of Trento, Italy, has conceived of the study, developed and analyzed the models.
Nico Stollenwerk, University of Trento, Italy, has conceived of the study, developed and analyzed the models.
- Nicole Cusimano, Post doctoral researcher of Mathematical and Theoretical Biology group (MTB). Nicole has produced the animated graphics for the SHARUCD model short predictions for the 5 dynamical variables, PCR tested positives, hopsitalizations, ICUs, deceased and recovered and the animated graphic for the combined growth rates.
- Marina Echeverría Ferrero, PhD Student (Predoc Severo Ochoa 2018) of Mathematical and Theoretical Biology group. Marina has collected the data and produced the bar chart race for the confirmed COVID-19 cases around the globe as well as the static plot for the global number of deceased cases.
- 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, 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. Retrieved from https://www.who.int/csr/don/12-january-2020-novel-coronavirus-china/en/
 World Health Organization. WHO announces COVID-19 outbreak a pandemic. Retrieved from http://www.euro.who.int/en/health-topics/health-emergencies/coronavirus-covid-19/news/news/2020/3/who-announces-covid-19-outbreak-a-pandemic
 World Health Organization. Coronavirus disease (COVID-2019) situation report 167 https://www.who.int/docs/default-source/coronaviruse/situation-reports/20200705-covid-19-sitrep-167.pdf?sfvrsn=17e7e3df_
 Maíra Aguiar, Nico Stollenwerk and BobW. Kooi. (2012). Modeling Infectious Diseases Dynamics: Dengue Fever, a Case Study. Chapter In book: Epidemiology Insights. DOI: 10.5772/31920
 FredBrauer. (2005). The Kermack–McKendrick epidemic model revisited. Mathematical Biosciences, 198(2), 119-131.
 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.
 Stollenwerk, N., \& Jansen, V. (2010). Population biology and criticality, ISBN 978-1-84816-401-7. (Imperial College Press, London). https://www.medrxiv.org/content/10.1101/2020.05.10.20086504v3
 Maíra Aguiar, Eduardo Millan Ortuondo, Joseba Bidaurrazaga Van-Dierdonck, Javier Mar, Nico Stollenwerk. (2020). Modeling COVID 19 in the Basque Country: from introduction to control measure response. Available at https://doi.org/10.1101/2020.05.10.20086504
 Maíra Aguiar, Joseba Bidaurrazaga Van-Dierdonck, Nico Stollenwerk. (2020). Reproduction ratio and growth rates: measures for an unfolding pandemic. Available at https://doi.org/10.1101/2020.05.18.20105528