Mechanistic Understanding of Pandemic Dynamics: A Multiscale Algorithmic Framework

1. Introduction

Pandemics have historically posed critical challenges to global public health, prompting extensive efforts across the mathematical and computational sciences to model, understand, and predict their progression. The COVID-19 pandemic, which began in late 2019, intensified this research, with early studies focusing on characterizing the initial outbreak dynamics [1,2,3], performing comparative analyses across countries [4,5,6], and examining spatial transmission patterns within national and regional settings [7,8]. As the pandemic progressed, attention focused toward quantifying the effectiveness of non-pharmaceutical interventions [9,10,11], mitigating healthcare system overload [12,13,14], and designing optimal exit strategies [15,16,17]. Throughout all stages, predictive modeling of epidemic trajectories and real-time estimation of the reproduction number remained central objectives [18,19,20].

This sustained research effort has produced a wide array of mathematical modeling approaches. Ordinary differential equation (ODE) compartmental models remain foundational [21,22,23,24], complemented by agent-based models [6,25,26] and meta-population formulations [27,28,29]. In parallel, data-driven approaches including machine learning [30,31], deep learning [32,33,34], Bayesian inference, and autoregressive time-series techniques [35,36,37] have been increasingly adopted. Systematic reviews of early COVID-19 forecasting efforts highlight both the breadth of methodologies employed and the substantial uncertainty associated with early predictions [38,39,40].

Table 1. Translational Glossary: Mapping Computational Singular Perturbation (CSP) terminology to biological and epidemiological interpretation.

Mathematical Term (CSP)	Biological / Epidemiological Translation	Practical Significance in This Study
Mode ($a_n{f}^n$)	Dynamical component	A component of the model, the action of which is characterized by an explosive or dissipative time scale. The mode that is characterized by the fastest time scale drives the initial surge, while all other modes sustain or attenuate the epidemic wave before stabilization occurs.
Explosive Time Scale (${\tau }_e$)	Characteristic expansion time	Represents the time frame of the action of a mode (explosive) that tends to drive the system away from equilibrium. A smaller value (e.g., ∼4 days during the 6th wave) indicates a more aggressive and rapidly spreading variant compared to larger values (e.g., ∼8 days during the 4th wave or ∼28 days during the 5th wave).
Dissipative Time Scale (${\tau }_d$)	Stabilization or damping rate	Represents the time frame of the action of a mode (dissipative) that tends to drive the system towards equilibrium. When a dissipative mode dominates, the outbreak transitions into decay and active cases begin to decline.
Slow Invariant Manifold (SIM)	Epidemic trajectory	The reduced-dimensional path followed by the epidemic once fast transient processes subside. It represents the established dynamical regime during the outbreak phase.
Amplitude (f)	Driver intensity	Provides a measure of the impact of a mode in driving the epidemic wave. In the case of an explosive mode, a high amplitude indicates dominant transmission dynamics, while a low amplitude signals increasing control or immunity effects.
Time Scale Participation Index (TPI)	Mechanism identifier	Quantifies the contribution of individual biological transitions (e.g., transmission, incubation, recovery) to the time scale. In the case of an explosive time scale, it determines the degree to which the outbreak is promoted by high transmission ($R^1$) or by rapid progression from exposure to infection ($R^2$), or it is obstructed by recovery ($R^4$).
Pointer (Po)	Population influence index	Identifies which population compartment most strongly influenced by a specific mode. For example, a high Pointer value for the exposed population during the 6th wave indicates that rapid progression from exposure to infection drives the surge.
Inflection Point	Turning point in epidemic acceleration	The moment when outbreak growth shifts from acceleration to deceleration. This precedes the peak of infected population and serves as an early warning indicator of the impending peak.

Table 1. Translational Glossary: Mapping Computational Singular Perturbation (CSP) terminology to biological and epidemiological interpretation.

