1 Introduction

Understanding how infectious diseases spread through populations requires an accurate characterization of social contact patterns, which form the foundation of transmission in compartmental epidemic models [34]. The social contact hypothesis [19] states that the number of secondary infections generated by an individual is proportional to their social contacts, with a pathogen-specific proportionality factor. Age-structured contact matrices, estimated by numerous studies [2, 14, 16, 20–23, 25–27], are widely used to represent the average number of interactions between individuals in different age groups, reflecting the heterogeneity of real-world social mixing [14, 16]. These matrices are essential for capturing how interventions, susceptibility, and behavioral factors influence disease dynamics. However, contact patterns are inherently uncertain and vary across time, geographic regions, and population demographics. During the COVID-19 pandemic, for instance, behavioral changes and non-pharmaceutical interventions (NPIs) such as mobility restrictions, lockdowns, and school closures caused substantial and often uneven shifts in contact rates [6, 15]. Furthermore, many contact matrices are derived from limited survey samples or extrapolated data [10], introducing additional uncertainty. These challenges can significantly impact model predictions and, consequently, public health decision-making.

Given these limitations, it is crucial not only to incorporate contact data into epidemic models but also to systematically assess how sensitive key epidemiological outcomes, such as the basic reproduction number ($\mathcal{R}_0$), hospitalizations, and mortality, are to variations in the underlying contact structure. Sensitivity analysis provides a structured way to do this. By identifying the contact matrix elements that most strongly influence model outcomes, such analysis enables researchers to prioritize data collection efforts and refine public health interventions [34].

Traditional global sensitivity techniques, such as Monte Carlo simulations, Latin Hypercube Sampling (LHS), and Partial Rank Correlation Coefficients (PRCC) [30–34], can capture uncertainty across parameter spaces but are computationally expensive and less informative about matrix-structured inputs like contact patterns. In contrast, analytical methods based on the Next Generation Matrix (NGM) provide more targeted insight into how changes in contact rates propagate through a model. Eigenvalue-based sensitivity analysis, in particular, allows one to assess the influence of contact perturbations on $\mathcal{R}_0$ by differentiating the NGM with respect to the contact matrix. For instance, Angeli et al. [4] applied formal perturbation techniques to the NGM to analyze the early spread of COVID-19 in Belgium. These approaches are especially useful in models where parameter interactions are more linear and tractable analytically, though they may be less effective in highly nonlinear or stochastic frameworks.

In this study, we extend the analytical approach proposed in [4] by developing a general sensitivity analysis framework centered on age-structured contact matrices. Unlike prior work, we construct a computational graph with the independent entries of the contact matrix as inputs, allowing us to explicitly trace how these elements influence key outcomes. Using the dominant eigenvalue of the NGM as a proxy for transmission potential, we can compute the gradient of $\mathcal{R}_0$ with respect to each contact element. We extend this framework beyond transmission analysis by linking the NGM to downstream clinical outcomes such as hospitalization, ICU admission, and mortality. Additionally, we demonstrate the flexibility of our framework by applying it to other epidemic models and contact data from the literature, showcasing its utility in various settings.

2 Methods

2.1 Framework for Sensitivity Analysis

The sensitivity analysis framework is designed to flexibly assess the impact of changes in contact patterns on disease transmission across various structured epidemic models. By updating key input files, users can adapt the framework to accommodate different social contact matrices, varying levels of model complexity, and customized model parameter values.

The framework’s primary input is a set of full and aggregated contact matrices that represent age-specific social interactions across various settings (e.g., homes, schools, workplaces, and other environments). These matrices are first scaled by the total population to produce a normalized total-contacts matrix, as defined in Eq. (9). Next, reciprocity is restored by averaging the total contacts between each pair of age groups (Eq. (7)), and the result is divided by the size of the contacting group to yield the full adjusted contact matrix, as described in Eq. (8).

The scaled contact matrix $C^{\mathrm{pop}}$, before these transformations, forms the input to a computational graph. This graph propagates the contact data through an age-structured transmission model to compute the next-generation matrix $\mathbf{K}$, from which the basic reproduction number $\mathcal{R}_0$ is obtained as its dominant eigenvalue. At this stage, automatic differentiation, implemented in PyTorch, is incorporated to compute the gradient of $\mathcal{R}_0$ with respect to the entries of $C^{\mathrm{pop}}$. These gradients form the foundation for all subsequent sensitivity calculations. We use automatic differentiation (AD) not to replace the closed-form eigenvalue sensitivity, but to propagate derivatives through the entire pipeline. For high-dimensional, multi-compartment models (e.g., the Röst model), manual chain-rule derivations are impractical; AD applies the exact chain rule on the computational graph, yielding tractable, reliable gradients, reducing algebraic errors, and making the framework readily extensible across models and outcomes.

