Modeling the initial phase of COVID-19 epidemic: The role of age and disease severity in the Basque Country, Spain.

Declared a pandemic by the World Health Organization (WHO), COVID-19 has spread rapidly around the globe. With eventually substantial global underestimation of infection, by the end of March 2022, more than 470 million cases were confirmed, counting more than 6.1 million deaths worldwide. COVID-19 symptoms range from mild (or no) symptoms to severe illness, with disease severity and death occurring according to a hierarchy of risks, with age and pre-existing health conditions enhancing risks of disease severity. In order to understand the dynamics of disease severity during the initial phase of the pandemic, we propose a modeling framework stratifying the studied population into two groups, older and younger, assuming different risks for severe disease manifestation. The deterministic and the stochastic models are parametrized using epidemiological data for the Basque Country population referring to confirmed cases, hospitalizations and deaths, from February to the end of March 2020. Using similar parameter values, both models were able to describe well the existing data. A detailed sensitivity analysis was performed to identify the key parameters influencing the transmission dynamics of COVID-19 in the population. We observed that the population younger than 60 years old of age would contribute more to the overall force of infection than the older population, as opposed to the already existing age-structured models, opening new ways to understand the effect of population age on disease severity during the COVID-19 pandemic. With mild/asymptomatic cases significantly influencing the disease spreading and control, our findings support the vaccination strategy prioritising the most vulnerable individuals to reduce hospitalization and deaths, as well as the non-pharmaceutical intervention measures to reduce disease transmission.

1 Introduction

More than two years have passed since COVID-19, a severe respiratory syndrome caused by a new coronavirus, was identified by the Chinese authorities in January 2020 [1]. Declared a global pandemic by the World Health Organization (WHO) in March 2020 [2], COVID-19 symptoms range from asymptomatic/mild to severe illness, with age and pre-existing health conditions increasing the likelihood of disease severity [3]. Vaccines against COVID-19 have been developed in record time and are now globally distributed [4, 5]. Although these vaccines are remarkably effective against severe disease, the so called sterilizing immunity, occurring when vaccinated individuals cannot transmit the virus, is still being evaluated.

Based on previous research experiences applied to other infectious diseases [6-14], and more recently applied to COVID-19 dynamics [15-19], the role of asymptomatic infections have been studied, showing that vaccine performance is driven by the ability of asymptomatic or mild disease cases transmitting the virus, with an eventual increase on the number of overall infections in a population [20, 21].

As an example of the pandemic’s impact in Europe, Spain has reported, by the end of March 2022, more than 11.5 million COVID-19 cases and over 100 thousand deaths [22, 23], with a significantly higher mortality rate for individuals older than 65 years of age [24, 25], in agreement with what was also observed in different European countries [26].

As the COVID-19 pandemic progressed, task forces have been created to assist public health managers and governments during the COVID-19 crisis, and research on mathematical modeling became critical to understand the epidemiological dynamics of COVID-19. Modeling studies to evaluate COVID-19 dynamics worldwide have been widely published. Using both, deterministic and stochastic approaches, models were developed to investigate disease spreading in different epidemiological contexts as well as the impact of the control measures so far implemented. Using the existing empirical data, these models have given insights on disease transmission rates, the effect of quarantine or use of facial masks, for example, with modeling assumptions statistically tested with the available empirical data [27-30].

Within the COVID-19 Basque Modeling Task Force (BMTF), a flexible stochastic framework was developed to describe the epidemics in terms of disease spreading and control in the Basque Country, Spain, giving projections on the national health system needs over time. The SHARUCD framework was parameterized and validated with epidemiological data continuously collected and provided by the Basque Health Department and the Basque Health Service (Osakidetza), and has been used, up to date, to monitor COVID-19 spreading and control over the course of the pandemic [15-21]. Model refinements and results on the evolution of the epidemics in the Basque Country are updated on a monthly basis and are publicly available as an online dashboard [5].

As a continuation of the BMTF efforts, we developed an age-stratified mathematical model framework to understand the epidemiological dynamics of COVID-19 introduction phase in the Basque Country. The models are calibrated with the available data referring to confirmed cases, hospitalizations and deaths, from February to the end of March 2020, in the Basque Country, prior to any intervention measure. After a careful data analysis, the population was divided into two groups, namely young and old. As opposed to the existing age structured models suggesting higher infection rate for individuals older than 60 years of age [31-34] than for younger individuals, our modeling assumption implies that while the risk for developing severe disease is higher for the older population, disease transmission is significantly driven by the mobile younger population.

