Minimum Action Path theory reveals the details of stochastic biochemical transitions out of oscillatory cellular states

TL;DR

This paper introduces a theoretical framework using Minimum Action Path theory to analyze stochastic transitions from oscillatory cellular states, revealing detailed escape paths and timing distributions in biochemical systems.

Contribution

It develops a method to identify the most probable escape paths and times from oscillatory states in biochemical networks using Freidlin-Wentzell action minimization.

Findings

Derived the stochastic spiral exit path from limit cycles.

Inferred escape time probability density functions.

Provided a detailed characterization of noise-driven state transitions.

Abstract

Cell state determination is the outcome of intrinsically stochastic biochemical reactions. Tran- sitions between such states are studied as noise-driven escape problems in the chemical species space. Escape can occur via multiple possible multidimensional paths, with probabilities depending non-locally on the noise. Here we characterize the escape from an oscillatory biochemical state by minimizing the Freidlin-Wentzell action, deriving from it the stochastic spiral exit path from the limit cycle. We also use the minimized action to infer the escape time probability density function.

Equations

$$\dot{r}=r(r^2-c^2)E(x,y),\dot{\theta }=\omega ,$$

$$\dot{x}={\rho }_x^{+}(x,y)-{\rho }_x^{-}(x,y),\dot{y}={\rho }_y^{+}(x,y)-{\rho }_y^{-}(x,y),$$

$${\rho }_x^{+}$$

$${\rho }_x^{-}$$

$${\rho }_y^{+}$$

$${\rho }_y^{-}$$

$$\dot{x}$$

$$\dot{y}$$

$$\mathcal{S}\equiv \mathcal{S}({\phi }^{\ast })=\underset{\phi }{\min }\mathcal{S}(\phi ).$$

$$\mathcal{S}({\phi }_{\tau })=\frac{1}{2}{\int }_0^{\tau }{\Vert {\dot{\phi }}_{\tau }(t)-f({\phi }_{\tau }(t))\Vert }_{g({\phi }_{\tau }(t))}^2\mathrm{d}t,$$

$$D(x,y)=\left(\begin{array}[]{cc}\rho_{x}^{+}+\rho_{x}^{-}&0\\ 0&\rho_{y}^{+}+\rho_{y}^{-}\end{array}\right).$$

$$T\simeq C{\mathrm{e}}^{\Omega \mathcal{S}}.$$

$$\mathcal{P}(t_n)=p(1-p{)}^{\frac{\omega t_n}{2\pi }-1}.$$

$$\mathcal{P}(t_n)=\frac{{\left(1-{\mathrm{e}}^{-\mathcal{S}\Omega }/C\right)}^{\frac{\omega t_n}{2\pi }}}{C{\mathrm{e}}^{\mathcal{S}\Omega }-1}.$$

$$\mathcal{P}(t)=\frac{1}{C}{\mathrm{e}}^{-\mathcal{S}\Omega +\frac{t}{C}\exp (-\mathcal{S}\Omega )}.$$

Full text

Introduction.

Cells are intrinsically noisy. Such stochasticity arises not only from the production and degradation of cellular components, but also from their mutual interaction or even the interaction with other cells. Nevertheless, some cellular processes require a precise deterministic output, and noise-suppression mechanisms are necessary within the cell Kepler and Elston (2001); Kaern et al. (2005); Maheshri and O’Shea (2007); Losick and Desplan (2008); Raj and van Oudenaarden (2008). On the other hand, since fluctuations are an intrinsic component of cellular dynamics, mechanisms are in place that cells exploit to improve its function Cai et al. (2008); Eldar and Elowitz (2010). For example, randomness can enhance the ability of cells to adapt and increase their fitness in random variable environments Kussell and Leibler (2005); Acar et al. (2008); Guerrero et al. (2015), or to sustain phenotypic variation Maheshri and O’Shea (2007); Losick and Desplan (2008); Raj and van Oudenaarden (2008); MacArthur et al. (2008); Balazsi et al. (2011).

Mean-field descriptions of biochemical processes can be analyzed using dynamical systems theory DiStefano (2015), where stable steady states, sustained oscillations or even transients of the ODEs correspond to different possible ceullar states Jaeger and Monk (2014); Verd et al. (2014). Relevant examples including sustained oscillations are circadian rhythms Goldbeter (1995); Leloup and Goldbeter (2000); Smolen et al. (2001), cAMP oscillations in Dictyostelium Goldbeter (1996), cell-cycle regulation Gerard and Goldbeter (2009, 2011, 2012), or patterns of bursting in neuronal activity Tsumoto et al. (2006); Ditlevsen and Greenwood (2013); Bodova et al. (2015); Morris and Lecar (1981); Tateno and Pakdaman (2004).

