9
$\begingroup$

I was wondering how to apply Markov steady-state chains to a simple case in genetics, eye colour.

Assuming the probabilities of being brown, green or blue eyed, given your parents were a certain colour combination are:

Br + Br - Br 75%  , Gr 18.75% , Bl 6.25%
Gr + Br - Br 50%  , Gr 37.5%  , Bl 12.5%
Bl + Br - Br 50%  , Gr 0%     , Bl 50%
Gr + Gr - Br 0.5% , Gr 74.75% , Bl 24.75%
Gr + Bl - Br 0%   , Gr 50%    , Bl 50%
Bl + Bl - Br 0%   , Gr 1%     , Bl 99%

Corresponding matrix:

{{3/4, 1/2, 1/2, 0.005, 0, 0}, {0.1875, 0.375, 0, 0.7475, 1/2, 1/
  100}, {0.0625, 0.125, 1/2, 0.2475, 1/2, 99/100}}

I want to find the steady state for the eye colour of the population.

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ If you've made any attempt at solving this yourself it would be a good idea to share that attempt. Providing data that's been formatted for Mathematica is also appreciated. $\endgroup$
    – C. E.
    Aug 28, 2013 at 1:43
  • $\begingroup$ This question requires an initial gentype frequency. some assumption or declaration of probabilities of genotypes mating and then use of Mendelian genetic rules to develop transition matrix for genotypes (as states). The eye colors can be inferred from the genotype-phenotype relationships after determining stationary(steady state) distribution for given initial state. The matrix you have, though rows add to 1, is not a transition matrix, eg entry (1,1) is not transition from state 1 to state 1 etc. $\endgroup$
    – ubpdqn
    Aug 28, 2013 at 4:20
  • $\begingroup$ I haven't made any attempts that led anywhere. I saw the table of values on twitter and immediately wondered what the steady-state would be. I was willing to assume that any one person can marry any other person on Earth to remove any regional differences there may be, but I was hoping that it was possible to somehow use the steady-state as the probability that person A is Br/Gr/Bl and same for person B and solve? I haven't studied Markov chains in that much detail and it seems like it is a lot more complicated than I had hoped. $\endgroup$
    – DannyBland
    Aug 28, 2013 at 8:14
  • $\begingroup$ I thought about setting: P(Brown)=0.75*P(Brown)*P(Brown)*... etc. It would be assuming that there is equal likelihood of a person with brown eyes having a child with a person with green, blue or brown eyes (so there is no bias from a particular eye colour to a particular eye colour) and that eye colour is spread the same universally. It would seem to lead to a solution but I don't have a maths package and there are too many terms for Wolfram Alpha. $\endgroup$
    – DannyBland
    Sep 9, 2013 at 17:03

4 Answers 4

Reset to default
4
$\begingroup$

@YayaManisya's answer can be acquired by numerical fixed point calculation:

