Pith. sign in

REVIEW 5 minor 54 references

Three gradient estimators make the Gillespie SSA usable for gradient-based parameter fitting, with complementary strengths and failure modes.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · grok-4.5

2026-07-13 13:57 UTC pith:VXBRTODF

load-bearing objection Solid methods paper that adapts three ML gradient estimators to Gillespie SSA, derives their variance scalings, and shows when GS-ST fails while SF stays usable.

arxiv 2604.02121 v2 pith:VXBRTODF submitted 2026-04-02 physics.comp-ph cond-mat.stat-mechcs.LGphysics.bio-phphysics.chem-ph

Gradient estimators for parameter inference in discrete stochastic kinetic models

Ludwig Burger , Annalena Kofler , Lukas Heinrich , Ulrich Gerland This is my paper
classification physics.comp-ph cond-mat.stat-mechcs.LGphysics.bio-phphysics.chem-ph
keywords stochastic kinetic modelsGillespie SSAgradient estimatorsparameter inferenceGumbel-Softmaxscore functionrepressilatordifferentiable simulators
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

Fitting parameters of stochastic chemical reaction networks usually means working without gradients, because the Gillespie algorithm samples discrete reaction events and waiting times that automatic differentiation cannot handle. This paper shows how three gradient estimators developed in machine learning—Gumbel-Softmax Straight-Through, Score Function, and Alternative Path—can be adapted to the Gillespie algorithm so that steady-state and time-dependent observables become differentiable with respect to rate constants. On a simple bimolecular association process the estimators recover correct gradients, but their variances scale differently with trajectory length and temperature-like hyperparameters; Gumbel-Softmax can explode exponentially in hard regimes while Score Function grows only linearly. The same estimators are then used to recover production and binding rates of a noisy three-gene repressilator from short oscillatory trajectories. Score Function succeeds on every trial; Gumbel-Softmax fails when variance diverges or when temperature is raised enough to introduce bias. The practical message is that gradient-based inference is now feasible for exact stochastic kinetic models, provided the user matches the estimator to the parameter regime.

Core claim

Gradient-based parameter inference can be effectively combined with the Gillespie stochastic simulation algorithm by adapting the Gumbel-Softmax Straight-Through, Score Function and Alternative Path estimators. Gumbel-Softmax generally produces low-variance gradients but can diverge in challenging regimes and thereby break inference, while the unbiased Score Function estimator remains robust; the three estimators therefore offer complementary advantages for both steady-state and time-dependent observables.

What carries the argument

Three Monte-Carlo gradient estimators (GS-ST, SF, AP) that convert the non-differentiable sampling of reaction channels and waiting times into unbiased or controllable-bias estimates of the gradient of an observable with respect to kinetic parameters, so that stochastic gradient descent can be run directly on Gillespie trajectories.

Load-bearing premise

The two biophysical models and the particular loss functions chosen are assumed to be representative enough that the observed variance-scaling regimes and the ranking of the three estimators will carry over to other reaction networks and observation schemes.

What would settle it

Apply the same three estimators, with identical trajectory budgets and temperatures, to a different network (for example a larger gene circuit or a non-oscillatory multi-species cascade) and check whether GS-ST still fails only when variance diverges and whether SF still recovers the ground-truth parameters on every trial.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

0 major / 5 minor

Summary. The manuscript adapts three machine-learning gradient estimators (Gumbel-Softmax Straight-Through, Score Function, and Alternative Path) to the Gillespie SSA so that gradients of steady-state and time-dependent observables can be obtained for parameter inference in discrete stochastic kinetic models. After defining the estimators for fixed-step trajectories and extending them to fixed-time trajectories (including waiting-time contributions), the authors characterize variance scaling analytically and numerically on the exactly solvable bimolecular association process, then demonstrate stochastic-gradient inference on the oscillatory repressilator over 50 random reference/initialization pairs. The central claim is that GS-ST is often low-variance but can diverge or become biased in hard regimes, while SF remains more robust (linear variance growth) and AP is generally higher-variance; the estimators are therefore complementary.

Significance. If the reported constructions and variance regimes hold, the work supplies a practical, reproducible route to gradient-based inference for a class of models that has historically relied on moment closures, approximate likelihoods, or likelihood-free methods. Strengths include exact Master-equation benchmarks for the association model (Appendix A), analytic variance scalings (Lyapunov exponent for GS-ST; linear score accumulation for SF; weight/path-difference interplay for AP), explicit SNR stopping criteria, and a public JAX implementation. Concurrent GS-ST/relaxation papers are acknowledged; the distinctive contribution is the systematic multi-estimator comparison that identifies concrete failure modes of GS-ST and the relative robustness of SF. The two-model scope is an ordinary methods caveat rather than a load-bearing flaw.

