An analysis of expected population states in the space of population states: 
the case study of three-element populations

Iwona Karcz-Duleba, Andrzej Cichon

(* 
the component integral (Subscript[I, 1])(p, r, s) 
*)
IntEx1[b_Real, p_Real,r_Real,s_Real]:=Chop[(Sqrt[(π/2)])(1/b)*Exp[-(1/2) (b^2) (((p-r))^2)]*Erfc[((b (2 s-p-r))/(Sqrt[6]))]+p π *(Erfc[((b (s-r))/(Sqrt[2]))]-Erf[((b (r-p))/(Sqrt[2]))]-(b/(Sqrt[π]))NIntegrate[Exp[-(b^2) ((((Subscript[x, 2])-r))^2)]* Erf[b ((Subscript[x, 2])-p)] *Erf[b ((Subscript[x, 2])-s)],{(Subscript[x, 2]),-∞,∞},AccuracyGoal-> 8,WorkingPrecision-> 9,
Method->{Automatic,"SymbolicProcessing"->0}])];


(*
(Subscript[x, 1]) coordinate of the expected state
Ex1[b, p, r, s, (Subscript[α, 1]), (Subscript[α, 2]), (Subscript[α, 3])]
*)
Ex1[b_Real,p_Real,r_Real,s_Real, a1_Real,a2_Real,a3_Real]:=
(3/(2 π))*(Plus @@ (Flatten[Outer[Times,{a1,a2,a3}, {a1,a2,a3},{a1,a2,a3}]] *
 Flatten[Outer[IntEx1,  {b},{p,r,s}, { p,r,s},{p,r,s}]])
)



(* 
the component integral (Subscript[I, 2])(p, r, s) 
*)
IntEx2[b_Real,p_Real,r_Real,s_Real]:=Chop[(Sqrt[(π/2)])(1/b)*(Exp[-(1/2) (b^2) (((r-s))^2)]-Exp[-(1/2) (b^2) (((r-p))^2)])+r π *(Erfc[((b (s-r))/(Sqrt[2]))]-Erf[((b (r-p))/(Sqrt[2]))])-b (Sqrt[π]) NIntegrate[(Subscript[x, 2])*Exp[-(b^2) ((((Subscript[x, 2])-r))^2)]* Erf[b ((Subscript[x, 2])-p)] *Erf[b ((Subscript[x, 2])-s)],{(Subscript[x, 2]),-∞, ∞},AccuracyGoal-> 8,WorkingPrecision-> 9,
Method->{Automatic,"SymbolicProcessing"->0}]];



Ex2[b_,p_,r_,s_, a1_,a2_,a3_]:=
(3/(2 π))*(Plus @@ (Flatten[Outer[Times,{a1,a2,a3}, {a1,a2,a3},{a1,a2,a3}]] *
 Flatten[Outer[IntEx2,  {b},{p,r,s}, { p,r,s},{p,r,s}]])
)



(* 
the component integral (Subscript[I, 3])(p, r, s) 
*)
IntEx3[b_Real,p_Real,r_Real,s_Real]:=Chop[
-(Sqrt[(π/2)])(1/b)*Exp[-(1/2) (b^2) (((r-s))^2)]*(Erfc[((b (r+s-2 p))/(Sqrt[6]))])+s π(Erfc[((b (s-r))/(Sqrt[2]))]-Erf[((b (r-p))/(Sqrt[2]))]-(b/(Sqrt[π]))NIntegrate[Exp[-(b^2) (((x-r))^2)]Erf[b(x-p)]*Erf[b(x-s)],{x, -∞, ∞},AccuracyGoal-> 8,WorkingPrecision->9,
Method->{Automatic,"SymbolicProcessing"->0}])];


Ex3[b_,p_,r_,s_, a1_,a2_,a3_]:=
(3/(2 π))*(Plus @@ (Flatten[Outer[Times,{a1,a2,a3}, {a1,a2,a3},{a1,a2,a3}]] *
 Flatten[Outer[IntEx3,  {b},{p,r,s}, { p,r,s},{p,r,s}]])
)

SetDirectory["C:\\"]; (* where output be directed  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*)

(* unimodal triangle fitness function *)
triangleuni[x_,A_,B_]:=Piecewise[{{x/A+1.0, -A≤ x<0},
{-x/B+ 1.0, 0≤ x≤  B}
}];

(* unimodal gaussian fitness function *)
gaussuni[x_,a1_,a2_]:=Piecewise[{{Exp[-a1 x^2], x≤ 0},
{Exp[-a2 x^2],x > 0}
}];


q[x_]:= Module[{a=5,
			a1=5, 
			a2=1,
			 h = 2,
			 A=0.5,
			 B=0.5}, 


(*
d = 1;
Global`funname = "gaussbi" <> ToString[d];
Global`params = "_"<>ToString[a1]<>"_"<>ToString[a2]<>"_"<>ToString[h];
Return[ Exp[-a1 x^2] + h Exp[-a2 (x-d)^2]  ];
*)



Global`funname = "gaussuni";
Global`params = "_"<>ToString[a1]<>"_"<>ToString[a2];
Return[gaussuni[x, a1,a2]];



];

simulation[σ_, s0_List, iter_Integer, saveToFile_]:=Module[{b,p,r,s, a1,a2,a3,x1i,x2i,x3i, ex1i,ex2i,ex3i,totalFitness,i, RES, filename, fits},

b=(1/((Sqrt[2])σ));

{x1i, x2i,x3i}={p,r,s} = Sort[s0, Greater];

RES = {};
AppendTo[RES, {x1i, x2i,x3i}];

Print["0.   [", x1i,", ", x2i,", ", x3i,"]"];


For[i=1, i <=iter, i++;

fits = Map[q, {x1i, x2i,x3i}];
totalFitness = Apply[Plus, fits ];
{a1,a2,a3} = fits/totalFitness;


ex1i=Ex1[b,x1i, x2i, x3i,a1,a2,a3];
ex2i=Ex2[b,x1i, x2i, x3i,a1,a2,a3];
ex3i=Ex3[b,x1i, x2i, x3i,a1,a2,a3];

{x1i, x2i,x3i}=Sort[{ ex1i,ex2i,ex3i}, Greater];


Print[i-1,".   [", x1i,", ", x2i,", ", x3i,"]"];

AppendTo[RES, { ex1i,ex2i,ex3i}];

];

(* export the results to a text file *)
If[saveToFile,
filename = "m3_" <> ToString[σ] <> "_" <> ToString[p] <> "_" <> ToString[r] <> "_" <> ToString[s] <> "_" <>funname<> params<>".txt";
Export[filename, RES, "List"];
];

];

σ = 0.2;
s0 = {0.5, 0.5, 0.5};
iterations =30;
saveToFile = False;

simulation[σ, s0, iterations, saveToFile]