Mathematical Term (CSP)	Biological / Epidemiological Translation	Practical Significance in This Study
Mode ($a_n{f}^n$)	Dynamical component	A component of the model, the action of which is characterized by an explosive or dissipative time scale. The mode that is characterized by the fastest time scale drives the initial surge, while all other modes sustain or attenuate the epidemic wave before stabilization occurs.
Explosive Time Scale (${\tau }_e$)	Characteristic expansion time	Represents the time frame of the action of a mode (explosive) that tends to drive the system away from equilibrium. A smaller value (e.g., ∼4 days during the 6th wave) indicates a more aggressive and rapidly spreading variant compared to larger values (e.g., ∼8 days during the 4th wave or ∼28 days during the 5th wave).
Dissipative Time Scale (${\tau }_d$)	Stabilization or damping rate	Represents the time frame of the action of a mode (dissipative) that tends to drive the system towards equilibrium. When a dissipative mode dominates, the outbreak transitions into decay and active cases begin to decline.
Slow Invariant Manifold (SIM)	Epidemic trajectory	The reduced-dimensional path followed by the epidemic once fast transient processes subside. It represents the established dynamical regime during the outbreak phase.
Amplitude (f)	Driver intensity	Provides a measure of the impact of a mode in driving the epidemic wave. In the case of an explosive mode, a high amplitude indicates dominant transmission dynamics, while a low amplitude signals increasing control or immunity effects.
Time Scale Participation Index (TPI)	Mechanism identifier	Quantifies the contribution of individual biological transitions (e.g., transmission, incubation, recovery) to the time scale. In the case of an explosive time scale, it determines the degree to which the outbreak is promoted by high transmission ($R^1$) or by rapid progression from exposure to infection ($R^2$), or it is obstructed by recovery ($R^4$).
Pointer (Po)	Population influence index	Identifies which population compartment most strongly influenced by a specific mode. For example, a high Pointer value for the exposed population during the 6th wave indicates that rapid progression from exposure to infection drives the surge.
Inflection Point	Turning point in epidemic acceleration	The moment when outbreak growth shifts from acceleration to deceleration. This precedes the peak of infected population and serves as an early warning indicator of the impending peak.

Among these approaches, compartmental epidemic models remain the most widely used due to their analytical tractability, interpretability, and ability to incorporate mechanistic understanding of disease transmission. The classical SIR model [41,42] and its extensions can capture key epidemiological characteristics [43,44], including presymptomatic transmission, asymptomatic infections, and disease-induced mortality [45,46]. For COVID-19, SEIRD-type models were shown particularly appropriate due to the presence of an incubation period and non-negligible fatality rates [47,48,49,50]. The SEIRD model has been extended to incorporate vaccination dynamics and intervention measures in order to improve epidemic prediction under evolving public-health conditions [51]. Comparative studies demonstrated that increasing model complexity does not necessarily lead to improved early predictive performance, especially under limited data availability [52].

Early forecasting based on compartmental models demonstrated both their utility and their limitations. While short-term predictions captured initial exponential growth, longer-term forecasts were often highly sensitive to parameter uncertainty and structural assumptions [53]. Parameter identifiability issues are particularly severe during the early outbreak phase, where sparse, delayed, and heterogeneous data lead to multiple parameter sets that fit observations equally well but diverge significantly in their predictions [54,55,56]. These challenges have been documented extensively for COVID-19 [38,57,58]. Hybrid approaches combining mechanistic models with data-driven methods have been proposed to mitigate these issues, though they remained sensitive to reporting delays, intervention changes, and calibration strategies [32,44,59,60,61].

A particularly influential line of work has focused on identifying early dynamical signatures of epidemic transitions. Change-point detection and Bayesian inference techniques have revealed how intervention measures alter transmission dynamics during the early exponential phase [40]. Early warning systems were designed to predict epidemic surges or peaks before they show up in official case data [62,63]. Nevertheless, many of these methods still rely heavily on sufficient observational data and may struggle precisely when early prediction is most critical.