When molecular populations are small, the mean-field framework is inaccurate and a stochastic description is required. This involves the formulation of the Master Equation (ME) describing the underlying multivariate biochemical birth-death process Anderson and Kurtz (2015). Unfortunately the ME is rarely solvable analytically, requiring the use of Monte Carlo methods (such as the Gillespie algorithm Gillespie (1977)). These numerical methods are often computationally costly and in usually infeasible Cao et al. (2005a, b). This is especially true in phenomena associated with rare fluctuations, as in the noise-induced escape from a basin of attraction. In these escape problems, approximations such as the Langevin equation van Kampen (2007); Gillespie (2000) or extreme event theory become necessary Touchette (2009); Perez-Carrasco et al. (2016). In spite of the importance of oscillatory phenomena in biology, most studies have tackled escape problems from point attractors, and a general theory of escape from stable limit cycles is lacking. In order to fill this gap, we consider a simplified oscillatory kinetic model and unveil the ability of the Minimum Action Path (MAP) method from large deviation theory Freidlin and Wentzell (1998); Weinan et al. (2004), which specifies the most probable path between attractors and the mean escape time by minimization of an action functional.

The model.

In order to study the limit cycle-fixed point transition we construct a tunable dynamical landscape and then derive the underlying kinetic reactions. This approach allows for a thorough analysis of the escape problem when changing key parameters, such as the angular velocity along the stable limit cycle or the distance between the fixed point and the stable limit cycle.

We construct this prototypical two-dimensional dynamical system, for species X and Y, such that there is a single stable fixed point at $(x_{c},y_{c})$, and a stable limit cycle of radius $c$ centered at the fixed point (see Fig.1). In order to determine the basins of attraction, we include an elliptical repulser, $E(x,y)\equiv 1-\frac{(x-x_{c}-x_{0})^{2}}{a^{2}}-\frac{(y-y_{c})^{2}}{b^{2}}=0$, separating the stable limit cycle from the fixed point. The evolution of the system is given by:

[TABLE]

where $r\equiv\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$ is the distance from the fixed point, and $\theta$ is the corresponding angular coordinate, $\tan\theta\equiv(y-y_{c})/(x-x_{c})$, and $\omega$ is the angular velocity. Trajectories are either attracted to the fixed point or to the limit cycle, with the exception of points lying on an unstable limit cycle, which remain on that orbit. It is important to note that the unstable limit cycle is not the same as the repulser.

Stochastic description.

To formulate a stochastic system whose mean-field behavior is given by Eqs. (1), we must derive a compatible set of biochemical reactions. Considering four reactions, each species being produced and degraded, the reaction rates are obtained by splitting the rhs of Eq. (1) into positive and negative contributions:

[TABLE]

where:

[TABLE]

Identifying $\rho_{x}^{\pm}$ and $\rho_{y}^{\pm}$ with the rates of the biochemical reactions, the deterministic system Eq. (2) corresponds to the macroscopic limit of the kinetic reaction set (see Fig. 1). Here, the system size, $\Omega$, relates the concentrations $x$ and $y$ with the numbers of molecules of each species, $X=x\Omega$ and $Y=y\Omega$. Reactions detailed in Fig. 1 describe a multivariate birth-death process that can be solved numerically using the Gillespie algorithm Gillespie (1977). As $\Omega$ grows, the intrinsic noise is reduced, recovering the mean-field limit (1) when $\Omega\rightarrow\infty$. For large (but finite) $\Omega$, the Master Equation can be approximated by the Chemical Langevin Equation (CLE) Gillespie (2000),

[TABLE]

where $\xi_{x}(t)$ and $\xi_{y}(t)$ are uncorrelated white Gaussian noises, of zero mean, and autocorrelation $\langle\xi_{x}(t)\xi_{x}(t^{\prime})\rangle=\langle\xi_{y}(t)\xi_{y}(t^{\prime})\rangle=\delta(t-t^{\prime})$. Within this formulation $\Omega$ contributes only to the stochastic terms of the CLE (4). Therefore tuning the value of $\Omega$ allows us to investigate the role of fluctuations in the transition between the stable limit cycle and the fixed point.

Minimum Action Path.

The intrinsic noise described in the previous section allows for transitions between the limit cycle and the fixed point. Such transitions can occur through many possible transient trajectories, $\varphi(x(t),y(t))$. Nevertheless, not all the transitions are equally probable. In particular, for reaction systems, unlikely transitions decay exponentially with $\Omega$, $\mathcal{P}\sim\mathrm{e}^{-\Omega\mathscr{S}(\varphi)}$ van Kampen (2007); Touchette (2009). Where the decay rate $\mathscr{S}(\varphi)$ is the so-called action of the transition. This means that for large enough values of $\Omega$, the stochastic transition will concentrate along the path, $\varphi^{*}$, which minimizes the action:

