Coding a simulation of multi-species population dynamics

Question

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}]

Answer 1

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