Motivated by the limitations of existing early prediction approaches, and particularly their reliance on extensive parameter calibration and large, high-quality datasets, this work proposes an alternative framework for early pandemic outbreak prediction grounded in time scale analysis of epidemic dynamics. Rather than focusing on solution profiles, which are often highly sensitive to parameter non-identifiability, the proposed approach exploits the intrinsic multiscale structure of compartmental epidemic models. Epidemic systems naturally exhibit processes evolving over multiple time scales, encompassing interactions between transmission, latency, recovery, and population turnover mechanisms [64,65]. By identifying and analyzing the time scales characterizing these interactions, it is possible to infer critical features of outbreak evolution, such as the driving paths, the critical population groups or regime shifts, using limited early-phase data. This time scale–based perspective provides a complementary and data-efficient pathway for early outbreak prediction, strengthening the theoretical foundations of epidemic forecasting and aligning with broader calls for system-oriented modeling frameworks capable of integrating complex epidemic dynamics [66,67].

Herein, the use of Computational Singular Perturbation (CSP) [68,69] is proposed as the core methodological tool for implementing this time scale–based framework in compartmental epidemic models. CSP is an algorithmic method of Geometrical Singular Perturbation Analysis (GSPA), designed for the analysis of multiscale dynamical systems that are driven by processes evolving over widely separated fast and slow time scales. In such systems, processes associated with fast time scales are rapidly equilibrated, leading to the formation of low-dimensional structures in the system’s tangent space, known as Slow Invariant Manifolds (SIMs) [70]. System trajectories subsequently evolve along these manifolds, governed by the remaining slow processes. CSP provides a complete and systematic set of diagnostic tools for identifying (i) the variables and processes associated with fast time scales that generate the SIM and (ii) the slow processes that dominate the long-term system dynamics [71].

Originally developed for the analysis of chemical kinetic mechanisms [72,73], CSP has since emerged as a general framework for multiscale dynamical systems, with applications spanning systems biology [74,75], pharmacokinetics [76,77], brain lactate metabolism [78,79], glycolysis [80] and population dynamics [81].

To demonstrate the systems-level understanding that can be acquired with the proposed CSP–based time scale analysis framework, we apply it to the case of the COVID-19 epidemic in Greece. The analysis will focus on the dynamics of the outbreak of the epidemic wave. Given prior knowledge regarding the existence of a latency period following infection with the SARS-CoV-2 virus [47,48], a SEIRD model is employed, incorporating exposed and deceased population groups. Specifically, we examine the 4^th, 5^th, and 6^th epidemic waves in Greece, which occurred in July 2021, October 2021, and December 2021, respectively. Each wave was driven by different underlying dynamics, reflecting changes in transmission intensity, population behavior, and intervention policies. Using CSP diagnostic tools, we focus on the mode that drives the outbreak and we identify the dominant populations and transition processes governing the outbreak phase of each epidemic wave, based solely on early-stage dynamics and limited observational data. This case study illustrates the method’s ability to disentangle evolving epidemiological mechanisms and highlights its potential as a robust early-warning tool for pandemic outbreak dynamics.

The overall methodological framework adopted in this study is summarized in Figure 1. The flowchart illustrates the sequential procedure from early outbreak data acquisition to CSP-based dynamical analysis and inflection point prediction.

2. Materials and Methods

2.1. The Mathematical Model

According to the SEIRD model, the population of susceptible individuals ($SP$) might become exposed ones ($EP$) upon transmission of the virus; i.e., individuals that are infected but not yet infectious. After the incubation period, the exposed individuals become infected ($IP$), so they are able to transmit the virus. The infected individual can then either become deceased ($DP$) or recovered ($RP$). The paths that connect the five population groups, according to the SEIRD model, are depicted schematically in Figure 2.

block = [draw=blue, line width=1pt,fill=gray!3, rectangle, minimum height=2em, minimum width=3em] sum = [draw, fill=white!3, stroke=white, circle, node distance=1cm] input = [coordinate] output = [coordinate] pinstyle = [pin edge=to-,thin,blue]

The governing equations of the SEIRD model are:

