Characterizing the transmission patterns of seasonal influenza in Italy: lessons from the last decade

Background Despite thousands of influenza cases annually recorded by surveillance systems around the globe, estimating the transmission patterns of seasonal influenza is challenging. Methods We develop an age-structured mathematical model to influenza transmission to analyze ten consecutive seasons (from 2010 to 2011 to 2019–2020) of influenza epidemiological and virological data reported to the Italian surveillance system. Results We estimate that 18.4–29.3% of influenza infections are detected by the surveillance system. Influenza infection attack rate varied between 12.7 and 30.5% and is generally larger for seasons characterized by the circulation of A/H3N2 and/or B types/subtypes. Individuals aged 14 years or less are the most affected age-segment of the population, with A viruses especially affecting children aged 0–4 years. For all influenza types/subtypes, the mean effective reproduction number is estimated to be generally in the range 1.09–1.33 (9 out of 10 seasons) and never exceeding 1.41. The age-specific susceptibility to infection appears to be a type/subtype-specific feature. Conclusions The results presented in this study provide insights on type/subtype-specific transmission patterns of seasonal influenza that could be instrumental to fine-tune immunization strategies and non-pharmaceutical interventions aimed at limiting seasonal influenza spread and burden. Supplementary Information The online version contains supplementary material available at 10.1186/s12889-021-12426-9.

Initial condition on the fraction of susceptible individuals at the beginning of each season are set by considering the observed age-specific vaccination coverage (see Tab. S1) and strain-specific vaccine effectiveness (VE) (see Tab S2). Tab S1. Seasonal influenza vaccine coverage by age group (m=month, y=year) for influenza vaccine in the period 2010-2020, as reported to the Italian National Institute of Health (1).

Virological data
The share of samples collected among ILI cases testing positive to each strain and total number of analyzed samples are obtained from data collected by the National Institute of Health (6) -data reported in Tab S3.
Tab S3. Circulating strains and relative share of samples testing positive to each strain over the period 2010-2020. We highlighted in blue the strains that above the 6% threshold that we used to conduct the Bayesian analysis. The shares in 2010-2011 season were approximated using data obtained from virological surveillance conducted in Lombardy region (7). In this study, we consider the circulation of a strain to be negligible during a season if the share of samples testing positive to that strain below 6%. All Influenza strains in all seasons are included in the analysis of the infection attack rate, while only influenza strains that are above a 6% threshold are considered in the Bayesian analysis.
The weekly number of swabs collected by the regional Reference Laboratories are shown in Fig.S1. . Weekly number of swabs analyzed by the regional reference laboratories plotted alongside the weekly ILI incidence.
The shares of ILI samples by age group a testing positive for strain s in season y, ! ( , ), are also obtained from the virological surveillance database of the National Institute of Health (6) -data reported in Tab S4.

