# Coding a simulation of multi-species population dynamics

I want to code several steps in a recursion equation (in my case, migration, followed by selection, followed by mating, etc in a population). I have seen this done using a lot of copying and pasting and Do loops. I am hoping someone has a better suggestion that makes use of Mathematica's recursion equation construction and memoization.

In step 1 I was able to code a pair of coupled recursion equations.

Step 1: Frequency of genotype i in population 1 (x1) and population 2 (x2) at time t depends on migration rate between the two populations and the frequency of genotype i in each at time t - 1.

parameters = {m->0.01}
x1[i_, t_] := x1[i, t] = (1 - m) x1[i, t - 1] + m x2[i, t - 1] /. parameters
x2[i_, t_] := x2[i, t] = (1 - m) x2[i, t - 1] + m x1[i, t - 1] /. parameters

x1[1,0] := 1
x2[1,0] := 0

ListPlot[Table[{{t, x1[1, t]}, {t, x2[1, t]}}, {t, 0, 1000}], PlotRange -> {0, 1}]


Now I want to add a step to my recursion equations so that selection also occurs in the same generation. Step 2 depends on x1[i, t] and x2[i, t] (as they are currently written) but I would not consider it time t + 1.

Not a reproducible example below because I don't know how to code this part

Step 2: Selection in population 1

x1Selection[1,t] = ((1 + s1) x1[1,t])/(1 + s1 x1[1,t])
x2Selection[1,t] = x2[1,t]


Step 3: Not shown

Step 4: Not shown

So I want to go through these steps and get the frequency of genotype i at time t, after migration, selection, etc have happened. Then repeat step 1 using x1Final[1, t] and x2Final[1, t] in place of x1[1, t - 1] and x2[1, t - 1], respectively.

Is there an elegant way to code a multistep recursion equation without copying and pasting the output of the previous step into the next step (because this gets considerably uglier in steps 3 and 4)? That is, without writing something like:

parameters = {m->0.01, s1->0.05}
x1[i_, t_] := x1[i, t] = ((1 + s1) ((1 - m) x1[i, t - 1] + m x2[i, t - 1]))/(1 + s1 ((1 - m) x1[i, t - 1] + m x2[i, t - 1])) /. parameters
x2[i_, t_] := x2[i, t] = (1 - m) x2[i, t - 1] + m x1[i, t - 1] /. parameters

x1[1,0] := 1
x2[1,0] := 0

ListPlot[Table[{{t, x1[1, t]}, {t, x2[1, t]}}, {t, 0, 1000}], PlotRange -> {0, 1}]

• You can use Composition to combine the different steps for your Genetic Algorithm and then use NestList to recursively simulate new generations. – Thies Heidecke Sep 17 '16 at 23:56
• I suggest omitting any attempt to memoise until you have fully debugged your code (and only implement it if essential for efficiency). – mikado Sep 18 '16 at 7:46
• @ThiesHeidecke I am reading through the documentation and having trouble understanding how to set up the function definitions to use with Composition. Would it be like x1Migration[i, t] = (1 - m) x1[i, t] + m x2[i, t] and then x1Selection[x1Migration[i,t]]:= ((1 + s1) x1Migration[i,t])/(1 + s1 x1Migration[i,t]) and then Composition[x1, x1Selection, x1Migration][i,t]? – biologyUser Sep 18 '16 at 18:13

I would drop the whole idea of using recursion equations and simply write a function that expresses how the next generation depends on the current one.

Like so.

With[{m = .01, s1 = .01},
nextGen[{x1_, x2_}] :=
Module[{x1tmp, x2tmp},
x1tmp = (1 - m) x1 + m x2;
x1tmp = (1 + s1) x1tmp/(1 + s1 x1tmp);
x2tmp = (1 - m) x2 + m x1;
{x1tmp, x2tmp}]]


The above function is easy to extend with further algorithmic steps (e.g., your steps 3 and 4, whatever they will be) and to the addition of species (e.g., {x1, x2, x3}) should that be needed. It also evaluates quite quickly.

With nextGen, the history over 1000 generations can be generated with

populations = Transpose @ NestList[nextGen, {1, 0}, 1000];
ListPlot[populations, PlotRange -> {0, 1}]


### Update

In response to the OP's comment, I give a vectorized version of nextGen.

With[{m = .01, s = {.0024, -.002}},
nextGen[pop_List] :=
Block[{sp},
sp[1] = (1 - m) pop[[1]] + m  pop[[2]];
sp[2] = (1 - m) pop[[2]] + m pop[[1]];
sp[1] = (1 + s[[1]]) sp[1]/(1 + s[[1]] sp[1]);
sp[2] = (1 + s[[2]]) sp[2]/(1 + s[[2]] sp[2]);
Table[sp[i], {i, Length[pop]}]]]


Nothing changes in the way the simulation is run.

populations = Transpose @ NestList[nextGen, {1, 0}, 1000];
ListPlot[populations, PlotRange -> {0, 1}]


The result is different because the 2nd species also undergoes selection and the selection parameters have been given different values (just to produce an different plot).

• Is there a way to do it without giving all the different species different names? I'll need to keep track of eight total. I was hoping to index them by population (k) and genes (pi and ti) to reduce the chance of making errors later. e.g. migration: xm[k_, pi_, ti_] := (1 - m) x[k, pi, ti] + m x[3 - k, pi, ti] /.parameters followed by selection: xs[k_,pi_,ti_] := ((1 + KroneckerDelta[ti, k] s[k]) xm[k, pi, ti])/(1 + s[k] (Sum[xm[k, pi, k], {pi, 1, 2}])) – biologyUser Sep 18 '16 at 17:26
• @user43111. What you bring up in your comment looks like a different question than the one you asked above. I don't like bait and switch. You should edit your question so it expresses your full, actual problem. Do not dribble it out a bit at a time -- especially not in a comment to an answer. – m_goldberg Sep 18 '16 at 20:11
• @biologyUser I'm a bit confused about what exactly you're trying to model and whether the equations make sense. Are you using the words population and species interchangeably, or are you envisioning multiple interacting species (or genotypes?) in a multi-patch system? – Chris K Sep 20 '16 at 13:19