$${\displaystyle \frac{d}{dt}}\left[\begin{array}{c}SP\\ EP\\ IP\\ RP\\ DP\end{array}\right]=\left[\begin{array}{c}-\beta IPSP\\ \beta IPSP-\sigma EP\\ \sigma EP-(\gamma +\mu )IP\\ \gamma IP\\ \mu IP\end{array}\right]$$

(1)

where $\beta$ is the transmission density-dependent ratio, $\sigma$ is the transition ratio from exposed to infected individuals, expressing the inverse of the incubation period of the disease, $\gamma$ is the recovery ratio, which also expresses the inverse of the infection period of the disease and $\mu$ is the fatality ratio. Note that $SP$, $EP$ and $IP$ are decoupled from $RP$ and $DP$. The SEIRD model in Eq. (1) can be cast in the form:

$${\displaystyle \frac{d\mathbf{y}}{dt}}={\mathbf{S}}_1{R}^1+{\mathbf{S}}_2{R}^2+{\mathbf{S}}_3{R}^3+{\mathbf{S}}_4{R}^4=\mathbf{g}\left(\mathbf{y}\right)$$

(2)

where $\mathbf{y}$ is the 5-dim. column state vector, the elements of which represent the number of individuals in the population groups; $\mathbf{y}={\left[\begin{array}{c}SP,EP,IP,RP,DP\end{array}\right]}^{\top }.$ The 5-dim. column vector field $\mathbf{g}\left(\mathbf{y}\right)$ incorporates the transition rates from a compartmental group to another. The column vectors ${\mathbf{S}}_k$ (k=1,4) denote the directions of the 5 transitions from a population group to another; ${\mathbf{S}}_1={\left[\begin{array}{c}-1,1,0,0,0\end{array}\right]}^{\top }$, ${\mathbf{S}}_2={\left[\begin{array}{c}0,-1,1,0,0\end{array}\right]}^{\top }$, ${\mathbf{S}}_3={\left[\begin{array}{c}0,0,-1,0,1\end{array}\right]}^{\top }$, ${\mathbf{S}}_4={\left[\begin{array}{c}0,0,-1,1,0\end{array}\right]}^{\top }$. The scalars $R^k$ (k=1,4) are the related to ${\mathbf{S}}_k$ rates of transition; $R^1=\beta IPSP$, $R^2=\sigma EP$, $R^3=\mu IP$ and $R^4=\gamma IP$, which denote the transmission, incubation, recovery and death rates, respectively. The model in Eq. (2) is subject to two conservation laws:

$$\begin{array}{c}\hfill {\displaystyle \frac{d(SP+EP+IP+RP+DP)}{dt}}=0{\displaystyle \frac{d(\gamma DP-\mu RP)}{dt}}=0\end{array}$$

(3)

2.2. Model Calibration

The model was calibrated against the daily cases of confirmed infected ($IP$), recovered ($RP$) and deceased ($DP$) individuals, as reported by the Greek government COVID-19 data repository [82]. The initial conditions of the susceptible ($SP$) and exposed ($EP$) individuals and the parameter values of $\beta$, $\sigma$, $\mu$ and $\gamma$ were estimated using the Genetic algorithm [83] provided by COPASI [84], on the basis of the reported data sets of $IP$, $RP$ and $DP$. Since this study focuses on the outbreak phase of three different epidemic waves, the model was calibrated for two periods of each wave, in order to demonstrate the robustness of the results. The starting point of both periods was the same (the start of the wave), while the end of these periods was different, but within the outbreak phase; see Section 3.1, Section 3.2 and Section 3.3 for details. The parameter values thus estimated for all six periods are displayed in Table 2.

2.3. Modeling Assumptions and Vaccination Effects

