Simulation-based inference using splitting schemes for partially observed diffusions in chemical reaction networks

We address the problem of simulation and parameter inference for chemical reaction networks described by the chemical Langevin equation, a stochastic differential equation (SDE) representation of the dynamics of the chemical species. This is challenging for two main reasons. First, the (multi-dimensional) SDEs cannot be explicitly solved and are driven by multiplicative and non-commutative noise, requiring the development of advanced numerical schemes for their approximation and simulation. Second, not all components of the SDEs are directly observed, as the available discrete-time data are typically incomplete and/or perturbed with measurement error. We tackle these issues via three contributions. First, we show that these models can be rewritten as perturbed conditionally Cox-Ingersoll-Ross-type SDEs, i.e., each coordinate, conditioned on all other coordinates being fixed, follows an SDE with linear drift and square root diffusion coefficient perturbed by additional Brownian motions. Second, for this class of SDEs, we develop a numerical splitting scheme that preserves structural properties of the model, such as oscillations, state space and invariant distributions, unlike the commonly used Euler-Maruyama scheme. Our numerical method is robust for large integration time steps. Third, we propose a sequential Monte Carlo approximate Bayesian computation algorithm incorporating "data-conditional" simulation and sequential learning of summary statistics, allowing inference for multidimensional partially observed systems, further developing previous results on fully observed systems based on the Euler-Maruyama scheme. We validate our approach on models of interest in chemical reaction networks, such as the stochastic Repressilator, Lotka-Volterra, and two-pool systems, demonstrating its effectiveness, in terms of both numerical and inferential accuracy, and reduced computational cost.