A detailed sensitivity analysis was performed to identify the key parameters influencing the transmission dynamics of COVID-19 in the population, opening new ways to understand the effect of age on disease severity during the pandemic. In terms of policy implications, our findings support the vaccination strategy prioritising the most vulnerable individuals [16], particularly to reduce hospitalization and deaths [21], as well as the non-pharmaceutical intervention measures that are still advised by the WHO to reduce disease transmission.

This paper is organized as follows. Section 2 presents the deterministic and the stochastic models formulation, followed by the model analysis. Section 3 is dedicated to data analysis, model calibration and parameter estimation. In Section 4 we present the models simulation and results, including a detailed sensitivity analysis for the parameters involved in reproduction number. We conclude this work with a discussion on the results obtained by both modeling approaches.

2 Materials and methods

Using age stratified data for COVID-19 incidences for tested positive cases, hospitalizations and deaths in the the Basque Country, this work is applied to the initial phase of the pandemic. Using statistical tools to analyse these data, we define as severe cases all hospitalized individuals, including the intensive care unit (ICU) admissions, for young (H1) and old patients (H2), reported from February 15 to March 25, 2020. It is important to mention that at the beginning of the pandemic, due to the testing capacity limitations, only patients with severe symptoms were tested using the PCR (Polymerase Chain Reaction) method.

2.1 The deterministic model

This model framework is a refinement of the model proposed by Srivastav et al. [27, 35]. For both age groups, young and old, susceptible individuals become exposed and infected E1(t) and E2(t), developing either mild/asymptomatic A1(t) and A2(t) or severe/hospitalized H1(t) and H2(t) disease. While mild/asymptomatic infections are assumed to recover, severe disease might evolve to death D. The parameter ϕ differentiates the disease transmission between hospitalized (H1 + H2) and mild/asymptomatic infections (A1 + A2), and the parameter ϵ is introduced to differentiate the infectivity of asymptomatic young individuals (A1) with respect to the baseline infectivity for the elderly individuals A2(t) in the Basque Country population of N = 2.6 million individuals.

The seriousness of symptoms from viral infections is often correlated with the amount of the virus in the body [36, 37]. Justified by the differences observed in viral load during the COVID-19 infection, lower for mild/asymptomatic and higher for severe/hospitalized cases, we assume ϵ < 1, indicating that young individuals have smaller infectivity than the elderly individuals. This assumption relies on the epidemiological observation of young individuals developing mild or no symptoms during the infection as opposed to the observation of severe symptoms occurring mostly in older ages, shaping the disease transmissibility pattern in a population. The parameter ϕ is a scaling factor used to differentiate the infectivity of mild/asymptomatic infections (ϕβ) with respect to the baseline infectivity of severe/hospitalized cases (β). The value of ϕ can be tuned to reflect different situations: a value of ϕ < 1 reflects the fact that severe cases have larger infectivity than mild cases (e.g., due to enhanced coughing and sneezing), while ϕ > 1 indicates that asymptomatic individuals and mild cases contribute more to the spread of the infection (e.g., due to their higher mobility and possibility of interaction) than the severe cases which are more likely to be detected and isolated [15]. Here, we assume ϕ > 1, with asymptomatic individuals contributing more to the force of infection than the hospitalized individuals [38, 39].

The total population N is divided into ten compartments, stratified into two age groups, young and old. Susceptible S1(t) and S2(t), Exposed E1(t) and E2(t), mild/Asymptomatic A1(t) and A2(t) or severe/Hospitalized H1(t) and H2(t) cases. Labels 1 and 2 refers to the young and to the old age populations respectively. Two extra classes to accommodate individuals from both age-groups are also considered. The deceased class D(t), for those who died from COVID-19, and finally the recovered class R(t), counting all individuals recovered from the disease.

For the mathematical modelling framework development, we make the following assumptions:

The total population N is constant.

The susceptible young individuals S1 become exposed to the infection E1 by contacting infectious individuals A1, A2 and H1, H2 at rates ϕβ and β, respectively.

The susceptible old individuals S2 become exposed to the infection E2 by contacting infectious individuals A1, A2 and H1, H2 at rates ϕβ and β, respectively.

