Scaling Epidemic Inference on Contact Networks: Theory and Algorithms

We sincerely appreciate your efforts to review our paper and provide valuable suggestions. Below we address each concern in detail.

• q1: Clarify relevance to NeurIPS.
• R1: Thank you for the comment regarding NeurIPS relevance. We respectfully point out that the NeurIPS 2025 Call for Papers explicitly welcomes research in probabilistic methods, scalable inference, optimization, theory, and machine learning for sciences—including health, network science, and the natural sciences. Our work addresses a core challenge in efficient, deterministic inference for structured infection stochastic processes on large networks, bridging algorithmic theory, scalable probabilistic modeling, and ML for science. We believe RAPID addresses a key bottleneck for scalable probabilistic inference under uncertainty on networks, which are central to many areas of modern ML.

Furthermore, probabilistic inference on graphs is a classic topic in machine learning as evidenced by Belief Propagation (BP), which was first published at AAAI’82 [1] and remains foundational in the field. We thus believe that our RAPID algorithm and its advances in scalable, deterministic inference will be of significant interest to the NeurIPS community.

• q2: Justify novelty beyond existing approximation methods.
• R2: As stated in our introduction, our contribution is twofold: (1) we derive a variance lower bound for Monte Carlo estimators on large contact networks, revealing how reproduction number, average out-degree, initial infection ratio, and simulation count affect estimator variance. This result is also empirically validated in Figure 1; (2) we develop RAPID, an efficient epidemic inference algorithm that innovatively replaces the residual in classical LocalPush algorithms with Jacobian-based perturbation terms, enabling asynchronous, selective updates instead of synchronous iterations. This residual-driven approach achieves 4-12x computational speedup while maintaining accuracy. We also establish worst-case complexity analysis and derive connections between propagation residuals and epidemic thresholds. Extensive experiments validate our method's effectiveness and efficiency. Reviewers iUqy and MyzN both recognized the novelty and practical impact of our algorithm.

Regarding the approximation methods mentioned, we have clarified their relationship to our algorithm in the preliminaries and related work. Belief Propagation (a variant of message passing) operates on factor graphs with sum-product updates for general probabilistic inference; in fact, our PID baseline (Algorithm 1) is equivalent to the BP implementation for epidemic diffusion on contact networks, as in [6]. Mean-field approximation assumes independence among variables to simplify joint distributions and forms the foundation of traditional MP/BP methods. Fokker–Planck equations describe the global density evolution of continuous systems via partial differential equations, and are unrelated to our discrete, network-based approach.

RAPID advances beyond epidemic BP by pioneering the use of a local push framework and Jacobian-based residuals for forward inference in nonlinear dynamics. This asynchronous, selective propagation strategy is not only novel but also highly efficient, making RAPID widely applicable to disease modeling, rumor spreading, and general network-based diffusion processes that require scalable, deterministic inference.

• q3: Use more realistic contact network data.
• R3: Thank you for the suggestion. Contact networks derived from mobility data providers like SafeGraph would indeed be valuable. However, we have also considered realistic settings by using hiv-Trans, a real-world HIV transmission network from the U.S. (1988-2001) collected by the National Addiction & HIV Data Archive Program [2]. Other datasets like email-EuAll and bitcoin-Alpha are commonly used for epidemic simulations [3,4], as viral spread shares similarities with rumor propagation and computer viruses. Following your suggestion, we added a hospital contact network derived from patient-provider interactions at Carilion Hospital, Virginia [5]. This dataset contains 11,810 nodes and 30,994 edges (undirected edges converted to bidirectional). Results compared against MC-1000 (for sufficiently small variance) are shown below, and we will incorporate these findings into Figure 2, Figure 4, and Tables 2-4.

• q4: Clarify when your approximations are accurate.
• R4: Our approximations are accurate when the contact network is sparse or locally tree-like. As mentioned in Remark 1 (Appendix), setting epidemic parameters $\beta$ and $\gamma$ to sufficiently small and proportionally reduced values ensures the linear approximation is accurate when modeling continuous infection processes. The primary source of error comes from the independent path assumption in Algorithm 1, which introduces inaccuracies due to non-independent path interactions (such as network loops). However, we emphasize that this is a standard and widely adopted assumption in network-based epidemic modeling. [6] shows Algorithm 1 is exact on tree graphs and provides rigorous upper bounds for infected individuals on arbitrary networks. [7] proves it typically yields accurate estimations on real sparse networks for a broad class of spreading dynamics. Therefore, RAPID's approximation errors are primarily related to graph sparsity and local structure: it achieves high accuracy on commonly encountered sparse networks or locally tree-like graphs, while on denser networks it may overestimate marginal infection probabilities—but still provides a meaningful worst-case estimate, which is valuable for robust epidemic assessment. As shown in Table 3, RAPID achieves highly accurate estimates with MAE of 0.0127 and 0.0103 on the sparsest networks hiv-Trans and email-EuAll, respectively.

