Obtaining generator matrix and first-passage time distribution for CTMC? 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:


*

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

*How can first-passage time distributions be obtained analytically, or via numerical estimates?
 A: 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.  
