Probability of one species reaching zero before the other in a Markov process on a 2d lattice

Question

$\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}$:

1. 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.

2. 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.

3. 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!

Answer 1

The probability $u(x_1,x_2)$ that the process hits $A_1$ before it hits $A_2$ starting at $(x_1, x_2)$ satisfies the boundary value problem $$ L u(x_1,x_2) = 0 \;, \quad u(x_1,0)=1 \;, \quad u(0,x_2)=0 \;, \quad (x_1,x_2) \in \mathbb{N}^2 \tag{$\star$} $$ where $L$ is the infinitesimal generator of the process defined as $$ L f(x_1,x_2) = \frac{x_1}{x_1+x_2} \Big( f(x_1+1,x_2) - 2 f(x_1, x_2) + f(x_1-1,x_2) \Big) \\ + \frac{x_2}{x_1+x_2} \Big( f(x_1,x_2+1) - 2 f(x_1, x_2) + f(x_1,x_2-1) \Big) \;. $$ It can be easily verified that the solution to ($\star$) is $u(x_1, x_2) = x_1 / (x_1 + x_2)$.

Answer 2

Here's an outline of an intuitive way to think about the answer (already given by Nawaf Bou-Rabee).

You can change the transition rates to $x_1\to x_1\pm 1$ at rate $x_1$ each, and $x_2\to x_2\pm 1$ at rate $x_2$ each. This doesn't change the distribution of the jump chain (it just rescales time).

Then you can interpret the system as a collection of $x_1$ type-$1$ individuals, and $x_2$ type-$2$ individuals. Each individual in the population gives birth (to another individual of the same type as itself) at rate $1$, and also dies at rate $1$.

Eventually the whole population dies out. Your question amounts to asking the probability that the last individual alive is of type $1$.

This last individual alive is a descendant of one of the time-$0$ individuals (possibly in the trivial sense of being one of the time-$0$ individuals itself). By symmetry, it is equally likely to be a descendant of any of the time-$0$ individuals.

So if the system starts at time $0$ with $x_1(0)$ type-$1$ individuals and $x_2(0)$ type-$2$ individuals, the probability that the last individual alive is of type $1$ is $x_1(0)/(x_1(0)+x_2(0))$.