The framework flexibly supports both simple and complex age-structured epidemic models for computing $\mathbf{K}$, making it adaptable to a variety of epidemiological settings. Mortality transition probabilities $\pi^{(\text{d})}_i$, defined in Eq. (4), can then be applied to transform $\mathbf{K}$, enabling sensitivity analyses focused specifically on the mortality outcome.

In the Eigenvector Calculator (see Fig. 1), the framework calculates the eigenvector associated with $\mathcal{R}_0$ according to Eq. (11). The NGM Gradient Calculator then uses Eq. (1) to generate $n_p$ sensitivity values in the Sensitivity Calculator, each capturing the impact of variations in specific elements of the contact matrix on $\mathcal{R}_0$. The sensitivity values are then summed using Eq. (2) to capture the total influence of each age group on $\mathcal{R}_0$. The results are presented as heat maps and charts, highlighting the interactions in the contact matrix that most significantly affect disease dynamics. The entire analysis process is automated in Python, using PyTorch for automatic differentiation and matplotlib for visualization.

2.2 Social Contact Matrices

To demonstrate the flexibility of our framework across different age-structured models, we use social contact matrices from two distinct sources representing two geographic contexts: the United Kingdom and Hungary. For the UK, we adopt contact data from Mossong et al. [14], which provides contact matrices for eight European countries based on 7,290 diary-based surveys. These diaries captured the characteristics of 97,904 contacts recorded over a single day, including details such as age, sex, location, duration, frequency, and whether physical contact occurred. For Hungary, we use synthetic contact matrices estimated by Prem et al. [16], who used Markov Chain Monte Carlo simulations and additional data sources to produce synthetic contact matrices for various countries based on contact data from Mossong et al. [14]. These matrices are stratified by setting—home, school, workplace, and other (e.g., public transit, leisure venues)—and reported in five-year age bands.

To ensure comparability across data sources, we apply a standard procedure to obtain an adjusted full contact matrix for each setting and country (details in Appendix A). The resulting matrices for Hungary are shown in Fig. 2, and the UK counterparts appear in the left panel of Fig. 3.

2.3 Sensitivity Measures

A precise characterization of disease transmission is crucial for developing reliable epidemic models. In age-structured models, transmission dynamics are governed by age-dependent parameters and a social contact matrix. We denote this contact matrix by $\bar{C} \in\mathbb{R}^{n \timesn}$, where $n$ represents the number of age groups—either $n = n_a$ for the full-resolution matrix or $n = n_g$ for its aggregated counterpart. The full-resolution matrix corresponds to the original age groups as provided in the empirical contact data sources (e.g., Prem et al. [16] or Mossong et al. [14]), while the aggregated version combines these groups into broader categories for comparison or to reduce model complexity. The entries of $\bar{C}$ are positive real numbers, since they represent the frequency of interactions between individuals across age groups. We compute the basic reproduction number $\mathcal{R}_0$ using the Next Generation Matrix $\mathbf{K}$ or the modified $\mathbf{K}$ based on the outcome as detailed in Section 2.4. The value of $\mathcal{R}_0$, defined as the dominant eigenvalue of $\mathbf{K}$, represents the expected number of secondary infections generated by a single infectious individual in a fully susceptible population.

To analyze how changes in contact rates influence $\mathcal{R}_0$, we first identify the elements of the contact matrix $\bar{C}$ that should be varied. Because $\bar{C}$ satisfies population-level reciprocity, each off-diagonal pair $(i,j)$ is determined by a single entry and the population sizes; hence, only the upper (or equivalently, lower) triangular part contains independent values. The total number of these unique contact values is then given by $n_p = \frac{n(n+1)}{2}.$ We collect these $n_p$ independent entries into a vector $\mathbf{c} \in\mathbb{R}^{n_p}$. Each element $c_p$ of $\mathbf{c}$ directly corresponds to one of the independent contact rates from the upper triangular part of $\bar{C}$, and thus represents a positive real value. The gradient of $\mathcal{R}_0$ is then computed with respect to the elements in $\mathbf{c}$.

