# Population exercise with Markov chains

I am totally stuck with this exercise on Markov chains. Maybe someone can help me :).

# Red and green bacteria

A growth medium at time $t = 0$ has 500 red bacteria and 500 green bacteria. Each hour, each bacterium divides in two. A colour-blind predator eats exactly 1000 bacteria per hour.

1. After a very long time, what is the probability distribution for the number $\alpha$ of red bacteria in the growth medium?
2. Roughly how long will it take to reach this final state?
3. Assume that the predator has a $1\%$ preference for green bacteria (implemented as you choose). Roughly how much will this change the final distribution?

Within the accuracy of this question, you may assume either that one bacterium reproduces and then one is eaten 1000 times per hour, or that at the end of each hour all the bacteria reproduce and then 1000 are consumed.

I want to get analytical solution. So, my attempt on former method, where steps are hours is based on binomial distribution. I think that probability for predator to consume number $n_r$ of red bacteria is: $$P_{\text{bin}} (n_r) = \frac{N!}{(N-n_r)!n_r!} (1-p)^{N-n_r} p^{n_r}$$ where $N=1000$ and $p$ is probability for predator to eat red bacteria which depends on current state. If I say that current number of red bacteria is $\alpha$, I can write
probability $p=\alpha/N$ and transition probability: $$P(\alpha+n|\alpha) = \frac{N!}{(N-\alpha+n)!\alpha-n!}(1-p)^{N-\alpha+n} p^{\alpha - n}$$ However when I try to simulate this with simple python program, where each step each bacteria duplicates itself and afterwards predator randomly eats half of the population I get different results. For example from simulation I get $P(501|500)=0.035$ and from analytical formula $P(501|500)=0.025$. Maybe someone can explain me why I get this difference and why binomial distribution is wrong choice.

After some benchmarking and inconsistencies I decided to switch to later method, where one bacterium reproduces and then one is eaten 1000 times per hour. At this case there are possible only three transitions from each state: \begin{eqnarray} P(\alpha + 1|\alpha) &=& P(\text{predator eats green})P(\text{red reproduces}) = \frac{N - \alpha}{N} \frac{\alpha}{N}\\ P(\alpha|\alpha) &=& P(\text{predator eats green})P(\text{green reproduces}) \\ && + P(\text{predator eats red})P(\text{red reproduces}) \\ &=& \frac{(N - \alpha)^2}{N^2} + \frac{\alpha^2}{N^2}\\ P(\alpha -1|\alpha) &=& P(\text{predator eats red})P(\text{green reproduces}) = \frac{\alpha}{N}\frac{N - \alpha}{N}\\ \end{eqnarray}

However further I don't know what to do next. Equilibrium distribution (state) is eigenvector $\vec\rho$ of transition matrix $$P_{ij} = P(i|j)$$ with eigenvalue 1. How I can find this vector. Writing eigenvalue problem in index notation gives me: $$\rho_\alpha = P(\alpha|\alpha-1)\rho_{\alpha - 1} + P(\alpha|\alpha)\rho_{\alpha} + P(\alpha|\alpha+1)\rho_{\alpha+1}$$ Or maybe I should use detailed balance condition: $$P(\alpha|\beta)\rho_{\beta} = P(\beta|\alpha)\rho_{\alpha}$$ Sorry for my english :).

• The idea here is to realize that eventually all reds will disappear, or all greens, and so the limiting probability distribution is extremal. – Michael Aug 14 '15 at 0:24
• Look at the Wright-Fisher model and the "time to fixation": en.wikipedia.org/wiki/Genetic_drift – user940 Aug 14 '15 at 13:55

The idea here is to realize that eventually either all reds disappear, or all greens, whichever comes first. So the limiting distribution is extremal.

Assume the subset of 1000 bacteria the predator eats every time slot is independent and uniformly distributed over all options. So, after each doubling step, each one of the 2000 bacteria has a probability of $1/2$ of being eaten.

One interesting aspect of this is that the expected number of reds is always the same, regardless of the number of reds we start with. To see this, let $N(t)$ be the number of reds (for time slots $t \in \{0, 1, 2, \ldots\}$), and assume $N(0)=n_0 \in \{0, 1, \ldots, 1000\}$ is a given constant. Then for each step $t\geq 0$ we get: $$E[N(t+1)|N(t)] = 2N(t)/2 = N(t)$$ So $E[N(t+1)]=E[N(t)]=n_0$ for all $t \in \{0, 1, 2, \ldots\}$.

This can be used to derive a simple expression for the probability that red eventually dies out, given a particular starting condition $n_0$.

• sorry, I don't get it. Why limiting distribution is extremal? How I can use fact that average number of red bacteria is constant to find stationary distribution? – skaydre zivri Aug 14 '15 at 9:29
• I said the limiting distribution is "extremal" because there are only two possibilities for the limiting outcome, and these two possibilities are polar opposites of each other. – Michael Aug 14 '15 at 16:51
• On expectations, we have $n_0 = E[N(t)] = \sum_{n=0}^{1000}nPr[N(t)=n]$, and we know something about $\lim_{t\rightarrow\infty}Pr[N(t)=2]$. – Michael Aug 14 '15 at 16:53
• However, to solve your first question, you do not need to use this fact about expectations (I just wrote that as an interesting side-note). Your first question can be more intuitively solved by realizing there are only two possible limiting outcomes, and reasoning out what their probabilities should be. – Michael Aug 14 '15 at 16:57