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 [3].

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 [5], 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 [5].

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 [9] 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.
• 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.

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” C_H, “asymptomatic” C_A, recovered C_R and ICU patients C_U. 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 x₁:=S/N , x₂:=H/N, x₃:=A/N , x₄:=R/N , x₅:=U/N , x₆:=CH/N, x₇:=CA/N , x₈:=CU/N and x₉:=D/N and x₁₀:=CR/N, hence state vector x := (x₁,…,x₁₀)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 C_H, C_A, C_R and C_U 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+e^x) . 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)=β₀ (1+ θ cos(ω(t+φ))) on top of the already implemented sigmoidal control function σ(x)=1/(1+e^x) from what we obtain the final format  β(t)=β₀ σ_–(x) + β₁ σ₊(x) σ_–(x₁) + β₂ σ₊(x₁ ) with σ_–=1/(1+e^x); σ₊=1/(1+e^-x); x=a(t-t_c) and x₁= a₁(t-t_c1).

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.