With i = 1, 2, for young and old respectively, exposed individuals E will develop mild/asymptomatic infection A with rate aη while the remaining individuals developing severe symptoms will be admitted to a hospital facility H with rate (1 − a)η.

While young and old asymptomatic individuals recover from COVID-19 infection (R) with rate α1 and α3 respectively, hospitalized young individuals will recover with rate α2 while hospitalized old individuals will recover with rate α4. Young and old hospitalized individuals will eventually die (D) with rate δ1 and δ2 respectively. The description of model framework parameters can be found in Table 1

Table 1

Description of model framework parameters.

Parameter	Description
β:	baseline COVID-19 transmission rate
ϕ:	scaling factor used to differentiate the infectivity of severe/hospitalized cases
ϵ:	scaling factor used to differentiate the infectivity of young and elderly mild/asymptomatic cases
δ1:	disease induced death rate for hospitalized young individuals
δ2:	disease induced death rate for hospitalized old individuals
η1:	hospitalization rate for young individuals
η2:	hospitalization rate for old individuals
α1:	recovery rate of asymptomatic young individuals
α3:	recovery rate of asymptomatic old individuals
α2:	recovery rate of hospitalized young individuals
α4:	recovery rate of asymptomatic old individuals
a:	Fraction of exposed population developing mild/asymptomatic disease
(1 − a):	Fraction of exposed population developing severe/hospitalized disease

The flow diagram for the disease related stages of our proposed model is shown in Fig 1, which translates into the following ODE system describing the temporal evolution of the number of individuals in each of the model compartments:

Fig 1

With ρ1(t) = βS1[ϕ(A1 + ϵA2) + (H1 + H2)] and ρ2(t) = βS2[ϕ(A1 + ϵA2) + (H1 + H2)], disease related stages are shown in orange color for young population and in light green for the old population.

Deceased and recovered population include both age groups and are shown in black and purple, respectively.

2.2 Existence of equilibrium points and the basic reproduction number (R0)

While the disease-free equilibrium of the system is given by the basic reproduction number R0 can be found by using the next generation matrix method [40], and is given by: The quantity is defined for the young group population and the quantity is defined for the old group population. The quantity R0 = R1 + R2 is the average number of secondary cases produced in a completely susceptible population by an index case, during the infectious period.

The calculation of the basic reproduction number R0 is shown in the S1 File. We can summarize our findings in the following theorems.

Theorem 1.1 If R0 < 1, the disease-free equilibrium of the system (1) is locally asymptotically stable, and if R0 > 1, the disease-free equilibrium E0 is unstable.

Next, we state globally asymptotically stability of disease-free equilibrium.

Theorem 1.2 If R0 < 1, the disease-free equilibrium of the system (1) is globally asymptotically stable whenever eigenvalue of the matrix F − V are having negative real parts, and if R0 > 1, the disease-free equilibrium E0 is unstable [41].

The proof of the global stability can be found in the S1 File.

2.3 The stochastic model

As all natural systems are prone to stochastic fluctuations, we extended our deterministic model, see Equation System 1, to the corresponding stochastic model. The derivation of the stochastic model and its analysis are important when populations are small, and hence with the dynamics being severely affected by small changes in the parameter values. Thus, for the initial phase of the COVID-19 outbreak, the stochastic model setup is the most appropriate modeling approach to be used for a local epidemiological evaluation.

The derivation of a stochastic differential equation (SDE) model is a diffusion approximation from the underlying state discrete Markov process [17, 41–45]. Let be a continuous random variable for where T denotes transpose of the matrix. Further, let ΔX = X(t + Δt) − X(t) = (ΔX1, ΔX2, ΔX3, ΔX4…) denotes the random vector for the change in random variables during time interval Δt. Here, we write the transition maps which define all possible changes between disease states in the SDE model. State changes and their probabilities are presented in Table 2, followed by the full SDE system 2.

Table 2

Possible changes of states and their probabilities.