[1] Pearl, Judea. "Reverend Bayes on inference engines: A distributed hierarchical approach." Proceedings of the Second AAAI Conference on Artificial Intelligence (AAAI-82), 1982, pp. 133–136.

[2] Morris, Martina, and Richard Rothenberg. "Hiv transmission network metastudy project: An archive of data from eight network studies, 1988--2001." (2011).

[3] He, Yinhan, et al. "Demystify Epidemic Containment in Directed Networks: Theory and Algorithms." Proceedings of the Eighteenth ACM International Conference on Web Search and Data Mining. 2025.

[4] Ru, Xiaolei, et al. "Inferring patient zero on temporal networks via graph neural networks." Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 37. No. 8. 2023.

[5] Adhikari, Bijaya, et al. "Fast and near-optimal monitoring for healthcare acquired infection outbreaks." PLoS computational biology 15.9 (2019): e1007284.

[6] Karrer, Brian, and Mark EJ Newman. "Message passing approach for general epidemic models." Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 82.1 (2010): 016101.

[7] Lokhov, Andrey Y., Marc Mézard, and Lenka Zdeborová. "Dynamic message-passing equations for models with unidirectional dynamics." Physical Review E 91.1 (2015): 012811.

We sincerely appreciate your efforts to review our paper and provide valuable suggestions. Below we address each concern in detail.

• w1: Purpose for the variance lower bound proved in Theorem 3.1 for Monte Carlo.
• R w1: The variance lower bound proved in Theorem 3.1 for Monte Carlo both motivates the development of efficient deterministic methods and highlights scenarios where MC simulations become computationally expensive—precisely the cases where RAPID offers the greatest advantage. Theorem 3.1 analyzes how Monte Carlo estimator variance depends on network structure and epidemic parameters $\beta,\gamma$, showing that under certain conditions, this variance can become very large. This high variance is intrinsic to the network itself, as discussed in Remark 1 and empirically validated on ER graphs in Figure 1. For example, small-world networks like hiv-Trans (sparse with small diameter) show higher MC simulation variance than other networks like email-EuAll under identical settings (also shown in R4), necessitating MC-1000 for stable, reliable results.

In contrast, PID and RAPID are deterministic methods that do not involve sampling, so they do not have estimator variance or confidence intervals. Their accuracy depends only on approximation error with respect to the true infection probabilities, not on stochastic variation.

• w2: The experiments are all somewhat "synthetic".
• R w2: Due to its theoretical nature, both classical [5,6] and modern work [2,3] in network-based epidemiology rely heavily on simulation-based validation. However, we have considered realistic epidemic settings by using hiv-Trans, a real-world HIV transmission network from the U.S. (1988-2001) collected by the National Addiction & HIV Data Archive Program [1]. Other datasets like email-EuAll and bitcoin-Alpha are commonly used for epidemic simulations [2,3], as viral spread shares similarities with rumor propagation and computer viruses propagation. Following your suggestion, we added a hospital contact network derived from patient-provider interactions at Carilion Hospital, Virginia [4]. This dataset contains 11,810 nodes and 30,994 edges (undirected edges converted to bidirectional). Results compared against MC-1000 (for sufficiently small variance) are shown below, and we will incorporate these findings into Figure 2, Figure 4, and Tables 2-4.

• w3: Not compared with Bayesian inference (e.g. approximate Bayesian computation (ABC) or MCMC).
• R w3: Bayesian methods (e.g., ABC, MCMC) solve the inverse problem: estimating epidemic parameters from observed time-series data (daily cases, deaths, etc.). Our RAPID method solves the inference problem: given known epidemic parameters and contact network structure, efficiently computing node-level infection probabilities at convergence. These are fundamentally different computational tasks with different inputs (time-series vs. network structure) and outputs (parameter estimates vs. node probabilities).

