Learning and Forecasting Multiscale Spatiotemporal Epidemic Dynamics with Lagrangian-Eulerian Graph Neural ODEs
Authors: Jinyi Duan, Qian Li, Biao Tang
Organizations: School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, 710049, People’s Republic of China . · School of Mathematics and Statistics, Xi’an University of Technology, Xi’an, 710054, People’s Republic of China . · Center for Intersection of Mathematics and Life Sciences, Xi’an Jiaotong University, Xi’an, 710049, People’s Republic of China .
Accurate medium- to long-range epidemic forecasting requires extrapolating high-dimensional, heterogeneous and nonstationary dynamics from partial observations. Yet hybrid models often leave implicit which observations inform which epidemic processes and how learned graph interactions modify system state. We introduce LE-EpiGNODE, which replaces generic message passing with process-specific, mechanism-guided neural operators: age-contact information drives Lagrangian infection generation, whereas mobility drives conservative Eulerian population transport within a structure-preserving universal differential equation. We evaluated LE-EpiGNODE across 14 controlled synthetic mechanisms spanning time-modulated, state-dependent, nonlinear and coupled dynamics, three complementary COVID-19 settings, 21 comparator models and five targeted ablations. Forecasting gains tended to become more pronounced at longer horizons; at the longest horizon, LE-EpiGNODE achieved the lowest WMAPE in all three real-world datasets, reducing error by 17.7-40.0% relative to the best baseline. Controlled ablations linked these gains to the Lagrangian-Eulerian decomposition and process-aligned priors, while synthetic ground truth showed that learned operators and induced state fluxes track known generating dynamics. These results support mechanism-organized graph neural dynamics as a route to reliable long-horizon extrapolation with directly testable, process-level mechanistic interpretation beyond trajectory accuracy.
Figures & tables
Figure 1: Conceptual motivation and framework of LE-EpiGNODE. (a) , Conceptual synthesis motivating LE-EpiGNODE. From the graph-learning side, flexible spatiotemporal message passing progresses toward continuous neural interactions but requires governing-dynamics guidance to define message semantics and admissible state updates. From the mechanistic side, compartmental models progress toward age- and spatially structured UDEs with learnable unresolved functions but require scalable graph neural learning to represent high-dimensional heterogeneous couplings. These reciprocal requirements converge on mechanism-guided graph neural operators embedded within structured governing dynamics. (b) , LE-EpiGNODE instantiates this convergence through a Lagrangian infection-generation operator and an Eulerian population-transport operator embedded in a structured UDE. Graph neural modules learn state- and time-dependent interaction strengths, while the epidemic mechanism specifies process semantics, admissible state updates, and conservation structure.
Figure 2: Fitting and forecasting trajectories in the synthetic scenarios. (a) , Lagrangian scenario, where infection generation is governed by age-structured transmission dynamics within regions. (b) , Eulerian scenario, where population transport couples regions while local infection generation and removal remain active. The vertical dashed line marks the forecast origin, and the gray-shaded region denotes the out-of-sample forecasting period. Colored envelopes show post-hoc normal-reference residual bands centered on the model predictions. All models are trained using only observations before the forecast origin.
Figure 3: Representative fitting and forecasting trajectories of LE-EpiGNODE and selected baselines on the real-world COVID-19 datasets. (a) , Ontario. (b) , Hubei. (c) , Spain. The vertical dashed line marks the forecast origin, and the gray-shaded region denotes the out-of-sample forecasting period. Models are trained using observations before the forecast origin. Ontario, Hubei, and Spain correspond to the Lagrangian, Eulerian, and multiscale propagation settings, respectively.
Figure 4: Integrated interaction fluxes for dynamical validation and process-level interpretation. (a) , LIF in the Lagrangian synthetic scenario. (b) , EIF in the Eulerian synthetic scenario. (c) , LIF in the multiscale synthetic scenario. (d) , EIF in the multiscale synthetic scenario. For synthetic systems, true and inferred fluxes are integrated over the full simulation horizon and min–max normalized for visualization. (e) , LIF in Ontario. (f) , EIF in Hubei. (g) , EIF in Spain. (h) , LIF in Spain. For real-world datasets, LIF is aggregated over regions and EIF over age groups to characterize age- and region-level interaction structure, respectively. Row and column normalization show the target distribution from each source and the source composition of each target. Diagonal conservation terms are excluded from EIF.
Model
1 week
4 weeks
12 weeks
20 weeks
Lagrangian
LE-EpiGNODE
44.38 / 5.67
51.09 / 5.84
58.66 / 8.45
63.30 / 12.45
SIR
270.78 / 34.57
299.29 / 34.21
238.84 / 34.42
186.39 / 36.68
MepoSIR
72.70 / 9.28
73.16 / 8.36
63.80 / 9.19
62.04 / 12.21
LSTM
70.16 / 8.96
119.28 / 13.63
220.36 / 31.75
373.05 / 73.40
GRU
50.56 / 6.45
87.05 / 9.95
215.85 / 31.10
383.14 / 75.39
Extended Data Table 1 — Forecasting performance on synthetic scenarios. Each cell reports MAE / WMAPE (%); the best and second-best results are shown in bold and underlined, respectively.
Ontario
Model
7 days
14 days
35 days
67 days
LE-EpiGNODE
128.57 / 3.43
165.71 / 4.49
190.60 / 7.80
248.66 / 16.20
MCM-SEIR
364.02 / 9.72
449.45 / 12.17
447.98 / 18.33
333.83 / 21.76
DCRNN
697.37 / 18.63
803.41 / 21.76
594.27 / 24.32
415.31 / 27.06
AGCRN
1255.42 / 33.53
1546.00 / 41.87
932.27 / 38.15
584.65 / 38.10
MTGNN
708.40 / 18.92
826.26 / 22.38
504.55 / 20.65
365.74 / 23.83
Extended Data Table 2 — Forecasting performance on the real-world COVID-19 datasets. Each cell reports MAE / WMAPE (%); the best and second-best results are shown in bold and underlined, respectively.
Hubei
Output
Training
7 days
14 days
26 days
OD mobility
0.0033 / 0.99
0.0120 / 28.75
0.0138 / 33.80
0.0161 / 40.15
Extended Data Table 3 — OD mobility fitting and forecasting errors. Each entry reports MAE / WMAPE (%) for the directed OD mobility values in Hubei and Spain. Errors are computed on each dataset’s original mobility scale; MAE values are not directly comparable across datasets.
Extended Data Fig. 1 — Fitting and forecasting trajectories across all coupled multiscale synthetic mechanisms. (a) , Multiscale scenarios combining the four Lagrangian infection-generation mechanisms (BMCM, TDCM, IACM and NMCM) with BRGM transport. (b) , Corresponding multiscale scenarios combining the same four Lagrangian mechanisms with CRGM transport. The vertical dashed line marks the forecast origin, and the gray-shaded region denotes the out-of-sample forecasting period. Colored envelopes show post-hoc normal-reference residual bands centered on the model predictions. All models are trained using only observations before the forecast origin.
Extended Data Fig. 2 — OD mobility extrapolation in Hubei. Observed and predicted trajectories are shown for all directed OD edges. The vertical dashed line marks the forecast origin. The auxiliary mobility predictor uses only OD observations available before the forecast origin. The overall R2 displayed in the plot summarizes the complete fitting and forecasting interval; forecasting MAE and WMAPE are reported separately in Extended Data Table Extended Data .
Extended Data Fig. 3 — OD mobility extrapolation in Spain. Observed and predicted trajectories are shown for all directed OD edges. The vertical dashed line marks the forecast origin. The auxiliary mobility predictor uses only OD observations available before the forecast origin. The overall R2 displayed in the plot summarizes the complete fitting and forecasting interval; forecasting MAE and WMAPE are reported separately in Extended Data Table Extended Data .
Mechanism
Dynamical class
Functional form
Lagrangian mechanisms
Baseline Mass-Action Contact Mechanism (BMCM)
Baseline mass-action
Fi→jL(t)=βCijkPik(t)Sjk(t)
Time-Decaying Contact Mechanism (TDCM)
Explicitly time-modulated
Fi→jL(t)=βe−rLtCijkPik(t)Sjk(t)
Infection-Avoidant Contact Mechanism (IACM)
State-dependent feedback
Fi→jL(t)=βe−αIIik(t)CijkPik(t)Sjk(t)
Nonlinear Mass-Action Contact Mechanism (NMCM)
Nonlinear state dependence
Fi→jL(t)=β[Iik(t)]qICijkPik(t)[Sjk(t)]pS
Eulerian mechanisms
Supplementary Table 1 : Synthetic dynamical mechanisms. Functional forms and dynamical classes used to generate the Lagrangian, Eulerian and multiscale synthetic systems. The eight multiscale conditions combine each Lagrangian mechanism with each Eulerian mechanism. Here, Nl and Nk are the fixed reference total populations of the source and destination regions, respectively, and N0 is the reference population used to scale the gravity terms.
Setting
Lagrangian
Eulerian
Multiscale
Graph size
M=1,A=4
M=4,A=1
M=4,A=4
Simulation horizon
0 – 200 days
0 – 200 days
0 – 200 days
Sampling interval
1 day
1 day
1 day
Transmission rate β
0.01
0.20
0.02
Removal rate γ
0.02
0.08
0.021
Mobility propensity η
–
[1]
[0.106,0.685,0.156,0.053]
Supplementary Table 2 : Synthetic scenario settings. Shared configurations used to generate the Lagrangian, Eulerian and multiscale systems.
Scenario
TDCM: rL
IACM: αI
NMCM: (pS,qI)
CRGM: rE
Lagrangian
0.01
0.10
(1.0,0.4)
–
Eulerian
–
–
–
0.30
Multiscale
0.0035
0.045
(1.0,0.2)
0.03
Supplementary Table 3 : Mechanism-specific parameters. Parameter values used by the explicitly time-modulated, state-dependent and nonlinear synthetic mechanisms.
Supplementary Fig. 1 : Synthetic trajectories in the isolated Lagrangian and Eulerian settings. Heatmaps show normalized infectious-state trajectories Ijk(t) generated by the controlled interaction mechanisms. Within each structural setting, population sizes, initial states, basic epidemiological parameters and observation rules are held fixed so that differences among trajectories primarily reflect the prescribed dynamical mechanisms.
Supplementary Fig. 2 : Multiscale ground-truth signatures of the coupled synthetic mechanisms. Normalized infection-intensity patterns are shown for the eight pairwise combinations of four Lagrangian mechanisms (BMCM, TDCM, IACM and NMCM) and two Eulerian mechanisms (BRGM and CRGM). The left block shows age–time patterns across four age groups, whereas the right block shows the corresponding city–time patterns across four cities. In both blocks, rows correspond to the four Lagrangian mechanisms and columns correspond to BRGM and CRGM. Comparisons across rows under the same Eulerian mechanism highlight Lagrangian mechanism-specific dynamical signatures, whereas comparisons between BRGM and CRGM within the same row highlight the additional effects associated with the Eulerian mobility mechanism. Together, these patterns illustrate the complementary age- and region-level signatures generated by the coupled multiscale mechanisms.
Setting
Ontario
Hubei
Spain
Epidemic setting
Age structured
Mobility driven
Coupled age–space
State structure
1 region × 4 age groups
3 regions × 1 age group
4 regions × 3 age groups
Epidemic observations
Age-stratified daily cases
City-level daily cases
Region–age daily cases
Process-aligned information
Age-contact matrix
Daily OD mobility; population and distance
Age-contact matrices; daily OD mobility; population and distance
Observation period
28 Nov. 2021–27 Feb. 2022
16 Jan. 2020–25 Feb. 2020
1 Mar. 2020–21 Jun. 2020
Training
25 days
15 days
20 days
Supplementary Table 4 : Data and evaluation settings for the three real-world COVID-19 datasets.
Input feature
Synthetic
Hubei
Spain
Source epidemic intensity
Yes
Yes
Yes
Destination epidemic intensity
Yes
Yes
Yes
Gravity prior
Yes
Yes
Yes
Normalized time
Yes
Yes
Yes
Weekly periodic encoding
No
No
Yes
Source region identity
No
No
Yes
Supplementary Table 5 : Mobility-predictor inputs. Input features used by the auxiliary mobility network in the synthetic, Hubei and Spain datasets.
Model
Role
Implementation in this study
SIR
Synthetic baseline
Independent nodewise SIR dynamics with nonnegative initial compartments and node-specific transmission and removal parameters fitted to the complete training trajectory.
MepoSIR
Synthetic baseline
Contact-coupled metapopulation SIR using the observed age-contact matrix, region-shared age-specific transmission parameters and node-specific removal parameters.
LSTM
Synthetic baseline
Stacked LSTM encoder with a direct 7-day regression head, trained on permitted node histories and rolled out autoregressively in 7-day blocks.
GRU
Synthetic baseline
Stacked GRU encoder with the same direct 7-day prediction and blockwise autoregressive rollout protocol as LSTM.
Neural ODE
Synthetic baseline
MLP vector field integrated by fixed-step RK4, with the historical window encoded as the initial latent state and decoded into a 7-day trajectory.
STGCN
Synthetic baseline
Temporal–Chebyshev graph convolution blocks using spatial supports constructed from permitted contact, distance and training-period mobility information.
Supplementary Table 6 : Baseline implementations. Experimental roles and study-specific configurations of the comparison models.
Supplementary Fig. 3 : Complete ground-truth and inferred FL(t) trajectories for BMCM and TDCM in the Lagrangian scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 4 : Complete ground-truth and inferred FL(t) trajectories for IACM and NMCM in the Lagrangian scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 5 : Complete ground-truth and inferred FE(t) trajectories for BRGM and CRGM in the Eulerian scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 6 : Complete ground-truth and inferred FL(t) trajectories for BMCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 7 : Complete ground-truth and inferred FE(t) trajectories for BMCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 8 : Complete ground-truth and inferred FL(t) trajectories for BMCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 9 : Complete ground-truth and inferred FE(t) trajectories for BMCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 10 : Complete ground-truth and inferred FL(t) trajectories for TDCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 11 : Complete ground-truth and inferred FE(t) trajectories for TDCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 12 : Complete ground-truth and inferred FL(t) trajectories for TDCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 13 : Complete ground-truth and inferred FE(t) trajectories for TDCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 14 : Complete ground-truth and inferred FL(t) trajectories for IACM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 15 : Complete ground-truth and inferred FE(t) trajectories for IACM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 16 : Complete ground-truth and inferred FL(t) trajectories for IACM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 17 : Complete ground-truth and inferred FE(t) trajectories for IACM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 18 : Complete ground-truth and inferred FL(t) trajectories for NMCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 19 : Complete ground-truth and inferred FE(t) trajectories for NMCM–BRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 20 : Complete ground-truth and inferred FL(t) trajectories for NMCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Supplementary Fig. 21 : Complete ground-truth and inferred FE(t) trajectories for NMCM–CRGM in the multiscale scenario. Blue solid lines denote the reference generating operators, and red dashed lines denote the inferred operators.
Epidemic forecasting predicts disease incidence, outbreaks, and healthcare burden across time and regions to guide surveillance and public-health decisions.
Epidemiological forecasting from surveillance data is a hard problem and hybridizing mechanistic compartmental models with neural models is a natural direction. The mechanistic structure helps keep trajectories epidemiologically plausible, while neural components can capture non-stationary, data-adaptive effects. In practice, however, many seemingly straightforward couplings fail under partial observability and continually shifting transmission dynamics driven by behavior, waning immunity, seasonality, and interventions. We catalog these failure modes and show that robust performance requires making non-stationarity explicit: we extract and extrapolate multi-scale structure from the observed infection series and use it as an interpretable control signal for a controlled neural ODE coupled to an epidemiological model. Concretely, we decompose infections into trend, seasonal, and residual components and use these signals to drive continuous-time latent dynamics while jointly forecasting and inferring time-varying transmission, recovery, and immunity-loss rates. Across early outbreak and multi-wave regimes, our approach attains the lowest RMSE on all five datasets (27-70% reduction over the strongest baseline), achieves the best peak detection accuracy, and recovers time-varying epidemiological rates within ground-truth ranges, without relying on auxiliary covariates.
Yiqi Su, Ray Lee, Jiaming Cui +1
Department of Computer Science, Virginia Tech, Alexandria, USA.
Accurate and reliable forecasting of epidemic incidences is critical for public health preparedness, yet it remains a challenging task due to complex nonlinear temporal dependencies and heterogeneous spatial interactions. Often, point forecasts generated by spatiotemporal models are unreliable in assigning uncertainty to future epidemic events. Probabilistic forecasting of epidemics is therefore crucial for providing the best or worst-case scenarios rather than a simple, often inaccurate, point estimate. We present deep spatiotemporal engression methods to generate accurate and reliable probabilistic forecasts on low-frequency epidemic datasets. The proposed methods act as distributional lenses, and out-of-sample probabilistic forecasts are generated by sampling from the trained models. Our frameworks encapsulate lightweight deep generative architectures, wherein uncertainty is quantified endogenously, driven by a pre-additive noise component during model construction. We establish geometric ergodicity and asymptotic stationarity of the spatiotemporal engression processes under mild assumptions on the network weights and pre-additive noise process. Comprehensive evaluations across six epidemiological datasets over three forecast horizons demonstrate that the proposal consistently outperforms several temporal and spatiotemporal benchmarks in both point and probabilistic forecasting. Additionally, we explore the explainability of the proposal to enhance the models' practical application for informed, timely public health interventions.
Rajdeep Pathak, Tanujit Chakraborty
SAFIR, Sorbonne University Abu Dhabi, United Arab Emirates · SCAI, Sorbonne Université, Paris, France
Spatio-temporal epidemic forecasting is critical for public health management, yet existing methods often struggle with insensitivity to weak epidemic signals, over-simplified spatial relations, and unstable parameter estimation. To address these challenges, we propose the Spatio-Temporal priOr-aware Epidemic Predictor (STOEP), a novel hybrid framework that integrates implicit spatio-temporal priors and explicit expert priors. STOEP consists of three key components: (1) Case-aware Adjacency Learning (CAL), which dynamically adjusts mobility-based regional dependencies using historical infection patterns; (2) Space-informed Parameter Estimating (SPE), which employs learnable spatial priors to amplify weak epidemic signals; and (3) Filter-based Mechanistic Forecasting (FMF), which uses an expert-guided adaptive thresholding strategy to regularize epidemic parameters. Extensive experiments on real-world COVID-19 and influenza datasets demonstrate that STOEP outperforms the best baseline by 11.1% in RMSE. The system has been deployed at a provincial CDC in China to facilitate downstream applications.
Sijie Ruan, Jinyu Li, Jia Wei +7
Beijing Institute of Technology, China · Beijing Institute of Technology, Beijing, China · Zhejiang Center for Disease Control and Prevention, China +6