#### 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_{1}:=S/N , x_{2}:=H/N, x_{3}:=A/N , x_{4}:=R/N , x_{5}:=U/N , x_{6}:=CH/N, x_{7}:=CA/N , x_{8}:=CU/N and x_{9}:=D/N and x_{10}:=CR/N, hence state vector x := (x_{1},...,x_{10})tr, giving the dynamics for the probabilities p(x,t) as

For a more detailed description of the SHARUCD model framework, see [7]. 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 η. The improved transitions are given in blue in the ODE system above.

#### 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 σ=1/(1+e^{x}), 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 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 the sigmoidal control function σ(x)=1/(1+ex) from what we obtain the final format

β(t)=β_{0} σ_{-}(x) + β_{1} σ_{+}(x) σ_{-}(x_{1}) + β_{2} σ_{+}(x_{1})

with σ_{-}=1/(1+e^{x}); σ_{+}=1/(1+e^{-x}); x=a(t-tc) and x_{1}= a1(t-tc1).

#### 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.

#### 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 I_{cum}(t) (PCR+ only), b) cumulative hospitalized cases C_{H}(t), c) cumulative ICU admissions C_{U}(t), and d) cumulative deceased cases D (t). Using data for “hospital discharges alive”, we match in e) the simulations for the cumulative recovered C_{R}(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 R_{0}. 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, I_{cum}, 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 I_{cum} 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**.

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 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).