Hybrid epidemic simulation framework coupling equation-based and individual-based models

Jaeyoung Kwak; Michael H. Lees; Chin Chun Ooi; Wentong Cai

physics.soc-phSep 28, 2026

Hybrid epidemic simulation framework coupling equation-based and individual-based models

Authors: Jaeyoung Kwak, Michael H. Lees, Chin Chun Ooi, Wentong Cai

Organizations: College of Computing and Data Science, Nanyang Technological University, Singapore · Informatics Institute, University of Amsterdam, Amsterdam, The Netherlands · Institute of High Performance Computing, Agency for Science, Technology and Research, Singapore

Abstract

Mass gathering events like concerts, sports matches, and festivals bring many people into close contact within a short period, creating localized bursts of infection that can shape epidemic outcomes across an entire city. To evaluate how these transient transmission events translate into broader urban impacts, we developed a simulation model linking event-scale contact dynamics with citywide commuting networks. Using Madrid, Spain, as a case study, we compared several types of gatherings and examined how their effects changed under different levels of disease transmissibility. We found that mass gatherings consistently amplified outbreak magnitude, accelerated progression, advanced district-level arrival times, and synchronized spatial spread. Remarkably, while the initial seed size generated at the event accounted for much of this acceleration, post-event transmission conditions provided complementary predictive signal regarding invasion timing. These findings demonstrate that mitigating transmission during mass gatherings can yield downstream public health benefits by delaying broader spatial spread. More generally, this multiscale framework offers a tool to evaluate how temporary, localized contact shifts produce longer-lasting consequences for urban populations.

Figures & tables

Fig 1: Distribution of active cases (exposed and infectious individuals) in the Madrid region under different mass gathering event (MGE) scenarios. Boxes indicate the interquartile range (IQR), red cross symbols $\times$ denote the median, and whiskers extend to 1.5×IQR. Results are averaged over 1000 independent stochastic simulations with basic reproduction number $R_{0}=1.5$ , event-specific transmission rate $\beta_{\text{event}}=0.5$ , and one initially infectious individual $I_{ss}=1$ .

Panel of Figure 2 — Fig 2: Temporal epidemic progression under different mass gathering event (MGE) scenarios. (a) Time evolution of active cases (exposed and infectious individuals) per 1000 persons in the Madrid region. (b) Summary of active case curves: peak time and peak incidence (active cases per 1000 persons). Different symbols indicate the mean values of peak time and peak incidence for different MGE scenarios. Rectangle width and height indicate the standard deviation of peak time and peak incidence, respectively. (c) Boxplots visualizing the distribution of final epidemic size. (d) The time required to reach the epidemic threshold of 1 cumulative case per 1000 persons across different MGE scenarios. Boxes indicate the interquartile range (IQR), red cross symbols $\times$ denote the median, and whiskers extend to 1.5×IQR. For all panels, the results are averaged over 1000 independent simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ .

Fig 3: Distribution of cumulative cases per 1000 persons. Red dashed vertical lines indicate the epidemic threshold of 1 cumulative case per 1000 individuals. Results are averaged over 1000 independent simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ .

Fig 4: Spatiotemporal epidemic progression across districts under the baseline and selected mass gathering event (MGE) scenarios. Rows correspond to the AMS_dance (top), Leipzig_3 (middle), and Baseline (bottom) scenarios, while columns show epidemic snapshots at Days 100, 150, and 200. Color intensity represents the mean number of active infections (E+I) in each district, with darker red shades indicating higher case counts. Results are averaged over 1000 independent simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ . Cross symbol $\times$ indicates the seed district where the MGE takes place.

Fig 5: Spatial distribution of mean infection arrival time across districts. Darker red shades indicate earlier epidemic arrival and lighter shades indicate later arrival. Results are averaged over 1000 independent stochastic simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ .

