Modelling COVID 19 in the Basque Country from introduction to control measure response

In March 2020, a multidisciplinary task force (so-called Basque Modelling Task Force, BMTF) was created to assist the Basque health managers and Government during the COVID-19 responses. BMTF is a modelling team, working on different approaches, including stochastic processes, statistical methods and artificial intelligence. Here we describe the efforts and challenges to develop a flexible modeling framework able to describe the dynamics observed for the tested positive cases, including the modelling development steps. The results obtained by a new stochastic SHARUCD model framework are presented. Our models differentiate mild and asymptomatic from severe infections prone to be hospitalized and were able to predict the course of the epidemic, providing important projections on the national health system’s necessities during the increased population demand on hospital admissions. Short and longer-term predictions were tested with good results adjusted to the available epidemiological data. We have shown that the partial lockdown measures were effective and enough to slow down disease transmission in the Basque Country. The growth rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \lambda $$\end{document}λ was calculated from the model and from the data and the implications for the reproduction ratio r are shown. The analysis of the growth rates from the data led to improved model versions describing after the exponential phase also the new information obtained during the phase of response to the control measures. This framework is now being used to monitor disease transmission while the country lockdown was gradually lifted, with insights to specific programs for a general policy of “social distancing” and home quarantining.

The Fokker-Planck equation gives a stochastic differential equation system with σ = 1/ √ N and Gaussian normal noise vector ε(t) = (ε x 1 (t), ..., ε x 10 (t)) tr as and using matrix square root from eigenvalue-eigenvector decomposition G 2 (x) = T ΛT −1 as G(x) = T √ ΛT tr to be numerically implemented easily, and much faster than the Gillespie algorithm for the master equation, when it comes to large population sizes N and longer runs of e.g. one year, as for a complete disease outbreak curve necessary. For possibly non-quadratic matrices B, expressing the covariance matrix as B tr B = G 2 , see [1], speeding up further the stochastic process simulations.
In mean field approximation we obtain explicitly from with the transitions w j (x)) specified in Eq. 4 in the main text. The deterministic version of the 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 in a complete form.
For a constant population size N , susceptible individuals (S) become infected with SARS- ratio of mild/asymptomatic infections contributing to force of infection isolated and would no longer transmit as much as mild/asymptomatic cases, often undetected and more mobile. Infected individuals may recover with a recovery rate γ, however, severe hospitalized individuals could be also admitted to the ICU facilities, with a rate ν, or die before being admitted to the ICU facilities, with disease induced death rate µ. Individuals admitted to the ICU facility could eventually recover with a recovery rate γ or die with disease induced death rate µ.
The model is calibrated using the empirical data for the Basque Country and the biological parameters are estimated and fixed as the model is able to describe the disease incidence during the exponential phase of the epidemic for each dynamical class. Parameter insecurities are quantified with likelihood functions. Model parameters and initial conditions are shown in Table 1, where β is the infection rate and φ is the ratio describing the asymptomatic/mild infections contribution to the force of infection. γ is the recovery rate, µ is the disease induced death rate and ν is the rate of hospitalized going to the ICU. η is the proportion of susceptible being infected, develop sever symptoms and being hospitalized whereas 1 − η is the proportion of susceptible becoming infected and developing mild disease or asymptomatic. ξ is the ratio of detected, via testing, mild/asymptomatic infect individuals. is the import rate needed to describe the introductory phase of the epidemics and for the present study, we assume to be much smaller than the other additive terms of the force of infection, given the strong observational insecurities on the data collected at the beginning of the outbreak.

The refined model
The original SHARUCD model was refined to synchronize ICU admissions to hospitalizations and positive tested infected, rather than to deceased and recovered. 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.
We change the transition into ICU admissions from the form like used in recovery γ and death µ, more to the one for distinction of hospitalized and asymptomatics η, to describe the results from the growth rate analysis of all data sets, see previous subsection. Hence we update the transitions just a bit from The other parameters values were kept the same as shown in Table 1.

Two parameter likelihood plots
We analyzed possible correlations between parameters by inspecting numerical two-parameter likelihood functions. Since the growth factor λ is mainly determined as λ ≈ β −γ we observe in  Figure S1: Two-parameter likelihood plot L(β, γ) for parameters the infection rate β and the recovery rate γ, showing fewer than expected correlations between these two parameters which essentially determin the common growth factor λ 1 . This well determined range of possible values for the parameters β and γ might be due to the available information from all 5 data sets, including the recovered in good quality.
many basic epidemiological models large correlations between especially these two parameters.
Large infectivity β and small recovery period γ −1 describe often data as well as small infectivity and long recovery period. However we observe in a first inspection of the numerical likelihood of these two parameters, obtained from considering all 5 paremter sets in comparison with the SHARUCD basic model, that the range of possible values for β and γ together, see in Fig. S1 the area lifted away from zero in L(β, γ), is quite restricted around the best values for each of the parameters individually, as they are shown in Fig. 3 a) (η, ξ), and c) for L(φ, ξ), the parameters which describe more internal effects of hospitalization ratio η versus asymptomatic/mild, and infectivity of mild/asymptomatic φ, and detection rate of mild/asymptmatic ξ. Largest correlations between parameters are observed between η and ξ, see Fig. part b). are considerable correlations between these two parametes detected, but they seem to be mild, and from the available data we can obtain good information also about the most likely values of these parameters. Finally, for the parameter combination of L(φ, ξ) we observe again a more restricted area of possible parameter sets, though also here some correlations are visible. In conclusion, it seems that the available data sets give quite some information about the possible parameter combinations to describe the system under investigation. Further refinements of these analyses might give further insight into the dynamics of the epidemic.
Due to initially small numbers of infected, hospitalized cases etc., we use initially the Gillespie algorithm for the state discrete Markov process described by the master equation [2,3]. The present two parameter likelihood plots become occasionally very computationally demanding, depending on system size N and numbers of transitions, especially for larger infectivity. hence the well tested approximation via the Fokker-Planck equation as stochastic differential equation system [4] speeds up the simulations, allowing for smoother likelihoods, but occasionally might give some errors in the small number situations. A carefull monitoring is therefor needed, see a first comparison of master equation likelihoods, using Eq. 3 from the main text, with stochastic differential equation likelihoods, using Eq. (6), with encouraging results Fig. S3, but need to be further investigated in future studies. However, the first results given here in Figs. S1 and S2 are already informative for the purpose of detecting eventual correlations between parameters, and encouraging for further intensive studies. The up to now detected correlations keep well within the expected ranges indicated from the one parameter numerical likelihood plots in Fig. 3 in the main text.