[TABLE]

For the $n$-dimensional stochastic differential equation $\dot{\varphi}=f(\varphi)+g(\varphi)\Omega^{-\frac{1}{2}}\xi(t)$, the action for any path $\varphi_{\tau}$ of duration $\tau$ is given by the Freidlin-Wentzell functional Freidlin and Wentzell (1998):

[TABLE]

where $f(\varphi_{\tau})$ is the deterministic field describing the dynamical system, given for our system by the rhs of Eq. (2). The multiplicative noise appears in the norm $\|\bullet\|_{g(\varphi_{\tau})}^{2}$, corresponding with the inner product $\left\langle\bullet,\left(g(\varphi_{\tau})g(\varphi_{\tau})^{\top})\right)^{-1}\bullet\right\rangle$, where $g(\varphi_{\tau})g^{\top}(\varphi_{\tau})\equiv D$ is the diffusion matrix. Here $D$ takes the form:

[TABLE]

Interestingly, the action and, consequently, the most probable path, are independent of $\Omega$. Additionally, the mean first passage time (MFPT) $T$ from one attractor to the other can be expressed as Freidlin and Wentzell (1998); Touchette (2009):

[TABLE]

In order to find numerically the path minimizing $\mathscr{S}(\varphi_{\tau})$, each path of duration $\tau$ was divided into a chain of $N$ segments with initial and final points in the relevant attractors. This reduces finding the optimal path to a minimization problem with $2(N-2)$ degrees of freedom. This was solved using the Broyden-Fletcher-Goldfarb-Shanno algorithm Fletcher (2000); GSL (2009), using the analytical expression for the gradient of the action in any of the $2(N-2)$ dimensions Perez-Carrasco et al. (2016).

Results.

To assess whether MAP theory can characterize the escape from a stable limit cycle, we have divided the analysis into two sections. First, we compare the MAP with paths obtained numerically from the Master Equation and the CLE. In the second section, we compare MAP predictions of the MFPT with those derived from CLE numerical solutions.

The MAP predicts average stochastic escape trajectories.

The MAP defines the most probable transient molecular concentrations during the escape from the stable limit cycle at low noise. Direct comparison of the MAP with trajectories obtained from numerical integration of the CLE or Gillespie simulations shows good agreement for $\Omega\geq 150$ (Fig. 2). This reveals that, as $\Omega$ increases, the stochastic escape trajectories converge to the MAP. Stochastic simulations for $\Omega<100$ reveal that, when the number of molecular species is low, oscillations have poor quality, and escape trajectories do not concentrate around a single path (data not shown). In addition, our simulations show that the MAP recapitulates changes in escape trajectory with the angular velocity $\omega$ (Fig. 2).

A more detailed comparison between stochastic simulations and MAP theory reveals that the prediction of the latter becomes less accurate close to the exit point from the cycle. The discrepancy originates in a fraction of trajectories following the limit cycle for a bit longer before starting the transition (Fig. 2). This results in the prediction of a smaller exit angle than the actual average exit angle (Fig. 3). To gain deeper understanding into the origin and magnitude of the discrepancy, we computed the angular probability distribution along the cycle, and the probability distribution of the escape angle from a thin annulus around the limit cycle (Fig. 3). The latter has been proposed in Hitczenko and Medvedev (2013) as a quantity that characterizes escape from a stable limit cycle in the low noise limit. Strikingly, our analysis shows that neither measure is as informative as the MAP regarding the escape angle. These results suggest that, even for a simple dynamical system, knowledge of the whole dynamical landscape is required to predict the exit angle from a stable limit cycle, since the purely local analysis around the stable limit cycle does not produce accurate predictions. In this respect, the MAP proves to be useful, since action minimization takes place along the whole escape trajectory.

Localized inaccuracies in the MAP prediction suggest a highly heterogeneous contribution to the action along the MAP. In order to study this, we have evaluated the density of the action along the MAP, i.e. the Lagrangian of the system. Results show that the density is highest in the middle of the MAP (Fig. 3), becoming negligible close to the stable and unstable limit cycles, where, in addition, the MAP is tangent to both limit cycles. This leads to the discrepancy observed in the exit angles. Note that the portion of the MAP inside the basin of attraction of the stable node does not contribute to the action functional since it corresponds with the deterministic trajectory ($\dot{\varphi}=f(\varphi(t))$).