Panel of Figure 6 — Fig 6: Epidemic arrival time (a) Mean epidemic arrival time across districts under different mass gathering event (MGE) scenarios. Each symbol represents the mean infection arrival time for each district. For each MGE scenario, districts are clustered based on the epidemic arrival time with arrival time threshold of 4 days. Different colors indicate different clusters, suggesting noticeable transitions in the spatiotemporal expansion of the epidemics. Results are averaged over 1000 independent stochastic simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ . (b) Boxplots visualizing the interval between mean epidemic arrival time across districts under different MGE scenarios. Boxes indicate the interquartile range (IQR), red cross symbols $\times$ denote the median, and whiskers extend to 1.5×IQR.

Fig 7: Epidemic arrival time shift across districts under different mass gathering event (MGE) scenarios. Color intensity represents the arrival time shift in each district, with darker blue grids indicating earlier epidemic arrival. Results are averaged over 1000 independent stochastic simulations with $R_{0}=1.5$ , $\beta_{\text{event}}=0.5$ , and $I_{ss}=1$ .

Fig 8: Seed–destination dependence of epidemic arrival time difference between AMS_dance and baseline scenarios. Arrival-time difference between AMS_dance and baseline simulations for each combination of initial epidemic seed district (rows) and destination district (columns). Values represent the arrival time difference $\Delta T=T_{\mathrm{MGE}}-T_{\mathrm{baseline}}$ in days; negative values indicate earlier epidemic arrival under the AMS_dance scenario. Results are averaged over 1000 independent stochastic simulations with $R_{0}=1.5$ , $\beta_{\mathrm{event}}=0.5$ , and $I_{ss}=1$ . Diagonal cells corresponding to the initial seed district are omitted.

Fig 9: Kaplan–Meier survival curves of district-level epidemic invasion under different mass gathering scenarios. Kaplan–Meier survival curves showing the probability that districts remain free of infectious individuals for the baseline (no-event) scenario and five representative mass gathering scenarios (AMS_dance, AMS_football, Leipzig_1, Leipzig_2, and Leipzig_3). Earlier declines in survival probability indicate more rapid epidemic invasion.

Fig 10: Hazard ratio estimates from sequential Cox proportional hazards models ( $R_{0}=1.5$ , $\beta_{\mathrm{event}}=0.5$ ). Forest plot showing hazard ratios and 95% confidence intervals estimated from three nested Cox proportional hazards models. Model M0 includes only the gathering-event indicator. Model M1 additionally adjusts for the logarithm of the initial epidemic seed size, while Model M2 further incorporates the average pre-invasion infectious pressure computed during the days preceding district invasion. The vertical dashed line indicates a hazard ratio of one (no association).

Panel of Figure 11 — Fig 11: Performance comparison of the sequential Cox proportional hazards models. (a) Model fit evaluated using the Akaike Information Criterion (AIC). Bars represent $\Delta$ AIC relative to the best-performing model, with smaller values indicating better model fit after accounting for model complexity. (b) Model discrimination evaluated using Harrell’s concordance index (C-index), which measures the model’s ability to correctly rank district-level epidemic arrival times. A C-index of 0.5 corresponds to random discrimination, whereas a value of 1.0 indicates perfect discrimination. Models M0–M2 correspond to the sequential Cox proportional hazards models described above, with Model M2 evaluated using 3-, 7-, and 14-day averaging windows for the pre-invasion force of infection.

Fig 12: Hazard ratio estimates across epidemiological parameter combinations for the AMS_dance scenario. Heatmaps showing hazard ratios estimated across combinations of the basic reproduction number ( $R_{0}$ ) and event transmission probability ( $\beta_{\mathrm{event}}$ ). (a) Hazard ratios estimated from the parameter-only Cox model (P0). (b) Hazard ratios after adjustment for the initial epidemic seed size (P1). (c) Difference between the two hazard-ratio estimates, illustrating how accounting for initial epidemic seed size changes the estimated association between epidemiological parameters and epidemic invasion timing. Hazard ratios are reported relative to the reference parameter combination ( $R_{0}=1.0$ , $\beta_{\mathrm{event}}=0.1$ ).

Seed Location	M0 HR	M1 HR	M1 seed-size HR
District 0	2.90 (2.70–3.11)	2.66 (2.47–2.86)	1.45 (1.43–1.47)
District 7	1.99 (1.85–2.14)	1.74 (1.61–1.87)	1.83 (1.80–1.86)
District 18	2.20 (2.05–2.36)	1.95 (1.81–2.10)	1.72 (1.70–1.75)