minor comments (5)
  1. The AP superlinear variance crossover for time-dependent observables is deferred to Supplementary Fig. 9 and Section III C; a one-sentence pointer in the main text (near Fig. 5) would make the ranking of estimators self-contained for readers who do not open the supplement.
  2. Notation for the one-hot / relaxed reaction indicator (X_s vs X_relaxed) and the stoichiometric update is introduced cleanly in Sec. II A but is reused without re-definition in the time-dependent section; a brief reminder would help.
  3. Fig. 1d error bars for low-τ GS-ST are stated to exceed the plotted range by many orders of magnitude; a log-scale inset or a note in the caption would make the divergence visually clearer.
  4. The concurrent arXiv works [34–36] are cited; a short comparative sentence on how the present multi-estimator variance analysis differs from those single-estimator demonstrations would strengthen the novelty paragraph in the Introduction.
  5. Appendix C (bias-driven failure at high τ) is important for the practical message; consider elevating a one-panel summary into the main Fig. 7 or 9 so that the bias–variance trade-off is not only in the appendix.

Circularity Check

0 steps flagged

No significant circularity: estimators are adapted from ML literature, validated against exact Master-equation solutions and independent Gillespie trajectories, and inference recovers held-out reference parameters without tautological fits.

full rationale

The paper adapts three established gradient estimators (GS-ST, SF, AP) to the Gillespie SSA, derives their forms for steady-state and time-dependent observables, and characterizes variance scaling analytically (Lyapunov exponent for GS-ST multiplicative accumulation; linear score accumulation for SF; weight/path-difference interplay for AP). All numerical results are generated from independent stochastic trajectories. For the bimolecular association model the reference distributions and gradients come from the exact null-space / matrix-exponential solution of the Chemical Master Equation (Appendix A), not from the estimators themselves. For the repressilator, 50 independent SGD runs recover randomly sampled ground-truth parameters from randomly initialized starting points using a fixed loss on log-mean copy numbers; failures of GS-ST are diagnosed by elevated variance or bias relative to the unbiased SF estimator, not by re-fitting the same data. No parameter is fitted to data and then re-presented as a prediction; no uniqueness theorem or ansatz is imported via self-citation; concurrent works are cited only for context. The derivation chain is therefore self-contained and externally falsifiable. Score 0 is the correct outcome.

Axiom & Free-Parameter Ledger

4 free parameters · 4 axioms · 0 invented entities

The work rests on the standard continuous-time Markov-chain formulation of chemical kinetics and on three published gradient estimators. No new physical entities are postulated; free parameters are ordinary hyper-parameters of the estimators and optimizers, not fitted physical constants that enter the scientific claim.

free parameters (4)
  • GS-ST temperature τ
    Controls bias-variance trade-off of the continuous relaxation; chosen by hand (0.02–0.75) and shown to determine success or failure of inference.
  • temporal soft-cutoff τ_time
    Softens the Heaviside that selects reactions before t_final; set to 2.5e-5 or 0.05 depending on the experiment.
  • SGD learning rate α
    Fixed at 0.1 for all repressilator runs; step-size choice affects convergence but is not derived.
  • number of trajectories per gradient estimate
    100–1000 trajectories; batch size directly sets reported variance and is chosen for computational convenience.
axioms (4)
  • domain assumption Gillespie SSA exactly samples the continuous-time Markov chain defined by the chemical master equation
    Standard foundation of all trajectory generation (Section II).
  • standard math Gumbel-max / softmax reparameterization yields a valid (biased) gradient estimator for categorical variables
    Taken from Jang et al. and Maddison et al.; used without re-derivation.
  • standard math Score-function (REINFORCE) estimator is unbiased for any discrete distribution
    Classic result of Williams; baseline subtraction preserves unbiasedness.
  • standard math Alternative-path estimator of Arya et al. is unbiased for discrete randomness
    Adopted from the cited NeurIPS 2022 paper.

pith-pipeline@v1.1.0-grok45 · 29071 in / 2416 out tokens · 28606 ms · 2026-07-13T13:57:40.377252+00:00 · methodology

0 comments
read the original abstract