It is noted that during the three epidemic waves considered here, vaccination coverage constituted and important factor. The vaccination campaign expanded during the first half of 2021; starting from older populations late January of 2021 and reaching 30-59 years old in April of that year [85]. Therefore, vaccination had no important effect during the 4th wave (July 2021), but had a stronger influence during the 5th and 6th waves (October-November 2021 and December 2021 - January 2022, respectively); although, due to the prevailing Omicron new variant during the 6th wave [86,87], the vaccine might have been not as effective as it was with the Delta variant that prevailed during the 4th and 5th waves [88]. The SEIRD framework employed here does not explicitly include dedicated compartments for vaccinated individuals or waning immunity (i.e., reinfections). This modeling choice was made to preserve structural parsimony and to focus on the dynamical properties of the outbreak phase rather than on detailed population stratification.

Importantly, the CSP-based time scale analysis relies on the local dynamical structure of the system, as encoded in the Jacobian matrix. Consequently, the net effects of vaccination, partial immunity, and reinfection are implicitly reflected in the coefficients of the model. Changes in susceptibility, reduced infectiousness, or accelerated recovery due to vaccination or reinfection manifest as modifications in the coefficients of the transition pathways, which in turn alter the characteristic dynamics. Therefore, while vaccination and reinfections are not modeled explicitly as a separate compartments, their aggregate dynamical impact is captured in the model.

2.4. Time Scale Analysis and CSP Tools

In the Computational Singular Perturbation (CSP) context, Eq. (2) is cast in the form:

$$\begin{array}{c}\hfill {\displaystyle \frac{d\mathbf{y}}{dt}}=\sum _{n=1}^5{\mathbf{a}}_n{f}^n{f}^n=\sum _{m=1}^4({\mathbf{b}}^n·{\mathbf{S}}_m)R^m\end{array}$$

(4)

where the column CSP vectors ${\mathbf{a}}_n$ (n=1,...,5) span the 5-dim. tangent space, the row vectors ${\mathbf{b}}^n$ are their dual (${\mathbf{b}}^i·{\mathbf{a}}_j={\delta }_j^i$) and the amplitudes $f^n$ are the projections of the vector field $\mathbf{g}\left(\mathbf{y}\right)$ on ${\mathbf{a}}_n$ (by proper adjustment of the signs of the CSP basis vectors) [68,69]. Each mode ${\mathbf{a}}_n{f}^n$ associates with a time scale ${\tau }_n$, ordered so that ${\tau }_1<{\tau }_2<\dots$. The amplitude $f^n$ provides a measure of the impact of the mode, while the time scale ${\tau }_n$ provides a measure of the time frame of its action. The time scale ${\tau }_n$ can be explosive or dissipative, depending on whether the components of the vector field responsible for its generation tend to move the trajectory away or towards equilibrium [89,90]. Due to the conservation laws in Eq. (3), $f^4=f^5=0$ and ${\tau }_4={\tau }_5=\infty$. As a result, only the first three modes will be retained in Eq. (4) and the last two will be neglected. According to CSP, the basis vectors ${\mathbf{a}}_n$ are computed via a refinement procedure that provides, order by order, the higher order corrections in the context of asymptotic analysis [91,92].

It will be shown next that during the outbreak phase of all three pandemic waves considered, the fastest time scale ${\tau }_1$ is of dissipative character (henceforth, ${\tau }_1={\tau }_d$ and $f^1=f^d$), while the other two time scales are of explosive one (henceforth, ${\tau }_2={\tau }_{e,f}$, ${\tau }_3={\tau }_{e,s}$ and $f^2=f^{e,f}$, $f^3=f^{e,s}$); where ${\tau }_d<{\tau }_{e,f}<{\tau }_{e,s}$. Due to the dissipative character of the fastest time scale ${\tau }_d$, the magnitude of the related amplitude $f^d$ is the outcome of significant cancellations among the additive terms in the rhs of its expression:

$$\begin{array}{c}\hfill f^d=({\mathbf{b}}^d·{\mathbf{S}}_1)R^1+({\mathbf{b}}^d·{\mathbf{S}}_2)R^2+({\mathbf{b}}^d·{\mathbf{S}}_3)R^3+({\mathbf{b}}^d·{\mathbf{S}}_4)R^4\end{array}$$

(5)

