Learning and Forecasting Multiscale Spatiotemporal Epidemic Dynamics with Lagrangian-Eulerian Graph Neural ODEs

Jinyi Duan; Qian Li; Biao Tang

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

R^{2}

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

R^{2}

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	$\displaystyle\mathcal{F}_{i\rightarrow j}^{L}(t)=\beta C_{ij}^{k}\frac{S_{j}^{k}(t)}{P_{i}^{k}(t)}$
Time-Decaying Contact Mechanism (TDCM)	Explicitly time-modulated	$\displaystyle\mathcal{F}_{i\rightarrow j}^{L}(t)=\beta e^{-r_{L}t}C_{ij}^{k}\frac{S_{j}^{k}(t)}{P_{i}^{k}(t)}$
Infection-Avoidant Contact Mechanism (IACM)	State-dependent feedback	$\displaystyle\mathcal{F}_{i\rightarrow j}^{L}(t)=\beta e^{-\alpha_{I}I_{i}^{k}(t)}C_{ij}^{k}\frac{S_{j}^{k}(t)}{P_{i}^{k}(t)}$
Nonlinear Mass-Action Contact Mechanism (NMCM)	Nonlinear state dependence	$\displaystyle\mathcal{F}_{i\rightarrow j}^{L}(t)=\beta[I_{i}^{k}(t)]^{q_{I}}C_{ij}^{k}\frac{[S_{j}^{k}(t)]^{p_{S}}}{P_{i}^{k}(t)}$
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,

N^{l}

and

N^{k}

are the fixed reference total populations of the source and destination regions, respectively, and

N_{0}

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 $\beta$	$0.01$	$0.20$	$0.02$
Removal rate $\gamma$	$0.02$	$0.08$	$0.021$
Mobility propensity $\eta$	–	$[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: $r_{L}$	IACM: $\alpha_{I}$	NMCM: $(p_{S},q_{I})$	CRGM: $r_{E}$
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

I_{j}^{k}(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 $\times$ 4 age groups	3 regions $\times$ 1 age group	4 regions $\times$ 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

\mathcal{F}^{L}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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

\mathcal{F}^{L}(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

\mathcal{F}^{E}(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.

Learning and Forecasting Multiscale Spatiotemporal Epidemic Dynamics with Lagrangian-Eulerian Graph Neural ODEs

Abstract

Figures & tables

Epidemic Forecasting

How (Not) to Hybridize Neural and Mechanistic Models for Epidemiological Forecasting

Deep Generative Spatiotemporal Engression for Probabilistic Forecasting of Epidemics

Prior Knowledge-enhanced Spatio-temporal Epidemic Forecasting

Abstract

Epidemic Forecasting

Explore similar work

How (Not) to Hybridize Neural and Mechanistic Models for Epidemiological Forecasting

Deep Generative Spatiotemporal Engression for Probabilistic Forecasting of Epidemics

Prior Knowledge-enhanced Spatio-temporal Epidemic Forecasting