$\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.

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

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


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)$.

  • $\begingroup$ For a more general system with different rates to those above, once you have the boundary values and infinitesimal generator, is there a strategy to find the probability $u(x_1,x_2)$? It sounds like, in this case, it was an intuitive guess. $\endgroup$ – Samuel Reid Nov 22 '18 at 23:09
  • $\begingroup$ Also, I'm wondering if you have a recommendation for a reference on this topic! $\endgroup$ – Samuel Reid Nov 22 '18 at 23:11
  • $\begingroup$ @SamuelReid To my knowledge, there are no general strategies to solve second-order, non-constant coefficient linear difference equations on unbounded spaces. A (very) nice resource are Andreas Eberle's notes on Markov processes: wt.iam.uni-bonn.de/fileadmin/WT/Inhalt/people/Andreas_Eberle/… See, in particular, Section 1.1 entitled "Potential Theory for Markov Chains" $\endgroup$ – Nawaf Bou-Rabee Nov 23 '18 at 17:02

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))$.

  • $\begingroup$ That's a very nice explanation. $\endgroup$ – Nawaf Bou-Rabee Nov 20 '18 at 1:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.