Clear["Global`*"]; iteration1Gen[populationPhenotypesWithGeneration_] := Block[{generationN = populationPhenotypesWithGeneration[[1]],populationPhenotypes = populationPhenotypesWithGeneration[[2]],NewGeneration,ParapatricFitness,PhilopatricFitness,fitness,survivors,deads,z,d,z1,d1,fertilityT}, If[Divisible[generationN, 500], PrintTemporary[sd]; PrintTemporary[generationN]; PutAppend[N[populationPhenotypes], stream] ,]; (* generation *) generationN++; (* philopatric and parapatric components of fitness *) fertilityT = ParallelTable[ Clip[ Sum[fertility[populationPhenotypes[[n,j,1]], populationPhenotypes[[n,k,1]]], {k,Delete[Table[m,{m,1,Nind}],j]}],{0.,Infinity} ] , {n, Nd}, {j, Nind}]; ParapatricFitness = ParallelTable[ fertilityT[[n,j]]*populationPhenotypes[[n,j,2]]*sd /(Nd-1), {n, Nd}, {j, Nind}] ; (* construct new generation from fitness and mutations *) NewGeneration = ParallelTable[ (* fitness wrt deme i *) PhilopatricFitness = Table[ fertilityT[[i,k]]*(1 - populationPhenotypes[[i,k,2]]), {k, Nind}]; fitness = ReplacePart[ ParapatricFitness, i -> PhilopatricFitness] ; (* find out the survivors dudes of deme i *) deads = {RandomChoice[ Range[Nind]]}; survivors = Complement[Range[Nind],deads]; Table[ Which[ MemberQ[survivors,j], {z, d} = populationPhenotypes[[i,j]], MemberQ[deads,j], {z1, d1} = RandomChoice[ Flatten[ fitness ] -> Flatten[populationPhenotypes, 1]]; event = RandomChoice[{1 - Pmutation, Pmutation} -> {0, 1}]; Which[ event == 0, {z, d} = {z1, d1}, event == 1, {z, d} = {z1, d1} + RandomVariate[BinormalDistribution[{0, 0}, {mutationStepZ,mutationStepD}, mutationCorr]]; z = Clip[z , {0.,1.}]; d = Clip[d , {0.,1.}] ]; ]; {z, d} , {j, Nind}], {i, Nd} ]; Prepend[{NewGeneration}, generationN] ]; (* Demographic parameters *) Nind = 8 (* number of individual per patch *); Nd = 1000 (* number of demes *); Ngenerations = 150000 (* Total number of generations to simulate *); (* fertility *) fertility[zi_,zj_] := 1. + b1 (zi + zj) + b2 (zi + zj)^2 - c1 (zi) - c2 (zi)^2; (* survival during dispersal *) sd; b1 = 6; b2 = -1.4; c1 = 4.56; c2 = -1.6; (* Mutation parameters *) Pmutation = 0.01 (* Probability of mutating *); mutationStepZ = .02 (* SD in mutation step *); mutationStepD = .02 (* SD in mutation step *); mutationCorr = 0. (* Correlation between mutations of two strategies *); (* Initial population *) ConvergentZ = N[ (c1 (-1 + Nind (-1 + sd)) + b1 (2 + Nind - sd - Nind sd))/(2 (2 b2 (-2 + Nind (-1 + sd) + sd) + c2 (1 + Nind - Nind sd)))]; ConvergentDisp = N[ 1/(1 + Nind - Nind sd)]; Do[ filename = "/Users/cmullon/Check_Wakano_N8_sd"<>ToString[sd]<>".csv"; stream = OpenWrite[filename, PageWidth -> Infinity]; InitialPopulation = Prepend[{ Table[{ConvergentZ, ConvergentDisp}, {i, Nd}, {j, Nind}]}, 0] ; PutAppend[InitialPopulation[[2]], stream]; Nest[iteration1Gen, InitialPopulation, Ngenerations]; Close[stream]; ,{sd,{0.04,0.2,0.5,0.8,0.91,0.92,0.93,0.94,0.94,0.96}} ] Quit[];