In usual escape problems, the path crosses from one basin of attraction to the other at the saddle point of the deterministic system. Here the boundary between basins of attraction is the unstable limit cycle, so the crossing point cannot be identified by a simple local stability analysis. This again shows the predictive power of the MAP approach.

Minimum action theory predicts MFPT for escape from the cycle.

Besides the optimal path, we are also interested in testing the ability of MAP theory to predict the MFPT to exit the basin of attraction. Eq. (8) shows that this can be achieved to logarithmic precision, up to a constant, $C$. When the basins of attraction are separated by a saddle point, $C$ can be determined by a Jacobian computed at the saddle Bouchet and Reygner (2016). However, in the current case, the separatrix is an unstable limit cycle and $C$ must be computed numerically by solving the CLE at low $\Omega$. It can then be used to predict MFPT for larger $\Omega$, where numerical integration of the CLE is computationally costly. Our results show that the minimum action theory allows us to capture the MFPT dependence on model parameters (Fig. 4). In our model, we observe an increase in the MFPT with $\omega$. In fact, $C$ also depends non-monotonically on $\omega$ (see Fig. 4). Nevertheless, as $\Omega$ grows, the contribution of the prefactor becomes less important ($\ln T\approx\Omega\mathcal{S}+\ln C$), and the minimum action dominates the escape time estimate.

In addition to the MFPT, we are interested in finding the probability distribution of escape times from the stable limit cycle. Assuming that escape is a rare event focused around a certain exit angle, the escape problem can be described as a Bernoulli process with low success probability $p$ taking place every period of the cycle $\tau=2\pi/\omega$, at times $t_{n}=2\pi n/\omega$. The probability of exiting at the $n-$th revolution follows the geometric distribution

[TABLE]

Using rare event theory, we can write the success probability as $p=\mathrm{e^{-\mathcal{S}\Omega}/C}$. The distribution Eq. (9) becomes

[TABLE]

Interestingly, since escape is rare, the probability $p$ will be very small and there will usually be many revolutions before the exit from the limit cycle occurs. In the limit $p\rightarrow 0$, the discrete geometric distribution (10) can be approximated by its continuum counterpart, the exponential distribution, which does not depend explicitly on the angular velocity,

[TABLE]

Comparing the distributions Eqs. (10) and (11), with the probability distribution of MFPT obtained over several CLE realizations, we obtained a good agreement (see Fig. 5). Surprisingly, even for realizations with a low average number of revolutions prior to escape, the resulting probability distribution is more similar to an exponential distribution than to a geometric one. This is true even for escapes that occur during the first revolution, suggesting that $\theta$ differs significantly from $\omega t$. A more accurate prediction would involve a convolution of geometric processes with the angular noise Berglund and Gentz (2014); Berglund (2016). However, for the parameters we used, the exponential distribution fits well independently of the average number of revolutions. Fitting the distribution (11) to the MPFTs from CLE realizations therefore provides an alternative method to compute the prefactor $C$.

Conclusions and perspectives.

We have shown that, within the rare event theory framework, escape problems from a stable limit cycle can be accurately characterized. The success of the method, in comparison with alternatives, relies on the fact that the Freidlin-Wentzell action is not a local property of the dynamical landscape but of the whole escape trajectory. For sufficiently large system sizes, we have shown that MAP theory accurately predicts escape trajectories and the escape time distribution. The method has also revealed properties of the escape trajectory, such as the tangent exit of the MAP from the stable limit cycle and tangent entry into the basin of attraction of the stable fixed point, as well as the dependence of the entry and exit points on the parameters of the system.

MAP predictions, whilst better than previously used methods, were less successful in determining the exit angle from the stable limit cycle. A deeper analysis of the exit angle could be carried out through a relaxation of the Laplace condition, i.e. the reduction to an integral along a single optimal path, by further exploring the distribution of suboptimal paths. Additionally, in the absence of a unique point ( the saddle) separating the basins of attraction, novel research into the calculation of the prefactor Bouchet and Reygner (2016) should be extended.

Eventually, the details of the dynamical landscape will be determined by concrete biological systems. For this reason, future work plans include the application of action minimization to stochastic excitable transitions in type II neurons, where rare event theory will prove its predictive power.

Acknowledgements.

RdC acknowledges AGAUR-Generalitat de Catalunya for funding under its doctoral scholarship program. RdC and TA are supported by the CERCA Program of the Generalitat de Catalunya, MINECO (grant MTM2015-71509-C2-1-R), and AGAUR (grant 2014SGR1307). TA acknowledges support from MINECO for funding awarded to the Barcelona Graduate School of Mathematics under the “María de Maeztu” program (grant MDM-2014-0445). RP-C, PG and KMP were supported by the Wellcome Trust (grant WT098325MA).