Table 1: Sequential Cox proportional hazards results for three selected seed districts with contrasting spatial and network characteristics under

R_{0}=1.5

\beta_{\mathrm{event}}=0.5

, and

I_{\mathrm{ss}}=1

. Model M0 includes only the gathering-event indicator, while Model M1 additionally adjusts for the logarithm of the initial epidemic seed size. Values are hazard ratios with 95% confidence intervals.

Fig 13: Robustness of survival-model estimates to the location of the initial epidemic seed. Results for three selected seed districts with contrasting spatial and network characteristics. (a) Gathering-event hazard ratios from Model M0 (unadjusted) and Model M1 (adjusted for log-transformed initial epidemic seed size), with 95% confidence intervals; the dashed vertical line indicates HR = 1. (b) Harrell’s C-index for M0 and M1, with connecting lines indicating the change in discrimination after adjustment for initial seed size. (c) Improvement in model fit after adjustment for initial seed size, expressed as $\Delta\mathrm{AIC}=\mathrm{AIC}_{M0}-\mathrm{AIC}_{M1}$ ; positive values indicate better AIC-based fit for M1.

Panel of Figure 15 — Fig 14: Conceptual interpretation of the association between mass gatherings, epidemic seeding, and invasion timing. (a) Model M0 estimates the overall association between mass gathering events and district-level epidemic invasion. Gathering events are associated with substantially earlier epidemic invasion (HR = 2.90). (b) Model M1 additionally adjusts for the logarithm of the initial epidemic seed size. The gathering-event hazard ratio decreases to HR = 2.66, while larger initial epidemic seeds are independently associated with earlier invasion (HR = 1.45 per unit increase in log-transformed seed size). The attenuation of the gathering-event estimate indicates that initial epidemic seed size accounts for part of the observed association, but the remaining association indicates that seed size alone does not account for the difference in invasion timing between gathering and no-event simulations.

Panel of Figure 16 — Fig 15: Schematic representation of the hybrid simulation logic. On event days, a subset of the population is recruited into an individual-based transmission model while the remaining population continues metapopulation susceptible-exposed-infected-recovered (SEIR) dynamics. Updated epidemiological states are reintegrated before advancing to the next time step.

Scenario	Location	$N_{\text{event}}$	$k_{\text{avg}}$	$k_{\text{std}}$	Refs.
Football match (AMS_football)	Amsterdam, NL	362	17.3	9.4	[ 27 ]
Dance event (AMS_dance)	Amsterdam, NL	1048	161.7	50.0	[ 27 ]
Indoor concert scenario 1 (Leipzig_1)	Leipzig, DE	1194	35.6	11.2	[ 20 ]
Indoor concertscenario 2 (Leipzig_2)	Leipzig, DE	1158	18.5	7.2	[ 20 ]
Indoor concertscenario 3 (Leipzig_3)	Leipzig, DE	1054	8.7	4.2	[ 20 ]

Table 2: Summary statistics for the individual-based model (IBM) datasets.

N_{\text{event}}

represents the total number of attendees. The mean nodal degree (

k_{\text{avg}}

) and its standard deviation (

k_{\text{std}}

) describe the distribution of unique contacts per person within the defined proximity threshold of 1.5 m.

Symbol	Description
$n$	Number of districts
$i$	Home-district index
$j$	Destination-district index
$\ell$	Generic district index
$N_{ij}$	Population whose home district is $i$ and current destination is $j$
$D_{ij}$	Probability of travel from district $i$ to district $j$ during the away phase

Table 3: Parameters and notation used in the stochastic metapopulation SEIR model

Symbol	Description
$i$	District index
$r$	Stochastic simulation replicate index
$T_{i}^{\mathrm{arr}}$	Epidemic arrival time in district $i$ , expressed in calendar days
$G$	Indicator for MGE versus baseline simulation
$Z$	Event-generated epidemic seed
$\Lambda_{i}(d)$	Daily transmission pressure in district $i$ on calendar day $d$

