Long-time analytic approximation of large stochastic oscillators: Simulation, analysis and inference

Abstract

In order to analyse large complex stochastic dynamical models such as those studied in systems biology there is currently a great need for both analytical tools and also algorithms for accurate and fast simulation and estimation. We present a new stochastic approximation of biological oscillators that addresses these needs. Our method, called phase-corrected LNA (pcLNA) overcomes the main limitations of the standard Linear Noise Approximation (LNA) to remain uniformly accurate for long times, still maintaining the speed and analytically tractability of the LNA. As part of this, we develop analytical expressions for key probability distributions and associated quantities, such as the Fisher Information Matrix and Kullback-Leibler divergence and we introduce a new approach to system-global sensitivity analysis. We also present algorithms for statistical inference and for long-term simulation of oscillating systems that are shown to be as accurate but much faster than leaping algorithms and algorithms for integration of diffusion equations. Stochastic versions of published models of the circadian clock and NF-κB system are used to illustrate our results.

TimeTeller: a New Tool for Precision Circadian Medicine and Cancer Prognosis

“…This stochastic model incorporates a system size parameter, Ω. We consider it has units L/nM and a value in the range 100-1000 based on typical cell size estimates (24). The individual trajectories were calculated in MATLAB, using an implementation of the Gillespie algorithm or the pc-LNA (24).…”

“…The 44 reaction rates for the stochastic Relogio model are shown in Note S12.1. In order to generate enough trajectories to estimate a distribution in a reasonable time the pc-LNA was used (24). = kd z4 P * C /C C , a(28) = kd z5 P C /C C , a(29) = kf z5 Cry C P er C , a(30) = kf z4 Cry C P er * C , a(31) = d z1 Cry C , a(32) = k p1 P er, a(33) = kd phz3 P er * C , a(34) = kph z2 P er C , a(35) = d z2 P er C , a(36) = d z3 P * C /C C , a(37) = d z4 P er * C , a(38) = d z5 P C /C C , a(39) = k p3 Reverb, a(40) = d z6 RV RB c , a(41) = k p4 Ror, a(42) = d z7 ROR c , a(43) = k p5 Bmal, a(44) = d z8 BM L c .…”

“…There are now several published algorithms which aim to estimate the time at which a transcriptomic 80 dataset was collected using the expression levels of the core clock genes (17,(23)(24)(25)(26)(27)(28). While these have 81 been used to study aspects of clock dysfunction by calculating a measure for a dataset as a whole and 82 comparing a diseased dataset with a healthy one (24,26), they are not designed to measure clock 83 dysfunction in single samples. The ideas from statistics and information theory mentioned above show 84 that the key quantities determining Θ are the covariance structure of and the way that moves 85 as time t advances (Note S4).…”

Multiplexing information flow through dynamic signalling systems

“…One can associate to the FIM I an s × s matrix s that satisfies I = s T s and certain optimality properties described in [17] (SI Sect. 2.5.2).…”

“…To do this we will introduce a quantity, called multiplexing capacity, that measures the ability of a noisy signalling system to multiplex a set of signals. The computation of this is underpinned by the pcLNA method [17] that allows fast and accurate computation of key information theoretic quantities, such as Kullback-Leibler divergences and the Fisher Information matrix [5], for stochastic dynamical systems.…”