Using Eq. (13) from Appendix B, each entry $S_{ij}$ quantifies the proportional contribution of each element $k_{ij}$ of $\mathbf{K}$ to $\mathcal{R}_0$ based on the left and right eigenvectors of $\mathbf{K}$. Following the low-level-parameters approach of [4], we treat contact entries as lower-level parameters and therefore compute the sensitivity of $\mathcal{R}_0$ with respect to $\mathbf{c}$ by applying the chain rule:

where

and

• $\frac{\partialk_{ij}}{\partialc_p}$ is a scalar that describes how the entry $k_{ij}$ of $\mathbf{K}$ changes with respect to the independent contact matrix elements.

The resulting gradient $\frac{\partial\mathcal{R}_0}{\partial\mathbf{c}} \in\mathbb{R}^{n_p}$ quantifies the sensitivity of $\mathcal{R}_0$ to changes in the $n_p$ independent elements of the contact matrix $\bar{C}$, representing the effect of perturbations in contact rates on disease transmission.

To summarize the overall impact of contact variations, we reconstruct the sensitivity information from Eq. (1) in matrix form, denoted by $\Sigma= [s_{ij}] \in\mathbb{R}^{n \timesn}$. Each entry $s_{ij}$ equals $\partial\mathcal{R}_0/\partialc_p$, where $p$ is the parameter index assigned to the upper-triangular contact pair $(i,j)$; the lower triangle is filled by mirroring by setting $s_{ij}=s_{ji}$. Thus, $s_{ij}$ quantifies the contribution of contacts between a susceptible individual in group $i$ and an infectious individual in group $j$. Diagonal entries capture within-group effects, while off-diagonal entries reflect cross-group contributions. Since each column $j$ of $\mathbf{K}$ represents the distribution of secondary infections produced by an infectious individual in group $j$, the cumulative transmission impact of group $j$ is naturally obtained by summing over its column of sensitivities:

2.4 Age-structured Epidemic Models

To illustrate our sensitivity analysis framework, we apply it to two age-structured epidemic models of varying complexity. Rather than constructing new models, we use two published compartmental models that incorporate age-specific contact matrices and transmission mechanisms: one by Pitman et al. (2012) [28], and the other by Röst et al. [17]. The objective of these case studies is to demonstrate that our eigenvector-based sensitivity framework can be applied consistently across models of different complexity, ranging from a relatively simple SEIRS-type system to a detailed multi-compartment COVID-19 model, and to highlight how demographic structures, contact matrices, and model assumptions influence the resulting sensitivity patterns. The Pitman et al. model serves as a classical influenza example for illustrating the methodology, whereas the Röst et al. model provides the opportunity to explore mortality-related outcomes in the context of COVID-19.

In both models, we compute the basic reproduction number $\mathcal{R}_0$ using the Next Generation Matrix (NGM) framework introduced by Diekmann et al. [8], where $\mathcal{R}_0$ is defined as the spectral radius of the NGM, $\mathbf{K}$. This approach involves linearizing the infection-related compartments of the system around the disease-free equilibrium. The resulting linear system is decomposed into two matrices: the transmission matrix $F$, which describes the rate at which new infections are produced, and the transition matrix $V$, which governs the movement between infectious compartments. We construct these matrices for each model in Sections 2.4.1 and 2.4.2.

2.4.1 Model 1: Pitman et al.

The first model is an SEIRS (Susceptible–Exposed–Infectious–Recovered–Susceptible) framework developed to simulate the transmission of influenza, incorporating age-specific contact patterns in England and Wales [28]. The model assumes uniform susceptibility across all age groups and sets the baseline reproduction number $\overline{\mathcal{R}}_0$ to 1.8. The system of governing equations is provided in Appendix D.1.

The infected state vector $X(t)$ comprises $n_s$ compartments, $E_i$ and $I_i$, for each age group $i = 1, \dots, n$. New infections enter only the $E_i$ compartments, while individuals in $I_i$ contribute to transmission but do not directly receive new inflow. The transmission matrix $F \in\mathbb{R}^{(n \cdotn_s) \times(n \cdotn_s)}$ is sparse, as new infections occur only in the $E_i$ compartments through contact with infectious individuals $I_j$. It is constructed as a block matrix with sub-blocks $F_{j,i}$, each representing transmission from infectious individuals in age group $j$ to susceptibles in age group $i$:

The progression matrix $V$ is block-diagonal, with each diagonal block $V_{i,i} \in\mathbb{R}^{n_s \timesn_s}$ defined by:

