**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?



Obtaining generator matrix and first-passage time distribution for CTMC?
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**](http://dx.doi.org/10.1137/060666457), SIAM Rev. 50, No. 2, 347-368 (2008). [ZBL1144.80011](https://zbmath.org/?q=an:1144.80011). 