Possible state change	Probability of state change
(ΔX)1 = (−1, 1, 0, 0, 0, 0, 0, 0, 0, 0)T Change when young susceptible meet infected individuals and move to the young exposed class	P1 = βX1[ϕ(X3 + ϵX7) + (X4 + X8)]Δt + O(Δt)
(ΔX)2 = (0, −1, 1, 0, 0, 0, 0, 0, 0, 0)T Change when fraction of young exposed become infectious and move to the young asymptomatic infected class	P2 = aη1 X2Δt + O(Δt)
(ΔX)3 = (0, −1, 0, 1, 0, 0, 0, 0, 0, 0)T Change when fraction of young exposed become infectious and move to the young hospitalized class	P3 = (1−a)η1 X2Δt + O(Δt)
(ΔX)4 = (0, 0, −1, 0, 0, 0, 0, 0, 1, 0)T Change when young asymptomatic infected recovers and move to the recovered class	P4 = α1 X3Δt + O(Δt)
(ΔX)5 = (0, 0, 0, −1, 0, 0, 0, 0, 0, 1)T Change when young hospitalized die and move to the deceased class	P5 = δ1 X4Δt + O(Δt)
(ΔX)6 = (0, 0, 0, −1, 0, 0, 0, 0, 1, 0)T Change when young hospitalized individuals recover and move the recovered class	P6 = α2 X4Δt + O(Δt)
(ΔX)7 = (0, 0, 0, 0, −1, 1, 0, 0, 0, 0)T Change when old susceptible meet infected individual and move to the old exposed class	P7 = βX1[ϕ(X3 + ϵX7) + (X4 + X8)]Δt + O(Δt)
(ΔX)8 = (0, 0, 0, 0, 0, −1, 1, 0, 0, 0)T Change when fraction of old exposed become infectious and move to the old asymptomatic infected class	P8 = aη2 X6Δt + O(Δt)
(ΔX)9 = (0, 0, 0, 0, 0, −1, 0, 1, 0, 0)T Change when fraction of old exposed become infectious and move to the old hospitalized class	P9 = (1−a)η2 X6Δt + O(Δt)
(ΔX)10 = (0, 0, 0, 0, 0, 0, −1, 0, 1, 0)T Change when old asymptomatic infected recovers and move to the recovered class	P10 = α3 X7Δt + O(Δt)
(ΔX)11 = (0, 0, 0, 0, 0, 0, 0, −1, 0, 1)T Change when old hospitalized die and move to the deceased class	P11 = δ2 X8Δt + O(Δt)
(ΔX)12 = (0, 0, 0, 0, 0, 0, 0, −1, 1, 0)T Change when old hospitalized individuals recover and move to the recovered class	P12 = α4 X8Δt + O(Δt)
(ΔX)13 = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0)T No change	P_13=1-∑i=112Δt+O(Δt)

The detailed derivation of the stochastic model can be found in the S1 File of this manuscript.

3 Data analysis and parameter estimation

3.1 Epidemiological data

Epidemiological data used in this study are provided by the Basque Health Department and the Basque Health Service (Osakidetza), continually collected with specific inclusion.

By March 4, 2022, around 600,000 cases were confirmed, with 32087 hospital admissions and 8788 deaths in the Basque Country. For the proposed model, the age stratification was decided after a careful data inspection and data fitting, followed by the parameter estimation.

We use the epidemiological data referring to the cumulative incidences of confirmed positive cases, hospitalizations, including ICU admissions, and deceased cases distributed by age groups available for the initial phase of the COVID-19 in the Basque Country, from February 15 to to March 25, 2020, as shown in Table 3.

Table 3

Cumulative disease cases by age in the Basque Country.

COVID-19 epidemiological data, from February 15 to March 25, 2020
Raw				Normalized by 105 people
age classes	positive cases	hospital admissions	deceased cases	positive cases	hospital admissions	deceased cases
0–9	19	3	0	10	2	0
10–19	34	5	0	17	3	0
20–29	188	34	1	97	18	1
30–39	388	118	2	146	45	1
40–49	600	255	4	168	71	1
50–59	796	393	6	230	118	2
60–69	714	518	20	263	191	8
70–79	638	622	44	316	308	22
80+	680	523	146	432	332	93

Note that during this period, testing capacity was limited and therefore the positive detected cases were restricted to symptomatic individuals and eventually to their close contacts during the process tracing and testing strategy.

3.2 Model calibration method