where

is the rate of progression from exposed to infectious, and

is the recovery rate. The negative diagonal entries capture the outflow from the exposed compartment (

) and from the infectious compartment (

).

To reduce the full $(n \cdotn_s)\times(n \cdotn_s)$ system to an $(n\timesn)$ NGM, we define the selection matrix $\mathcal{E} \in\mathbb{R}^{n \times(n \cdotn_s)}$. The $i$-th row of $\mathcal{E}$ selects the newly infected state $E_i$ from the local pair $[E_i,\, I_i]$ (i.e., it has the block $[1\ 0]$ at the columns of age group $i$) and zeros elsewhere; $I_i$ is not selected. The NGM is then computed as:

2.4.2 Model 2: Röst et al.

The Röst et al. model (section D.2) includes $n = 16$ age groups, each with $n_s = 9$ infectious compartments: $L^{(1)}, L^{(2)}, I^{(p)}, I^{(a,1)}, I^{(a,2)}, I^{(a,3)}, I^{(s,1)}, I^{(s,2)}, I^{(s,3)}$. The infected state vector $X(t)$ therefore contains all of these compartments for each age group. Consequently, the full transmission matrix has dimension $F \in\mathbb{R}^{(n \cdotn_s) \times(n \cdotn_s)}$. Each entry in this matrix reflects contributions governed by the contact rate $\overline{c}_{j,i}$, the age-specific susceptibility $\sigma_i$, and the transmission rate $\beta$, as described in the original formulation [17]. The matrix is composed of $n \timesn$ blocks $F_{j,i}$, where each block has shape $n_s \timesn_s$. Importantly, only the first row of each block is non-zero, since all new infections enter the first latent compartment $L^{(1)}_i$. The nonzero entries capture contributions from pre-symptomatic ($I^{(p)}_j$), asymptomatic ($I^{(a,m)}_j$), and symptomatic ($I^{(s,m)}_j$) individuals with reduced transmissibility for asymptomatic cases ($\inf^{(a)} = 0.5$).

The progression matrix $V \in\mathbb{R}^{(n \cdotn_s) \times(n \cdotn_s)}$ is block-diagonal, with each diagonal block $V_{i,i} \in\mathbb{R}^{n_s \timesn_s}$ encoding transitions within the infectious compartments of age group $i$, including progression from latent to pre-symptomatic, branching into asymptomatic or symptomatic paths, and recovery. To reduce the full system to an $(n \timesn)$ NGM, we define a selection matrix $\mathcal{E} \in\mathbb{R}^{n \times(n \cdotn_s)}$. The $i$-th row of $\mathcal{E}$ selects the newly infected compartment $L^{(1)}_i$ from the set of infectious states for age group $i$. The reduced NGM is then computed as in Eq. (3).

As seen in Eq. (15), this model is sufficiently complex to include various compartments for hospitalization and death. While it supports analysis of health outcomes such as hospitalization ($I^{(h)}$), ICU admission ($I^{(c)}$), and death ($D$), our focus here is solely on mortality. Accordingly, we consider only the probability of death when modifying $\mathbf{K}$.

The age-specific probability of death is derived from the model’s compartmental parameters and reflects the conditional pathway from presymptomatic infection to a fatal outcome:

where

is the probability that an individual in age group

progresses from the presymptomatic compartment

to death.

To explore mortality-related sensitivities in the epidemic process, we define a mortality-weighted matrix $\mathbf{M} = [m_{ij}]$ by scaling the entries of $\mathbf{K}$ with the age-specific probabilities of death:

Unlike the transmission analysis in which

is a well-defined epidemiological target, the mortality-weighted analysis does not introduce a new explicit target quantity; rather, it yields a conceptual target that captures the effectiveness of lethality (i.e., transmission weighted by fatality risk). Importantly, this measure should not be interpreted as the sensitivity of the cumulative mortality trajectory, but as a heuristic indicator of which contact pathways most strongly influence mortality-weighted transmission. This construction comes with important limitations. By scaling the NGM in this way, deaths are effectively assumed to occur instantaneously after infection, overlooking temporal delays between infection, hospitalization, and fatal outcomes. Therefore,

should not be regarded as a new NGM, and its spectral radius is not a mortality-specific reproduction number. Instead, the role of

is heuristic: it reweights transmission pathways according to their mortality consequences, and the associated left and right eigenvectors highlight which age-specific contacts contribute most to fatal outcomes. In this way, the sensitivity framework of Eq. (1) is adapted to emphasize mortality-related effects rather than transmission itself.