Estimation of ILI reporting rate
We consider, as a starting point, levels of immunity of the Italian population in different age classes before and after the A/H1N1 influenza pandemic in 2009-2010 (8,9). The two sero-epidemiological studies are based on 1,152 and 1,436 leftover sera collected among individuals of the general population before and right after the end of the 2009 influenza pandemic, respectively.
Using these two datasets we estimate the infection attack rates of the 2009 influenza pandemic for two age groups: 0-14y and 15y+ (denoted as "#$% and $&' , respectively), as follows: where: • a denotes the age class; • POS . /012 and POS . /34 are the shares of samples testing positive to A/H1N1 A/California/07/2009 (A/H1N1pdm09) in age group a before and after the 2009-10 pandemic, respectively.
The attack rate detected by the Italian surveillance system in the same age groups (denoted as "#$% 565 and $&' 565 for the 0-14y and 15+y age groups, respectively), can be computed as follows: ;<$ where: is the the share of samples collected among ILI cases testing positive to A/H1N1pdm09; • a denotes the age class; • w denotes the surveillance week; • . ( ) is the incidence of ILI cases reported to the Italian epidemiological influenza surveillance system (10) during the pandemic for age group a on week w.
Note that during the 2009-2010 influenza season, 96.4% of samples collected among ILI cases testing positive for influenza infection, tested positive to A/H1N1pdm09 infection (6).
The reporting rates of the surveillance system for the age group 0-14 years and 15+ years, (denoted as "#$% and $&' , respectively) can be estimated by the following equation: Essentially, "#$% and $&' represent the share of influenza infections that are detected by the Italian influenza surveillance system. In the rest of the analysis, we assume these reporting rates are the same for other types and subtypes of influenza virus and for each season. Although this is a strong assumption, we find similar estimates of the reporting rate for other surveillance systems estimated for different influenza seasons after the 2009 influenza pandemic (11,12).

Estimation of age-specific infection attack rates
For each analyzed influenza season (from 2010-11 to 2019-20), the age-specific infection attack rates are calculated by combining the reporting rates mentioned above with the observed ILI incidence in the different age classes (0-4, 5-14, 15-64, and 65+ years) and virological data by age. In particular, we use the following equation: where: • a is the age class; • s is the strain; • w is the surveillance week; • y is the season; • ! ( ; ) is the incidence of ILI cases reported to the Italian surveillance system (10) for age group a on week w of season y; • 1 . (y) is the share of ILI samples for age group a testing positive for strain s in season y; • ra is the reporting rate for age group a (see previous section).

Age-structured influenza transmission model
Influenza transmission for the three considered strains is simulated through a deterministic nonstationary age-structured SEIR model (9) stratified in 85-years age classes and based on the assumption of heterogeneous mixing by age. The epidemiological transitions for each individual's age are described by the following system of ordinary differential equations: where: • t represents the time; • S(a,t), E(a,t), I(a,t) and R(a,t) represent the number of susceptible, latent, infectious and removed individual of age a at time t, respectively;

• N(a,t)=S(a,t)E(a,t)+I(a,t)+R(a,t) represents the total population of age a at time t;
• is the transmission rate; • ! denotes the susceptibility to infection of individuals of age a relative to the susceptibility to infection of individuals aged 0-4 years that is set to 1; • C is contact matrix by age, where each element !;!A represents the average number of contact between individuals of age a and individuals of age a'; • 1/ is the average duration of the latent period; • 1/ is the average duration of the infectious period. We initialized the model with 10 infectious individuals (and zero latent individual) at the beginning of each season, distributed randomly across different ages. The initial numbers of susceptible and removed individuals account for the observed age-specific vaccination coverage and strain-specific vaccine effectiveness and are computed as follows: where: • (0) and (0) are respectively the initial numbers of susceptible and removed individuals in the population (irrespectively of the age) at the beginning of each season; • (0) is the total population at the beginning of each season; • ( , 0) is the coverage for individuals of age a at the beginning of each season (see Sec. S1); • ( ) represents the strain-specific vaccine effectiveness (see Sec. S1); • ( , ′) is the Dirac delta function, equal to 1 if a = a' and 0 otherwise.
The contact matrix for the Italian population is taken from (13). The average duration of the latent period (1/ ) is set to 1.5 days (14) and the infectious period (1/ ) is set to 1.2 days in such a way that the resulting generation (15) is 2.7 days, in agreement with influenza literature (16).

Bayesian analysis
We calibrate the influenza transmission model (see Sec. 5) separately for each season and circulating strain. For each season and strain, we have four free parameters (the transmission rate, , and the susceptibilities to infection $,...9 of individuals in age classes 5-14, 15-64 and >64 years, respectively) whose posterior distributions are estimated through a Bayesian approach. In particular, model calibration is carried out by means of a differential evolution Markov chain Monte Carlo sampling (17) applied to the binomial likelihood of the season-and strain-specific influenza infection attack rates * ! ( ). In particular, for each season y, the likelihood of the transmission rate and susceptibility to infection by age group given the infection attack rates by age is defined as: where: • a runs over the four age classes for which the infection attack rates is known (i.e., 0-4, 5-14, 15-64, and 65+ years); • n(a) is the number of individuals in the age class a tested for a specific strain during seasonspecific virological studies; • * ! ( ) is the age-specific infection attack rates for strain s during season y (see Sec. S4); • ( ; , ! ) is the estimated infection attack rate in age class a as obtained by model simulation with a specific transmission rate β and susceptibility to infection $,...9 .
For each influenza season and strain, we run five chains of 50,000 iterations using different starting points for β and $,...9 . For β and $,...9 we use the same prior distribution: a flat distribution between 0 and 1,000. MCMC convergence and the length of the burn-in period are assessed by checking via visual inspection of the trace plots associated with the different chains, i.e. the sequence of accepted parameter values are approximately characterized by a constant width and average, therefore proving good mixing of the parameters.
As an example, we show the trace plots associated with one chain for the 2010-11 season and A/H1N1pdm09 strain is shown in Fig. S3 (panels A-D), along with the estimated infection attack rate by age (panel E) and posterior distributions of Reff and susceptibility to infection by age group (panels F and G). A Trace plot of the susceptibility to infection of the 5-14y age group relative to the 0-4y age group. B As A, but for the 15-64y age group. C As A, but for the 65+y age group. D As A, but for the transmission rate. E Observed mean age-specific infection attack rates (see Sec.

Sensitivity analyses
We conducted two sensitivity analyses to account for possible variability in the generation time. We assumed respectively a generation time of 1.7 days and 3.7 days, chosen according to the work of Vink et al. (18). Based on these values, we therefore computed the incubation and infectious periods as the same fractions of the generation time used in the baseline analysis, respectively 55% and 45%. Tab.S5 and Fig.S4 show estimates on the effective reproduction number and the age-specific susceptibility to infection obtained considering a generation time of 1.7 days. Similarly, Tab.S6 and Fig. S5 show the same estimates obtained considering a generation time of 3.7 days.  S5. Estimated posterior distributions of the susceptibility to infection by age group (mean and 95%CI) relative to the 0-4 years age group (for which the susceptibility to infection is set to the reference value of 1). Only types/subtypes that accounted for more than 15% of the positive samples are considered. The values reported above the vertical lines in the right panels represent the 97.5% percentile of the distribution, when the value exceed the limit of the vertical axis.