Stochastic Kinetic Models: Gradient-Based Optimisation Accurately Recovers Rates

Stochastic kinetic models are crucial for understanding systems in biology, chemistry and physics where randomness and small numbers matter. Researchers Francesco Mottes, Qian-Ze Zhu, and Michael P. Brenner, all from the School of Engineering and Applied Sciences at Harvard University, have tackled a significant challenge in this field: optimising these models when traditional methods fail due to their inherent non-differentiability. Their new approach, utilising straight-through Gumbel-Softmax estimation, allows for efficient gradient-based optimisation whilst maintaining the accuracy of exact stochastic simulations, a breakthrough that enables systematic inference and design across a wide range of applications. The team demonstrate its power by accurately recovering kinetic rates in gene expression and characterising trade-offs in stochastic thermodynamics, paving the way for rational design in complex systems governed by continuous-time Markov dynamics.

Both can be written as deterministic functions of parameter-independent random draws, enabling gradients to flow through the simulation while preserving exact discrete dynamics (Fig0.1, right). For waiting times, this reparameterization is standard.

Kinetic rates recovered from gene expression data often

The team demonstrated this capability across 25 parameter sets, successfully recovering the true values of kon and ktx, parameters defining promoter switching and mRNA degradation, with high fidelity. Results demonstrate that optimized models accurately reproduce target moments, as evidenced by data presented in Figure 0.2b, and achieve Pearson correlations of 0.68-0.74 in less challenging regimes, surpassing previous automatic differentiation studies. Full results are detailed in Figure S1, confirming the method’s effectiveness even with ill-conditioned loss landscapes featuring flat ridges. Experiments revealed that despite the problem’s analytical minimum, the loss landscape presented significant optimization challenges due to low sensitivity directions.
Achieving full convergence necessitated deviating from standard hyperparameter settings, a testament to the complexity of the optimization process. The research team tackled this by reformulating the moment matching problem, effectively preconditioning the landscape by reparameterizing to fit burst frequency (konkoff/(kon + koff)) and mean burst size (ktx/ koff). However, they demonstrated that their method could directly solve the original ill-conditioned problem without requiring such reparameterization, highlighting its inherent robustness. Data shows that inference from full steady-state RNA copy-number distributions, containing substantially more information than moments alone, allowed for the unique determination of additional kinetic parameters.

With three free parameters (kon, ktx, kmdeg), the optimization problem admitted a unique minimum, though the loss landscape exhibited near-degeneracies where coupled changes in promoter switching and mRNA degradation produced nearly indistinguishable distributions. Scientists quantified distributional mismatch using the 1-Wasserstein (Earth Mover’s) distance, calculating the L1 distance between cumulative distribution functions, achieving a more robust metric than cross-entropy or KL divergence. Measurements confirm that accurately resolving the distribution, particularly in the tails, required a large number of simulations, but backpropagating through all of them would have exceeded memory constraints. To overcome this, the team computed distributional statistics from a large pool of forward-only simulations combined with a smaller set of gradient-tracked trajectories, decoupling sample size from memory cost and reducing gradient variance. Across the 25 synthetic parameter sets, the procedure accurately recovered ktx, while errors in kon and kmdeg exhibited anticorrelation, reflecting a near-degenerate direction in the loss landscape, a characteristic of “sloppy parameter” structure in kinetic models. Importantly, the optimized models faithfully reproduced the full target distribution over several orders of magnitude in probability, as shown in Figure 0.3d and detailed in Figures S5, S6, demonstrating the efficacy of ST-GS gradients in solving challenging distribution-level inverse problems.

Gumbel-Softmax for Stochastic System Optimisation offers a differentiable

Their method employs straight-through Gumbel-Softmax estimation, enabling differentiation through exact stochastic simulations by approximating gradients with a continuous relaxation applied solely during the backward pass. This technique maintains the accuracy of stochastic modelling while facilitating optimisation, a significant advancement for fields reliant on these complex systems. Researchers successfully demonstrated the robustness of their method in two key applications: inferring parameters in stochastic gene expression and performing inverse design in stochastic thermodynamics. Furthermore, it characterised Pareto-optimal trade-offs between non-equilibrium currents and entropy production, showcasing its utility in thermodynamic design.

The authors acknowledge a bias introduced by their gradient estimator, specifically a shrinkage factor dependent on the number of gradient-tracked versus baseline simulations. However, they demonstrate that this bias is mitigated in the limit of a large number of baseline simulations, and that adaptive optimisers can effectively absorb the effect. By efficiently differentiating through exact stochastic simulations, the method offers a powerful tool for understanding and manipulating complex biological, chemical, and physical processes. Future research could explore extending this approach to even more complex systems or investigating methods to further reduce the remaining bias in the gradient estimation, potentially enhancing the precision and efficiency of parameter inference and design.