The choice between infection-based and mortality-weighted sensitivities ultimately depends on the epidemiological context and intervention goals: reducing transmission is critical when epidemic control and eradication are priorities, while focusing on mortality-weighted measures is more relevant in settings where severe outcomes dominate the public health burden. This trade-off is analogous to vaccination strategies, where prioritizing younger high-contact groups curbs spread, while targeting older groups reduces deaths.

3 Results

In this section, we present the results of our framework for quantifying the sensitivity of transmission dynamics to age-specific contact patterns. The contact matrices serve as key inputs to the models and are the primary focus of the sensitivity analysis. As a case study, we apply the method to the influenza SEIR model of Pitman et al. (Section D.1) using age-structured contact data from the UK. To validate the framework in a different epidemiological context, we repeat the analysis using Röst et al. model described in Section D.2, using the Hungarian contact data shown in Fig. 2. In addition to the full-resolution analysis, we also investigate the effect of contact matrix aggregation on the sensitivity results using the Röst model. Model parameter values and assumptions are taken directly from the original studies [17, 28].

3.1 Age-structured Influenza SEIR model by Pitman et al

This SEIR model provides a straightforward example to demonstrate the framework, using age-specific mixing patterns from England and Wales (Fig. 3, left panel). The model quantifies the sensitivity of the basic reproduction number $\mathcal{R}_0$ to contact rates, although it does not allow exploration beyond transmission-related outputs.

The middle panel of Fig. 3 shows the sensitivity of $\mathcal{R}_0$ to pairwise age group interactions. Notably, contacts among individuals aged 5–19 have the highest influence, suggesting that transmission dynamics are particularly sensitive to interactions in this age range. The right panel presents age-group-level sensitivity scores summed across contact pairs. It reveals that both younger individuals (ages 5–19) and middle-aged adults (ages 35–49) make the largest contributions to variations in $\mathcal{R}_0$.

3.2 Sensitivity Analysis Using the Röst et al. COVID-19 Model

To assess the applicability of our framework to more complex epidemiological settings, we apply it to the age-structured COVID-19 model developed by Röst et al. [17], which was calibrated for the Hungarian population. This model features 15 compartments, incorporating Erlang-distributed latent and infectious periods, distinct hospitalization stages, and fatality tracking. It uses synthetic age-specific contact matrices for Hungary, derived from Prem et al. [16] and illustrated in Fig. 2. Model parameters, including age-dependent recovery, transition rates, and susceptibility, are taken from the original study. The basic reproduction number is fixed at $\mathcal{R}_0 = 1.8$.

We compute the sensitivity of two key epidemiological outcomes, $\mathcal{R}_0$ and cumulative mortality, with respect to perturbations in contact rates. Here, cumulative mortality refers to the total number of deaths predicted by the model over the full course of the epidemic, aggregated across all age groups. Results are summarized in Fig. 4. The left panel shows that contacts among individuals aged 20–54 exert the strongest influence on $\mathcal{R}_0$, highlighting their central role in transmission dynamics. However, mortality (right panel) sensitivities shift toward contacts involving older adults, particularly those aged 60 and above.

To better understand which age groups most influence model outcomes, we compute age-group-level sensitivities by summing the pairwise sensitivity values associated with each age group (as defined in Eq. (2)). These results are shown in Fig. 5. Age-group sensitivity for $\mathcal{R}_0$ peaks in the 20–54 age range, indicating that improved estimation of contact patterns and transmission parameters in this demographic could substantially enhance model fidelity. Conversely, the sensitivities for mortality are concentrated among individuals aged 70 and above, underscoring the critical importance of accurate clinical parameterization and data collection for these older age groups. These insights can help prioritize efforts to refine age-specific contact matrices and outcome probabilities in the most impactful segments of the population, thereby improving the overall reliability and policy relevance of age-structured epidemic models.

3.2.1 Sensitivity Analysis Using the Model with Aggregated Age Groups

To assess the influence of age aggregation on sensitivity results, we applied the framework to a contact matrix aggregated from 16 to 7 age groups, as described in Eq. (10). The resulting sensitivities for transmission and clinical outcomes are presented in Fig. 6. Compared to the full-resolution analysis, the sensitivity for $\mathcal{R}_0$ now shifts slightly toward a broader age group, 15–59 years, indicating that aggregation blends the influence of younger and middle-aged groups. For the clinical outcomes of mortality, the highest sensitivities remain concentrated in the oldest age categories. However, with aggregation, the peak shifts from 70+ in the full model (Fig. 5) to 60+ in the coarser model.