Stochastic kinetic models are ubiquitous in physics, yet inferring their parameters from experimental data remains challenging. For deterministic models, parameter inference often relies on gradients, which can be obtained efficiently through automatic differentiation (AD). However, AD cannot be applied directly to the Gillespie stochastic simulation algorithm (SSA), since sampling from a discrete set of reactions introduces non-differentiable operations. In this work, we adopt three gradient estimators from machine learning for the Gillespie SSA: the Gumbel-Softmax Straight-Through (GS-ST) estimator, the Score Function estimator, and the Alternative Path estimator. We use the estimators to evaluate gradients of steady-state and time-dependent observables, and compare their performance in representative biophysical systems with relaxation dynamics (bimolecular association) and oscillatory dynamics (repressilator). We find that the GS-ST estimator generally yields well-behaved gradient estimates, but exhibits diverging variance in challenging parameter regimes, which can cause parameter inference to fail. In these cases, other estimators provide more robust, lower variance gradients. Our results demonstrate that gradient-based parameter inference can be effectively combined with the Gillespie SSA, with different estimators offering complementary advantages.

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

54 extracted references · 3 linked inside Pith

  1. [1]

    Brehmer, F

    J. Brehmer, F. Kling, I. Espejo, and K. Cranmer, Mad- Miner: Machine learning-based inference for particle physics, Comput. Softw. Big Sci.4, 3 (2020)

    2020

  2. [2]

    A. I. Khan, M. M. Billah, C. Ying, J. Liu, and P. Dutta, Bayesian method for parameter estimation in transient heat transfer problem, Int. J. Heat Mass Transf.166, 120746 (2021)

    2021

  3. [3]

    Carteret al., A benchmark JWST near-infrared spec- trum for the exoplanet WASP-39 b, Nat

    A. Carteret al., A benchmark JWST near-infrared spec- trum for the exoplanet WASP-39 b, Nat. Astron.8, 1008 (2024)

    2024

  4. [4]

    Kofler, M

    A. Kofler, M. Dax, S. R. Green, J. Wildberger, N. Gupte, J. H. Macke, J. Gair, A. Buonanno, and B. Sch¨ olkopf, Flexible gravitational-wave parameter estimation with transformers (2025), arXiv:2512.02968

    Pith/arXiv arXiv 2025

  5. [5]

    Gavrikov, others, and JUNO-Italia Consortium, Simulation-based inference for precision neutrino physics through neural Monte Carlo tuning, Commun

    A. Gavrikov, others, and JUNO-Italia Consortium, Simulation-based inference for precision neutrino physics through neural Monte Carlo tuning, Commun. Phys.9, 63 (2026)

    2026

  6. [6]

    Mitra and W

    E. Mitra and W. Hlavacek, Parameter estimation and un- certainty quantification for systems biology models, Curr. Opin. Syst. Biol.18, 9 (2020)

    2020

  7. [7]

    Fr¨ ohlich, B

    F. Fr¨ ohlich, B. Kaltenbacher, F. Theis, and J. Hasenauer, Scalable parameter estimation for genome-scale biochem- ical reaction networks, PLoS Comp. Biol.13, e1005331 (2017)

    2017

  8. [8]

    Stapor, F

    P. Stapor, F. Fr¨ ohlich, and J. Hasenauer, Optimization and profile calculation of ODE models using second or- der adjoint sensitivity analysis, Bioinformatics34, i151 (2018)

    2018

  9. [9]

    Frank, Automatic differentiation and the optimization of differential equation models in biology, Front

    S. Frank, Automatic differentiation and the optimization of differential equation models in biology, Front. Ecol. Evol.10, 1010278 (2022)

    2022

  10. [10]

    Arkin, J

    A. Arkin, J. Ross, and H. McAdams, Stochastic ki- netic analysis of developmental pathway bifurcation in phageλ-infected Escherichia coli cells, Genetics149, 1633 (1998)

    1998

  11. [11]

    Schnoerr, G

    D. Schnoerr, G. Sanguinetti, and R. Grima, Approxima- tion and inference methods for stochastic biochemical ki- netics – a tutorial review, J. Phys. A: Math. Theor.50, 093001 (2017)

    2017

  12. [12]

    Gillespie, Exact stochastic simulation of coupled chemical reactions, J

    D. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem.81, 2340 (1977)

    1977

  13. [13]

    Komorowski, M

    M. Komorowski, M. Costa, D. Rand, and M. Stumpf, Sensitivity, robustness, and identifiability in stochastic chemical kinetics models, Proc. Natl. Acad. Sci. USA 108, 8645 (2011)

    2011

  14. [14]

    Ruess, H

    J. Ruess, H. Koeppl, and C. Zechner, Sensitivity estima- tion for stochastic models of biochemical reaction net- works in the presence of extrinsic variability, J. Chem. Phys.146, 124122 (2017)

    2017

  15. [15]

    Zechner, J

    C. Zechner, J. Ruess, P. Krenn, S. Pelet, M. Peter, J. Lygeros, and H. Koeppl, Moment-based inference pre- dicts bimodality in transient gene expression, Proc. Natl. Acad. Sci. USA109, 8340 (2012)

    2012

  16. [16]

    Schnoerr, G

    D. Schnoerr, G. Sanguinetti, and R. Grima, Comparison of different moment-closure approximations for stochas- tic chemical kinetics, J. Chem. Phys.143, 185101 (2015)

    2015

  17. [17]

    L¨ uck and V

    A. L¨ uck and V. Wolf, Generalized method of moments for estimating parameters of stochastic reaction networks, BMC Syst. Biol.10, 98 (2016)

    2016

  18. [18]

    Fr¨ ohlich, P

    F. Fr¨ ohlich, P. Thomas, A. Kazeroonian, F. Theis, R. Grima, and J. Hasenauer, Inference for stochastic chemical kinetics using moment equations and system size expansion, PLoS Comp. Biol.12, e1005030 (2016)

    2016

  19. [19]

    Komorowski, B

    M. Komorowski, B. Finkenst¨ adt, C. Harper, and D. Rand, Bayesian inference of biochemical kinetic pa- rameters using the linear noise approximation, BMC Bioinformatics10, 343 (2009)

    2009

  20. [20]

    Golightly and D

    A. Golightly and D. Wilkinson, Bayesian parameter in- ference for stochastic biochemical network models using particle markov chain monte carlo, Interface Focus1, 807 (2011)

    2011

  21. [21]

    T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. Stumpf, Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems, J. R. Soc. Interface6, 187 (2009)

    2009

  22. [22]

    Warne, R

    D. Warne, R. Baker, and M. Simpson, Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art, J. R. Soc. Interface.16, 20180943 (2019)

    2019

  23. [23]

    Mohamed, M

    S. Mohamed, M. Rosca, M. Figurnov, and A. Mnih, Monte carlo gradient estimation in machine learning, J. Mach. Learn. Res.21, 1 (2020)

    2020

  24. [24]

    E. Jang, S. Gu, and B. Poole, Categorical reparameteri- zation with gumbel-softmax, inInternational Conference on Learning Representations(2017)

    2017

  25. [25]

    C. J. Maddison, A. Mnih, and Y. W. Teh, The concrete distribution: A continuous relaxation of discrete random variables, inInternational Conference on Learning Rep- resentations(2017)

    2017

  26. [26]

    Williams, Simple statistical gradient-following algo- rithms for connectionist reinforcement learning, Mach

    R. Williams, Simple statistical gradient-following algo- rithms for connectionist reinforcement learning, Mach. Learn.8, 229 (1992)

    1992

  27. [27]

    G. Arya, M. Schauer, F. Sch¨ afer, and C. V. Rackauckas, Automatic differentiation of programs with discrete ran- domness, inConference on Neural Information Process- ing Systems(2022)

    2022

  28. [28]

    Kagan and L

    M. Kagan and L. Heinrich, Branches of a tree: Taking derivatives of programs with discrete and branching ran- domness in high energy physics (2023), arXiv:2308.16680

    Pith/arXiv arXiv 2023

  29. [29]

    Heinrich and T

    L. Heinrich and T. Magorsch, Differentiable quantum- trajectory simulation of Lindblad dynamics for QGP transport-coefficient inference (2026), arXiv:2601.14399

    arXiv 2026

  30. [30]

    G. Arya, R. Seyer, F. Sch¨ afer, K. Chandra, A. K. Lew, M. Huot, V. Mansinghka, J. Ragan-Kelley, C. V. Rackauckas, and M. Schauer, Differentiating metropolis- hastings to optimize intractable densities, inICML Work- shop on Differentiable Almost Everything(2023)

    2023

  31. [31]

    Elowitz and S

    M. Elowitz and S. Leibler, A synthetic oscillatory net- work of transcriptional regulators, Nature403, 335 (2000)

    2000

  32. [32]

    Loinger and O

    A. Loinger and O. Biham, Stochastic simulations of the repressilator circuit, Phys. Rev. E76, 051917 (2007)

    2007

  33. [33]

    Potvin-Trottier, N

    L. Potvin-Trottier, N. Lord, G. Vinnicombe, and J. Paulsson, Synchronous long-term oscillations in a syn- thetic gene circuit, Nature538, 514 (2016)

    2016

  34. [34]

    Mottes, Q.-Z

    F. Mottes, Q.-Z. Zhu, and M. Brenner, Gradient-based optimization of exact stochastic kinetic models (2026), arXiv:2601.14183

    arXiv 2026

  35. [35]

    J. M. G. Vilar and L. Saiz, Exact discrete stochastic simulation with deep-learning-scale gradient optimiza- tion (2026), arXiv:2602.19775. 19

    arXiv 2026

  36. [36]

    Rijal and P

    K. Rijal and P. Mehta, A differentiable gillespie algo- rithm for simulating chemical kinetics, parameter esti- mation, and designing synthetic biological circuits, eLife 14, RP103877 (2025)

    2025

  37. [37]

    Rao and M

    R. Rao and M. Esposito, Nonequilibrium thermodynam- ics of chemical reaction networks: Wisdom from stochas- tic thermodynamics, Phys. Rev. X6, 041064 (2016)

    2016

  38. [38]

    D. F. Anderson, Lecture notes on stochastic processes with applications in biology (2017)

    2017

  39. [39]

    E. J. Gumbel,Statistical Theory of Extreme Values and Some Practical Applications: A Series of Lectures(U.S. Government Printing Office, 1954)

    1954

  40. [40]

    Harris and S

    D. Harris and S. Harris,Digital Design and Computer Architecture, 2nd ed. (Elsevier, 2012)

    2012

  41. [41]

    Zheng and A

    A. Zheng and A. Casari,Feature Engineering for Machine Learning: Principles and Techniques for Data Scientists, 1st ed. (O’Reilly Media, Inc., Sebastopol, CA, 2018)

    2018

  42. [42]

    D. J. Rezende, S. Mohamed, and D. Wierstra, Stochastic backpropagation and approximate inference in deep gen- erative models, inInternational Conference on Machine Learning(2014)

    2014

  43. [43]

    Kingma and M

    D. Kingma and M. Welling, Auto-encoding variational bayes, inInternational Conference on Learning Repre- sentations(2014)

    2014

  44. [44]

    Schulman, N

    J. Schulman, N. Heess, T. Weber, and P. Abbeel, Gra- dient estimation using stochastic computation graphs (2016), arXiv:1506.05254

    Pith/arXiv arXiv 2016

  45. [45]

    Rubner, C

    Y. Rubner, C. Tomasi, and L. Guibas, The Earth Mover’s Distance as a metric for image retrieval, Int. J. Comput. Vis.40, 99 (2000)

    2000

  46. [46]

    M. B. Paulus, C. J. Maddison, and A. Krause, Rao- Blackwellizing the Straight-Through Gumbel-Softmax gradient estimator, inInternational Conference on Learning Representations(2021)

    2021

  47. [47]

    Greensmith, P

    E. Greensmith, P. L. Bartlett, and J. Baxter, Variance reduction techniques for gradient estimates in reinforce- ment learning, J. Mach. Learn. Res.5, 1471–1530 (2004)

    2004

  48. [48]

    Rosenberger, T

    J. Rosenberger, T. G¨ oppel, P. Kudella, D. Braund, U. Gerland, and B. Altaner, Self-assembly of informa- tional polymers by templated ligation, Phys. Rev. X11, 031055 (2021)

    2021

  49. [49]

    Burger and U

    L. Burger and U. Gerland, Toward stable replication of genomic information in pools of RNA molecules, eLife 14, RP104043 (2025)

    2025

  50. [50]

    Harth-Kitzerow, T

    J. Harth-Kitzerow, T. G¨ oppel, L. Burger, T. Enßlin, and U. Gerland, Sequence motif dynamics in RNA pools, Phys. Rev. E113, 024407 (2026)

    2026

  51. [51]

    Duane, A

    S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, Hybrid monte carlo, Phys. Lett. B195, 216 (1987)

    1987

  52. [52]

    Neal, MCMC using hamiltonian dynamics, inHand- book of Markov Chain Monte Carlo(Chapman and Hall/CRC, 2011)

    R. Neal, MCMC using hamiltonian dynamics, inHand- book of Markov Chain Monte Carlo(Chapman and Hall/CRC, 2011)

    2011

  53. [53]

    Bradbury, R

    J. Bradbury, R. Frostig, P. Hawkins, M. J. John- son, Y. Katariya, C. Leary, D. Maclaurin, G. Nec- ula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2026)

    2026

  54. [54]

    Bois and M

    J. Bois and M. Elowitz, The repressilator enables self- sustaining oscillations (2020), Lecture notes ”Design Principles of Genetic Circuits”, California Institute of Technology

    2020