FixedPoint[
 Apply[# / Total@# & @*
   Function[{br, gr, bl},
    {br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
     {{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
      {.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}}]],
 {1, 1, 1} / 3]

{0.0218527, 0.154522, 0.823625}

Here we repeatedly apply a function which computes a vector of pairwise probabilities for parents, and multiply that with the corresponding offspring probability matrix, and in the end normalizes the sum of probabilities to 1, mostly to avoid numerical errors destroying the result. FixedPoint nests these calls until the result doesn't change (much).

This result doesn't seem to vary in this case depending on the starting population.

Change in population eye color can be seen here for one starting population:

NestList[
   Apply[# / Total@# & @*
     Function[{br, gr, bl},
      {br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
       {{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
        {.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}}]],
   {0.9, 0.075, 0.025}, 100] // Transpose // 
 ListLinePlot[#, PlotStyle -> {Brown, Green, Blue}] &

enter image description here

Also, the solution can be found symbolically (especially now when we are convinced there's just one steady state):

Solve[0 <= br <= 1 && 0 <= gr <= 1 && 0 <= bl <= 1 && 
   br + gr + bl == 1 &&
   {br, gr, bl} == {br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
     Rationalize@{{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
       {.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}},
  {br, gr, bl}] // FullSimplify // ToRadicals

$$\text{br}\to \frac{2 \left(-342-46790 \sqrt{2}+5 \sqrt{4438369759-409761860 \sqrt{2}}\right)}{22339241},\\\text{ gr}\to \frac{10 \left(-93580-342 \sqrt{2}+5 \sqrt{2 \left(4438369759-409761860 \sqrt{2}\right)}\right)}{22339241} ,\\\text{bl}\to \frac{46551450+194000 \sqrt{2}-22339241 \sqrt{\frac{87264648203600}{499041 688456081}+\frac{9394337092000 \sqrt{2}}{499041688456081}}}{44678 482}$$

These values correspond to the numerical values achieved with FixedPoint.

Improve this answer
$\endgroup$
3
$\begingroup$

Here's an easy and probably useless implementation. Imagine nine states in a Markov chain. Each state is an ordered pair, drawn from $(Br,\ Gr,\ Bl)$ that describes the eye color of a person and their mate. If you assume that a person chooses a mate with one of the three eye colors with equal probability, then I think the following MarkovProcess tracks a single path through an infinite family tree (or the 100 generations I simulate below):

proc = DiscreteMarkovProcess[Table[1/9, {9}], 
    Flatten@Table[Table[p, {3}], {p, #}]/3 & /@ {
       {0.75, 0.1875, 0.0625},
       {0.5, 0.375, 0.125},
       {0.5, 0, 0.5},
       {0.5, 0.375, 0.125},
       {0.005, 0.7475, 0.2475},
       {0, 0.5, 0.5},
       {0.5, 0, 0.5},
       {0, 0.5, 0.5},
       {0, 0.01, 0.99}}];
ListPlot[RandomFunction[proc, {0, 100}],
    Filling -> Axis, TicksStyle -> 16, ImageSize -> 500,
    Ticks -> {Automatic, Flatten[Table[{3 (i - 1) + j, 
        StringJoin[{"(", {"Br", "Gr", "Bl"}[[i]], ",", {"Br", "Gr", "Bl"}[[j]], ")"}]},
        {i, 3}, {j, 3}], 1]}
]

enter image description here

Hardly useful! But maybe interesting as a simple example of a Markov Chain simulation? Of course, there's no steady state here....

Improve this answer
$\endgroup$
3
$\begingroup$

This is a simulation and perhaps will motivate better answers or can be adapted. In the following:

  • the starting population is (1/3,1/3,1/3) blue:brown:green

  • each generation produces one child from one coupling with probability of color as per OP

  • no deaths

Again this is just to a toy example to illustrate evolution over time of a population with the above assumptions. The gif is a downsampling of vis/@sim[100,1000]:

mat = Transpose[{{3/4, 1/2, 1/2, 0.005, 0, 0}, {0.1875, 0.375, 0, 
     0.7475, 1/2, 1/100}, {0.0625, 0.125, 1/2, 0.2475, 1/2, 99/100}}];
mnpar = {{"Blue", "Blue"} -> mat[[6]], {"Blue", "Brown"} -> 
    mat[[3]], {"Blue", "Green"} -> mat[[5]], {"Brown", "Brown"} -> 
    mat[[1]], {"Brown", "Green"} -> mat[[2]], {"Green", "Green"} -> 
    mat[[4]]};
fun[1, 0, 0] := "Blue";
fun[0, 1, 0] := "Brown";
fun[0, 0, 1] := "Green";
sim[n_, num_] := 
 Module[{rv = 
    RandomVariate[MultinomialDistribution[1, {1/3, 1/3, 1/3}], n]},
  Map[fun,
   NestList[{RandomVariate[
        MultinomialDistribution[
         1, (Sort[fun @@@ RandomSample[#, 2]] /. mnpar)]]}~Join~# &, 
    rv, num], {2}]
  ]
vis[lst_] := Module[{col = {"Blue", "Brown", "Green"}, tly},
  tly = {#, Count[lst, #]} & /@ col;
  BarChart[#2/Total[tly[[All, 2]]] & @@@ tly, ChartLabels -> col, 
   PlotLabel -> Row[{"Population:", Total[tly[[All, 2]]]}], 
   LabelingFunction -> Above, PlotRange -> {0, 0.5}, Frame -> True, 
   ChartStyle -> {Blue, Brown, Green}]]

enter image description here

Improve this answer
$\endgroup$
2
$\begingroup$

Actually you just need to assign unknowns to the steady state probabilities and sum them up accordingly.

Let $P(W)$, $P(G)$ and $P(B)$ represent the probabilities of being brown, green and blue-eyed, and thus

$P(W) + P(G) + P(B) = 1$ ,

assuming that you get no other eye colours in children due to illness etc.

Then you will get

$\begin{eqnarray}P(W) = &&0.75[P(W)]^2 + [P(W)][P(G)] \\ &&+[P(W)][P(B)] + 0.005[P(G)]^2\end{eqnarray}$

$\begin{eqnarray}P(G) = &&0.1875[P(W)]^2 + 0.75[P(W)][P(G)] + 0.7475[P(G)]^2 \\ &&+ [P(G)][P(B)] + 0.01[P(B)]^2\end{eqnarray}$

$\begin{eqnarray}P(B) = &&0.0625[P(W)]^2 + 0.25[P(W)][P(G)] + 0.2475[P(G)]^2\\ &&+ [P(W)][P(B)] + [P(G)][P(B)] + 0.99[P(B)]^2\end{eqnarray}$

From here on you can use normal calculation method to find the answer. Simplify the equations first, and after that you can use any calculator that enable square root calculation to find the final numerical answers.

I get

$\begin{eqnarray} P(W) &&= 1 / (1 + 2.5 * \sqrt{199 + 20 * \sqrt{2}} + 5 * \sqrt{2})\\ &&= 0.0218527 \\ &&= 2.2\% \text{ of the population,}\end{eqnarray}$

$\begin{eqnarray}P(G) &&= P(W) * 5 * \sqrt{2}\\ &&= 0.1545219 \\ &&= 15.5\% \text{ of the population,}\end{eqnarray}$

$\begin{eqnarray}P(B) &&= 1 - P(W) - P(G) \\ &&= 82.3\% \text{ of the population.}\end{eqnarray}$

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ This site is about the Mathematica software. Can you complement you answer with valid MMA code to carry out the calculations? $\endgroup$
    – corey979
    Sep 13, 2018 at 17:42
  • $\begingroup$ Welcome to Mathematica.SE! I have edited your answer using MathJax to make it more readable. Please check, if I introduced any errors. $\endgroup$
    – Johu
    Sep 13, 2018 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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