• q1: The Monte Carlo simulations have stochasticity coming from random sampling of the initial infection node set.
• R1: No, the stochasticity in MC simulations comes from random transmission and recovery events during propagation, not from sampling of the initial infected nodes. In simulations/mc.py, both run_sir_simulation and run_monte_carlo_simulations use a fixed set of initial infected node indices as inputs. The test_sir_simulation function, which samples initial seeds randomly, is only used for testing and is not called by the main code. The initial infected nodes are selected once in main.py (lines 207–208) and are kept the same across all methods (MC, PID, RAPID, centrality measures) to ensure a fair and direct comparison. This follows standard experimental protocol in epidemic modeling and matches our problem definition in the preliminaries ("given initial infected node set $\mathcal{I}_0$").

Additionally, as shown in Appendix D.4, we tested different distributions for initial infected nodes, and the conclusions remained consistent, further demonstrating the efficiency and effectiveness of our RAPID algorithm.

• q2: Can this framework be transferred into the Bayesian inference setting.
• R2: Our method and Bayesian inference address different epidemic modeling problems. Our approach focuses on efficiently computing node-level marginal state distributions given known parameters and network structure. In contrast, Bayesian modeling estimates epidemic parameters from macroscopic observations like aggregate infection cases. When complete contact networks are available (e.g., from sensor data or mobility data), our RAPID method can accelerate forward simulations for methods like ABC, enabling faster parameter updates through efficient simulation-observation comparisons.

• q3: Why is RAPID's accuracy higher than PID.
• R3: This is an insightful question. As [7] notes in section "Dynamic Belief Propagation", when contact graphs contain many short loops, nodes within these loops cyclically influence each other's states over time, reducing the accuracy of the original PID algorithm. Our RAPID algorithm improves upon this through asynchronous max-heap iterations that accelerate convergence and selective updates that avoid some short loops, resulting in higher accuracy than PID.

• q4: How accurate are 50-run Monte Carlo simulations.
• R4: This is an important methodological concern. We provide the standard deviation of global mean infection probability for MC-50 across datasets (see table below; for hiv-Trans and carilion-Hospital, we use MC-1000 for higher reliability). For computational efficiency, we use MC-50 as baseline on large networks, and report the standard deviation across independent simulations as error bars in Table 3 and Table 4, ensuring that stochasticity is fully accounted for and avoiding bias from single-simulation outliers. These standard deviations and error bars (e.g. 1.27±0.47 for RAPID on hiv-Trans in Table 3) are consistently one to two orders of magnitude smaller than the effect sizes (e.g., MAE), confirming that MC-50 is sufficient for stable inference in practice. When propagation variance is higher (as in hiv-Trans and carilion-Hospital), we switch to MC-1000 to ensure robustness.

*High-variance datasets use MC-1000 as baseline.

• q5: Paper Formatting Concerns.
• R5: Thank you for pointing out. We will delete the repeated "variance of" in Line 113 and repeated "state" in Line 76. We will also change the title of Table 2 to "Dataset Statistics".

[1] Morris, Martina, and Richard Rothenberg. "Hiv transmission network metastudy project: An archive of data from eight network studies, 1988--2001." (2011).

[2] He, Yinhan, et al. "Demystify Epidemic Containment in Directed Networks: Theory and Algorithms." Proceedings of the Eighteenth ACM International Conference on Web Search and Data Mining. 2025.

[3] Ru, Xiaolei, et al. "Inferring patient zero on temporal networks via graph neural networks." Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 37. No. 8. 2023.

[4] Adhikari, Bijaya, et al. "Fast and near-optimal monitoring for healthcare acquired infection outbreaks." PLoS computational biology 15.9 (2019): e1007284.

[5] Prakash, B. Aditya, et al. "Threshold conditions for arbitrary cascade models on arbitrary networks." Knowledge and information systems 33.3 (2012): 549-575.

[6]Karrer, Brian, and Mark EJ Newman. "Message passing approach for general epidemic models." Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 82.1 (2010): 016101.

[7] Lokhov, Andrey Y., Marc Mézard, and Lenka Zdeborová. "Dynamic message-passing equations for models with unidirectional dynamics." Physical Review E 91.1 (2015): 012811.

Dear Reviewer,

I hope this message finds you well. As the discussion period is nearing its end with less than one day remaining, I wanted to ensure we have addressed all your concerns satisfactorily. If there are any additional points or feedback you'd like us to consider, please let us know. Your insights are valuable to us, and we're eager to address any remaining issues to improve our work.

Thank you for your time and effort in reviewing our paper.