Table 4: Parameters and notation used in the presented statistical analysis of district-level epidemic invasion timing

Fig 16: Mobility matrix heatmap of inter-district commuting patterns in Madrid. Columns represent the fraction of residents from each district who commute to other districts, while rows indicate the corresponding inflow of commuters into each destination district. Darker red colors denote higher commuting volumes between district pairs. The pronounced diagonal elements indicate that the majority of residents remain within their home districts rather than commuting elsewhere on a daily basis.

Fig 18: Mean epidemic arrival time for the baseline scenario as a function of seed and destination districts. Each cell represents the mean arrival time over 1,000 simulation runs for the corresponding seed–destination pair. Rows indicate the seed district and columns indicate the destination district. The diagonal is masked because the destination coincides with the initial epidemic seed. Results correspond to $R_{0}=1.5$ and $\beta_{\mathrm{event}}=0.5$ .

Fig 19: Mean epidemic arrival time for the AMS dance scenario as a function of seed and destination districts. Each cell represents the mean arrival time over 1,000 simulation runs for the corresponding seed–destination pair. Rows indicate the seed district and columns indicate the destination district. The diagonal is masked because the destination coincides with the initial epidemic seed. Results correspond to $R_{0}=1.5$ and $\beta_{\mathrm{event}}=0.5$ .

Explore similar work

A Hybrid ABM-PDE Framework for Real-World Infectious Disease Simulations
Apr 11, 2025Kristina Kehrer, Tim O. F. ConradAgent-Based Model Epidemiological Models
Zuse Institute Berlin Berlin, 14195
An Infectious Disease Spread Simulation Based on Large Language Model Decision Making
Jun 4, 2026Yonchanok Khaokaew, Ruochen Kong, Andreas Zufle +6Epidemiological Models Large Language Model Decisions
Computer Science and Engineering Faculty of Engineering The University of New South Wales Sydney, NSW, Australia · King Mongkut’s University of Technology North Bangkok (KMUTNB) · Department of Computer Science Emory University Atlanta, Georgia, USA +5
Beyond Time Series: Spatial Reasoning for Epidemic Forecasting via Multimodal Learning
Jun 20, 2026Diana Guadalupe Gomez, Chenwei Wu, Zhiyi Wang +2Epidemiological Models Spatio-Temporal Forecasting
University of Michigan Ann Arbor, MI, USA

Apr 11, 2025, cs.MAJ/K move · Enter open · S save

A Hybrid ABM-PDE Framework for Real-World Infectious Disease Simulations

Kristina Kehrer, Tim O. F. Conrad

Zuse Institute Berlin Berlin, 14195

Agent-Based Model Epidemiological Models

Open paper AskPDF arXiv

This paper presents a hybrid modeling approach that couples an Agent-Based Model (ABM) with a partial differential equation (PDE) model in an epidemic setting to simulate the spatial spread of infectious diseases using a compartmental structure with seven health states. The goal is to reduce the computational complexity of a full-ABM by introducing a coupled ABM-PDE model that offers significantly faster simulations while maintaining comparable accuracy. Our results demonstrate that the hybrid model not only reduces the overall simulation runtime (defined as the number of runs required for stable results multiplied by the duration of a single run) but also achieves smaller errors across both 25% and 100% population samples. The coupling mechanism ensures consistency at the model interface: agents crossing from the ABM into the PDE domain are removed and represented as density contributions, while surplus density in the PDE domain is used to generate agents with plausible trajectories derived from mobile phone data. We evaluate the hybrid model using real-world mobility and infection data for the Berlin-Brandenburg region in Germany, showing that it captures the core epidemiological dynamics while enabling efficient large-scale simulations. These results demonstrate that the proposed ABM-PDE framework provides a robust and computationally efficient alternative to full-scale agent-based simulations, making it suitable for realistic epidemic modeling and scenario analysis.

Abstract

Explore similar work

A Hybrid ABM-PDE Framework for Real-World Infectious Disease Simulations

An Infectious Disease Spread Simulation Based on Large Language Model Decision Making

Beyond Time Series: Spatial Reasoning for Epidemic Forecasting via Multimodal Learning