Action functional gradient descent algorithm for estimating escape paths in stochastic chemical reaction networks

We first derive the Hamilton-Jacobi theory underlying continuous-time Markov processes, and then we use the construction to develop a variational algorithm for estimating escape (least improbable or first passage) paths for a generic stochastic chemical reaction network that exhibits multiple fixed points. The design of our algorithm is such that it is independent of the underlying dimensionality of the system, the discretization control parameters are updated toward the continuum limit, and there is an easy-to-calculate measure for the correctness of its solution. We consider several applications of the algorithm and verify them against computationally expensive means such as the shooting method and stochastic simulation. While we employ theoretical techniques from mathematical physics, numerical optimization and chemical reaction network theory, we hope that our work finds practical applications with an inter-disciplinary audience including chemists, biologists, optimal control theorists and game theorists.

Figure 1

Diagrammatic representation of 1-Schlogl, 2-Schlogl, 3-Schlögl, and 4-Schlögl model CRNs from top left to bottom right. Species, complexes, stoichiometry, and transitions are represented with solid circles, empty circles, solid lines, and dashed arrows, respectively.

Figure 2

Diagrammatic proof of the Lyapunov property of the nonequilibrium potential with respect to mass action kinetics.

Figure 3

Escape paths and heteroclinic network for the 2-Schlögl model.

Figure 4

Comparing the algorithm against Gillespie stochastic simulation for the 2-Schlögl model (Sec. 4b). Scaled log-improbability of the stationary distribution $\pi$ is defined to be $(-1/V)log\left(\pi \right)$.

Figure 5

Lifting a curve in configuration space to a trajectory in phase space by assigning a momentum and time difference to each segment.

Figure 6

Functional gradient descent for 2-Schlögl model without filtering. Note that the low-frequency waves first amplify and then concentrate at higher frequencies near a central, badly-performing place on the trajectory.

Figure 7

The AFGD algorithm in a nutshell. For the CRN Hamiltonian, when the end points are fixed points of the system then $h=0$ must be chosen. For a detailed description, see Sec. 3b.

Figure 8

Diagrammatic representation of the reaction network (left) and a proxy-heteroclinic network using relaxation trajectories for the Selkov model.

Figure 9

Comparing escape paths obtained from AFGD algorithm against the shooting method (left panel) and plotting the nonequilibrium potential along the escape (right panel). The data is taken from Ref. [32], and can be compared against Figs. 9 and 15 of the same reference. (For details on the differences between the two methods, see the last paragraph in Sec. 3c.)

Figure 10

Descent progress, summary and proof of convergence for the Selkov model after annealing to a trajectory with 4000 points. For initial descent with 2500 points, see Fig. 20.

Figure 11

$H(p,q)=0$ submanifold in the tangent and cotangent space for the 1D Schlögl model.

Figure 12

Fixed points labeled by repelling directions (left) and proxy-heteroclinic network using relaxation trajectories (right) for the 2-Schlögl model.

Figure 13

Extracting functional gradient from phase-space trajectory and smoothing it (top panel). Smoothed gradient and power spectrum at selected iterations (bottom panel) for the 2-Schlögl model.

Figure 14

Progress of integrated Hamilton's equations of motion against the AFGD algorithm at selected iterations for the 2-Schlögl model.

Figure 15

Evolving distribution by $\mathrm{\Delta }t.$

Figure 16

Initial conditions for the AFGD algorithm. (a) 2-Schlögl model and (b) Selkov model.

Figure 17

Phase space trajectory obtained by lifting the initial conditions shown in Fig. 16. (a) 2-Schlögl model initial condition phase space trajectory and (b) Selkov model initial condition phase space trajectory.

Figure 18

Extracting functional gradient from phase-space trajectory and smoothing it. In all panels, the unfiltered and filtered curves are represented by dashed and solid lines, respectively. (a) 2-Schlögl model original and filtered functional gradient at the first iteration and (b) Selkov model original and filtered functional gradient at the first iteration.

Figure 19

Descent progress and summary for the 2-Schlögl model for a trajectory with 500 points.

Figure 20

Descent progress and summary for the Selkov model for a trajectory with 2500 points.

Figure 21

Integrating Hamilton's equations of motion against AFGD algorithm output for the 6D Schlögl model after 300 iterations. The convention is the same as used in Fig. 14, and is omitted for neatness.