Using MATLAB software, parameter estimation was performed using nonlinear least square method [47]. In detail, we search for the set of parameters that minimizes the sum of squared differences between the observed data and the corresponding model solution denoted by (f(t, Θ) The Root Mean Square Error (RMSE) values for the deterministic and stochastic models are calculate using the following formula, where t are the time points at which the time series data are observed, and n is the number of data points available for parameter inference. Hence, the model solution f(t, Θ) yields the best-fit to the time series data .

3.3 Raw data and model fitting

These raw data distribution by age groups are shown in Fig 2.

Fig 2

From February 15 to March 25, 2020, raw data distribution for (a) total positive cases, (b) Hospital admission including ICU cases and (c) deceased cases.

During the initial phase of the pandemic, a strong correlation of positive cases and severe disease leading to hospitalizations is observed, see y-axis of Fig 2a) and 2b). Increased age appears to be a strong risk factor for developing severe illness with COVID-19 infections, however, by looking at the raw data referring to the hospital admissions, this consideration is not very clear, with similar high hospitalization rates for individuals younger than 50 years of age and individuals older than 70 years of age. Nevertheless, when looking at the deceased cases, it is indeed observed that older adults have higher risk of severe outcomes. With potential underlying health conditions [3], most of deaths occurred in those older than 70 year of age, see Fig 2c).

Aiming to understand the role of population age heterogeneity on disease transmission and severe outcomes in the absence of vaccines and other non-pharmaceutical interventions, we consider the information obtained from the raw hospitalization data for the initial phase of the pandemic. While the young group includes individuals between 0–39 years of age, the old group considers the remaining individuals in the population older than 40 years of age. Models are calibrated with the data and the parameters reflecting the differences in disease transmission by age group are estimated.

The available data referring to cumulative hospital admission cases, for the young and the old age groups, are matched with both models, deterministic and stochastic. Fig 3 shows the models fitting to the empirical raw data. In this data matching scenario, the RMSE values for the deterministic and stochastic models are 0.55 and 0.47 respectively, indicating that the stochastic modeling approach explains better the existing data. The scaling factor parameter used to differentiate the infectivity of severe/hospitalized cases ϕ, and scaling factor parameter used to differentiate the infectivity of young and elderly mild/asymptomatic cases ϵ, were estimated to be for the young group, and for the old group. The other parameter values are fixed as suggested in [17]. Referring to the raw data, the used parameter values for the data fitting are listed in Table 4.

Fig 3

On the left hand side, the deterministic model curve (blue line) and on the right hand side, the stochastic model realizations (in blue), fitting the cumulative empirical data referring to hospital admissions (red dots).

In (a) and in (b) data matching with model simulations for the young (0–39 years of age) age group. In (c) and (d) data matching with model simulations for the old (40 years and older) age group.

Table 4

Parameters values used for model calibration.

Parameter	Normalized data values (fitting)	Raw data values (fitting)	Ref.
β:	0.15	0.15	[17]
ϕ (young):	1.5 [1–2]	1.2 [1–2]	fitted
ϵ (young):	0.3 [0–1]	0.25 [0–1]	fitted
ϕ (old):	1.3 [1–2]	1.55 [1–2]	fitted
ϵ (old):	0.2 [0–1]	0.4 [0–1]	fitted
δ1:	0.003 [0.001–0.004]	0.0012 [0.001–0.004]	fitted
δ2:	0.04 [0.02–0.05]	0.025 [0.02–0.05]	fitted
η1:	0.035 [0.0–0.5]	0.035 [0.0–0.5]	[17]
η2:	0.03 [0.0–0.05]	0.03 [0.0–0.05]	[17]
α1:	0.02 [0.0–0.09]	0.02 [0.0–0.09]	[17]
α3:	0.05 [0.0–0.09]	0.05 [0.0–0.09]	[17]
α2:	0.01 [0.0–0.09]	0.01 [0.0–0.09]	[17]
α4:	0.03 [0.0–0.09]	0.03 [0.0–0.09]	[17]
a:	0.02	0.02	[17]

3.4 Normalized data and model fitting

The normalized raw data relative to the population size for each age class in the Basque Country is shown Table 3, with its visual age distribution shown in Fig 4.

Fig 4

From February 15 to March 25, 2020, normalized data distribution by age group.

The data is presented as confirmed cases per 100000 people. In (a-b) total positive cases, (c-d) hospitalized cases and (e-f) deceased cases.

Similarly to what was observed with the raw data, positive cases are increasing with age. The large majority of the deceased cases have been reported for the 80 years and older population group, confirming the strong correlation of severe disease outcome and age. Nevertheless, the normalization of the raw data shows clearer an increase of hospitalization rates for older age classes, allowing us to modify our modeling age stratification definition for young and old groups.

We summarize the distribution of disease cases using box plots to represent the deviation in the reported cases by age, see Fig 5, with the median being the measure of central tendency of the underlying distribution of the data as shown on Table 3.