Best regards, The authors of Submission 17683

We sincerely appreciate your efforts to review our paper and provide valuable suggestions. Below we address each concern in detail.

• w1: The error inherent to their theoretical simplifications is not well explained.
• R w1: The approximation in Theorem 3.1 mainly relies on the sparse graph assumption, while the linearization approximation in Theorem 3.2 is justified by Appendix Remark 1. The error inherent to our RAPID algorithm comes from two main sources: (1) the linearization approximation error when computing the propagation residual, and (2) the independent path assumption used during both the preheating step and state updates, which matches the error sources in Algorithm 1. As discussed in Remark 1 (Appendix), the linearization error can be made negligible by proportionally reducing the epidemic parameters $\beta$ and $\gamma$ to sufficiently small values when modeling continuous infection processes.

The primary source of RAPID's approximation error is the independent path assumption. While this assumption is standard in network epidemiology analysis, it introduces errors from non-independent path interactions (e.g., network loops). As shown by [1], Algorithm 1 is exact on tree graphs and provides rigorous upper bounds for the number of infected individuals on arbitrary networks. [2] demonstrates that this approach typically yields accurate estimates on real sparse networks for a broad class of spreading dynamics. Therefore, RAPID’s approximation error is primarily influenced by graph sparsity and local structure: it achieves high accuracy on commonly encountered sparse or locally tree-like networks, while on denser networks it may overestimate marginal infection probabilities—but still provides a meaningful worst-case estimate, which is valuable for robust epidemic assessment. As shown in Table 3, RAPID achieves highly accurate results with MAE of 0.0127 and 0.0103 on the sparsest networks (hiv-Trans and email-EuAll), respectively.

• q1: The most predictable weaknesses of the linear approximation used in RAPID.
• R1: As discussed in Appendix Remark 1, epidemic parameters $\beta,\gamma$ need to be sufficiently small for the linear approximation to be accurate. If the parameters are large or the network is dense, RAPID may underestimate propagation residuals by neglecting higher-order effects, leading to missed activations for some nodes. The table below (bitcoin-Alpha) empirically demonstrates that smaller $\beta$ values result in lower MAE: | $\beta$ | $\gamma$ | MAE ($\times 10^{-2}$) | | ------- | -------- | ---------------------- | | 0.2 | 0.4 | 6.86±0.08 | | 0.1 | 0.2 | 5.72±0.10 | | 0.05 | 0.1 | 4.25±0.06 | | 0.01 | 0.02 | 3.01±0.11 |

• q2: How to choose $\epsilon$.
• R2: We have provided sensitivity analysis for both preheating steps $p$ and residual threshold $\epsilon$ in Appendix Section D.5 and Figure 10. For $\epsilon$ specifically, our analysis shows that large values skip informative updates while very small values have diminishing influence since residual scales are approximatedly bounded by $\bar{k}\beta$. In our experiments, $\epsilon=10^{-3}$ works well across datasets. Generally, $\epsilon$ should be 1-2 orders of magnitude smaller than $\bar{k}\beta$ to capture significant state changes while maintaining efficiency.

• q3: Is it only possible for nodes to be activated once in the complexity analysis?
• R3: No, in our algorithm, nodes can be activated and updated multiple times during asynchronous propagation. In our worst-case complexity analysis, we account for up to $\mathcal{O}(N)$ node activations. Here we present statistics on activation frequency versus node count for the bitcoin-Alpha dataset: | activate times | 0 | 1 | 2 | 3 | 4 | 5 | | -------------- | ---- | ---- | ---- | ---- | ---- | ---- | | Node count | 35 | 3730 | 12 | 2 | 3 | 1 |

We can see from the table that most nodes are activated once, making RAPID highly efficient.

• q4: $\ell_{ij}$ under equation 5.
• R4: Sorry for the confusion. $\ell_{ij}$ represents the path length between nodes $i$ and $j$. We will add this to Table 1 for clarity.

• q5: Content arangement.
• R5: We will move Remark 1 to precede the linearization analysis (i.e., starting from line 157 in the main text). Thank you for the feedback!

[1] Karrer, Brian, and Mark EJ Newman. "Message passing approach for general epidemic models." Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 82.1 (2010): 016101.

[2] Lokhov, Andrey Y., Marc Mézard, and Lenka Zdeborová. "Dynamic message-passing equations for models with unidirectional dynamics." Physical Review E 91.1 (2015): 012811.