i.e., $f^d/max\left|({\mathbf{b}}^d·{\mathbf{S}}_k)R^k\right|\ll 1$, where the magnitude of this ratio is proportional to the ratio ${\tau }_d/{\tau }_{e,f}$ [69,91]. As it will be shown next, when the pandemic waves will be analyzed, these cancellations tend to keep the amplitude $f^d$ small, so that the process is mainly driven by the dominant of the two slower explosive modes; see Appendix A for details. Interested in leading order accuracy, the CSP basis vectors can be represented by the right eigenvectors of the Jacobian $\mathbf{J}\left(\mathbf{y}\right)$ of $\mathbf{g}\left(\mathbf{y}\right)$; i.e., ${\mathbf{a}}_n={\alpha }_n$ and ${\mathbf{b}}^n={\beta }^n$, where ${\alpha }_n$ and ${\beta }^n$ are the right and left eigenvectors of $\mathbf{J}\left(\mathbf{y}\right)$, satisfying ${\beta }^i·{\alpha }_j={\delta }_j^i$ [68]. The CSP formulation in the form of Eq. (4) allows the development of various diagnostics tools, among them the Pointer (Po) and the Timescale Participation Index (TPI).

Each population group is associated to each CSP mode with varying intensity. The degree to which each group associates to the n-th CSP mode ($n=1,2,3$) is quantified by the Pointer (Po) metric:

$${\mathbf{D}}^n=diag\left[{\alpha }_n{\beta }^n\right]=\left[{\beta }_1^n{\alpha }_n^1,{\beta }_2^n{\alpha }_n^2,{\beta }_3^n{\alpha }_n^3,{\beta }_4^n{\alpha }_n^4,{\beta }_5^n{\alpha }_n^5\right]$$

(6)

where, due to the orthogonality of the CSP basis vectors, the following relation holds ${\sum }_{i=1}^5{\beta }_i^n{\alpha }_n^i=1$ [93]. The components of ${\mathbf{D}}^n$ assess the way the four populations will adjust when the system is perturbed along the direction of the basis vector ${\alpha }_n$ [94].

The time scales will be approximated by the inverse of the eigenvalues of the Jacobian $\mathbf{J}\left(\mathbf{y}\right)$; i.e., ${\tau }_n={\left|{\lambda }_n\right|}^{-1}$ ($n=1,2,3$) [89]. This formulation allows the introduction of the Timescale Participation Index (TPI) as:

$$J_k^n={\displaystyle \frac{c_k^n}{{\sum }_{i=1}^4\left|c_i^n\right|}}{c}_k^n={\beta }^n\nabla \left({\mathbf{S}}_k{R}^k\right){\alpha }_n{\lambda }_n=c_1^n+\dots +c_4^n$$

(7)

where $k=1,4$, $n=1,3$ and by definition ${\sum }_{k=1}^4\left|J_k^n\right|=1$ [95,96]. $c_k^n$ denotes the contribution of the k-th path to the n-th time scale and can be either positive or negative. $J_k^n$ identifies the paths in the SEIRD model that contribute significantly to the generation of the timescale ${\tau }_n$, while negative (positive) $J_k^n$ implies that the k-th path contributes to a dissipative (explosive) character of the n-th timescale ${\tau }_n$.

In the case of a complex pair of eigenvectors, where ${\alpha }_k={\alpha }_r+i{\alpha }_i$ and ${\alpha }_{k+1}={\alpha }_r-i{\alpha }_i$ (${\lambda }_k={\lambda }_r+i{\lambda }_i$ and ${\lambda }_{k+1}={\lambda }_r-i{\lambda }_i$), the related pair of basis vectors are defined as ${\mathbf{a}}_k={\alpha }_r$ and ${\mathbf{a}}_{k+1}={\alpha }_i$. In that case, the quantities ($f^k$, ${\mathbf{D}}^k$) will refer to the real part (${\alpha }_r$) and the quantities ($f^{k+1}$, ${\mathbf{D}}^{k+1}$) will refer to the imaginary part (${\alpha }_i$) [95].