$\textbf{Background}$: Say we've got a two-variable system of stochastic chemical reactions, with quantities $\vec{x}(t) = (x_1(t),x_2(t)) \in \mathbb{N}^2$ evolving according to the following system, with rates equal to their current frequencies:

$x_1 \xrightarrow{\frac{x_1}{x_1+x_2}} x_1+1$ , $\; \; \; \; x_1 \xrightarrow{\frac{x_1}{x_1+x_2}} x_1-1$,

$x_2 \xrightarrow{\frac{x_2}{x_1+x_2}} x_2+1$ , $\; \; \; \; x_2 \xrightarrow{\frac{x_2}{x_1+x_2}} x_2-1$

With probability $1$, the system will eventually reach the point $(0,0)$. Also, define $$ A_1 = \{(x_1,0) \; | \; x_1>0\}$$ and $$ A_2 = \{(0,x_2) \; | \; x_2>0\}$$ Then for a system with non-zero initial conditions $x_1(0),x_2(0)>0$, the event of hitting $A_1$ and the event of hitting $A_2$ are mutually exclusive. Furthermore, $P(\text{hit } A_1) = 1-P(\text{hit }A_2)$.

$\textbf{Question}$: Assuming $x_1(0)>0 $ and $x_2(0)>0$, what's the probability that $x_2$ will reach $0$ before $x_1$? In other words, what is the hitting probability of $A_1$?

$\textbf{Ideas and comments}$:

- Since we are considering probability, the timing of the events is irrelevant, so we can easily reformulate the question to be a Markov Chain with state space $S=\mathbb{N}^2$ and transition probabilities equivalent to the normalized rates. Then we can write out the system of equations for the hitting probabilities $h^{A_1}_{\vec{x}}$ as per standard theory.

$$ \begin{cases} h_{\vec{x}}^{A_1} = 1 \; \; \; \; \text{for} \; \; \; \vec{x} \in A_1 \\ h_{\vec{x}}^{A_1} = \sum_{\vec{y} \in S} P(\vec{x},\vec{y}) h_{\vec{y}}^{A_1} \; \; \; \; \text{for} \; \; \; \vec{x} \not \in A_1 \end{cases} $$ Unfortunately, this is an infinite system of equations, and I don't know how to go about attempting to solve it.

One could run a simulation of such a system via the Gillespie algorithm given an initial point $x(0) = (x_1(0),x_2(0))$, so for a specific starting point, one could do a large set of simulations and get a good approximation. But this will be very slow. I expect that there are more clever ways to numerically solve this.

A very similar question, with the same question but with rates equal to $1/4$ (and absorbing boundary conditions), is solved in the 1940 paper entitled "Random Paths in Two and Three Dimensions", by W.H. McCrea et al. I'm unclear but doubtful that the same type of method can be applied here.

Thank you for any insight!