Fig 5

Box plots for (a) total positive cases, (b) hospitalized cases and (c) deceased cases.

Horizontal lines denote lower quartile, median and upper quartile, with dots showing outliers.

Fig 5a) shows similar median values for individuals of 30 years and older, suggesting that they are more likely to develop symptoms than individuals at younger ages. In respect to the hospitalizations, see Fig 5b), the median values are similar for the individuals older than 50 years of age, suggesting that infections within these age groups are likely to be more severe requiring hospitalizations than for the younger ages, with individuals older than 80 years of age more likely to die from COVID-19 infection than any other age class, see Fig 5c). For these data, the age distribution assumption is now modified, considering individuals between 0–69 years of age as part of the young group and individuals older than 70 years of age as part of the old age group.

The cumulative empirical data for both age groups are matched with the deterministic system 1 and stochastic system 2 model simulations, see Fig 6.

Fig 6

On the left hand side, the deterministic model curve (blue line) and on the right hand side, the stochastic model realizations (in blue), fitting the cumulative empirical data (red dots).

In (a-b) the hospitalizations for the young group (0–69 years), in (c-d) the cumulative hospitalizations for the old group (70 years and older) and in (e-f) overall deceased cases.

The estimated values for the scaling factors used to differentiate the infectivity within the population are slightly smaller than the values obtained with the raw data. With for the young population, and for the old population, the disease induced death rate for hospitalized young and old groups, δ, are also estimated. Referring to the normalized data, the model parameters used for fitting the data are shown in Table 4. This parameter set will be used in the further sections of this manuscript.

The RMSE values were calculated to be 0.35 and 0.2 for the deterministic and stochastic models respectively. With lower values than the values obtained by fitting the raw data, again, the stochastic model has a better fitting (with a lower RMSE value than the deterministic model), confirming that the stochastic approach explains better the existing normalized data.

4 Results

4.1 Sensitivity analysis

A detailed sensitivity analysis is performed to determine how the parameter values variation will affect the reproduction number (R0) of the system. These results are of use to guide public health authorities during a disease outbreak.

In order to detect which are the parameters with higher impact on the R0 measure, with effects to increase or to decrease its value and consequently to define which parameters are to be targeted by intervention measures, we use the the normalized forward sensitivity method index of a variable to a parameter [48]. The normalized forward sensitivity index of R0 is defined using partial derivatives, showing the variation of the variable with respect to a given parameter p, as follows While the magnitude of the R0 measure increases as the values of β, a, ϕ, and ϵ parameters increase (positive indices), an inverse relation with the R0 value is observed for the δ1, α1, α2, α3, α4, and δ2 parameters, with negative indices, i.e., as the parameter values increase, the magnitude of R0 decreases. The sensitivity index of R0 for the parameter β is 1, meaning that R0 increases or decreases with the same percentage as the parameter β varies, see Fig 7.

Fig 7

Normalized forward sensitivity indices of R0.

Complementary to the forward sensitivity method index analysis above, we use the spline regression method to fit 10000 points for a range of each parameter value. The quantification of the parameter variation effect on the R0 value is shown in Fig 8, confirming that the increase of the transmission rate β, the fraction of asymptomatic individuals a, and the scaling factors differentiating the disease transmission withing the population, ϕ and ϵ, values affects significantly the behaviour of the R0 measure.

Fig 8

Spline regression method to quantify the effect of the model parameter variation on R0 behaviour.

4.2 Model simulations: An exploratory analysis

In this section, we explore different parameter combinations for the disease infectivity factors ϕ and ϵ and for the disease induced mortality rate δ that are able to explain the exponential phase of the COVID-19 epidemic in the Basque Country. For both, the deterministic and stochastic models, the assumed biological parameters for COVID-19 dynamics were estimated for the normalized data, see Table 4. While for the deterministic model simulations we have used the function ode45 in MATLAB, for the stochastic model simulations we have obtained 100 realizations using the Euler-Maruyama approach.

As an exploratory exercise to understand the impact of the key parameters on disease severity dynamics during the initial phase of the pandemics, numerical simulations are performed to describe the available data in the Basque Country, from February 15 to March 25, 2020, see Fig 10a). This is a dynamic work. While the present analysis focus on the introductory phase of the pathogen in the Basque Country, the evaluation of the effect of the imposed control measures will be carried out later.