4 Discussion

This study introduces a generalizable sensitivity analysis framework that integrates age-structured social contact data into compartmental epidemic models to assess the impact of contact heterogeneity on disease outcomes. The core objective is to quantify how variations in contact patterns between age groups influence key epidemiological quantities, such as the basic reproduction number $\mathcal{R}_0$, hospitalizations, ICU admissions, and mortality.

A distinguishing feature of our approach is its ability to decompose sensitivities across the full contact matrix and sum them to identify age groups with the most influential contact patterns. Building on ideas from [4], we incorporated cumulative sensitivity metrics to highlight which age groups contribute most substantially to model outputs, guiding targeted data collection and intervention design. To address the known issue that arbitrary age bin definitions can distort sensitivity analyses, we further extended the framework to support analyses over both fine- and coarse-grained age structures.

A significant implication of this framework is its ability to inform data collection. Because contact matrices are difficult to estimate accurately, especially in underrepresented age groups, our method helps prioritize empirical efforts by pinpointing which contact interactions carry the greatest influence or uncertainty. This focus enables more efficient and impactful refinement of contact data used in infectious disease modeling. At the same time, two important interpretive limitations should be noted. First, contact rates between age groups are not truly independent quantities. In practice, a change in contacts between groups $i$ and $j$ will typically influence other entries of the contact matrix, reflecting constraints from household composition, social networks, and cultural factors. Our framework, like most perturbation analyses, treats contacts as independent for tractability. The resulting sensitivities should therefore be understood as hypothetical responses to isolated perturbations rather than literal predictions of feasible behavioral change. To mitigate this issue, we emphasize cumulative sensitivities at the age-group level, which provide more robust insights than pair-wise sensitivities.

Second, the sensitivity patterns depend on the choice of age partition. Our results with aggregated contact matrices demonstrate that grouping can shift the apparent drivers of transmission and mortality (e.g., transmission sensitivities shifting from 20–54 to 15–59, mortality from 70+ to 60+). This dependence is inherent to age-structured modeling and highlights the importance of carefully selecting age bins. In practice, partitions should balance epidemiological relevance (e.g., school-age, working-age, elderly), policy needs, and data availability. Sensitivity results should therefore be interpreted in light of the chosen age structure, and complementary analyses using alternative partitions can help assess robustness.

The framework offers several key advantages. First, its deterministic nature provides reproducibility and computational efficiency, making it well-suited for rapid policy analysis. Second, the method is outcome-flexible, enabling sensitivity assessments not only for transmission metrics like $\mathcal{R}_0$, but also for downstream epidemiological outcomes such as hospitalization, ICU occupancy, or mortality, within the same model infrastructure. This versatility supports localized, age-specific public health strategies across diverse model types and settings. Despite these strengths, the framework also has other limitations. It is based on point estimates and does not currently incorporate statistical uncertainty, which constrains its interpretability in probabilistic terms. The underlying sensitivity analysis assumes local linearity around a baseline scenario and may not fully capture dynamics in highly nonlinear or unstable regimes. Another limitation is the assumption of static contact matrices, which restricts the framework’s ability to account for temporal variation in contact behavior due to evolving policies, interventions, or social norms. Nonetheless, the framework is readily extendable to accommodate time-varying contact matrices. Provided that contact data are available across distinct time intervals, the sensitivity analysis can be applied independently to each temporal slice over specified periods as in [4]. Future work will also aim to embed uncertainty quantification techniques, such as bootstrapping or Bayesian inference, to produce confidence intervals for sensitivity estimates.