The authors' responses were very helpful! Especially in relation to high accuracy in sparse or local tree-like structures. I have updated my score contingent on a more extensive discussion of the practical implications of theoretical simplifications being included in the paper.

We sincerely appreciate your efforts to review our paper and provide valuable suggestions. Below we address each concern in detail.

• w1: The model is limited to non-reinfection dynamics.
• R w1: Our method targets non-reinfection epidemic models (SIR and variants), which represent a well-established and widely important problem setting in epidemic modeling. This framework naturally applies to many practical scenarios: acute outbreaks (SARS, MERS), high-fatality diseases (Ebola, Marburg), vaccine-preventable diseases (measles, MPox), and cases where immunity duration exceeds analysis time frames. Extensive prior work in network epidemiology focuses on this same model class [1,2], reflecting its fundamental importance.

Our contribution achieves significant computational improvements for these established and practically relevant models. We plan to extend our algorithm to reinfection diseases in future work.

• q1: Insight or empirical justification for the preheating phase.
• R1: The preheating phase is designed to ensure that the linear approximation error remains small when the epidemic parameter $\beta$ is not sufficiently small. Theoretically, if both $\beta$ and $\gamma$ are proportionally reduced to very small values according to Remark 1 in the Appendix, the linear approximation becomes highly accurate and preheating is unnecessary. However, in our experiments with $\beta=1/18$ and moderately dense networks like email-Enron ($\bar{k}=10.02$), the product $\beta\bar{k}$ is not much less than 1, resulting in non-negligible linearization error. When further scaling of parameters is not practical, the preheating phase guarantees that $\max_{v\in\mathcal{V}} P_I^v\leq \delta$, so that $\beta\bar{k}\delta\ll1$ and the linear approximation remains valid. Preheating iterates rapidly as only a small fraction of nodes are initially affected. Therefore, we recommend choosing the number of preheating steps $p$ based on a target $\delta$. Our sensitivity analysis in Appendix D.5 explores $p$ values from 0 to 50 on email-Enron, effectively including both ablation ($p=0$) and varying preheating, providing a more comprehensive understanding than a simple ablation study.

• q2: Bitcoin-Alpha data is not found under /data/.
• R2: Sorry for the confusion. We uploaded bitcoin-Alpha under the original filename "soc-sign-bitcoinalpha". We did not upload soc-Pokec due to file size limits, but it can be downloaded from Stanford's SNAP Dataset (social networks section). Our code includes the necessary preprocessing logic for this dataset.

• q3: Where is Network centrality measures and how it works.
• R3: Figure 2 shows the Kendall-tau rank correlation and runtime trade-off between RAPID and baselines including network centrality measures against MC-50. We use node centrality rankings (e.g., degree, betweenness) as proxies for infection probability rankings, as high-centrality nodes are commonly considered high-risk in epidemic studies before [3,4].

While centrality measures show some alignment with true infection probabilities in certain contexts, no single centrality measure consistently ranks nodes by infection probability. This limitation becomes more pronounced at higher initial infection ratios $\alpha$, as shown in Appendix D.3.

• q4: Why MC-50 showed in Table 4 but not in Table 3.
• R4: In our experimental setup, MC-50 serves as the ground truth for infection probabilities. Therefore, Table 3 shows MAE relative to this baseline, while Table 4 compares absolute runtime performance including MC-50.

[1] Karrer, Brian, and Mark EJ Newman. "Message passing approach for general epidemic models." Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 82.1 (2010): 016101.

[2] Lokhov, Andrey Y., Marc Mézard, and Lenka Zdeborová. "Dynamic message-passing equations for models with unidirectional dynamics." Physical Review E 91.1 (2015): 012811.

[3] Bell, David C., John S. Atkinson, and Jerry W. Carlson. "Centrality measures for disease transmission networks." Social networks 21.1 (1999): 1-21.

[4] Cohen, Reuven, Shlomo Havlin, and Daniel Ben-Avraham. "Efficient immunization strategies for computer networks and populations." Physical review letters 91.24 (2003): 247901.

Dear Reviewer,

I hope this message finds you well. As the discussion period is nearing its end with less than one day remaining, I wanted to ensure we have addressed all your concerns satisfactorily. If there are any additional points or feedback you'd like us to consider, please let us know. Your insights are valuable to us, and we're eager to address any remaining issues to improve our work.

Thank you for your time and effort in reviewing our paper.

Best regards, The authors of Submission 17683