Epidemiological data used in this study are provided by the Basque Health Department and the Basque Health Service (Osakidetza), continually collected with specific inclusion and exclusion criteria. We use the following incidence and cumulative data for RT-PCR (reverse transcriptase-polymerase chain reaction), see Fig 9. While the incidence data are shown in Fig 9a), the cumulative data used here refer to the overall hospital admissions, including the ICU cases, are shown in orange and the decease cases in black in Fig 9b).

Fig 9

COVID-19 epidemiological data in the Basque Country.

In (a) the cumulative hospital admissions and deceased cases. In (b) incidences for disease cases referring to hospitalizations including ICU and deaths.

Within the timeline of first wave of the pandemic in the Basque Country, the black line shows the date of the partial lockdown implementation, followed by the full lockdown, see red line. The light blue line shows the last data point used in this study, March 25, 2020, ten days after the partial lockdown was implemented, when the exponential growth of disease cases decelerates into a growth close to zero towards a linear phase [19].

To investigate the possible dynamics of hospitalizations for the young (H1) and for the old (H2) groups, as well as the dynamics for the overall deceased cases when no control measure would have been implemented in the Basque Country, a 100 days simulation time is shown, from February 15 to May 25, 2020, covering the post-lockdown period. The effects of different parameter combinations of the scaling factors of disease transmission and the disease induced mortality rates are shown in Fig 10. For the hospital admission cases dynamics, we evaluate the effect of the scaling factors affecting the disease transmission individually. By fixing ϕ = 1.5 as estimated from the normalized data, we vary the value of the ϵ parameter, see Fig 10a) and 10b), while in Fig 10c and 10d) we fixed ϵ = 1.3, varying the value of the ϕ parameter. The same experiment was performed for the deceased cases, see Fig 10e), varying the combination of the disease induced mortality δ, always assuming δ1 < δ2.

Fig 10

Deterministic model simulations.

Without any control measure, the epidemic would follow its course with a massive number of hospitalizations and deaths within the first 100 days of the pandemic. While a qualitatively similar dynamical behavior is observed when varying those key parameters, with an increase on the number of disease cases as the parameter value increases, the scaling factor ϕ, differentiating the transmission between the mild/asymptomatic and the severe/hospitalized individuals, appears to affect significantly the older population, eventually reaching its maximum towards stationary, much faster than the dynamics in the young population.

This effect is also observed for the overall infection cases (A1 + H1 + A2 + H2), and for overall hospitalizations (H1 + H2), see Fig 11. Using both modeling approaches, deterministic and stochastic, our results have shown that disease cases would have eventually reached stationarity after 100 days if no control measure was implemented.

Fig 11

By varying the infectivity scaling factors ϕ and ϵ, the dynamics of the overall disease cases (A1 + H1 + A2 + H2), and the dynamics of the overall hospitalization (H1 + H2) are plotted for 100 and 300 days respectively, using the following parameter set: β = 0.15, δ1 = 0.003, δ2 = 0.04, η1 = 0.035, η2 = 0.03, α1 = 0.02, α2 = 0.01, α3 = 0.05, α4 = 0.03 and a = 0.02.

In (a) and (c) the deterministic model simulations and in (b) and (d) 100 stochastic realizations.

5 Discussion

Declared a pandemic by the World Health Organization (WHO) in March 2020 [2], the collective behavior of societies has been significantly affected by the extreme measures implemented to control disease transmission. As the COVID-19 pandemic progressed, research on mathematical modeling became imperative and very influential to understand the epidemiological dynamics of disease spreading and control under different scenarios. The hypothesis of a new pathogen able to cause a very severe disease with an extremely high transmission rate were gradually adjusted overtime. It is now accepted that COVID-19 disease severity and death occur according to a hierarchy of risks, with age and pre-existing health conditions enhancing risks of disease severity.

In this paper, a mathematical model framework for COVID-19 transmission is proposed. Applied to the first wave of COVID-19 epidemic in Basque country, Spain, we stratify the population into young and old groups, after a detailed data analysis for the available epidemiological data referring to confirmed positive cases, hospitalization and deceased cases. The deterministic and the stochastic models are analyzed and results are compared.

For the deterministic approach, we calculate the disease-free equilibrium and the basic reproduction number (R0). We show that disease-free equilibrium is global asymptotically stable. A detailed sensitivity analysis is performed to identify the key parameters influencing the basic reproduction number, and hence, regulating the transmission dynamics of COVID-19.