On behalf of the RRF-2.3.1-21-2022-00006 project, we are grateful for the possibility of using HUN-REN Cloud (see [35], https://science-cloud.hu/) which helped us achieve the results published in this paper.

Declarations

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

• Funding

  This project was supported by the National Laboratory for Health Security Program (RRF-2.3.1-21-2022-00006) and the NKFIH KKP 129877, both funded by the Ministry of Innovation and Technology of Hungary through the National Research, Development, and Innovation Fund. The views and opinions expressed in this work are solely those of the author(s) and do not necessarily reflect those of the National Laboratory for Health Security Program or the NKFIH KKP 129877 project. Neither entity assumes responsibility for these views.

• Conflict of interest/Competing interests

  The authors declare that they have no competing interests.

• Ethics approval and consent to participate

  Not applicable.

• Consent for publication

  Not applicable.

• Data and Code Availability

  All the code and data required to reproduce our results are available in the repository [29].

• Author contribution

  Evans Kiptoo Korir: Contributed to methodology, conceptualization, drafting the original manuscript, reviewing and editing, software development, and visualization.

  Zsolt Vizi: Responsible for data curation, conceptual framework, formal analysis, methodology, project management, software development, visualization, manuscript review and editing, and supervision.

A Social Contact Matrices

We use Hungarian contact data estimated in [16] to describe the contact patterns applied in this study. The contact matrices for different settings—home ($\widetilde{C}^{\text{H}}$), school ($\widetilde{C}^{\text{S}}$), work ($\widetilde{C}^{\text{W}}$), and other environments ($\widetilde{C}^{\text{O}}$)—are combined to form the full contact matrix:

This full contact matrix

consists of elements

, representing the average number of daily contacts reported by individuals in age group

with those in age group

. Due to sampling variability, the matrix is often asymmetric and does not satisfy reciprocity. In theory, contact matrices should satisfy the reciprocity condition

where

and

are the population sizes of age groups

and

. To enforce reciprocity, we apply the adjustment method of Knipl et al. [11]. First, we construct the total contact matrix

This matrix is generally not symmetric. We then enforce symmetry at the level of total contacts by averaging:

By construction,

is symmetric, and when rescaled to per-capita form (Eq. (8)), the resulting matrix satisfies reciprocity. To obtain the per-capita contact matrix used in transmission modeling, we divide by the size of the contacting age group:

The resulting

specifies, for each individual in age group

, the average number of contacts with individuals in age group

. For sensitivity analysis, we also introduce a population-scaled version

which is not required in the classical definition of contact matrices, but is used in our computation graph. This normalized form ensures that contact entries are scaled consistently with the total population, serving as the foundation for gradient-based sensitivity calculations.

The adjusted contact matrices $\overline{C} \in\mathbb{R}^{n_a \timesn_a}$, computed via Eq. (8), are visualized in Fig. 2. Because of the reciprocity condition, the computational graph processes only the upper (or equivalently, lower) triangular entries of $C^{\mathrm{pop}}$, which amounts to $n_p$ unique values. These values sufficiently describe all pairwise contact interactions and serve as inputs for sensitivity analysis.

To examine how the level of age resolution influences model sensitivity, we also construct aggregated versions of the contact matrix from the full $\widetilde{C} \in\mathbb{R}^{n_a \timesn_a}$, following the aggregation approach described in [13]. Let $G_i$ denote the set of original age indices grouped into the $i$-th new age bin, with $i = 1, \dots, n_g$ and $n_g < n_a$. The population in the new age bin is defined as $\mathbf{N}_i = \sum_{m \inG_i} N_m$. The aggregated contact matrix $\mathbf{\widetilde{C}} \in\mathbb{R}^{n_g \timesn_g}$ is then computed by taking a population-weighted average of contact rates between all combinations of original age groups:

This aggregated matrix is subsequently processed using the same reciprocity-restoring and normalization procedures as the full contact matrix, yielding the final per-capita aggregated matrix

used in sensitivity analysis.

B General Sensitivity Measures

For a given matrix $\mathbf{K}$, the eigenvalue $\lambda$ is associated with the left and right eigenvectors, which provide crucial insight into the relationship between $\mathbf{K}$ and $\lambda$. These eigenvectors satisfy the following equations:

Since

is a strictly positive matrix, its dominant eigenvalue

is real, positive, and algebraically simple. It admits a corresponding right eigenvector

and left eigenvector

, both with strictly positive real components, satisfying:

Thus, the leading eigenvalue

governs the asymptotic growth rate of the system. Following the methodology outlined in [4, 36], we introduce a perturbation to

in Eq. (11):

Expanding both sides and omitting second-order terms:

Subtracting

from both sides yields:

Rearranging terms:

To isolate

, we project both sides onto the left eigenvector

:

Using the identity

, we find that:

which implies the last term vanishes. Thus:

Now, consider a perturbation that affects only one entry of

denoted by

, then:

Taking the limit as

yields

the sensitivity of

to the NGM entry

. Following [4], we obtain the corresponding proportional sensitivity,

In our framework, the entries of

depend on contact data and epidemiological rates, so we treat contact entries as low-level parameters and obtain derivatives via the chain rule (Eq.( 1)). Consequently, contact-level derivatives do not admit a closed form like (12) and are computed through the chain rule.

C Automatic Differentiation to Calculate Gradients

In our framework, a computational graph is the sequence of mathematical operations that maps the independent entries of the population-scaled contact matrix $C^{\mathrm{pop}}$ to the epidemiological outcomes of interest, such as the basic reproduction number $\mathcal{R}_0$. Each node of the graph represents an algebraic operation (e.g., reciprocity restoration, scaling, matrix multiplication, or eigenvalue computation), and the directed edges represent the flow of intermediate results.

The inputs to this graph are the $n_a \timesn_a$ elements of the fully adjusted and population-scaled contact matrix, $C^{\mathrm{pop}} \in\mathbb{R}^{n_a \timesn_a}$, as defined in Eq. (9). This matrix, which encodes the total contact rates between age groups normalized by the overall population, forms the foundation of the computational graph. It serves as the starting point for automatic differentiation, enabling gradient-based sensitivity analysis of transmission dynamics and outcome measures.

In the first phase, we treat only the upper-triangular (including diagonal) entries of $C^{\mathrm{pop}}$ as independent parameters and complete the matrix by reflecting them to the lower triangle. The resulting symmetry arises from the reciprocity-restoring averaging at the total-contacts level (Eq. (7)), not from assuming per-capita symmetry in the raw data. Subsequently, this matrix is converted into the adjusted contact matrix $\overline{C}$, as defined in Eq. (8). The matrix $\overline{C} \in\mathbb{R}^{n_a \timesn_a}$ serves as the contact input to the transmission model.

The computation graph then proceeds by feeding $\overline{C}$ into an age-structured compartmental model to compute the next-generation matrix $\mathbf{K}$. Depending on the outcome of interest, $\mathbf{K}$ may be further transformed (e.g., by incorporating hospitalization, ICU, or mortality probabilities), and its dominant eigenvalue is extracted to compute $\mathcal{R}_0$. Automatic differentiation is applied at this stage using the PyTorch library. We perform backpropagation through the entire graph to compute the gradient of $\mathcal{R}_0$ with respect to the independent entries of $C^{\mathrm{pop}}$. These gradients quantify the sensitivity of epidemic outcomes to variations in specific age-specific contact rates, forming the basis of cumulative sensitivity analyses.

D Epidemic models

D.1 Pitman et al. model

The SEIR model in Eq. (14) describes the progression of individuals through susceptible ($S$), exposed ($E$), infectious ($I$), and recovered ($R$) states, with key parameters: rate of latent infection $\alpha$, and recovery rate $\gamma$ as adapted from [28]. The corresponding system of ordinary differential equations (ODEs) is given by:

The model incorporates age-specific contact patterns through the matrix

, where contacts are based on pre-pandemic estimates. The parameter sets used for both the full-resolution and aggregated models are summarized in Table 1.

D.2 Röst et al. model

In this COVID-19 model, the population is divided into multiple compartments, each representing a distinct stage of infection. These compartments include susceptible individuals ($S$), who are at risk of contracting the virus, and latent individuals ($L$), who have been infected but are not yet symptomatic. As the infection progresses, individuals transition into the presymptomatic stage ($I^{(p)}$) before the onset of noticeable symptoms. Some individuals may remain asymptomatic ($I^{(a)}$), carrying the virus without exhibiting symptoms. Others will develop symptoms of varying severity, classifying them as symptomatic individuals ($I^{(s)}$). Those who experience severe symptoms and require hospitalization are categorized as hospitalized individuals ($I^{(h)}$). Among these, the most critically ill patients who require intensive care are placed in the intensive care unit (ICU) compartment ($I^{(c)}$). Finally, individuals who recover from the infection and gain immunity are classified as recovered ($R$), while those who succumb to the disease are placed in the deceased ($D$) compartment. The compartments and their transitions are depicted in Fig. 7, and the complete set of model equations is provided in (15).

The following set of ordinary differential equations governs the dynamics of the system:

In this model, individuals are divided into age groups indexed by

, where

denotes the probability of transmission upon contact. Each parameter indexed by

corresponds to a specific age group. The model leverages the contact matrix shown in Fig. 2 to compute sensitivity values relevant to transmission dynamics. A complete list of model parameters used for both the full-resolution and aggregated versions of the model is provided in the accompanying Zenodo repository [29].

References