Two-Age-Structured COVID-19 Epidemic Model: Estimation of Virulence Parameters to Interpret Effects of National and Regional Feedback Interventions and Vaccination

: The COVID-19 epidemic has recently led in Italy to the implementation of different external strategies in order to limit the spread of the disease in response to its transmission rate: strict national lockdown rules, followed ﬁrst by a weakening of the social distancing and contact reduction feedback interventions and ﬁnally the implementation of coordinated intermittent regional actions, up to the application, in this last context, of an age-stratiﬁed vaccine prioritization strategy. This paper originally aims at identifying, starting from the available age-structured real data at the national level during the speciﬁc aforementioned scenarios, external-scenario-dependent sets of virulence parameters for a two-age-structured COVID-19 epidemic compartmental model, in order to provide an interpretation of how each external scenario modiﬁes the age-dependent patterns of social contacts and the spread of COVID-19.


Introduction
COVID-19 (SARS-CoV-2) is at the root of the recent economic and public health crisis worldwide.Since it was reported in December 2019 in China, the virus quickly took pandemic proportions throughout six continents and over 210 countries.Over 100 countries declared lockdowns and curfews, with an estimated global economic loss of one trillion US dollars in 2020 (see [1] and references therein).By October 2020, over 36 million people were definitely reported to be infected with COVID-19 and more than one million people had died from virus-related complications.
COVID-19 causes respiratory disease: common symptoms include fever, dry cough, fatigue, shortness of breath and loss of smell or taste, with possible complications including pneumonia and acute respiratory distress syndrome up to severe respiratory failures, septic shocks and death [2].
A huge amount of effort has been spent with the aim of finding novel methods for mathematical modeling and control of epidemics [3].Mathematical models can in fact accurately portray the epidemic's dynamic spread [3][4][5][6][7][8][9].Eight stages of infection, namely susceptible (S), infected (I), diagnosed (D), ailing (A), recognized (R), threatened (T), healed (H) and extinct (E), are presented in [8] to illustrate how restrictive social-distancing measures have to include a combination with widespread testing and contact tracing.A control-oriented SIR model that stresses the effects of delays and compares the outcomes of different containment policies is proposed in [5], whereas stochastic transmission models are considered in [3,9].A method for the joint optimal lockdown and release design in a pandemic is proposed in [4] and then applied in a realistic simulation scenario based on the data of COVID-19's evolution in Italy.A dynamical model specifically designed for COVID-19 is used in [6] to describe the epidemic evolution in Italy, with different kinds of control actions (social, political, and medical) being explicitly modeled.A parameter-varying modification of the SIRD model is finally proposed in [7] for describing and predicting the behavior of the COVID-19 contagion in Italy through identification of model parameters, written as linear combinations of basis functions.
An important feature of COVID-19 is its highly non-uniform attack of different age strata of society [2]: the infection fatality ratio for individuals older than 80 is likely significantly higher than the infection fatality ratio for individuals younger than 50 years.In particular, [10], which fit an age-structured mathematical model to epidemic data from China, Italy, Japan, Singapore, Canada and South Korea, shows that susceptibility to infection in individuals under 20 years of age is approximately half that in adults aged over 20 years, and that clinical symptoms manifest in 21% of infections in 10-to 19-year-olds, rising to 69% of infections in people aged over 70 years.On the other hand, ref. [11] estimates an overall infection fatality rate of 1.29%, as well as large differences by age, with a low infection fatality rate of 0.05% for those under 60 years old and a substantially higher 4.25% for people above 60 years of age; even if only 10% of the population were infected, the infection fatality rate would not rise above 0.2% for people under 60.
Since social contacts are influenced by age structure of the population and the frequency of contacts across the population [12], the spread of the disease relies on contact patterns among different subjects in the infected population.Epidemic models that are stratified by age are therefore particularly relevant when the hospital load and fatalities related to COVID-19 are to be estimated, with age-dependent patterns of social contacts being incorporated.Such models, which take into account the mechanism of its transmission, including the (possibly heterogeneous) pattern of mixing among the population, the susceptibility within the population, the virulence of the infection, the probability of transmission per contact, and the changes in behavior in the affected population in response to an epidemic [13], constitute useful tools to understand and characterize the complex transmission dynamics acting among different groups of the population (see [13]), as well as to identify and predict the effects of different age-stratified intervention strategies in slowing the spread.They concurrently (i) provide valuable information for public-health policy makers, (ii) avoid health systems saturation, and (iii) mitigate the impact on costs.
In this respect, despite containing simplifying assumptions, common variants of SIR-type models, including age-dependent substructures, are of great help in characterizing epidemics (see [14] and references therein).The reader is referred to the very recent [1,12,[14][15][16][17][18] for age-structured modeling of the COVID-19 epidemic and related age-dependent analyses.The reader is also referred to [19] for the latest study adopting a model with an age-dependent pre-pandemic contact matrix, reflecting the goal of a return to pre-pandemic routines once a vaccine is available, to compare five age-stratified vaccine prioritization strategies.A modified age-structured SIR model-based on known patterns of social contact and distancing measures within Washington, USA-is presented in [20]: population age-distribution has a significant effect on disease spread and mortality rate and contributes to the efficacy of age-specific contact and treatment measures.On the other hand, ref. [21] shows that if transmission rates' return to normal in the future and the epidemic ends only when population immunity is sufficient to survive reintroduction of infection, then age-targeted mitigations can still achieve a large mortality reduction.
This paper, unlike any other approach in the literature, aims to illustrate how six different diseases transmission scenarios and concurrently adopted social distancing and feedback interventions-including an age-stratified vaccine prioritization strategy-modify the age-dependent patterns of social contacts and the spread of COVID-19 disease.To this purpose, we use a (deterministic) two-age-structured COVID-19 epidemic compartmental model, in which two-age-classes (lower than 60 years old and not lower than 60 years old) are adopted, in order to comply with the adopted vaccination strategy while avoiding issues related to the lack of identifiability.In innovation is that identification of virulence parameters within the two groups is performed during the different phases.Real data, indeed, are taken from the Italian context, in which the implementation of the following subsequent-in-time different strategies has been carried out over time in response to different disease transmission scenarios: (i) a strict national lockdown rule (scenario a), as necessary in the first place (in the presence of a relatively high estimate of the disease transmission rate) to remove social contacts in workplaces, schools, markets and other public areas; (ii) a weakened feedback social distancing and contact reduction intervention (in the presence of a relatively low estimate of the disease transmission rate), which is composed of a weakened lockdown phase (scenario b), a low distancing phase (scenario c), a low distancing + workplace/school-contacts re-activation phase (scenario d), with a progressive release of the population back to their daily routine appearing; (iii) a coordinated intermittent regional action (scenario e)-in the presence of a newly alarming increase in the estimated disease transmission rate-where social distancing measures are put in place or relaxed independently by each region according to the ratio between hospitalized individuals and the total capacity of the health system in that region; and (iv) direct mRNA-vaccination of subjects-especially the elderly-(scenario f ) at highest risk for severe outcomes, along with Vaxzevria-vaccination of young subjects belonging to crucial occupational categories, to indirectly protect subjects at highest risk for severe outcomes.
In particular, the recent [22] is at the foundation of the intermittent intervention, which is inspired by the regionalism as an integral part of the Italian constitution and involves regional feedback strategies, where each of the twenty regions strengthens or weakens local mitigating actions, namely social distancing, inflow-outflow control, as a function of the saturation of their hospital capacity.It is worth noticing that the regional action was actually intermittent in a strict sense, owing to the presence, over time, of Italian regions with actually weakened social distancing.
The resulting framework originally exploited in the paper thus enriches the one proposed in [19], while turning out to be useful to compare impacts of national or regional intervention contexts and age-stratified vaccine prioritization strategies on the virulence parameters (within the age-dependent groups) of the age-stratified model (though the effects of seasonality on COVID-19 remain unsettled, they might be represented by the estimation results.).The results of this paper thus move in the direction of understanding the impact of human contact networks and human behavior on the spread of infectious diseases (no prediction purpose is declared), while assessing the implications of this for the planning of public health policy.In fact, even though the simplest mathematical models assume that the population mixes homogeneously, such an assumption is often only sufficient to obtain general insights, with the pattern of contacts between different age groups playing an essential role in determining the spread of disease.Indeed, in a real epidemic, the behavioral changes will not only reduce the number of contacts and intensity but will even change the structure of the contact network.

The Model
The (deterministic) compartmental model used in this paper is reported in this section.It is a natural extension of the classical SIR model [23].The compartments are subdivided into two different age groups: group of subjects with age lower than 60; group of subjects with age not lower than 60.There are many reasons for the choice to split the population into these two groups.The first one, which is more conceptual, is that this division more or less coincides with a division of active and retired populations, in which one may expect different patterns of social interactions, leading to different transmission dynamics.The second one is in accordance with the recent literature, where empirical estimates based on population-level data show a sharp difference in fatality rates between young and old people and firmly rule out overall fatality ratios below 0.5% in populations with more than 30% being over 60 years old [11].The third reason relies on the fact that this choice is compatible with the vaccination strategy adopted in scenario f.Finally, a closer look at the Italian data reveals that this choice is the one that divides the number of COVID-19 cases in the most uniform way, allowing a better exploitation of the data information and reducing the bias due to non-uniform cardinality division between the two classes.The resulting model thus reads: where: t is the time, measured in days; S i , I i and C i , i = y, o are the numbers of susceptible, infected and reported cases for the two age classes, respectively; and N(t) is the number of persons who are not quarantined, hospitalized or dead at time t.The parameters v ij , i, j = 1, 2 represent the virulence of the virus among the different age classes, while 1/τ i , i = 1, 2 is the average time for disease identification and γ = 0.07 is the rate of asymptomatic infected who recover without being reported [22,24].Note that the model is not autonomous, since N(t) cannot be reconstructed from the state variables.In other words, such a model does not count the number of quarantined, hospitalized, recovered or dead, and thus it needs this time-series in order to be simulated.Even though this model might be certainly extended so it can used for prediction, this would be out of the scope of this paper.As aforementioned, the goal here is to understand the demographical habits (how people interconnected) starting from the epidemiological data while estimating model parameters to overview people's reactions under the six different scenarios.

Estimation of Model Parameters
To fit the model, we first detected the six scenarios of Section 1 that characterized the COVID-19 pandemic in Italy: a.
From t a 0 = 9 March 2020 to t a e = 28 April 2020: strict national lockdown rule in which social contacts in workplaces, schools, markets and other public areas are removed; b.
From t b 0 = 7 May 2020 to t b e = 3 June 2020: weakened feedback social distancing and contact reduction intervention, with a slow release of the population back to their daily routine appearing (especially the elderly, as a psychological toll due to the suffered isolation); c.
From t c 0 = 9 June 2020 to t c e = 8 September 2020: low feedback social distancing and contact reduction intervention, due to a low ratio between hospitalized individuals and the total capacity of the national health system; d.
From t d 0 = 15 September 2020 to t d e = 27 October 2020: low feedback social distancing and contact reduction intervention, with social contacts in workplaces and schools being re-activated; e.
From t e 0 = 7 November 2020 to t e e = 29 December 2020: coordinated intermittent regional action, where social contacts in schools is decreased at national level and social distancing measures are put in place or relaxed independently by each region according to the ratio between hospitalized individuals and the total capacity of the health system in that region; and f. from t Real data about the pandemic on each of these scenarios are taken from the official Ministerial website https://www.epicentro.iss.it/coronavirus/aggiornamenti(accessed on 1 June 2021).In particular, in the data centre, it is possible to find: After noting that we fit the model to the real data by minimizing the squared relative error between the measured data of the detected cases and the ones predicted from the model.More precisely, starting from the time windows that identify scenario a, i.e., t ∈ [t i 0 , t i e ], i = a, we compute the trajectory of the model that starts from the initial condition for a particular set of the parameters that we are identifying, which are We then compute as cost the relative error between the predicted and the real new daily cases, i.e., We repeat the same procedure for the following scenarios, i.e., i = b, . . ., f , but, this time, to guarantee the continuity of the identified solution, we impose that the initial condition of the current scenario (parameters I i t0y and I i t0o ) is different from the final condition of the previous one of at most 10%.We then tune the 48 parameters (10 of them are constrained) in order to minimize the sum of the six cost functions: through the fmincon routine in Matlab©.The total number of data we use for our fitting procedure over the considered six scenarios is 128 [22,10,26,14,18,38, respectively], with the problem of estimation needing at least 48 data points (one for each of the parameter we are estimating, eight per scenario).A practical identifiability analysis [25][26][27] of the parameters around the estimation point confirms that the values we obtained with this procedure can be locally determined from the data we used (the local minimum we have found has no directions on which the cost function does not significantly increase with respect to the parameter variations).This allows us to provide the picture of the agedependent patterns of social contacts and the spread of COVID-19 disease in the Italian context, which is depicted in Figure 1.The estimated parameters (in the different scenarios i ∈ {a, . . ., f }) v i kl , k, l ∈ {1, 2} are reported in Table 1, while the estimated parameters τ i l , l ∈ {1, 2}, and I i t0y , I i t0o appear in Table 2.The modulus and phase of the complex numbers are reported in Table 3 (in the different scenarios a -f ), with i here representing the imaginary unit satisfying i 2 = −1.Notice that a non-zero real number λ i 1 or λ i 2 thus represents the case in which only intra-group virulences appear.Conversely, a purely imaginary number λ i 1 or λ i 2 represents the case in which only cross-group virulences appear.On the other hand, the alternative option for coordinates in the complex plane constituted by the polar coordinate system-which uses the distance of a point from the origin and the angle subtended between the positive real axis and the line segment connecting the origin and the point in a counterclockwise way-allows us to directly visualize the complex number λ i 1 or λ i 2 through modulus and phase.Zero phase therefore means a real number, whereas π/2 phase means a purely imaginary number (according to the previously reported interpretation), with the modulus 21 constituting a measure of both the contributions coming from the involved virulence parameters.

Discussion
The following comments are in order.They provide a meaningful interpretation of the estimation results in accordance with Tables 1-3.

•
All the estimates corresponding to the different scenarios, including the estimated I i t0y , I i t0o (initial young subjects infected; initial old subjects infected), allow the estimated profile to satisfactorily reproduce the actual one along the different scenarios, as shown by Figure 1.

•
Comments for estimated 1/τ i 1 (average time for disease identification in young subjects).This average time takes homogeneous values: it varies from 3 to 7 days weeks over all the scenarios, with about 7 days passing for scenarios b and c.Actually, after the lockdown period and the related concerns, young subjects paid much less attention to their symptoms (recall that scenarios b and c cover a period starting from 7 May 2020 up to 8 September 2020).In addition, recall that young subjects have a higher probability of being asymptomatic (or even weakly symptomatic), while old subjects have a lower probability of being asymptomatic.Asymptomatic subjects usually continue their social interactions, infecting many people before recognizing that they are sick, and are then isolated.

•
Comments for estimated 1/τ i 2 (average time for disease identification in old subjects).This average time varies from 1 to 5 days, with less than 3 days occurring in scenarios cf in which the elderly paid a higher level of attention to symptoms, as a psychological toll due to the suffered isolation in scenarios a-b.

-
During scenario a, a very small intra-elder virulence appears due to the strict national lockdown rule, with an increase during scenario b, due to the weakened feedback social distancing and contact reduction intervention.

-
During scenarios c and d, a larger increase in the intra-elder virulence occurs, during summer holidays (as a consequence of the juvenile-elder virulence of scenario b) and owing to the re-activation of contacts in workplaces and schools.
Recall that school closures during epidemics and pandemics aim to decrease transmission among children.They seemingly have whole-population effects, whenever children are major contributors to community transmission rates.

-
During scenarios e and f, a decrease in the intra-elder virulence is exhibited (when compared to scenarios c and d), as a consequence of an imposed decrease in social contacts in schools and in the direct mRNA vaccination of subjects (the elderly) at highest risk for severe outcomes, in spite of a re-activation of social contacts in schools and in Christmas-related activities.Notice that the intermittent intervention of scenarios e and f, in which each of the twenty regions strengthens or weakens local mitigating actions as a function of the saturation of their hospital capacity, has been largely lighter than the lockdown intervention of scenario a, leading to the possibility of reinvigorating economy and mitigating costs due to the epidemic's spread.
-Large intra-juvenile virulence (>0.39) is exhibited in scenarios a, c, e and f, i.e., during the strict lockdown (with the virus circulating within families), as well as during summer holidays and after the first days of November, whereas small values accordingly appear in scenarios b and d, in which more attention was paid by young subjects after the perceived social alarms coming after the end of the strict lockdown and the end of summer vacations.-Large juvenile-elder virulence (>1.36) is exhibited in scenarios d-f, after 15 September 2020, owing to the (typically Italian) juvenile-elder contacts coming from school re-activation, with the smallest value actually occurring in scenario e, in which social contacts in schools are decreased at national level and social distancing measures are put in place or relaxed independently by each region according to the ratio between hospitalized individuals and the total capacity of the health system in that region.Nevertheless, a large phase of λ b 1 (with a rather small modulus of λ b 1 ) is exhibited in scenario b, owing to a weakened feedback social distancing and contact reduction intervention after the strict lockdown.

-
The elder-juvenile virulence appears to be relatively small (about zero) in all the scenarios, except for scenario a, in which the virus circulated within families (see also the phase of λ a 2 ).

-
The sum of the two λ 1 and λ 2 phases is small only in scenario c, i.e., during holiday vacations, in which a sort of decoupling between the two age classes appeared.
On the other hand, once the estimates of the model parameters have been obtained (Tables 1 and 2), the values of the reproduction number R i t [m] associated with model (1) in each scenario i [m stands for model-based computation], as average number of new infections caused by an infected person, can be computed through the formula (adapted from [24]) where σ 1 (•) denotes the biggest among the moduli of the eigenvalues of the matrix argument.The resulting mean reproduction numbers R i i [m] over the scenarios i ∈ {a, . . ., f } read: 1.2, 0.7, 0.7, 2.8, 0.9, 1.1.They are compatible -excepting for a relatively slight overestimate of scenarios f -with the maximum likelihood values of the national reproduction number in Figure 2, computed from raw data through the EpiEstim toolbox [28] (a gamma distribution is used as prior distribution, with shape α = 1.87,rate β = 0.28, whereas the R t -maximum likelihood value is taken as (α − 1)/β regarding the analytically expressed posterior distributions), showing that model (1) -though relying on just two age-based subgroups and neglecting, for instance, gender-based different behaviours -is able to catch the main epidemic features along the considered scenarios.
Finally, estimating τ 1 and τ 2 in the six scenarios leads to the possibility of identifying the number of infected cases, in accordance to coming from the last two equations in (1).This is an advantageous feature of our approach.Looking at Table 2: τ −1 1 (young age class) is close to 3 in scenarios a, e, f, whereas it is larger than 5 in the least juvenile-action-restrictive scenarios b-d; τ −1 2 (old age class) is close to 1.5 in the most-decoupled or vaccine-characterized scenarios c, e, f, whereas it is close to 5 in the scenarios a-b in the middle of the pandemic wave.

Conclusions
Starting from the available age-structured real data at national level (from 9 March 2020 up to 12 May 2021), the parameters of a two-age-structured COVID-19 epidemic compartmental model-with the same two-age-classes definition adopted in the implemented vaccination strategy and more or less coinciding with a division of active and retired popu-lations, as well as with a division characterized by relatively large differences in fatality rates and by uniformity in subgroups cardinality-have been identified in the different scenarios reported in Section 3: the ways in which external scenarios have modified the age-dependent patterns of social contacts and the spread of COVID-19 disease has been assessed.In particular, an epidemiological model for Covid-19 has been developed, which considers the epidemic within the younger age group and older age group separately.Such a model provides insight, at national level, in the different evolution of the epidemic within these two interacting age groups, while simultaneously evaluating, through the estimation of model parameters along time, the impact of changes in social distancing measures and vaccination due to varying external strategies.The results of the present study exhibit some specific limits at their root, such as the inclusion of just two age-dependent subgroups or the absence of gender differences in the subgroups, which might be further investigated in future studies (see [29] for related results).In this respect, it is worth noticing that such a modeling choice is, however, motivated by the fact that adopting more complex SIR variants may fall into unidentifiability problems owing to insufficient data in the details of the many involved compartments, or because of their overly complex structure whose different features cannot be caught in the initial fast-increasing phases.Nevertheless, our study also possesses points of strength.It certainly gives a deep interpretative insight into the time-varying action of parameters within a well-defined COVID epidemics model structure, while providing useful information regarding not only the number of undetected infected cases but also the effects of strategic actions and behaviours.This is rather meaningful, especially on the eve of possible occurrences of variant-based epidemic waves in response to an exact duration of the immunity for the vaccines that is, at this moment, uncertain.The exact structure of the contact patterns in the general population is, in fact, still unknown to a large extent and merits specific research efforts, especially when model-based prediction has to be performed in the presence of political choices that change rules governing social distancing.Such political choices, when they vary along time and affect interactions between subgroups, influence variations of the internal parameters of deterministic models that should be suitably taken into consideration in the related context.

f 0 = 5
January 2021 to t f e = 15 May 2021: direct mRNA-vaccination of subjects (the elderly) at highest risk for severe outcomes and indirect protection through Vaxzevria vaccination of young subjects belonging to crucial occupational categories.

Figure 1 .
Figure 1.Data fitting for the compartmental model: actual and estimated cumulative profiles for young subjects infected and old subjects infected (logarithmic scale).

•
The cumulative detected cases on a weekly scale C(t) divided by age (so it is possible to compute C y (t) and C o (t)); • The number of recovered people (not divided by age) R(t).

Table 3 .
National reproduction number R t within the considered time windows.Each shaded portion of the plane corresponds to a specific scenario.The mean value Ri t (among the values corresponding to our sampling) within each time window i ∈ {a, . . ., f } is reported at the bottom of each shaded region.Modulus and phase of the complex numbers λ i 1 = v i 11 + iv i 12 and λ 2 = v i 22 + iv i 21 in the different scenarios i ∈ {a, . . ., f }.