Further, the deterministic model was extended to its stochastic counterpart. The stochastic differential equation (SDE) model is derived from a diffusion process. Simulations were obtained by the Euler-Maruyama method. Model derivation is shown in the S1 File.

Both models were able to fit well the empirical data using similar parameter value range, with the stochastic model always presenting a better result. A detailed sensitivity analysis was performed allowing us to identify the key parameters affecting the disease dynamics.

An exploratory analysis to understand the impact of those key parameters on disease severity dynamics during the initial phase of the pandemics, from February 15 to March 25, 2020, was performed. Numerical simulations have demonstrated that differences in infectivity from severe/hospitalized cases and mild/asymptomatic cases are the most important factors influencing the disease spreading in the population and without any control measure, the epidemic would have followed its course with a massive number of hospitalizations and deaths within the first 100 days of the pandemic.

These results are of use to guide public health authorities on disease control. The sensitivity analysis results shown in Figs 7 and 8 give insights on how to control the disease outbreak, suggesting possible ways of action for an effective containment of the disease transmission towards its elimination by limiting the increase of parameters with positive indices. On the other hand, by increasing the parameters with positive indices, such as providing treatment for a fast recovery or decreasing mortality, for example.

The numerical simulations have shown that without the lockdown, disease cases would increase continuously with severe cases eventually reaching its maximum numbers towards the herd immunity scenario, i.e, when a large portion of the population become immune to the disease. This behaviour was observed to affect the old group population much faster than the young population. Minimizing the scaling transmission factors, ϕ and ϵ, via social distancing or vaccination, for example, would significantly reduce the disease burden in the population. Therefore, in terms of policy implications, our findings support the vaccination strategy prioritising the most vulnerable individuals to reduce hospitalization and deaths, as well as the non-pharmaceutical intervention measures, e.g social distancing and use of masks, that are still advised by the public health authorities, to reduce disease transmission.

This is a dynamic work. While the present analysis has focused on the initial phase of the COVID-19 epidemic in the Basque Country, it is important to mention that the evaluation of the effect of the imposed lockdown and other control measures is ongoing. As continuation of this work, the models are under refinement, using this framework as baseline to describe the progression of COVID-19 epidemics in the Basque Country and to understand the impact of lockdown implementation and the increased of testing capacity over time. As our model is able to describe the available data, see Fig 12, we will be also able to measure the impact of mild/asymptomatic cases on disease spreading and control, including non-sterilizing vaccine performance [21].

Fig 12

The following parameter set: ϕ = 1.4, ϵ = 0.25, δ1 = 0.003, δ2 = 0.04, η1 = 0.035, η2 = 0.03, α1 = 0.02, α2 = 0.01, α3 = 0.05, α4 = 0.03 and a = 0.02, the deterministic dynamics for the overall hospitalizations (H1 + H2) is shown with and without control.

Cumulative data on overall hospitalizations are shown in blue. The simulation plotted as red line includes a control function (β(t) = β0 σ−(x(t)) + β1 σ+(x(t)), with a standard sigmoid function , see [17]) which is able to describe the empirical data, while the green line shows the solution without any control. The black line shows the last data point used in this study, March 25, 2020, ten days after the partial lockdown was implemented. By that date, the exponential growth of disease cases decelerates into a growth close to zero towards a linear phase. The full lockdown started on March 31, 2020 (black dashed line).

Supporting information includes computation of the basic reproduction number (R0) (S1 Appendix), proof of theorem 1.2 (S2 Appendix), stochastic process of the SIRS model (S3 Appendix), stochastic process of the proposed deterministic model 1 (S4 Appendix).

(PDF)

Click here for additional data file.

8 Jun 2022

13 Jun 2022

A pdf file with the detailed answer to the reviewers is uploaded.

Submitted filename: reviewers_June13_Aguiar.pdf

Click here for additional data file.

17 Jun 2022

Modeling the initial phase of COVID-19 epidemic: The role of age and disease severity in the Basque Country, Spain

PONE-D-22-11062R1

Dear Dr. Aguiar,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Constantinos Siettos, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

22 Jun 2022

PONE-D-22-11062R1

Modeling the initial phase of COVID-19 epidemic: The role of age and disease severity in the Basque Country, Spain

Dear Dr. Aguiar:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Professor Constantinos Siettos

Academic Editor

PLOS ONE