Obtaining generator matrix and first-passage time distribution for CTMC?

Question

Setup:

I have a model of a biological process described by two ODEs as follows: $$\dot{X_1} = (\beta_1-d-1)X_1 + 2X_1^2 - X_1^3 + dX_2$$ $$\dot{X_2} = (\beta_2-d-1)X_2 + 2X_2^2 - X_2^3 + dX_1$$

I want to analyze the stochastic version of this system using an appropriate underlying mechanistic process. My choice of representation is a chemical reaction network as follows:

$$ X_1 \overset{\beta_1}{\rightharpoonup} 2X_1 $$ $$ X_1 \overset{1}{\rightharpoonup} \emptyset $$ $$ 2X_1 \overset{4}{\rightharpoonup} 3X_1 $$ $$ 3X_1 \overset{6}{\rightharpoonup} 2X_1 $$ $$ X_2 \overset{d}{\rightharpoonup} X_1 $$ $$ X_2 \overset{\beta_2}{\rightharpoonup} 2X_2 $$ $$ X_2 \overset{1}{\rightharpoonup} \emptyset $$ $$ 2X_2 \overset{4}{\rightharpoonup} 3X_2 $$ $$ 3X_2 \overset{6}{\rightharpoonup} 2X_2 $$ $$ X_1 \overset{d}{\rightharpoonup} X_2 $$

Following the procedure in Section 5.3.6 of Edward Allen's Modeling with Ito Stochastic Differential Equations, we can formulate a system of SDEs for the above model using the chemical reaction network. This allows for a noise vector that is derived from first principles, i.e. not tagged on in an ad-hoc manner to account for observed phenomenology.

I've been working with numerical simulations of this system for a while now. I've also surveyed a ton of literature for tools to derive analytical results. However, analytical progress is very slow (due to the cubic nonlinearities within a multi-dimensional system).

---

Questions:

1. Is there a way to obtain the infinitesimal generator matrix for the continuous-time Markov chain associated with this stochastic process? If so, how?

2. How can first-passage time distributions be obtained analytically, or via numerical estimates?

Answer 1

The infinitesimal generator $\mathscr{A}$ corresponding to the OP's chemical reaction network can defined by its action on a function $f: \mathbb{Z}^2_{\ge 0} \to \mathbb{R}$ as follows
$$ \mathscr{A}f(x) = \sum_{\ell} a_{\ell}(x) ( f(x+\nu_{\ell}) - f(x) ) $$ where we introduced the propensity functions and reaction channels defined respectively as $$ a_{\ell}(x) = \begin{cases} \beta_1 x_1 & \ell=1 \\ x_1 & \ell=2 \\ 2 x_1 (x_1-1) & \ell=3 \\ x_1 (x_1-1) (x_1-2) & \ell=4 \\ d x_2 & \ell=5 \\ \beta_2 x_2 & \ell=6 \\ x_2 & \ell=7 \\ 2 x_2 (x_2-1) & \ell=8 \\ x_2 (x_2-1) (x_2-2) & \ell=9 \\ d x_1 & \ell=10 \end{cases} \qquad \nu_{\ell} = \begin{cases} (1,0) & \ell=1 \\ (-1,0) & \ell=2 \\ (1,0) & \ell=3 \\ (-1,0) & \ell=4 \\ (1,-1) & \ell=5 \\ (0,1) & \ell=6 \\ (0,-1) & \ell=7 \\ (0,1) & \ell=8 \\ (0,-1) & \ell=9 \\ (-1,1) & \ell=10 \end{cases} $$ Here $x_1$ and $x_2$ are the first and second components of $x$, respectively. The propensities were constructed using the procedure detailed on page 352 of the following paper. Note that all of the reactions in the OP's network are either first-order, dimerizations, or trimerizations.

Higham, Desmond J., Modeling and simulating chemical reactions, SIAM Rev. 50, No. 2, 347-368 (2008). ZBL1144.80011.