Epidemiological Model

Deterministic Dynamics

Differential Equations

Symmetry and scaling assumption of transmission rates

βSub={β[1,1]->(X/κ),β[1,2]->(Y/κ),β[2,1]->(Y/κ),β[2,2]->(X/κ)};

Differential Equations

((ds[h,t])/dt)=f1[h,t]

f1[h_,t_]:=b(s[h,t]+Sum[i[h,l,t],{l,1,2}])(1-((Sum[s[k,t]+Sum[i[k,l,t],{l,1,2}],{k,1,2}])/κ))-δ s[h,t]-Sum[β[h,p] s[h,t]Sum[i[k,p,t],{k,1,2}],{p,1,2}]+γ Sum[i[h,l,t],{l,1,2}]/.βSub

((dsi[h,p,t])/dt)=f2[h,p,t]

f2[h_,p_,t_]:=β[h,p]s[h,t]Sum[i[k,p,t],{k,1,2}]-(δ+α) i[h,p,t]-γ i[h,p,t]/.βSub

Selection coefficients

σH=((f1[1,t]+f2[1,1,t]+f2[1,2,t])/(s[1,t]+i[1,1,t]+i[1,2,t]))-((f1[2,t]+f2[2,1,t]+f2[2,2,t])/(s[2,t]+i[2,1,t]+i[2,2,t]))//FullSimplify

σP=Simplify[((f2[1,1,t]+f2[2,1,t])/(i[1,1,t]+i[2,1,t]))-((f2[1,2,t]+f2[2,2,t])/(i[1,2,t]+i[2,2,t]))]

Transformation of Variables

DerivSub=Flatten[Join[Table[D[s[h,t],t]->f1[h,t],{h,1,2}],Table[D[i[h,p,t],t]->f2[h,p,t],{h,1,2},{p,1,2}]]];

VarTrans=Flatten[Solve[{nH[t]==Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}],nP[t]==Sum[i[h,p,t],{h,1,2},{p,1,2}],pH[t]==((s[1,t]+i[1,1,t]+i[1,2,t])/(nH[t])),pP[t]==((i[1,1,t]+i[2,1,t])/(nP[t])),ρH[t]==((i[1,1,t]+i[1,2,t])/(nP[t])),Δ[t]==i[1,1,t]-nP[t]ρH[t]pP[t]},Flatten[Join[Table[s[h,t],{h,1,2}],Table[i[h,p,t],{h,1,2},{p,1,2}]]]]];

g1[t_]:=D[Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}],t]/.DerivSub/.VarTrans//Simplify

g2[t_]:=D[Sum[i[h,p,t],{h,1,2},{p,1,2}],t]/.DerivSub/.VarTrans//Simplify

g3[t_]:=D[((s[1,t]+i[1,1,t]+i[1,2,t])/(Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}])),t]/.DerivSub/.VarTrans//Simplify

g4[t_]:=D[((i[1,1,t]+i[2,1,t])/(Sum[i[h,p,t],{h,1,2},{p,1,2}])),t]/.DerivSub/.VarTrans//Simplify

g5[t_]:=D[((i[1,1,t]+i[1,2,t])/(Sum[i[h,p,t],{h,1,2},{p,1,2}])),t]/.DerivSub/.VarTrans//Simplify

g6[t_]:=D[i[1,1,t]-(((i[1,1,t]+i[1,2,t])(i[1,1,t]+i[2,1,t]))/(Sum[i[h,p,t],{h,1,2},{p,1,2}])),t]/.DerivSub/.VarTrans//FullSimplify

equListTrans={g1[t],g2[t],g3[t],g4[t],g5[t],g6[t]};

varsTrans={nH[t],nP[t],pH[t],pP[t],ρH[t],Δ[t]};

Equilibria

equiSolTrans=Solve[equListTrans=={0,0,0,0,0,0},varsTrans]//Simplify;

equiSolTrans//Length

Equilibria 1-3 are disease extinction

Equilibria 4,5,6,7,10,11,12,13 are fixation equilibria

Equilibria  8 and 9 however are endemic and polymorphic.

equiSolTrans[[8]]

varsTrans

gx[e_]:=ToExpression["g"<>ToString[e]<>"[t]"]

Table[{v1,v2},{v1,1,6},{v2,1,6}]//MatrixForm;

jMtrxTrans=Table[D[gx[v1],varsTrans[[v2]]],{v1,1,6},{v2,1,6}];

testPars={b->1,X->0.5,Y->0.2,α->0.1,δ->0.1,γ->0,κ->1};

equiSolTrans[[8]]/.testPars

esys=Simplify[jMtrxTrans/.equiSolTrans[[8]]]/.testPars//Eigensystem//Chop

esys[[1]]

esys[[2]]//MatrixForm

equiSolTrans[[8]]

equTemp={nH[t]->(((b X+b Y-X α-Y α-X δ-Y δ+AA) κ)/(2 b (X+Y))),nP[t]->(((b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+AA) κ)/(2 b (X+Y))),pH[t]->(1/2),pP[t]->(1/2),ρH[t]->(1/2),Δ[t]->(((X-Y) (b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+AA) κ)/(8 b (((X+Y))^2)))};

CPTrans=Factor[CharacteristicPolynomial[jMtrxTrans/.equTemp,λ]]

Collect[CPTrans[[6]],λ]

Collect[CPTrans[[7]],λ,Simplify]

Equilibria

The list of differential equations and variables

equList=Join[Table[f1[h,t],{h,1,2}],Table[f2[h,p,t],{h,1,2},{p,1,2}]]//Flatten;

vars=Join[Table[s[h,t],{h,1,2}],Table[i[h,p,t],{h,1,2},{p,1,2}]]//Flatten;

Jacobian Matrix

jMtrx=Join[Table[{D[f1[h1,t],s[1,t]],D[f1[h1,t],s[2,t]],D[f1[h1,t],i[1,1,t]],D[f1[h1,t],i[1,2,t]],D[f1[h1,t],i[2,1,t]],D[f1[h1,t],i[2,2,t]]},{h1,1,2}],Flatten[Table[{D[f2[h1,p1,t],s[1,t]],D[f2[h1,p1,t],s[2,t]],D[f2[h1,p1,t],i[1,1,t]],D[f2[h1,p1,t],i[1,2,t]],D[f2[h1,p1,t],i[2,1,t]],D[f2[h1,p1,t],i[2,2,t]]},{h1,1,2},{p1,1,2}],1]]//Simplify;

Extinction equilibrium

equList/.{s[1,t]->0,s[2,t]->0,i[1,1,t]->0,i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0}

jMtrx/.{s[1,t]->0,s[2,t]->0,i[1,1,t]->0,i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0}//Eigenvalues

Disease Extinction equilibrium

assump={(*s[1,t]->0,s[2,t]->0,*)i[1,1,t]->0,i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0};

temp=equList/.{(*s[1,t]->0,s[2,t]->0,*)i[1,1,t]->0,i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0}

equilSol=Solve[temp[[{1,2}]]=={0,0},{s[1,t],s[2,t]}]

Stability

λs=Eigenvalues[jMtrx/.assump]//Simplify

Eigenvalue 3: Escape from extinction

For eigenvalue 3 this reduces to:

((b κ-δ κ)/κ)-2 b ((((((b-δ) κ)/b)))/κ)//Simplify

Eigenvalue 4: Equilibrium “manifold”

For eigenvalue 4 this reduces to:

((b κ-δ κ-b((((b-δ) κ)/b)) )/κ)//Simplify

Eigenvalue 5 and 6: R0

For eigenvalue 5

Let φ be the proportion of the population that is host type 1.

temp=λs[[5]]/.{s[1,t]->ψ ((((b-δ) κ)/b)),s[2,t]->(1-ψ) ((((b-δ) κ)/b))}//FullSimplify

Reduce[{temp<0,0<δ<b,0<Y<X,0<κ,0<α,0<γ,0<ψ<1},X,FullSimplify]

Looking at the condition: X<((-b α-b γ-b δ+b Y ψ-Y δ ψ)/(-b+δ+b ψ-δ ψ))

Collect[((-b α-b γ-b δ+b Y ψ-Y δ ψ)/(-b+δ+b ψ-δ ψ)),Y,Simplify]

Hence this condition can be rewritten as:

(1-ψ)X+ψ Y<((b (α+γ+δ))/((b-δ)))=((κ(α+ γ+ δ))/sEqu)

For eigenvalue 6

temp=λs[[6]]/.{s[1,t]->ψ ((((b-δ) κ)/b)),s[2,t]->(1-ψ) ((((b-δ) κ)/b))}//FullSimplify

Reduce[{temp<0,0<δ<b,0<Y<X,0<κ,0<α,0<γ,0<ψ<1},X,FullSimplify]

Collect[((-b Y+b α+b γ+b δ+Y δ+b Y ψ-Y δ ψ)/(b ψ-δ ψ)),Y,Simplify]

Hence this condition can be rewritten as:

ψ X+(1-ψ) Y<((b (α+γ+δ))/((b-δ)))=((κ(α+ γ+ δ))/sEqu)

Single Host Single Pathogen (Matching)

Solving for equilibrium

assump={(*s[1,t]->0,*)s[2,t]->0,(*i[1,1,t]->0,*)i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0};

temp=equList/.{(*s[1,t]->0,*)s[2,t]->0,(*i[1,1,t]->0,*)i[1,2,t]->0,i[2,1,t]->0,i[2,2,t]->0}

s1sol=Flatten[Solve[temp[[3]]==0,s[1,t]]]

Factor[temp/.s1sol][[1]]

This is a square in i[1,1,t].  Solving for a variable may be useful instead

bSol=Flatten[FullSimplify[Solve[Factor[temp/.s1sol][[1]]==0,b]]]

i11sol=FullSimplify[Solve[Factor[temp/.s1sol][[1]]==0,i[1,1,t]],Assumptions->{0<Y<X,0<δ<b,0<α,0<γ,0<κ}]

equ1=Join[s1sol/.i11sol[[1]],i11sol[[1]],assump]

equ2=Join[s1sol/.i11sol[[1]],i11sol[[1]],assump]

Biological Validity

Reduce[{0<(i[1,1,t]/.i11sol[[1]]),0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//FullSimplify

Reduce[{0<(i[1,1,t]/.i11sol[[2]]),0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//Simplify

Stability

jMtrx[[{1,4},{1,4}]]/.{s[2,t]->0,i[1,1,t]->0,i[2,1,t]->0,i[2,2,t]->0}

equilibraMM={s[1,t]->(((α+γ+δ) κ)/Y),i[1,1,t]->(1/(2 b Y))(-Y (α+δ)+b (Y-2 (α+γ+δ))+(Sqrt[Y ((b^2) Y+Y (((α+δ))^2)+b (-2 Y (α+δ)+4 α (α+γ+δ)))])) κ}

jMtrx[[{1,3},{1,3}]]

temp=jMtrx/.assump/.s1sol/.bSol//Simplify;

jmtrxInteral=jMtrx[[{1,3},{1,3}]]/.equ1//Simplify

Simplify[jmtrxInteral//Eigenvalues]/.{b->1,X->6,δ->0.5,α->1.5,γ->0.1}

Simplify[jmtrxInteral//Eigenvalues]/.{b->1,X->7,δ->0.5,α->1.5,γ->0.1}

Simplify[jmtrxInteral//Eigenvectors]/.{b->1,X->7,Y->3.5,δ->0.05,α->0.1,γ->0.1}

Collect[CharacteristicPolynomial[jmtrxInteral,λ],λ,Simplify]

CP=CharacteristicPolynomial[temp,λ];

crit=Numerator[Factor[CP/.λ->0]]

Reduce[{(crit[[5]]/.i11sol[[1]])>0,0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//FullSimplify

temp2=Numerator[Factor[CP]]

Solve[temp2[[3]]==0,λ]//Simplify

temp3=Solve[temp2[[4]]==0,λ]//Simplify

Reduce[{(λ/.temp3[[1]])<0,0<Y<X,0<δ<b,0<α,0<γ,0<κ,0<i[1,1,t]}]

Reduce[{(λ/.temp3[[2]])<0,0<Y<X,0<δ<b,0<α,0<γ,0<κ,0<i[1,1,t]}]

Single Host Single Pathogen (Mismatching)

Solving for equilibrium

assump={(*s[1,t]->0,*)s[2,t]->0,i[1,1,t]->0,(*i[1,2,t]->0,*)i[2,1,t]->0,i[2,2,t]->0};

temp=equList/.assump

s1sol=Flatten[Solve[temp[[4]]==0,s[1,t]]]

Factor[temp/.s1sol][[1]]

This is a square in i[1,1,t].  Solving for a variable may be useful instead

bSol=Flatten[FullSimplify[Solve[Factor[temp/.s1sol][[1]]==0,b]]]

i12sol=FullSimplify[Solve[Factor[temp/.s1sol][[1]]==0,i[1,2,t]],Assumptions->{0<Y<X,0<δ<b,0<α,0<γ,0<κ}]

equ1=Join[s1sol/.i12sol[[1]],i12sol[[1]],assump]

equ2=Join[s1sol/.i12sol[[1]],i12sol[[1]],assump]

Biological Validity

Reduce[{0<(i[1,2,t]/.i12sol[[1]]),0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//FullSimplify

Reduce[{0<(i[1,2,t]/.i12sol[[2]]),0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//Simplify

Stability

jMtrx[[{1,4},{1,4}]]/.{s[2,t]->0,i[1,1,t]->0,i[2,1,t]->0,i[2,2,t]->0}

equilibraMM={s[1,t]->(((α+γ+δ) κ)/Y),i[1,1,t]->(1/(2 b Y))(-Y (α+δ)+b (Y-2 (α+γ+δ))+(Sqrt[Y ((b^2) Y+Y (((α+δ))^2)+b (-2 Y (α+δ)+4 α (α+γ+δ)))])) κ}

jMtrx[[{1,3},{1,3}]]

temp=jMtrx/.assump/.s1sol/.bSol//Simplify;

jmtrxInteral=jMtrx[[{1,3},{1,3}]]/.equ1//Simplify

Simplify[jmtrxInteral//Eigenvalues]/.{b->1,X->6,δ->0.5,α->1.5,γ->0.1}

Simplify[jmtrxInteral//Eigenvalues]/.{b->1,X->7,δ->0.5,α->1.5,γ->0.1}

Simplify[jmtrxInteral//Eigenvectors]/.{b->1,X->7,Y->3.5,δ->0.05,α->0.1,γ->0.1}

Collect[CharacteristicPolynomial[jmtrxInteral,λ],λ,Simplify]

CP=CharacteristicPolynomial[temp,λ];

crit=Numerator[Factor[CP/.λ->0]]

Reduce[{(crit[[5]]/.i11sol[[1]])>0,0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//FullSimplify

temp2=Numerator[Factor[CP]]

Solve[temp2[[3]]==0,λ]//Simplify

temp3=Solve[temp2[[4]]==0,λ]//Simplify

Reduce[{(λ/.temp3[[1]])<0,0<Y<X,0<δ<b,0<α,0<γ,0<κ,0<i[1,1,t]}]

Reduce[{(λ/.temp3[[2]])<0,0<Y<X,0<δ<b,0<α,0<γ,0<κ,0<i[1,1,t]}]

Single host two pathogens- No solutions

temp=equList/.{(*s[1,t]->0,*)s[2,t]->0,(*i[1,1,t]->0,*)(*i[1,2,t]->0,*)i[2,1,t]->0,i[2,2,t]->0}//Factor

s1sol=Flatten[Solve[temp[[3]]==0,s[1,t]]]

temp[[4]]/.s1sol//Factor

Implies i[1,2,t]=0 which is not allowed

Two host one pathogen- No solution

temp=equList/.{(*s[1,t]->0,*)(*s[2,t]->0,*)(*i[1,1,t]->0,*)i[1,2,t]->0,(*i[2,1,t]->0,*)i[2,2,t]->0}//Factor

i11sol=Flatten[Simplify[Solve[temp[[3]]==0,i[1,1,t]]]]

Factor[temp[[5]]/.i11sol]

s2sol=Flatten[Simplify[Solve[Factor[temp[[5]]/.i11sol]==0,s[2,t]]]]

All the variables are squares so solving for a parameter.

bSol=Flatten[Solve[Numerator[Factor[temp/.i11sol/.s2sol][[2]]]==0,b]]

Factor[temp/.i11sol/.s2sol/.bSol]

Implies i[2,1,t]=0 or s[1,t]=0 which are not allowed.

All variables

Solving for Equilibrium

i22sol=Flatten[Simplify[Solve[Factor[equList[[6]]]==0,i[2,2,t]]]]

s2sol=Simplify[Flatten[Solve[Factor[equList[[5]]/.i22sol]==0,s[2,t]]]]

i11sol=Flatten[Simplify[Solve[Numerator[Factor[equList[[3]]/.i22sol/.s2sol]]==0,i[1,1,t]]]]

s1sol=Simplify[Flatten[Solve[Numerator[Factor[equList[[4]]/.i22sol/.s2sol/.i11sol]]==0,s[1,t]]]]

All variables are squares so solving for a parameter

bSol=Flatten[Simplify[Solve[Numerator[Factor[equList[[2]]/.i22sol/.s2sol/.i11sol/.s1sol]]==0,b]]]

i12sol=Flatten[Solve[Numerator[Factor[equList[[1]]/.i22sol/.s2sol/.i11sol/.s1sol/.bSol]]==0,i[1,2,t]]]

i21sol=Simplify[Solve[b==(b/.bSol/.i12sol),i[2,1,t]],Assumptions->{0<Y<X,0<δ<b,0<α,0<γ,0<κ}]

equ1=Join[i21sol[[1]],i12sol/.i21sol[[1]],s1sol/.i12sol/.i21sol[[1]],i11sol/.s1sol/.i12sol/.i21sol[[1]],s2sol/.i11sol/.s1sol/.i12sol/.i21sol[[1]],i22sol/.s2sol/.i11sol/.s1sol/.i12sol/.i21sol[[1]]]//Simplify;

equ2=Join[i21sol[[2]],i12sol/.i21sol[[2]],s1sol/.i12sol/.i21sol[[2]],i11sol/.s1sol/.i12sol/.i21sol[[2]],s2sol/.i11sol/.s1sol/.i12sol/.i21sol[[2]],i22sol/.s2sol/.i11sol/.s1sol/.i12sol/.i21sol[[2]]]//Simplify;

test={b->1,(*κ->3000,*)δ->0.5,α->0.4,X->6,Y->5,γ->0.75};

s[1,t]+s[2,t]+i[1,2,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]/.equ1/.test

Solve[(s[1,t]+s[2,t]+i[1,2,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]/.equ1/.test)==500,κ]

Solve[(i[1,2,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]/.equ1/.test)==250,κ]

equ1//FullSimplify

equEnd={i[2,1,t]->(1/(4 b (((X+Y))^2)))Y (b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+√((X+Y) ((b^2) (X+Y)+(X+Y) (((α+δ))^2)-2 b (X (α+δ)+Y (α+δ)-4 α (α+γ+δ))))) κ,i[1,2,t]->(1/(4 b (((X+Y))^2)))Y (b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+√((X+Y) ((b^2) (X+Y)+(X+Y) (((α+δ))^2)-2 b (X (α+δ)+Y (α+δ)-4 α (α+γ+δ))))) κ,s[1,t]->(((α+γ+δ) κ)/(X+Y)),i[1,1,t]->(1/(4 b (((X+Y))^2)))X (b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+√((X+Y) ((b^2) (X+Y)+(X+Y) (((α+δ))^2)-2 b (X (α+δ)+Y (α+δ)-4 α (α+γ+δ))))) κ,s[2,t]->(((α+γ+δ) κ)/(X+Y)),i[2,2,t]->(1/(4 b (((X+Y))^2)))X (b X+b Y-4 b α-X α-Y α-4 b γ-4 b δ-X δ-Y δ+√((X+Y) ((b^2) (X+Y)+(X+Y) (((α+δ))^2)-2 b (X (α+δ)+Y (α+δ)-4 α (α+γ+δ))))) κ};

Biological Validity

temp=Reduce[{0<(s[1,t]/.equ1),0<(s[2,t]/.equ1),0<(i[1,1,t]/.equ1),0<(i[1,2,t]/.equ1),0<(i[2,1,t]/.equ1),0<(i[2,2,t]/.equ1),0<Y<X,0<δ<b,0<α,0<γ,0<κ},Y]//FullSimplify

Let ν be the average transmission rate (ν=((X+Y)/2) hence X=2ν-Y).  Given that Y<X we also have that Y<ν

Reduce[Join[{0<(s[1,t]/.equ1),0<(s[2,t]/.equ1),0<(i[1,1,t]/.equ1),0<(i[1,2,t]/.equ1),0<(i[2,1,t]/.equ1),0<(i[2,2,t]/.equ1),0<Y<X,0<δ<b,0<α,0<γ,0<κ}/.{X->2 ν-Y},{Y<ν}]]//FullSimplify

Reduce[{0<(i[1,1,t]/.equ2),0<(i[1,2,t]/.equ2),0<(i[2,1,t]/.equ2),0<(i[2,2,t]/.equ2),0<Y<X,0<δ<b,0<α,0<γ,0<κ}]//Simplify

Allele Frequency

Pathogen AF

Simplify[((i[1,1,t]+i[2,1,t])/(i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]))/.equ2]

Host AF

Simplify[((s[1,t]+i[1,1,t]+i[2,1,t])/(s[1,t]+s[2,t]+i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]))/.equ2]

Stability

βSub={β[1,1]->(X/κ),β[1,2]->(Y/κ),β[2,1]->(Y/κ),β[2,2]->(X/κ)};

An alternative set of conditions that satisfy the endemic equilibrium is:

equSub={s[1,t]->(H/2),s[2,t]->(H/2),i[1,1,t]-> (IX/2),i[1,2,t]-> (IY/2),i[2,1,t]-> (IY/2),i[2,2,t]->(IX/2)};

bequ=Flatten[Solve[(f1[1,t]/.βSub/.equSub)==0,b]]

Xequ=Flatten[Solve[(f2[1,1,t]/.βSub/.equSub)==0,X]]//Simplify

Yequ=Flatten[Solve[(f2[1,2,t]/.βSub/.equSub)==0,Y]]//Simplify

jMtrx1=Simplify[jMtrx/.equSub];
CP=CharacteristicPolynomial[Simplify[jMtrx1],λ]//Factor//Numerator;

The characteristic polynomial consists of 3 factors.  The first one is boring

Playing with the second polynomial

temp=Collect[CP[[2]]/.bequ/.Xequ/.Yequ//Simplify,λ,FullSimplify]

coeffs=CoefficientList[(temp/(Coefficient[temp,(λ^2)])),λ](*/.bequ*)//FullSimplify

Using the R-H conditions

Reduce[{coeffs[[1]]>0,coeffs[[2]]>0,0<α,0<γ,0<δ,0<κ,0<H,0<IY<IX,H+IX+IY<κ},κ]//FullSimplify

This polynomial always has all negative roots. (Given biologically valid)

Playing with the third polynomial

temp=Collect[CP[[3]]/.bequ/.Xequ/.Yequ,λ,FullSimplify]

coeffs=CoefficientList[(temp/(Coefficient[temp,(λ^3)])),λ](*/.bequ*)//FullSimplify

Using the R-H conditions

Reduce[{Simplify[coeffs[[1]]]>0,Simplify[coeffs[[3]]]>0,Simplify[coeffs[[3]]coeffs[[2]]]>Simplify[coeffs[[1]]],0<α,0<γ,0<δ<b,0<κ,0<H,0<IY<IX},κ]

The roots of this polynomial are always negative. (Given biologically valid)

Hence the conditions for stability are equivalent to those for biological validity.

Special case of α->∞, X->∞,Y->∞,γ->0: Stability of the system as delay gets small.

sSub=Flatten[Solve[equList[[{3,5}]]=={0,0},{s[1,t],s[2,t]}]]

temp=Flatten[Solve[equList[[{3,4,5,6}]]=={0,0,0,0},{s[1,t],s[2,t],i[1,1,t],i[2,2,t]}]]

equList[[{1,2}]]/.temp/.{i[1,2,t]->iTot-i[2,1,t]}//Simplify

i21sol=Flatten[Solve[%[[1]]==0,i[2,1,t]]]

lastequ=Factor[equList[[2]]/.temp /.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol]

We seek a limit such that infections last an infinitesimally short time but there remains selection.  This choice won’t give a finite number of infected individuals in the limit:

Collect[Numerator[lastequ]/.{α->Α f,X->x f,Y->y f,(*δ->0,*)γ->γ f},f,Factor]

as the leading order term is -(f^5) iTot y (((x+y))^3) Α κ and requires iTot goes to zero in the limit as f goes to infinity.  Possibilities that appear to work and that retain selection are:

Collect[Numerator[lastequ]/.{α->Α,X->x f,Y->y f,(*δ->0,*)γ->γ},f,Factor]

Solve[(b iTot x+b iTot y-b y κ+y Α κ+y δ κ)==0,iTot]//Simplify

{s[1,t],s[2,t],i[1,1,t],i[1,2,t],i[2,1,t],i[2,2,t]}/.temp/.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol/.iTot->((Y (b-α-δ) κ)/(b (X+Y)))//Simplify

where both clearance and selection are constant (is this really allowed??  not clear that the above solution is permissible as X and Y get large as some of the terms).

Alternatively, if clearance also gets large:

Collect[Numerator[lastequ]/.{α->Α,X->x f,Y->y f,(*δ->0,*)γ->γ f},f,Factor]

iTotsol=Solve[ (-b (iTot^2) (x^4)-4 b (iTot^2) (x^3) y-6 b (iTot^2) (x^2) (y^2)-4 b (iTot^2) x (y^3)-b (iTot^2) (y^4)+b iTot (x^3) y κ+3 b iTot (x^2) (y^2) κ+3 b iTot x (y^3) κ+b iTot (y^4) κ-iTot (x^3) y Α κ-3 iTot (x^2) (y^2) Α κ-3 iTot x (y^3) Α κ-iTot (y^4) Α κ-4 b iTot (x^2) y γ κ-8 b iTot x (y^2) γ κ-4 b iTot (y^3) γ κ-iTot (x^3) y δ κ-3 iTot (x^2) (y^2) δ κ-3 iTot x (y^3) δ κ-iTot (y^4) δ κ+2 b x (y^2) γ (κ^2)+2 b (y^3) γ (κ^2)-4 b (y^2) (γ^2) (κ^2)-2 x (y^2) γ δ (κ^2)-2 (y^3) γ δ (κ^2))==0,iTot]//Simplify

which is the expected disease equilibrium but with γ assumed large relative to α and δ.

jMtrxTemp=Simplify[jMtrx/.temp/.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol];

CharacteristicPolynomial[jMtrxTemp/.{α->Α,X->x f,Y->y f,(*δ->0,*)γ->γ f},λ]

Factor[CharacteristicPolynomial[jMtrxTemp/.{α->Α,X->x f,Y->y f,(*δ->0,*)γ->γ f},λ]];

%232

Numerator[%232][[2]]

Numerator[%232][[3]]

Coefficient[%287,(f^7)]//Factor

Solve[%==0,λ]

Collect[((-2 b iTot (x^2)-4 b iTot x y-2 b iTot (y^2)+b x y κ+b (y^2) κ-x y Α κ-(y^2) Α κ-4 b y γ κ-x y δ κ-(y^2) δ κ)/(y (x+y) κ)),{iTot,b},Factor]

Collect[((-2 b iTot (x^2)-4 b iTot x y-2 b iTot (y^2)+b x y κ+b (y^2) κ-x y Α κ-(y^2) Α κ-4 b y γ κ-x y δ κ-(y^2) δ κ)/(y (x+y) κ)),{b, κ},Factor]

Solve[%==0,iTot]//Simplify

(-b (iTot^2) (x^4)-4 b (iTot^2) (x^3) y-6 b (iTot^2) (x^2) (y^2)-4 b (iTot^2) x (y^3)-b (iTot^2) (y^4)+b iTot (x^3) y κ+3 b iTot (x^2) (y^2) κ+3 b iTot x (y^3) κ+b iTot (y^4) κ-iTot (x^3) y Α κ-3 iTot (x^2) (y^2) Α κ-3 iTot x (y^3) Α κ-iTot (y^4) Α κ-4 b iTot (x^2) y γ κ-8 b iTot x (y^2) γ κ-4 b iTot (y^3) γ κ-iTot (x^3) y δ κ-3 iTot (x^2) (y^2) δ κ-3 iTot x (y^3) δ κ-iTot (y^4) δ κ+2 b x (y^2) γ (κ^2)+2 b (y^3) γ (κ^2)-4 b (y^2) (γ^2) (κ^2)-2 x (y^2) γ δ (κ^2)-2 (y^3) γ δ (κ^2))/.iTot->((y (b (x+y-4 γ)-(x+y) (Α+δ)) κ)/(2 b (((x+y))^2)))//Factor

Collect[CoefficientList[Numerator[%232][[3]],λ][[-1]],f]

temp3=Factor[(CoefficientList[Numerator[%232][[3]],f][[-1]]/.iTotsol[[1]])]

Solve[temp3==0,λ]

Depletion of Susceptible.

Simplify[equEnd/.α->0,{0<Y<X,0<δ<b}]

FullSimplify[jMtrx/.equEnd/.α->0,{0<Y<X,0<δ<b}]

FullSimplify[Eigenvalues[%336],{0<Y<X,0<δ<b}]

Simplify[Eigensystem[%336],{0<Y<X,0<δ<b}]

The second and fourth eigenvalues determine susceptible host population size.  They both become more negative as b->∞ and δ->0.

Eqtemp=Flatten[Solve[equList[[{3,4,5,6}]]=={0,0,0,0},{s[1,t],s[2,t],i[1,1,t],i[2,2,t]}]]

equList[[{1,2}]]/.Eqtemp/.{i[1,2,t]->iTot-i[2,1,t]}//Simplify

i21sol=Flatten[Solve[%344[[1]]==0,i[2,1,t]]]//Simplify

lastequ=Factor[equList[[2]]/.Eqtemp /.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol]

Collect[Numerator[lastequ]/.{b->b f,δ->(δ/f)},f,Factor]

Solve[(iTot (X^2)+2 iTot X Y+iTot (Y^2)+2 Y α κ+2 Y γ κ)==0,iTot]

iTotBV=Solve[(iTot (X^2)+2 iTot X Y+iTot (Y^2)-X Y κ-(Y^2) κ+2 Y α κ+2 Y γ κ)==0,iTot]//Simplify//Flatten

CP=CharacteristicPolynomial[jMtrx/.Eqtemp/.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol/.iTotBV/.{b->b f,δ->(δ/f)}//Simplify,λ]

CP2=Numerator[Factor[CP]]

temp=Collect[CoefficientList[CP2[[3]],f][[-1]],λ,Simplify]

Solve[(temp/.{b->10000,X->0.5,Y->0.2,α->0.1,δ->(1/10000),γ->0.1})==0,λ]

Collect[Numerator[lastequ]/.{b->b f},f,Factor]

Solve[(iTot (X^2)+2 iTot X Y+iTot (Y^2)+2 Y α κ+2 Y γ κ+2 Y δ κ) ==0,iTot]

iTotBV=Solve[(iTot (X^2)+2 iTot X Y+iTot (Y^2)-X Y κ-(Y^2) κ+2 Y α κ+2 Y γ κ+2 Y δ κ)==0,iTot]//Simplify

iTotBV

CP=CharacteristicPolynomial[jMtrx/.Eqtemp/.{i[1,2,t]->iTot-i[2,1,t]}/.i21sol/.iTotBV/.{b->b f,δ->(δ/f)}//Simplify,λ];

CoefficientList[Numerator[CP],f]

CoefficientList[Numerator[CP],f]//Length

temp=CoefficientList[Numerator[CP],f][[-1]]

Playing with the eigenvalue limits numerically

Time Lag: What if X,Y and γ get big

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1,X->0.5f,Y->0.2f,α->0.1,δ->0.1,γ->0.1f}/.{f->1}]//Chop,Re[#1]>Re[#2]&]

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1,X->0.5f,Y->0.2f,α->0.1,δ->0.1,γ->0.1f}/.{f->10000}]//Chop,Re[#1]>Re[#2]&]

Like we predict above this becomes more stable and not complex

{{λ->((-2 b iTot (x^2)-4 b iTot x y-2 b iTot (y^2)+b x y κ+b (y^2) κ-x y Α κ-(y^2) Α κ-4 b y γ κ-x y δ κ-(y^2) δ κ)/(y (x+y) κ))},{λ->((-b iTot (x^2)-2 b iTot x y-b iTot (y^2)+b x y κ+b (y^2) κ-x y Α κ-(y^2) Α κ-2 b y γ κ-x y δ κ-(y^2) δ κ)/(y (x+y) κ))}}/.iTotsol[[1]]/.{Α->α,x->X,y->Y}/.{b->1,X->0.5,Y->0.2,α->0.1,δ->0.1,γ->0.1}

Time Lag 2: What if X,Y, γ, and α get big

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1,X->0.5f,Y->0.2f,α->0.1f,δ->0.1,γ->0.1f}/.{f->1}]//Chop,Re[#1]>Re[#2]&]

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1,X->0.5f,Y->0.2f,α->0.1f,δ->0.1,γ->0.1f}/.{f->10000}]//Chop,Re[#1]>Re[#2]&]

equEnd/.{b->1,X->0.5f,Y->0.2f,α->0.1f,δ->0.1,γ->0.1f}/.{f->10000}

The infection goes away and the equilibrium is more stable.

Depletion of susceptible

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1 f,X->0.5,Y->0.2,α->0.1,δ->(0.1/f),γ->0.1}/.{f->1}]//Chop,Re[#1]>Re[#2]&]

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1 f,X->0.5,Y->0.2,α->0.1,δ->(0.1/f),γ->0.1}/.{f->10000}]//Chop,Re[#1]>Re[#2]&]

equEnd/.{b->1 f,X->0.5,Y->0.2,α->0.1,δ->(0.1/f),γ->0.1}/.{f->100000}

It doesn’t matter what γ is.

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1 f,X->0.5,Y->0.2,α->0.1,δ->(0.1/f),γ->0}/.{f->1}]//Chop,Re[#1]>Re[#2]&]

Sort[Eigenvalues[jMtrx/.equEnd/.{b->1 f,X->0.5,Y->0.2,α->0.1,δ->(0.1/f),γ->0}/.{f->10000}]//Chop,Re[#1]>Re[#2]&]

Eq. population sizes, and generation times

HtotEqu=Sum[s[h,t]/.equEnd,{h,1,2}]+Sum[i[h,p,t]/.equEnd,{h,1,2},{p,1,2}]//FullSimplify

PtotEqu=Sum[(i[h,p,t])/.equEnd,{h,1,2},{p,1,2}]//FullSimplify

δHEqu=(Sum[δ s[h,t],{h,1,2}]+Sum[(δ+α)i[k,p,t],{h,1,2},{k,1,2},{p,1,2}])/.βSub/.equEnd//FullSimplify

δPEqu=(Sum[(δ+α+γ) i[h,p,t],{h,1,2},{p,1,2}])/.βSub/.equEnd//FullSimplify

Stability

Leading Eigenvalue

The general expression for the R0 as a function of the fraction of matching hosts (ϕ)

R0[ϕ_]:=(((X ϕ+Y(1-ϕ))(b-δ))/(b(α+γ+δ)));

The numerical leading eigenvalue of this system is:

λMax[pars_]:=Sort[Chop[Eigenvalues[jMtrx/.equEnd/.pars]],Re[#1]>Re[#2]&][[1]]

Clear[λMaxRe]
λMaxRe[pars_]:=λMaxRe[pars]=If[(R0[(1/2)]/.pars)<1||(X<Y)/.pars,0,Re[Sort[Chop[Eigenvalues[jMtrx/.equEnd/.pars]],Re[#1]>Re[#2]&]][[1]]]

Clear[λMaxIm]
λMaxIm[pars_]:=λMaxIm[pars]=If[(R0[(1/2)]/.pars)<1||(X<Y)/.pars,0,Im[Sort[Chop[Eigenvalues[jMtrx/.equEnd/.pars]],Re[#1]>Re[#2]&]][[1]]]

Stability Region Plots

Critical Lines

The value of R0(ϕ).  When ϕ=1 the plot is a becomes a vertical line and hence not a valid function. It is approximated below by (Subscript[R, 0])(0.999)

R0Plot[δ_,γ_,α_,ϕ_,range_,color_,filling_,fstyle_]:=If[ϕ<1,Plot[-((x δ ϕ+b (α+γ+δ-x ϕ))/((b-δ) (-1+ϕ)))/.{b->1},{x,0,range},PlotRange->{{0,range},{0,range}},PlotStyle->Directive[color,Thickness[0.01]],Filling->filling,FillingStyle->fstyle]]

A plot that excludes the area where X<Y

XYPlot[range_]:=Plot[x,{x,0,range},PlotRange->{{0,range},{0,range}},PlotStyle->Directive[Gray,Thick],Filling->Top,FillingStyle->Black,AspectRatio->1]

Plotting the critical point where ||Im((Subscript[λ, max]))||=0

cyclePlot[δin_,γin_,αin_,range_,color_]:=ContourPlot[λMaxIm[{b->1,δ->δin,γ->γin,α->αin,X->x,Y->y}]==((-1)/350),{x,((2 (αin+ δin+γin))/(1-δin))+(1/100),range},{y,0,x-(1/3)},Contours->{((-1)/350)},ContourStyle->Directive[Thick,color],ContourShading->None,PerformanceGoal->"Quality",PlotRange->{{0,range},{0,range}}]

cyclePlot[(1/2),0,(1/50),3,Black]

Plots

leg=LineLegend[{Directive[Thickness[0.02],Darker[Green]],Directive[Thickness[0.02],Red],Directive[Thickness[0.02],Darker[Purple]],Black},{Style["(Subscript[R, o])(0)=1",Black,16],Style["(Subscript[R, o])(1/2)=1",Black,16],Style["(Subscript[R, o])(1)=1",Black,16],Style["Im[(Subscript[λ, max])]=0",Black,14]},Background->White,LegendFunction->"Frame",LegendMargins->Automatic]

R0Region[range_,δin_,γin_,αin_]:=Show[R0Plot[δin,γin,αin,0,range,Darker[Green],Top,ColorData["GreenPinkTones"][0.2]],R0Plot[δin,γin,αin,0,range,Darker[Green],Bottom,ColorData["GreenPinkTones"][0.4]],R0Plot[δin,γin,αin,0.5,range,Red,Bottom,ColorData["GreenPinkTones"][0.6]],R0Plot[δin,γin,αin,0.99,range,Purple,Bottom,ColorData["GreenPinkTones"][0.8]],XYPlot[range],cyclePlot[δin,γin,αin,range,Black],AspectRatio->1,Frame->True,FrameTicks->None,(*FrameLabel->{{"Mismatching Transmission Y",None},{"Matching Transmission X",None}},*)LabelStyle->Directive[Bold,Black, 16],ImageSize->Large,Epilog->{Text[Style["Not Valid:\n X<Y",18,White],Scaled[{0.2,0.7}]],
Text[Style["Disease \n always spreads",18,Black],Scaled[{0.85,0.7}]],
Text[Style["Coevolution facilitated \n disease spread ",18,Black],Scaled[{0.62,0.27}]],Text[Style["Coevolution \n driven \n disease free",18,Black],Scaled[{0.46,0.08}]],Text[Style["Always \n disease free",18,Black],Scaled[{0.2,0.08}]],Text[Style["No Cycles",14,Black,Italic],Scaled[{0.9,0.15}]],Inset[leg,Scaled[{0.15,0.4}]]}]

R0Region[3,0.5,0,0.02]

Export[NotebookDirectory[]<>"Plots/Region.pdf",R0Region[3,0.5,0,0.02],ImageResolution->500];

Neutral Expectation

Given a sample path (s[1,t], s[2,t], i[1,1,t], i[1,2,t] i[2,1,t],i[2,2,t]) we can derive a differential for the heterozygosity at a neutral (unlinked) locus. We do so by considering the effect of each the 6 possible events (susceptible birth, infected birth, susceptible death, infected death, transmission, and recovery) on the allele frequency at the neutral locus and hence on the heterozygosity at the neutral locus.  Since the population is structured into susceptible and infected compartments we calculate the population level heterozygosity by first calculating the dynamics within each component (Hs[t] and Hi[t]) as well as between components Hsi[t].

H[t]=((((Ns/(Ns+Ni))))^2)Hi[t]+2((Ns/(Ns+Ni)))((Ni/(Ns+Ni)))Hsi[t]+((((Ni/(Ns+Ni))))^2)Hi[t]
Ns=s[1,t]+s[2,t]

Letting the rate of susceptible births be λS, infected births λI, susceptible deaths μS, infected deaths μI, and transmission events be B, the resulting differential equations are given by:

((dHs[t])/dt)=λS ()+λI()+μS()+B ()
((dHi[t])/dt)=μI()+B ()
((dHsi[t])/dt)=λI()+B ()

ratesSub={λS->b(s[1,t]+s[2,t])(1-((s[1,t]+s[2,t]+i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t])/κ)),λI->b(i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t])(1-((s[1,t]+s[2,t]+i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t])/κ)),μS->δ(s[1,t]+s[2,t]),μI->(δ+α)(i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t]),Β-> (X/κ)(s[1,t](i[1,1,t]+i[1,2,t])+s[2,t](i[2,1,t]+i[2,2,t]))+(Y/κ)(s[2,t](i[1,1,t]+i[1,2,t])+s[1,t](i[2,1,t]+i[2,2,t])),Γ->γ(i[1,1,t]+i[1,2,t]+i[2,1,t]+i[2,2,t])};

Effect of individual events

Heterozygosity within the susceptible population

Susceptible Birth

temp=(pS 2pS2 qS2/.qS2->1-pS2/.pS2->((Ns pS+1)/((Ns+1))))+((1-pS)2pS2 qS2/.qS2->1-pS2/.pS2->((Ns pS)/((Ns+1))))-2 pS(1-pS)//Simplify;
Hsλs=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}]

Infected Birth

temp=(pI 2pS2 qS2/.qS2->1-pS2/.pS2->((Ns pS+1)/((Ns+1))))+((1-pI)2pS2 qS2/.qS2->1-pS2/.pS2->((Ns pS)/((Ns+1))))-2 pS(1-pS)//FullSimplify;
Hsλi=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//FullSimplify

Susceptible Death

temp=(pS 2 pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS-1)/((Ns-1))) )+((1-pS)2pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS)/((Ns-1))))-2 pS(1-pS)//FullSimplify;
Hsμs=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}]

Infected Death- No effect

Transmission

temp=(pS 2 pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS-1)/((Ns-1))) )+((1-pS)2pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS)/((Ns-1))))-2 pS(1-pS)//FullSimplify;
HsΒ=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}]

Recovery

temp=(pI 2 pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS+1)/((Ns+1))) )+((1-pI)2pS2 qS2/.qS2-> 1-pS2/.pS2->((Ns pS)/((Ns+1))))-2 pS(1-pS)//FullSimplify;
HsΓ=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//FullSimplify

Heterozygosity within the infected population

Susceptible Birth - No effect

Infected Birth- No effect

Susceptible Death- No effect

Infected Death

temp=(pI 2 pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI-1)/((Ni-1))) )+((1-pI)2pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI)/((Ni-1))))-2 pI(1-pI)//FullSimplify;
Hiμi=FullSimplify[Expand[temp]/.{(pI^2)-> ((2pI-Hi[t])/2)}]

Transmission

temp=(pS 2pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI+1)/((Ni+1))))+((1-pS)2pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI)/((Ni+1))))-2 pI(1-pI)//FullSimplify;
HiΒ=FullSimplify[Expand[temp]/.{(pI^2)-> ((2pI-Hi[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//Simplify

Recovery

temp=(pI 2 pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI-1)/((Ni-1))) )+((1-pI)2pI2 qI2/.qI2-> 1-pI2/.pI2->((Ni pI)/((Ni-1))))-2 pI(1-pI)//FullSimplify;
HiΓ=FullSimplify[Expand[temp]/.{(pI^2)-> ((2pI-Hi[t])/2)}]

Cross Heterozygosity

Susceptible Birth- No effect

Infected Birth

temp=(pI (pS2(1-pI)+pI(1-pS2))/.pS2->((Ns pS+1)/((Ns+1))))+((1-pI) (pS2(1-pI)+pI(1-pS2))/.pS2->((Ns pS)/((Ns+1))))-pS(1-pI)-pI(1-pS)//Simplify;
Hsiλi=FullSimplify[Expand[temp]/.{(pI^2)-> ((2pI-Hi[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//Simplify

Susceptible Death- No effect

Infected Death- No effect

Transmission

temp=(pS (pS2(1-pI2)+pI2(1-pS2))/.pS2->((Ns pS-1)/((Ns-1)))/.pI2->((Ni pI+1)/((Ni+1))))+((1-pS) (pS2(1-pI2)+pI2(1-pS2))/.pS2->((Ns pS)/((Ns-1)))/.pI2->((Ni pI)/((Ni+1))))-pS(1-pI)-pI(1-pS)//Simplify;
HsiΒ=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}/.{(pI^2)-> ((2pI-Hi[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//FullSimplify

Recovery

temp=(pI (pS2(1-pI2)+pI2(1-pS2))/.pS2->((Ns pS+1)/((Ns+1)))/.pI2->((Ni pI-1)/((Ni-1))))+((1-pI) (pS2(1-pI2)+pI2(1-pS2))/.pS2->((Ns pS)/((Ns+1)))/.pI2->((Ni pI)/((Ni-1))))-pS(1-pI)-pI(1-pS)//Simplify;
HsiΓ=FullSimplify[Expand[temp]/.{(pS^2)-> ((2pS-Hs[t])/2)}/.{(pI^2)-> ((2pI-Hi[t])/2)}]/.{pI->((-Hsi[t]+pS)/(-1+2 pS))}//FullSimplify

Differential Equations

odes={D[Hs[t],t]==λs Hsλs+ λi Hsλi+μs Hsμs+Β HsΒ+Γ HsΓ,D[Hi[t],t]== μi Hiμi+Β HiΒ+Γ HiΓ,D[Hsi[t],t]==λi Hsiλi+ Β HsiΒ+Γ HsiΓ,Hs[0]==0.5,Hi[0]==0.5,Hsi[0]==0.5}(*/.{λs->b sTot[t](1-((sTot[t]+iTot[t])/κ)),λi->b iTot[t](1-((sTot[t]+iTot[t])/κ)),μs->δ sTot[t],μi->(δ+α) iTot[t],Β->β sTot[t]iTot[t]}*)

{D[Hs[t],t]==λs Hsλs+ λi Hsλi+μs Hsμs+Β HsΒ+Γ HsΓ,D[Hi[t],t]== μi Hiμi+Β HiΒ+Γ HiΓ,D[Hsi[t],t]==λi Hsiλi+ Β HsiΒ+Γ HsiΓ,Hs[0]==0.5,Hi[0]==0.5,Hsi[0]==0.5}/.{λs->b sTot[t](1-((sTot[t]+iTot[t])/κ)),λi->b iTot[t](1-((sTot[t]+iTot[t])/κ)),μs->δ sTot[t],μi->(δ+α) iTot[t],Β->β sTot[t]iTot[t],Γ-> γ iTot[t]}

Simulations

λ[s1_,s2_,i1_,i2_]:=Max[0,b (1-((s1+s2+i1+i2)/κ))];
μS[s1_,s2_,i1_,i2_]:=δ ;
μI[s1_,s2_,i1_,i2_]:=(δ+α);
Β[s1_,s2_,i1_,i2_]:=β(i1+i2)
Γ[s1_,s2_,i1_,i2_]:=γ

TotalRate[s1_,s2_,i1_,i2_,pars_]:=λ[s1,s2,i1,i2](s1+s2+i1+i2)+μS[s1,s2,i1,i2](s1+s2)+μI[s1,s2,i1,i2](i1+i2)+Β[s1,s2,i1,i2](s1+s2)+Γ[s1,s2,i1,i2](i1+i2)/.pars
samp[s1_,s2_,i1_,i2_,pars_]:=Block[{out,Δt},
Δt=RandomVariate[ExponentialDistribution[TotalRate[s1,s2,i1,i2,pars]]];
RandomChoice[{λ[s1,s2,i1,i2](s2+i2)/.pars,λ[s1,s2,i1,i2](s1+i1)/.pars,μS[s1,s2,i1,i2]s2/.pars,μS[s1,s2,i1,i2]s1/.pars,μI[s1,s2,i1,i2]i2/.pars,μI[s1,s2,i1,i2]i1/.pars,Β[s1,s2,i1,i2]s2/.pars,Β[s1,s2,i1,i2]s1/.pars,Γ[s1,s2,i1,i2]i2/.pars,Γ[s1,s2,i1,i2]i1/.pars}-> 
{{Δt,0,1,0,0},{Δt,1,0,0,0},{Δt,0,-1,0,0},{Δt,-1,0,0,0},{Δt,0,0,0,-1},{Δt,0,0,-1,0},{Δt,0,-1,0,1},{Δt,-1,0,1,0},{Δt,0,1,0,-1},{Δt,1,0,-1,0}}]]

sim[s0_,i0_,nRep_,pars_]:=Module[{out,OUT,e,r},OUT={};
For[r=1,r≤nRep,r++,
out={{0,Round[(s0/2)],Round[(s0/2)],Round[(i0/2)],Round[(i0/2)]}};
While[out[[-1,1]]<tMax/.pars,
If[Total[out[[-1,{2,3,4,5}]]]==0,AppendTo[out,{tMax/.pars,0,0}],
AppendTo[out,out[[-1]]+samp[out[[-1,2]],out[[-1,3]],out[[-1,4]],out[[-1,5]],pars]]
]
];
AppendTo[OUT,out]
];
OUT
]

Pars={b->1,δ->0.1,α->0.1,β->0.01,γ->0.1,κ->500,tMax->1};

ExampleSims=sim[100,30,2000,Pars];
Clear[NDyn,AFDyn]
NDyn[vec_]:=NDyn[vec]=Map[Map[{#[[1]],Total[#[[vec]]]}&,#]&,ExampleSims];
AFDyn[vec1_,vec2_]:=AFDyn[vec1,vec2]=Map[Map[{#[[1]],((Total[#[[vec1]]])/(Total[#[[vec1]]]+Total[#[[vec2]]]))//N}&,#]&,ExampleSims];

Clear[AFDynTab]
AFDynTab[npts_,AFDyn_,pars_]:=AFDynTab[npts,AFDyn,pars]=Block[{f,e,r,out,OUT},OUT={};
For[r=1,r≤Length[AFDyn],r++,
out={{0,0.5}};
f=Interpolation[AFDyn[[r]],InterpolationOrder->0];
For[e=1,e≤npts,e++,
AppendTo[out,{e ((tMax/.pars)/npts),f[e ((tMax/.pars)/npts)]}];
];
AppendTo[OUT,out];
];
OUT//N]

GV[npts_,AFDyn_,pars_]:=Map[Map[{#[[1]],2#[[2]](1-#[[2]])}&,#]&, AFDynTab[npts,AFDyn,pars]]

Clear[GVExp];(*Clear[finalnum];*)
GVExp[npts_,NDyn_,pars_]:=GVExp[npts,NDyn,pars]=Block[{r,sTot,iTot,nSol,HsSol,HiSol,HsiSol,out},out={};
For[r=1,r≤Min[50,Length[NDyn]](*Approximation for speed purposes*),r++,
sTot=Interpolation[Map[{#[[1]],#[[2]]+#[[3]]}&,ExampleSims[[r,;;]]],InterpolationOrder->0];
iTot=Interpolation[Map[{#[[1]],#[[4]]+#[[5]]}&,ExampleSims[[r,;;]]],InterpolationOrder->0];
nSol=Flatten[NDSolve[{Hs'[t]==((γ (-(1+2 Ns) Hs[t]+2 Ns Hsi[t]) iTot[t])/(((1+Ns))^2))-((δ Hs[t] sTot[t])/(((-1+Ns))^2))-((β Hs[t] iTot[t] sTot[t])/(((-1+Ns))^2))+((b (-(1+2 Ns) Hs[t]+2 Ns Hsi[t]) iTot[t] (1-((iTot[t]+sTot[t])/κ)))/(((1+Ns))^2))-((b Hs[t] sTot[t] (1-((iTot[t]+sTot[t])/κ)))/(((1+Ns))^2)),Hi'[t]==-((γ Hi[t] iTot[t])/(((-1+Ni))^2))-(((α+δ) Hi[t] iTot[t])/(((-1+Ni))^2))+((β (-(1+2 Ni) Hi[t]+2 Ni Hsi[t]) iTot[t] sTot[t])/(((1+Ni))^2)),Hsi'[t]==((γ (Ni Hi[t]+Hsi[t]-Ni Hsi[t]) iTot[t])/((-1+Ni) (1+Ns)))+((β (Ns Hs[t]+Hsi[t]-Ns Hsi[t]) iTot[t] sTot[t])/((1+Ni) (-1+Ns)))+((b (Hi[t]-Hsi[t]) iTot[t] (1-((iTot[t]+sTot[t])/κ)))/(1+Ns)),Hs[0]==0.5,Hi[0]==0.5,Hsi[0]==0.5}/.{Ns->sTot[t],Ni->iTot[t]}/.pars,{Hs[t],Hi[t],Hsi[t]},{t,0,tMax/.pars}]];
Clear[HsSol,HiSol,HsiSol];
HsSol[τ_]:=HsSol[τ]=Hs[t]/.nSol/.{t->τ};
HiSol[τ_]:=HiSol[τ]=Hi[t]/.nSol/.{t->τ};
HsiSol[τ_]:=HsiSol[τ]=Hsi[t]/.nSol/.{t->τ};
AppendTo[out,Table[{τ,(((((sTot[τ])/(sTot[τ]+iTot[τ]))))^2)HsSol[τ]+(((((iTot[τ])/(sTot[τ]+iTot[τ]))))^2)HiSol[τ]+2(((sTot[τ])/(sTot[τ]+iTot[τ])))(((iTot[τ])/(sTot[τ]+iTot[τ])))HsiSol[τ](*((1/2) (1-(Sqrt[1-2 HiSol[τ]]) (Sqrt[1-2 HsSol[τ]])))*)},{τ,0,tMax/.pars, ((tMax/.pars)/npts)}]];
];
out
]

Show[ListLinePlot[{Mean[GV[40,AFDyn[{2,4},{3,5}],Pars]],Mean[GV[40,AFDyn[{2,4},{3,5}],Pars]]+2((StandardDeviation[GV[40,AFDyn[{2,4},{3,5}],Pars]])/(Sqrt[Length[GV[40,AFDyn[{2,4},{3,5}],Pars]]])),Mean[GV[40,AFDyn[{2,4},{3,5}],Pars]]-2((StandardDeviation[GV[40,AFDyn[{2,4},{3,5}],Pars]])/(Sqrt[Length[GV[40,AFDyn[{2,4},{3,5}],Pars]]]))},PlotStyle->{Black,Gray,Gray},Filling->2->{3}],ListLinePlot[{Mean[GVExp[40,NDyn[{2,3,4,5}],Pars]],Mean[GVExp[40,NDyn[{2,3,4,5}],Pars]]+2((StandardDeviation[GVExp[40,NDyn[{2,3,4,5}],Pars]])/(Sqrt[Length[GVExp[40,NDyn[{2,3,4,5}],Pars]]])),Mean[GVExp[40,NDyn[{2,3,4,5}],Pars]]-2((StandardDeviation[GVExp[40,NDyn[{2,3,4,5}],Pars]])/(Sqrt[Length[GVExp[40,NDyn[{2,3,4,5}],Pars]]]))},PlotStyle->{Blue,LightBlue,LightBlue},Filling->2->{3}],Frame->True,FrameLabel->{"Time","Expected Heterozygosity"}]

Show[ListLinePlot[GV[40,AFDyn[{2,4},{3,5}],Pars],PlotRange->All,PlotStyle->Directive[Gray,Opacity[0.2]]],
ListLinePlot[GVExp[40,NDyn[{2,3,4,5}],Pars],PlotStyle->Directive[Blue, Opacity[0.1]]],ListLinePlot[GVExp[40,NDyn[{2,3,4,5}],Pars]//Mean,PlotStyle->Darker[Blue]],ListLinePlot[Mean[GV[40,AFDyn[{2,4},{3,5}],Pars]],PlotStyle->Black]]

Numerical Integration

ODEs[pars_]:=Flatten[Join[Table[D[s[h,t],t]==f1[h,t],{h,1,2}],Table[D[i[h,p,t],t]==f2[h,p,t],{h,1,2},{p,1,2}]]]/.pars

Inits[pars_]:=Flatten[Join[Table[s[h,0]==Sum[s[k,t],{k,1,2}]*If[h==1,0.8,0.2],{h,1,2}],Table[i[h,p,0]==Sum[i[k,l,t],{k,1,2},{l,1,2}]*If[h==1,0.8,0.2]*If[p==1,0.8,0.2],{h,1,2},{p,1,2}]]]/.equEnd/.pars

vars=Flatten[Join[Table[s[h,t],{h,1,2}],Table[ i[h,p,t],{h,1,2},{p,1,2}]]];

Clear[nsol]
nsol[pars_]:=nsol[pars]=Flatten[NDSolve[Join[ODEs[pars],Inits[pars]],vars,{t,0,1000}]]
pHsol[pars_,τ_]:=((s[1,t]+Sum[i[1,p,t],{p,1,2}])/(Sum[s[h,t]+Sum[i[h,p,t],{p,1,2}],{h,1,2}]))/.nsol[pars]/.{t->τ}
pPsol[pars_,τ_]:=((Sum[i[h,1,t],{h,1,2}])/(Sum[i[h,p,t],{h,1,2},{p,1,2}]))/.nsol[pars]/.{t->τ}

DetAFPlot[file_,run_,rep_,color_]:=Show[ListLinePlot[Table[{pHsol[simPars[file,run,rep],τ],pPsol[simPars[file,run,rep],τ]},{τ,0,500}],PlotStyle->color]/.Line[x_]:>{Arrowheads[Table[.05,{r,{2}}]],Arrow[x]},
ListPlot[{{pHsol[simPars[file,run,rep],0],pPsol[simPars[file,run,rep],0]}},PlotStyle->Red,PlotMarkers->{●,14}],PlotRange->{{0,1},{0,1}},AspectRatio->1,ImageSize->Medium,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,16],Axes->False]

DetAFPlot[fileE,1,3,Black]
Export[NotebookDirectory[]<>"Plots/Epi_Det_AF.pdf",%,ImageResolution->300];

DetHetPlot[file_,run_,rep_,color_]:=Show[ListLinePlot[Table[{τ,2pHsol[simPars[file,run,rep],τ](1-pHsol[simPars[file,run,rep],τ])},{τ,0,500}],PlotStyle->color,PlotRange->{0,0.5}],ListPlot[{{0,2pHsol[simPars[file,run,rep],0](1-pHsol[simPars[file,run,rep],0])}},PlotStyle->Red,PlotMarkers->{●,14}],ImageSize->Medium,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,16],Axes->False]

DetHetPlot[fileE,1,3,Black]
Export[NotebookDirectory[]<>"Plots/Epi_Det_Het.pdf",%,ImageResolution->300];

Ensemble Moment Dynamics

Substituitons

Moment Substitutions

first=Join[Table[δS[h]->μδS[h],{h,1,2}],Table[δI[k,p]->μδI[k,p],{k,1,2},{p,1,2}]]//Flatten;

second=Block[{out},out={};
Table[AppendTo[out,Solve[CδSS[h1,h2]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])]/.δS[h1]δS[h2]->X,X][[1]]/.first/.X-> δS[h1]δS[h2]],{h1,1,2},{h2,h1,2}];
Table[AppendTo[out,Solve[CδSI[h1,p1,k1]==Expand[(δS[h1]-μδS[h1])(δI[p1,k1]-μδI[p1,k1])]/.δS[h1]δI[p1,k1]->X,X][[1]]/.first/.X->δS[h1]δI[p1,k1]],{h1,1,2},{p1,1,2},{k1,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2,Solve[CδII[p1,k1,p2,k2]==Expand[(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])]/.δI[p1,k1]δI[p2,k2]->X,X][[1]]/.first/.X->δI[p1,k1]δI[p2,k2],{}]],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}];
Flatten[out]
];

third=Block[{out},out={};
Table[AppendTo[out,Solve[CδSSS[h1,h2,h3]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])(δS[h3]-μδS[h3])]/.δS[h1]δS[h2]δS[h3]->X,X][[1]]/.second/.first/.X->δS[h1]δS[h2]δS[h3]],{h1,1,2},{h2,h1,2},{h3,h2,2}];
Table[AppendTo[out,Solve[CδSSI[h1,h2,p1,k1]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])(δI[p1,k1]-μδI[p1,k1])]/.δS[h1]δS[h2]δI[p1,k1]->X,X][[1]]/.second/.first/.X->δS[h1]δS[h2]δI[p1,k1]],{h1,1,2},{h2,h1,2},{p1,1,2},{k1,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2,Solve[CδSII[h1,p1,k1,p2,k2]==Expand[(δS[h1]-μδS[h1])(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])]/.δS[h1]δI[p1,k1]δI[p2,k2]->X,X][[1]]/.second/.first/.X->δS[h1]δI[p1,k1]δI[p2,k2],{}]],{h1,1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2≤2p3+k3,Solve[CδIII[p1,k1,p2,k2,p3,k3]==Expand[(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])(δI[p3,k3]-μδI[p3,k3])]/.δI[p1,k1]δI[p2,k2]δI[p3,k3]->X,X][[1]]/.second/.first/.X->δI[p1,k1]δI[p2,k2]δI[p3,k3],{}]],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2}];
Flatten[out]
];

fourth=Block[{out},out={};
Table[AppendTo[out,Solve[CδSSSS[h1,h2,h3,h4]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])(δS[h3]-μδS[h3])(δS[h4]-μδS[h4])]/.δS[h1]δS[h2]δS[h3]δS[h4]->X,X][[1]]/.second/.first/.X->δS[h1]δS[h2]δS[h3]δS[h4]],{h1,1,2},{h2,h1,2},{h3,h2,2},{h4,h3,2}];
Table[AppendTo[out,Solve[CδSSSI[h1,h2,h3,p1,k1]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])(δS[h3]-μδS[h3])(δI[p1,k1]-μδI[p1,k1])]/.δS[h1]δS[h2]δS[h3]δI[p1,k1]->X,X][[1]]/.second/.first/.X->δS[h1]δS[h2]δS[h3]δI[p1,k1]],{h1,1,2},{h2,h1,2},{h3,h2,2},{p1,1,2},{k1,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2,Solve[CδSSII[h1,h2,p1,k1,p2,k2]==Expand[(δS[h1]-μδS[h1])(δS[h2]-μδS[h2])(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])]/.δS[h1]δS[h2]δI[p1,k1]δI[p2,k2]->X,X][[1]]/.second/.first/.X->δS[h1]δS[h2]δI[p1,k1]δI[p2,k2],{}]],{h1,1,2},{h2,h1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2≤2p3+k3,
Solve[CδSIII[h1,p1,k1,p2,k2,p3,k3]==Expand[(δS[h1]-μδS[h1])(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])(δI[p3,k3]-μδI[p3,k3])]/.δS[h1]δI[p1,k1]δI[p2,k2]δI[p3,k3]->X,X][[1]]/.second/.first/.X->δS[h1]δI[p1,k1]δI[p2,k2]δI[p3,k3],{}]],{h1,1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2}];
Table[AppendTo[out,If[2p1+k1≤2p2+k2≤2p3+k3≤2p4+k4,Solve[CδIIII[p1,k1,p2,k2,p3,k3,p4,k4]==Expand[(δI[p1,k1]-μδI[p1,k1])(δI[p2,k2]-μδI[p2,k2])(δI[p3,k3]-μδI[p3,k3])(δI[p4,k4]-μδI[p4,k4])]/.δI[p1,k1]δI[p2,k2]δI[p3,k3]δI[p4,k4]->X,X][[1]]/.second/.first/.X->δI[p1,k1]δI[p2,k2]δI[p3,k3]δI[p4,k4],{}]],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2},{p4,p3,2},{k4,1,2}];
Flatten[out]
];

Simplification

βSub={β[1,1]->(X/κ),β[1,2]->(Y/κ),β[2,1]->(Y/κ),β[2,2]->(X/κ)};

Let δH[i]=H[i]-(OverHat[H[i]]) and δP[i]=P-(OverHat[P[i]]).  We will assume that the deviations from equilibrium, δH and δP, are small.

Note that the covariance between the deviations from equilibrium are equal to the the covariance in the original variables.

eSub=Join[Table[s[h]->s[h,t]+ϵ δS[h],{h,1,2}],Table[i[h,p]->i[h,p,t]+ϵ δI[h,p],{h,1,2},{p,1,2}]]/.equEnd//Flatten//Simplify;

μeSubRev=Flatten[Join[Table[μδS[h]->μS[h,t]-s[h,t],{h,1,2}],Table[μδI[h,p]->μI[h,p,t]-i[h,p,t],{h,1,2},{p,1,2}],Table[CδSS[h1,h2]-> CSS[h1,h2,t],{h1,1,2},{h2,h1,2}],Table[CδSI[h1,p1,k1]-> CSI[h1,p1,k1,t],{h1,1,2},{p1,1,2},{k1,1,2}],Table[If[2p1+k1≤2p2+k2,CδII[p1,k1,p2,k2]-> CII[p1,k1,p2,k2,t],{}],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}],Table[CδSSS[h1,h2,h3]-> CSSS[h1,h2,h3,t],{h1,1,2},{h2,h1,2},{h3,h2,2}],Table[CδSSI[h1,h2,p1,k1]-> CSSI[h1,h2,p1,k1,t],{h1,1,2},{h2,h1,2},{p1,1,2},{k1,1,2}],Table[If[2p1+k1≤2p2+k2,CδSII[h1,p1,k1,p2,k2]-> CSII[h1,p1,k1,p2,k2,t],{}],{h1,1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}],Table[If[2p1+k1≤2p2+k2≤2p3+k3,CδIII[p1,k1,p2,k2,p3,k3]-> CIII[p1,k1,p2,k2,p3,k3,t],{}],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2}]]]/.equEnd//Flatten;

Moment Closure

assump3=Flatten[Join[Table[CδSSSS[h1,h2,h3,h4]->0,{h1,1,2},{h2,h1,2},{h3,h2,2},{h4,h3,2}],Table[CδSSSI[h1,h2,h3,p1,k1]->0,{h1,1,2},{h2,h1,2},{h3,h2,2},{p1,1,2},{k1,1,2}],Table[If[2p1+k1≤2p2+k2,CδSSII[h1,h2,p1,k1,p2,k2]->0,{}],{h1,1,2},{h2,h1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}],Table[If[2p1+k1≤2p2+k2≤2p3+k3,CδSIII[h1,p1,k1,p2,k2,p3,k3]->0,{}],{h1,1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2}],Table[If[2p1+k1≤2p2+k2≤2p3+k3≤2p4+k4,CδIIII[p1,k1,p2,k2,p3,k3,p4,k4]->0,{}],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2},{p4,p3,2},{k4,1,2}]]];

assump2=Flatten[Join[Table[CδSSS[h1,h2,h3]->0,{h1,1,2},{h2,h1,2},{h3,h2,2}],Table[CδSSI[h1,h2,p1,k1]->0,{h1,1,2},{h2,h1,2},{p1,1,2},{k1,1,2}],
Table[If[2p1+k1≤2p2+k2,CδSII[h1,p1,k1,p2,k2]->0,{}],{h1,1,2},{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}],
Table[If[2p1+k1≤2p2+k2≤2p3+k3,CδIII[p1,k1,p2,k2,p3,k3]->0,{}],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2},{p3,p2,2},{k3,1,2}],assump3]];

assump1=Join[Table[CδSS[h1,h2]->0,{h1,1,2},{h2,h1,2}],
Table[CδSI[h1,p1,k1]->0,{h1,1,2},{p1,1,2},{k1,1,2}],
Table[If[2p1+k1≤2p2+k2,CδII[p1,k1,p2,k2]->0,{}],{p1,1,2},{k1,1,2},{p2,p1,2},{k2,1,2}],assump2]//Flatten;

Rates f[x] and change in x

fEpi=Join[(*Host birth*)Table[b(s[h]+Sum[i[h,l],{l,1,2}])(1-((Sum[s[k]+Sum[i[k,l],{l,1,2}],{k,1,2}])/κ)),{h,1,2}],
(*Susceptible Host death*)Table[δ s[h],{h,1,2}],
(*Infection*)Table[β[h,p] s[h]Sum[i[k,p],{k,1,2}],{h,1,2},{p,1,2}],
(*Infected Host Death*)Table[(δ+α) i[h,p],{h,1,2},{p,1,2}],
(*Pathogen Clearance*)Table[γ i[h,p],{h,1,2},{p,1,2}]
]/.βSub//Flatten;

ΔS[h1_]:=Join[(*Host birth*)Table[If[h==h1,1,0],{h,1,2}],
(*Susceptible Host death*)Table[If[h==h1,-1,0],{h,1,2}],
(*Infection*)Table[If[h==h1,-1,0],{h,1,2},{p,1,2}],
(*Infected Host Death*)Table[0,{h,1,2},{p,1,2}],
(*Pathogen Clearance*)Table[If[h==h1,1,0],{h,1,2},{p,1,2}]
]/.βSub//Flatten

ΔI[h1_,p1_]:=Join[(*Host birth*)Table[0,{h,1,2}],
(*Susceptible Host death*)Table[0,{h,1,2}],
(*Infection*)Table[If[h==h1&&p==p1,1,0],{h,1,2},{p,1,2}],
(*Infected Host Death*)Table[If[h==h1&&p==p1,-1,0],{h,1,2},{p,1,2}],
(*Pathogen Clearance*)Table[If[h==h1&&p==p1,-1,0],{h,1,2},{p,1,2}]
]/.βSub//Flatten

Differential Equations

Means

Clear[dμSdt]
dμSdt[h_,t_]:=dμSdt[h,t]=Expand[Normal[Series[fEpi.ΔS[h]/.eSub,{ϵ,0,2}]]]/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev/.{τ->t}

Clear[dμIdt]
dμIdt[h_,p_,t_]:=dμIdt[h,p,t]=Expand[Normal[Series[fEpi.ΔI[h,p]/.eSub,{ϵ,0,2}]]]/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev/.{τ->t}

Variances

Clear[dSSdt]
dSSdt[h1_,h2_,t_]:=dSSdt[h1,h2,t]=Expand[Normal[Series[Simplify[fEpi.Simplify[(s[h1]+ΔS[h1])(s[h2]+ΔS[h2])-s[h1]s[h2]]/.eSub],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev/.{τ->t}

Clear[dSIdt]
dSIdt[h1_,p1_,k1_,t_]:=dSIdt[h1,p1,k1,t]=Expand[Normal[Series[Simplify[fEpi.Simplify[(s[h1]+ΔS[h1])(i[p1,k1]+ΔI[p1,k1])-s[h1]i[p1,k1]]/.eSub],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev/.{τ->t}

Clear[dIIdt]
dIIdt[p1_,k1_,p2_,k2_,t_]:=dIIdt[p1,k1,p2,k2,t]=Expand[Normal[Series[Simplify[fEpi.Simplify[(i[p1,k1]+ΔI[p1,k1])(i[p2,k2]+ΔI[p2,k2])-i[p1,k1]i[p2,k2]]/.eSub],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev/.{τ->t}

Genetic Variance

HEM=Expand[Normal[Series[2 (((s[1]+Sum[i[1,p],{p,1,2}])/(Sum[s[h]+Sum[i[h,p],{p,1,2}],{h,1,2}])))(1-((s[1]+Sum[i[1,p],{p,1,2}])/(Sum[s[h]+Sum[i[h,p],{p,1,2}],{h,1,2}])))/.eSub,{ϵ,0,2}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump2/.μeSubRev

Individual-Based Simulations

Choosing Parameters

R0[(1/2)]//Simplify

Solve[(((X+Y) (b-δ))/(2 b (A+δ)))==1,A]/.b->1//FullSimplify

fPstar=((Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd//Simplify;

fHstar=((Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd//Simplify;

Randomly sampling parameter space

Fixing the time scale assume b=1

Stability

Solve[(R0[(1/2)]/.{α->A-γ}/.b->1)==1,A]//FullSimplify

If 0<Y<X<1 then α+γ+δ<1

Clear[randParsEpi]
randParsEpi[npts_]:=randParsEpi[npts]=Block[{ct,fHr,fPr,Xr,Yr,δr,αr,γr,br,temp,temp2,out},out={};
For[ct=1,ct≤npts,ct++,
br=1;
Xr=RandomReal[{0,1}];Yr=RandomReal[{0,Xr}];
δr=RandomReal[{0.1,1}];
γr=RandomReal[{0.1,(1/2) (Xr+Yr-(2+Xr+Yr) δr)}];
αr=RandomReal[{0,(1/2) (Xr+Yr-(2+Xr+Yr) δr)-γr}];
temp2={X->Xr,Y->Yr,α->αr,γ->γr,b->br,δ->δr};
temp={fH->((Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd/.temp2,fP->((Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd/.temp2,X->Xr,Y->Yr,α->αr,γ->γr,b->br,δ->δr};
If[(R0[(1/2)]/.temp)>1&&0.1<δr<br,
(*Print[temp];*)
AppendTo[out,temp],ct--
];
];
out
]

γ=0;

Clear[randPars2Epi]
randPars2Epi[npts_]:=randPars2Epi[npts]=Block[{ct,fHr,fPr,Xr,Yr,δr,αr,γr,br,temp2,temp,out},out={};
For[ct=1,ct≤npts,ct++,
br=1;
Xr=RandomReal[{0,1}];Yr=RandomReal[{0,Xr}];
δr=RandomReal[{0.1,1}];
γr=0;
αr=RandomReal[{0,(1/2) (Xr+Yr-(2+Xr+Yr) δr)-γr}];
temp2={X->Xr,Y->Yr,α->αr,γ->γr,b->br,δ->δr};
temp={fH->((Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd/.temp2,fP->((Sum[i[h,p,t],{h,1,2},{p,1,2}])/κ)/.equEnd/.temp2,X->Xr,Y->Yr,α->αr,γ->γr,b->br,δ->δr};
If[(R0[(1/2)]/.temp)>1&&0.1<δr<br,
(*Print[temp];*)
AppendTo[out,temp],ct--
];
];
out
]

Epidemiological Model

κSol[pars_,rep_]:=κ/.Flatten[Solve[(Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.pars)== (10^(2+0.25*(rep-1))),κ],1]
nRep=7;
nDemes=1000;
gmax=250;
(*sampFreq=1;*)

epiFile[pars_]:=Block[{out,rep},
out="Input file for Ecological simulations

number of reps (nRep): "<> ToString[nRep]<>"
number of demes (nDemes): "<> ToString[nDemes]<>"
maximum number of generations (gmax): "<> ToString[gmax]<>"
host birth rate (b): "<> ToString[b/.pars]<>"
matching infection rate (X): "<> ToString[X/.pars]<>"
mismatching infection rate (Y): "<> ToString[Y/.pars]<>"
host natural death rate (delta): "<> ToString[δ/.pars]<>"
death probability (alpha): "<> ToString[α/.pars]<>"
pathogen death rate (gamma): "<>ToString[γ/.pars]<>"
carrying capacity (kappa): ";
For[rep=1,rep≤nRep,rep++,out=out<>ToString[κSol[pars,rep]]<>" "];
out
]

Do[Export[NotebookDirectory[]<>"/sims_11_1/inputEpi_"<>ToString[i]<>".txt",epiFile[randPars2Epi[100][[i]]]],{i,1,100}]

Importing Data

Parameters and Data

(*Random Sampling*)
fileE=NotebookDirectory[]<>"sims_11_1/epi_11_4/runEpi_";
fileE2=NotebookDirectory[]<>"sims_11_1/epi_11_4/";
nRuns=100;

(*γ=0*)
(*fileE=NotebookDirectory[]<>"sims_11_1/epi_11_6/runEpi_";
fileE2=NotebookDirectory[]<>"sims_11_1/epi_11_6/";*)

Clear[pars]
pars[file_,run_]:=pars[file,run]=Import[file<>ToString[run]<>"/out.csv"];

pars[fileE,1]

Clear[nReps,nDemes,gmax]
Clear[X,Y,δ,α,γ,κ]
nReps[file_,run_]:=pars[file,run][[2,2]];
nDemes[file_,run_]:=pars[file,run][[3,2]];
gmax[file_,run_]:=pars[file,run][[4,2]];
b[file_,run_]:=pars[file,run][[5,2]]
X[file_,run_]:=pars[file,run][[6,2]]
Y[file_,run_]:=pars[file,run][[7,2]]
δ[file_,run_]:=pars[file,run][[8,2]]
α[file_,run_]:=pars[file,run][[9,2]]
γ[file_,run_]:=pars[file,run][[10,2]]
κ[file_,run_,rep_]:=pars[file,run][[11,rep+1]]

simPars[file_,run_,rep_]:={b->b[file,run],X->X[file,run],Y->Y[file,run],δ->δ[file,run],α->α[file,run],γ->γ[file,run],gmax->gmax[file,run],κ->κ[file,run,rep]}

Clear[time,state,stateN]
time[file_,run_,rep_]:=time[file,run,rep]=Import[file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/time.csv"][[;;,;;-2]]
state[file_,run_,rep_]:=state[file,run,rep]=Table[Import[file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/state_"<>ToString[v]<>".csv"][[;;,;;-2]],{v,1,6}]//N
stateN[file_,run_,rep_]:=stateN[file,run,rep]=Table[Import[file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/state_"<>ToString[v]<>".csv"][[;;,;;-2]],{v,7,12}]//N

fileExist[file_,run_,rep_,t_]:=FileExistsQ[file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/state_"<>ToString[t]<>".csv"]

Clear[MoranNeut]
MoranNeut[file_,run_,rep_]:=MoranNeut[file,run,rep]=Table[0.5(((1-2/(Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}])))^g),{g,1,gmax[file,run]}]/.equEnd/.simPars[file,run,rep]

Processing

Bootstrap Functions

BootstrapMean[data_,nBoot_]:=BootstrapMean[data,nBoot]=Block[{out,b},out={};
For[b=1,b≤nBoot,b++,
AppendTo[out,Mean[RandomChoice[data,Length[data]]]]
];
out]

Bootstrap[data_,dummy_]:=Bootstrap[data,dummy]=RandomChoice[data,Length[data]]

Equilibrium Population Size

Clear[nHEqu]
nHEqu[file_,run_,rep_]:=nHEqu[file,run,rep]=Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep]

Clear[nPEqu]
nPEqu[file_,run_,rep_]:=nPEqu[file,run,rep]=Sum[i[h,p,t],{p,1,2},{h,1,2}]/.equEnd/.simPars[file,run,rep]

Allele Frequency

Clear[pHMtrx]
pHMtrx[file_,run_,rep_]:=pHMtrx[file,run,rep]=((state[file,run,rep][[1]]+state[file,run,rep][[3]]+state[file,run,rep][[4]])/(Total[state[file,run,rep]]))

Clear[pHNMtrx]
pHNMtrx[file_,run_,rep_]:=pHNMtrx[file,run,rep]=((stateN[file,run,rep][[1]]+stateN[file,run,rep][[3]]+stateN[file,run,rep][[4]])/(Total[stateN[file,run,rep]]))

Heterozygosity

Clear[HetMtrx,HetMtrx2]
HetMtrx[file_,run_,rep_]:=HetMtrx[file,run,rep]=2pHMtrx[file,run,rep](1-pHMtrx[file,run,rep])//Chop
HetMtrx2[file_,run_,rep_]:=HetMtrx2[file,run,rep]=Table[HetMtrx[file,run,rep][[;;,g]]/.Indeterminate->0,{g,1,250}]//Transpose

Clear[HetNMtrx,HetNMtrx2]
HetNMtrx[file_,run_,rep_]:=HetNMtrx[file,run,rep]=2pHNMtrx[file,run,rep](1-pHNMtrx[file,run,rep])//Chop
HetNMtrx2[file_,run_,rep_]:=HetNMtrx2[file,run,rep]=Table[HetNMtrx[file,run,rep][[;;,g]]/.Indeterminate->MoranNeut[file,run,rep][[g]],{g,1,250}]//Transpose

Heterozygosity at generation τ

Clear[HetCompMtrx2]
HetCompMtrx2[file_,rep_,τ_]:=HetCompMtrx2[file,rep,τ]=Import[fileE2<>"/HetMtrx/HetMtrx_"<>ToString[rep]<>"_"<>ToString[τ]<>".csv"]

Harmonic Mean Population Size

Clear[NeMtrx]
NeMtrx[rep_]:=NeMtrx[rep]=Import[fileE2<>"/NeMtrx/NeMtrx_"<>ToString[rep]<>".csv"]

Clear[NeMtrx2]
NeMtrx2[rep_,g_]:=NeMtrx2[rep,g]=Import[fileE2<>"/NeMtrx/NeMtrx_"<>ToString[rep]<>"_"<>ToString[g]<>".csv"]

Processing- Run only once

Preprocessing- RUN ONLY ONCE!

Clear[HetCompMtrx]
HetCompMtrx[file_,rep_,τ_]:=HetCompMtrx[file,rep,τ]=Table[{Mean[HetNMtrx2[file,run,rep][[;;,τ]]],StandardDeviation[BootstrapMean[HetNMtrx2[file,run,rep][[;;,τ]],1000]],Mean[HetMtrx2[file,run,rep][[;;,τ]]],StandardDeviation[BootstrapMean[HetMtrx2[file,run,rep][[;;,τ]],1000]]},{run,1,100}]

(*Run This Only Once*)
(*CreateDirectory[fileE2<>"/HetMtrx"];*)
Table[Export[fileE2<>"/HetMtrx/HetMtrx_"<>ToString[rep]<>"_"<>ToString[τ]<>".csv",HetCompMtrx[fileE,rep,τ]],{rep,3,7},{τ,10,250,10}];

Harmonic Mean Population Size

state[fileE,1,3]//Dimensions

Clear[HarmonicMeanNe]
HarmonicMeanNe[file_,run_,rep_,g_]:=HarmonicMeanNe[file,run,rep]=Block[{temp},temp=Total[state[file,run,rep][[;;,;;,;;g]]];
Map[If[Min[#]==0,0,(g/(Total[((#)^(-1))]))]&,temp]]

NeExp[gen_]:=Do[Export[fileE2<>"NeMtrx/NeMtrx_"<>ToString[rep]<>"_"<>ToString[gen]<>".csv",Table[HarmonicMeanNe[fileE,run,rep,gen],{run,1,100}]];Print[rep],{rep,1,7}];

(*Run This Only Once*)
(*CreateDirectory[fileE2<>"NeMtrx"];*)
NeExp[10]

Plotting

Dynamics of a single rep

dynPlot[file_,run_,rep_]:=Show[
(*Neutral Simulations*)
ListLinePlot[{Mean[HetNMtrx[file,run,rep]],Mean[HetNMtrx[file,run,rep]]+2((StandardDeviation[HetNMtrx[file,run,rep]])/(Sqrt[nDemes[file,run]])),Mean[HetNMtrx[file,run,rep]]-2((StandardDeviation[HetNMtrx[file,run,rep]])/(Sqrt[nDemes[file,run]]))},PlotStyle->{Black,Gray,Gray},Filling->2->{3},FillingStyle->Directive[Gray,Opacity[0.2]]],
(*Coevolutionary Simulations*)
(*ListLinePlot[Table[Mean[Bootstrap[Htable[file,run,rep],bs]],{bs,0,100}],PlotStyle->Directive[Pink,Opacity[0.1]]],*)ListLinePlot[{Mean[HetMtrx[file,run,rep]],Mean[HetMtrx[file,run,rep]]+2((StandardDeviation[HetMtrx[file,run,rep]])/(Sqrt[nDemes[file,run]])),Mean[HetMtrx[file,run,rep]]-2((StandardDeviation[HetMtrx[file,run,rep]])/(Sqrt[nDemes[file,run]]))},PlotStyle->{Red,Pink,Pink},Filling->2->{3},FillingStyle->Directive[Pink,Opacity[0.2]]],
(*Moran Neutral Expectation*)
ListLinePlot[MoranNeut[file,run,rep],PlotStyle->Purple](*,PlotRange->{{0,5},{0.48,0.5}}*)]

dynPlot[fileE,8,4]

Mean Heterozygosity

Mean Smoothing

Project data onto the X=Y line.

Clear[dataNew]
dataNew[data_]:=dataNew[data]=Transpose[{Cos[Abs[ArcTan[((data[[;;,2]])/(data[[;;,1]]))]-(π/4)]]*Sqrt[Total[Transpose[(data^2)]]],Sin[ArcTan[((data[[;;,2]])/(data[[;;,1]]))]-(π/4)]*Sqrt[Total[Transpose[(data^2)]]]}]

A gaussian weighting of the data  near the point xNew along the X=Y line

Clear[weights]
weights[xNew_,αW_,data_]:=weights[xNew,αW,data]=Exp[-((((dataNew[data][[;;,1]]-xNew))^2)/(2 ((αW)^2)))]//Chop;

Find the weighted mean and standard deviation and standard error.

Clear[meanNew]
meanNew[xNew_,αW_,data_]:=meanNew[xNew,αW,data]=Total[weights[xNew,αW,data]dataNew[data][[;;,2]]]/Total[weights[xNew,αW,data]]

SDNew[xNew_,αW_,data_]:=Sqrt[Total[weights[xNew,αW,data]*(((dataNew[data][[;;,2]]-meanNew[xNew,αW,data]))^2)]/Total[weights[xNew,αW,data]]]

SENew[xNew_,αW_,data_]:=((SDNew[xNew,αW,data])/(Sqrt[Total[weights[xNew,αW,data]]]))

The rolling average fit along the X+Y line (rotated plot)

Clear[RollFit]
RollFit[αW_,npts_,data_]:=RollFit[αW,npts,data]=Table[{x,meanNew[x,αW,data]},{x,0.001,Sqrt[2*(0.5^2)]-0.001,((Sqrt[2*(0.5^2)])/npts)}];

Projecting back into the original coordinates.

RollFitOld[αW_,npts_,data_]:=Transpose[{Cos[ArcTan[((RollFit[αW,npts,data][[;;,2]])/(RollFit[αW,npts,data][[;;,1]]))]+(π/4)]Sqrt[((RollFit[αW,npts,data][[;;,2]])^2)+((RollFit[αW,npts,data][[;;,1]])^2)],Sin[ArcTan[((RollFit[αW,npts,data][[;;,2]])/(RollFit[αW,npts,data][[;;,1]]))]+(π/4)]Sqrt[((RollFit[αW,npts,data][[;;,2]])^2)+((RollFit[αW,npts,data][[;;,1]])^2)]}]

WF Diploid Over-dominance as a comparison

transition probability of i copies to j copies of the allele

w[i_,n_,s_]:=((i/n))(1-s)+(((n-i))/n)-((((n-i))/n)(1-s)+(i/n))//Simplify

npSol=Flatten[NDSolve[{D[p[t],t]==(w[p[t]n,n,s])p[t](1-p[t])/.{n->100,s->0.1},p[0]==0.1},p[t],{t,0,100}]]

p[i_,j_,n_,s_]:=If[i==0,If[j==0,1,0],If[i==n,If[j==n,1,0],Binomial[ n,j]((((((i/n) (1+w[i,n,s]))/((i/n) (1+w[i,n,s])+((n-i)/n) (1+w[n-i,n,s])))))^j)(((((((n-i)/n) (1+w[n-i,n,s]))/((i/n) (1+w[i,n,s])+((n-i)/n) (1+w[n-i,n,s])))))^(n-j))]]

p0Vec[n_]:=Table[If[i==(n/2),1,0],{i,0,n}]

Clear[tMtrx]
tMtrx[n_,s_]:=tMtrx[n,s]=Table[p[i,j,n,s],{j,0,n},{i,0,n}]

Clear[pOut]
pOut[n_,s_,t_]:=pOut[n,s,t]=MatrixPower[tMtrx[n,s],t].p0Vec[n]

hOver[n_,s_,t_]:=Table[2(i/n)(((n-i)/n)),{i,0,n}].pOut[n,s,t]

hNeut[n_,t_]:=(1/2) Exp[((-t)/n)]

temp=Table[{hNeut[Round[(10^npower),2],250],hOver[Round[(10^npower),2],0.01,250]},{npower,2,2.75,0.25}]

HetOver=Join[{{0,0}},temp,{{0.5,0.5}}]//N;

Show[ListLinePlot[HetOver,PlotStyle->Black],ListPlot[Map[{#}&,HetOver],PlotStyle->Join[{Black},Table[ ColorData["CMYKColors"][((rep-1)/6)],{rep,1,4}],{Black}]]]

Moran Diploid Over-dominance as a comparison

transition probability of i copies to j copies (out of n total copies) of the allele

w[i_,n_,s_]:=((i/n))(1-s)+(((n-i))/n)-((((n-i))/n)(1-s)+(i/n))//Simplify

wA[i_,n_,s_]:=((((i/n))(1-s)+(((n-i))/n))/(((((i/n)))^2)(1-s)+2(i/n)(((n-i))/n)+((((((n-i))/n)))^2)(1-s)))

wa[i_,n_,s_]:=(((((n-i))/n)(1-s)+(i/n))/(((((i/n)))^2)(1-s)+2(i/n)(((n-i))/n)+((((((n-i))/n)))^2)(1-s)))

npSol=Flatten[NDSolve[{D[p[t],t]==(w[p[t]n,n,s])p[t](1-p[t])/.{n->100,s->0.1},p[0]==0.1},p[t],{t,0,100}]]

npSol=Flatten[NDSolve[{D[p[t],t]==(wA[p[t]n,n,s]-1)p[t]/.{n->100,s->0.1},p[0]==0.1},p[t],{t,0,100}]]

p[i_,j_,n_,s_]:=If[i==j,1-wA[i,n,s](i/n)(((n-i))/n)-wa[i,n,s](((n-i))/n)(i/n),If[j==i+1,wA[i,n,s](i/n)(((n-i))/n),If[j==i-1,wa[i,n,s](((n-i))/n)(i/n),0]]]

p0Vec[n_]:=Table[If[i==(n/2),1,0],{i,0,n}]

Clear[tMtrx]
tMtrx[n_,s_]:=tMtrx[n,s]=Table[p[i,j,n,s],{j,0,n},{i,0,n}]

Clear[pOut]
pOut[n_,s_,t_]:=pOut[n,s,t]=MatrixPower[tMtrx[n,s],t n].p0Vec[n]

hOver[n_,s_,t_]:=Table[2(i/n)(((n-i)/n)),{i,0,n}].pOut[n,s,t]

hNeut[n_,t_]:=(1/2) (((1-(2/n)))^t)//N

Clear[HetOver]
HetOver[s_,τ_]:=HetOver[s,τ]=Block[{temp},
temp=Table[{hNeut[Round[(10^npower),2],250],hOver[Round[(10^npower),2],s,τ]},{npower,2,3.25,0.25}];
Join[{{0,0}},temp,{{0.5,0.5}}]//N]

ListLinePlot[{HetOver[0.01,250],HetOver[0.05,250],HetOver[0.1,250]},PlotStyle->{Gray,Darker[Gray],Black}]

Variation within versus across replicates

?ErrorListPlot

Needs["ErrorBarPlots`"]

VarVisual[rep_]:=Block[{temp},temp=Sort[Map[{#[[3]],1.96#[[4]]}&,HetCompMtrx2[fileE,rep,250]],#1[[1]]<#2[[1]]&];
ErrorListPlot[temp,PlotRange->{Min[Map[#[[1]]-#[[2]]&,temp]]0.95,Max[Map[#[[1]]+#[[2]]&,temp]]*1.05},Frame->True,FrameTicksStyle->Directive[Black,12],ImageSize->Medium,PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]]

Show[Table[VarVisual[rep],{rep,1,7}],PlotRange->All]
Export[NotebookDirectory[]<>"Plots/Epi_Variation.pdf",%,ImageResolution->300];

VarWB=Table[{Round[Mean[HetCompMtrx2[fileE,rep,250][[;;,4]]],0.0001],Round[StandardDeviation[HetCompMtrx2[fileE,rep,250][[;;,3]]],0.001],Round[StandardDeviation[HetCompMtrx2[fileE,rep,250][[;;,3]]]/Mean[HetCompMtrx2[fileV,rep,250][[;;,4]]],0.01]},{rep,1,7}]

Clear[temp]

ErrorListPlot2[file_,rep_,τ_]:=Block[{temp},

temp=Sort[HetCompMtrx2[file,rep,τ],#1[[3]]<#2[[3]]&];
Show[ListLinePlot[Table[{{r,temp[[r,3]]-2temp[[r,4]]},{r,temp[[r,3]]+2temp[[r,4]]}},{r,1,100}],PlotStyle->Directive[Lighter[ColorData["CMYKColors"][((rep-1)/6)]]],PlotRange->{0,0.5}],
ListPlot[Table[{r,temp[[r,3]]},{r,1,100}],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)],PlotMarkers->{"●",4}]
]
]

Show[Table[ErrorListPlot2[fileE,rep,250],{rep,1,7}],PlotRange->All,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,16],ImageSize->Medium]
Export[NotebookDirectory[]<>"Plots/Epi_Variation.pdf",%,ImageResolution->300];

Export[fileE2<>"VarWB.csv",VarWB]

Comparison to the simulated neutral expectation

meanPlot[file_,τ_]:=Block[{sigData,meanData},
meanData=Flatten[Table[HetCompMtrx2[file,rep,τ][[;;,{1,3}]],{rep,1,7}],1];
sigData[rep_]:=sigData[rep]=Select[HetCompMtrx2[file,rep,τ],(#[[3]]-2 #[[4]])>(#[[1]]+2#[[2]])||(#[[3]]+2 #[[4]])<(#[[1]]-2#[[2]])&];
Show[
Table[
ListPlot[{sigData[rep][[;;,{1,3}]],Complement[HetCompMtrx2[file,rep,τ],sigData[rep]][[;;,{1,3}]]},PlotStyle-> ColorData["CMYKColors"][((rep-1)/6)],PlotMarkers->{●,○}],{rep,1,7}],
ListLinePlot[Table[RollFitOld[0.1,40,Bootstrap[meanData,x]],{x,1,100}],PlotStyle->Directive[Lighter[Blue],Opacity[0.1]]],
ListLinePlot[RollFitOld[0.1,40,meanData],PlotStyle->Blue],Plot[x,{x,0,0.5},PlotStyle->{Black,Thick,Dashing[0.02]}],
(*Compariosn to Overdominance*)
ListLinePlot[{HetOver[0.01,250],HetOver[0.05,250],HetOver[0.1,250]},PlotStyle->{Gray,Darker[Gray],Black}],
(*Plot Options*)
PlotRange->{{0,0.5},{0,0.5}},AspectRatio->1,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,16],Axes->False]]

meanPlot[fileE,250]
Export[NotebookDirectory[]<>"Plots/Epi_MeanHet_Rand2.pdf",%,ImageResolution->300];

meanData[file_,τ_]:=Flatten[Table[HetCompMtrx2[file,rep,τ][[;;,{1,3}]],{rep,1,7}],1];

meanExport[file_,τ_]:=Transpose[Flatten[Join[{Transpose[RollFitOld[0.1,50,meanData[file,τ]]]},Table[Transpose[RollFitOld[0.1,50,Bootstrap[meanData[file,τ],x]]],{x,1,100}]],1]]

meanExport[fileE,250]

Export[fileE2<>"MeanHet.csv",meanExport[fileE,250]];

Explaining Variation

Perturbations

γδPlot[file_,τ_]:=Block[{data,plot,min,max,lm},
data[rep_]:=Transpose[{Table[(δ[file,run]Sum[s[h,t],{h,1,2}]-γ[file,run]Sum[i[h,p,t],{h,1,2},{p,1,2}])/.equEnd/.simPars[file,run,rep],{run,1,nRuns}],HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min=Min[data[1][[;;,1]]];max=Max[data[1][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y][x];
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervalus"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min,max},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
Show[Table[plot[rep],{rep,1,7}],PlotRange->All,AxesOrigin->{0,0},Frame->True,FrameLabel->{"γ-δ","ΔH"},ImageSize->Large]
]

γδPlot[fileE,100]

λMax

λMaxPlot[file_,τ_]:=Block[{data,plot,min,max,lm},
data[rep_]:=Transpose[{Table[-λMaxRe[simPars[file,run,rep]],{run,1,nRuns}],HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min=Min[data[1][[;;,1]]];max=Max[data[1][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y][x];
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min,max},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
Table[R2[rep]=lm[data[rep],"RSquared"],{rep,1,7}];
Table[Slope[rep]={lm[data[rep],"BestFitParameters"][[2]],lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]-lm[data[rep],"BestFitParameters"][[2]]},{rep,1,7}];
Show[Table[plot[rep],{rep,1,7}],PlotRange->All,AxesOrigin->{0,0},Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,16],ImageSize->Medium]
]

λMaxPlot[fileE,250]
Export[NotebookDirectory[]<>"Plots/Epi_lambda.pdf",%,ImageResolution->300];

Table[Join[Round[{R2[rep]},0.01],Round[Slope[rep],0.001]],{rep,1,7}]
Export[NotebookDirectory[]<>"RSquaredTemp.csv",%];

Harmonic Mean Population Size

Does harmonic mean population size at time t explain ΔH at time t

NePlot[file_,τ_]:=Block[{eqTab,data,plot,min,max,lm},
eqTab[rep_]:=Table[Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
data[rep_]:=Transpose[{((Mean[Transpose[NeMtrx[rep]]]-eqTab[rep])/(eqTab[rep])),HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min[rep_]:=Min[data[rep][[;;,1]]];max[rep_]:=Max[data[rep][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y][x];
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],(*Print[lm[data[rep],"RSquared"]];*)Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
(*Storing Fit Info*)
FitTab=Table[Flatten[{StandardDeviation[data[rep][[;;,1]]],lm[data[rep],"BestFitParameters"],lm[data[rep],"ParameterConfidenceIntervals"]}],{rep,1,7}];
(*Making Plot*)
Show[Table[plot[rep],{rep,1,7}],Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[lm[data[rep],"RSquared"],0.001]],12,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.2,1.2-0.06 rep}]],{rep,5,7}],PlotRange->All,AxesOrigin->{0,0},Frame->True,(*FrameLabel->{"Ne-(OverHat[N])","ΔH"},*)FrameTicksStyle->Directive[Black,12], ImageSize->Medium]
]

NePlot[fileE,250]
Export[NotebookDirectory[]<>"Plots/Epi_Ne.pdf",%,ImageResolution->300];

Does the harmonic mean at time g explain ΔH at time t

NePlot2[file_,τ_,g_]:=Block[{eqTab,data,plot,min,max,lm},
eqTab[rep_]:=Table[Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
data[rep_]:=Transpose[{((Mean[Transpose[NeMtrx2[rep,g]]]-eqTab[rep])/(eqTab[rep])),HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min[rep_]:=Min[data[rep][[;;,1]]];max[rep_]:=Max[data[rep][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y,ConfidenceLevel->.99][x];
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],(*Print[lm[data[rep],"RSquared"]];*)Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
(*Storing Fit Info*)
FitTab2[g]=Table[Flatten[{StandardDeviation[data[rep][[;;,1]]],lm[data[rep],"BestFitParameters"],lm[data[rep],"ParameterConfidenceIntervals"],lm[data[rep],"RSquared"]}],{rep,1,7}];
(*Plotting*)
Show[Table[plot[rep],{rep,1,7}],(*Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[lm[data[rep],"RSquared"],0.001]],12,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.2,1.2-0.06 rep}]],{rep,5,7}],*)PlotRange->All,AxesOrigin->{0,0},Frame->True,(*FrameLabel->{"Ne-(OverHat[N])","ΔH"},*)FrameTicksStyle->Directive[Black,12], ImageSize->Medium]
]

NePlot2[fileE,250,10]
(*Export[NotebookDirectory[]<>"Plots/Epi_Ne_50.pdf",%,ImageResolution->300];*)

Round[FitTab2[10],0.001]//MatrixForm

NePlot2[fileE,250,50]

Round[FitTab2[50],0.001]//MatrixForm

NePlot2[fileE,250,250]
(*Export[NotebookDirectory[]<>"Plots/Epi_Ne_50.pdf",%,ImageResolution->300];*)

Round[FitTab2[250],0.001]//MatrixForm

Does the effective population size at generation g explain ΔH at generation g?

NePlot3[file_,g_]:=Block[{eqTab,data,plot,min,max,lm},
eqTab[rep_]:=Table[Sum[s[h,t],{h,1,2}]+Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
data[rep_]:=Transpose[{((Mean[Transpose[NeMtrx2[rep,g]]]-eqTab[rep])/(eqTab[rep])),HetCompMtrx2[file,rep,g][[;;,3]]-HetCompMtrx2[file,rep,g][[;;,1]]}];
min[rep_]:=Min[data[rep][[;;,1]]];max[rep_]:=Max[data[rep][[;;,1]]];
lm[data_,x_]:=lm[data,x]=LinearModelFit[data,y,y,ConfidenceLevel->.99][x];
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],(*Print[lm[data[rep],"RSquared"]];*)Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
(*Storing Fit Info*)
FitTab3[g]=Table[Flatten[{StandardDeviation[data[rep][[;;,1]]],lm[data[rep],"BestFitParameters"],lm[data[rep],"ParameterConfidenceIntervals"],lm[data[rep],"RSquared"]}],{rep,1,7}];
(*Plotting*)
Show[Table[plot[rep],{rep,1,7}],(*Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[lm[data[rep],"RSquared"],0.001]],12,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.2,1.2-0.06 rep}]],{rep,5,7}],*)PlotRange->All,AxesOrigin->{0,0},Frame->True,(*FrameLabel->{"Ne-(OverHat[N])","ΔH"},*)FrameTicksStyle->Directive[Black,12], ImageSize->Medium]
]

NePlot3[fileE,10]

FitTab3[10]

Relative Turnover rate

turnoverPlot[file_,τ_]:=Block[{δTab,γTab,βTab,data,plot,min,max,lm},
δTab[rep_]:=Table[δ Sum[s[h,t],{h,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
γTab[rep_]:=Table[γ Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
βTab[rep_]:=Table[Sum[β[h,p]s[h,t]i[k,p,t],{h,1,2},{p,1,2},{k,1,2}]/.equEnd/.βSub/.simPars[file,run,rep],{run,1,100}];
data[rep_]:=data[rep]=Transpose[{δTab[rep]-γTab[rep]-βTab[rep],HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min[rep_]:=Min[data[rep][[;;,1]]];max[rep_]:=Max[data[rep][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y][x];
(*lm[data_,x_]:=lm[data,x]=NonlinearModelFit[data,{y0+((-L)/(1+Exp[k(y-x0)])),0<L,0<k},{y0,L,k,x0},y][x];*)
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],(*Print[lm[data[rep],"RSquared"]];*)Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
(*plot[rep_]:=Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]];*)
Show[Table[plot[rep],{rep,1,7}],(*Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[lm[data[rep],"RSquared"],0.001]],12,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.2,1.2-0.06 rep}]],{rep,5,7}],*)PlotRange->All,AxesOrigin->{0,0},Frame->True,(*FrameLabel->{"Ne-(OverHat[N])","ΔH"},*)FrameTicksStyle->Directive[Black,12], ImageSize->Medium]
]

turnoverPlot[fileE,250]

turnoverPlot2[file_,τ_]:=Block[{δTab,γTab,βTab,data,plot,min,max,lm},
δTab[rep_]:=Table[δ Sum[s[h,t],{h,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
γTab[rep_]:=Table[γ Sum[i[h,p,t],{h,1,2},{p,1,2}]/.equEnd/.simPars[file,run,rep],{run,1,100}];
βTab[rep_]:=Table[Sum[β[h,p]s[h,t]i[k,p,t],{h,1,2},{p,1,2},{k,1,2}]/.equEnd/.βSub/.simPars[file,run,rep],{run,1,100}];
data[rep_]:=data[rep]=Transpose[{δTab[rep]-γTab[rep],HetCompMtrx2[file,rep,τ][[;;,3]]-HetCompMtrx2[file,rep,τ][[;;,1]]}];
min[rep_]:=Min[data[rep][[;;,1]]];max[rep_]:=Max[data[rep][[;;,1]]];
lm[data_,x_]:=LinearModelFit[data,y,y][x];
(*lm[data_,x_]:=lm[data,x]=NonlinearModelFit[data,y0+((-L)/(1+Exp[k(y-x0)])),{y0,0<L,0<k,x0},y][x];*)
plot[rep_]:=If[Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,1]]]==Sign[lm[data[rep],"ParameterConfidenceIntervals"][[2,2]]],(*Print[lm[data[rep],"RSquared"]];*)Show[ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],Plot[lm[data[rep],ξ],{ξ,min[rep],max[rep]},PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]],ListPlot[data[rep],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]]
];
Show[Table[plot[rep],{rep,1,7}],(*Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[lm[data[rep],"RSquared"],0.001]],12,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.2,1.2-0.06 rep}]],{rep,5,7}],*)PlotRange->All,AxesOrigin->{0,0},Frame->True,(*FrameLabel->{"Ne-(OverHat[N])","ΔH"},*)FrameTicksStyle->Directive[Black,12], ImageSize->Medium]
]

turnoverPlot2[fileE,250]

Allele Frequency Dynamics

Clear[smallest,largest]

(*The really small element of rep=3*)
smallest[rep_]:=Ordering[HetCompMtrx2[fileE,rep,250][[;;,3]],1]

largest[4]

(*The really high element of rep=3*)
largest[rep_]:=Ordering[HetCompMtrx2[fileE,rep,250][[;;,3]],-1]

AFPlot[file_,run_,rep_,nd_]:=Block[{h1,h2,p1,p2,coev,dirSel,fixed,extinct,all},
{h1,h2,p1,p2}={SparseArray[Total[state[file,run,rep][[{1,3,4},1;;nd,;;]]]],SparseArray[Total[state[file,run,rep][[{2,5,6},1;;nd,;;]]]],SparseArray[Total[state[file,run,rep][[{3,5},1;;nd,;;]]]],SparseArray[Total[state[file,run,rep][[{4,6},1;;nd,;;]]]]};
all=Flatten[Table[{deme,g},{deme,1,nd},{g,1,250}],1];
coev=Intersection[h1["NonzeroPositions"],h2["NonzeroPositions"],p1["NonzeroPositions"],p2["NonzeroPositions"]];
fixed=Union[Complement[all,h1["NonzeroPositions"]],Complement[all,h2["NonzeroPositions"]]];
extinct=Complement[Intersection[Complement[all,p1["NonzeroPositions"]],Complement[all,p2["NonzeroPositions"]]],fixed];
dirSel=Complement[Union[Complement[p1["NonzeroPositions"],p2["NonzeroPositions"]],Complement[p2["NonzeroPositions"],p1["NonzeroPositions"]]],Union[fixed,extinct]];
coev=Table[Select[coev,#[[1]]==d&],{d,1,nd}][[;;,;;,2]];dirSel=Table[Select[dirSel,#[[1]]==d&],{d,1,nd}][[;;,;;,2]];fixed=Table[Select[fixed,#[[1]]==d&],{d,1,nd}][[;;,;;,2]];
extinct=Table[Select[extinct,#[[1]]==d&],{d,1,nd}][[;;,;;,2]];

Show[ListLinePlot[Table[pHMtrx[file,run,rep][[d,coev[[d]]]],{d,1,nd}],PlotStyle->Directive[Gray,Opacity[0.2]]],ListLinePlot[Table[If[Length[dirSel[[d]]]>0,Transpose[{dirSel[[d]],pHMtrx[file,run,rep][[d,dirSel[[d]]]]}],Missing[]],{d,1,nd}],PlotStyle->Directive[Red,Opacity[0.3]]],ListLinePlot[Table[If[Length[extinct[[d]]]>0,Transpose[{extinct[[d]],pHMtrx[file,run,rep][[d,extinct[[d]]]]}],Missing[]],{d,1,nd}],PlotStyle->Directive[Green,Opacity[0.3]]],ListLinePlot[Table[If[Length[fixed[[d]]]>0,Transpose[{fixed[[d]],pHMtrx[file,run,rep][[d,fixed[[d]]]]}],Missing[]],{d,1,nd}],PlotStyle->Directive[Blue,Opacity[0.2]]],Epilog->Text[Style["(Subscript[λ, Max])="<>ToString[Round[λMaxRe[simPars[file,run,rep]],0.001]],14,Black],Scaled[{0.8,0.8}]],PlotRange->{0,1},Axes->False,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,12],ImageSize->Medium]
]

Clear[h1,h2,p1,p2]

smallest[4]

AFPlot[fileE,1,4,100]
Export[NotebookDirectory[]<>"Plots/Epi_AF_1.pdf",%,ImageResolution->300];

largest[4]

AFPlot[fileE,72,4,100]
Export[NotebookDirectory[]<>"Plots/Epi_AF_2.pdf",%,ImageResolution->300];

Fixation Order

Clear[hpState]
hpState[file_,run_,rep_]:=hpState[file,run,rep]={Total[state[file,run,rep][[{1,3,4}]]],Total[state[file,run,rep][[{2,5,6}]]],Total[state[file,run,rep][[{3,5}]]],Total[state[file,run,rep][[{4,6}]]]}

hostFixTime[file_,run_,rep_]:=Table[FirstPosition[hpState[file,run,rep][[{1,2},d,;;]],0.],{d,1,1000}]

pathFixTime[file_,run_,rep_]:=Table[FirstPosition[hpState[file,run,rep][[{3,4},d,;;]],0.],{d,1,1000}]

pathExtinctTime[file_,run_,rep_]:=Table[FirstPosition[Total[hpState[file,run,rep][[{3,4},d,;;]]],0.],{d,1,1000}]

Clear[fixOrder]
fixOrder[file_,rep_]:=fixOrder[file,rep]=Block[{temp1,temp2,temp3,temp4,temp5,run,out},out={};
For[run=1,run≤100,run++,
(*All Data Points*)
temp1=Transpose[{hostFixTime[file,run,rep],pathFixTime[file,run,rep],pathExtinctTime[file,run,rep]}];
(*Host Fixed*)
temp2=DeleteCases[temp1,_?(#[[1]]==Missing["NotFound"]&)];
(*Host Fixed and Pathogen Fixed*)
temp3=DeleteCases[temp2,_?(#[[2]]==Missing["NotFound"]&)];
(*Pathogen Fixed (Pathogen may have gone extict) then Host fixed*)
temp4=Select[temp3,#[[1,2]]>#[[2,2]]&];
(*Pathogen Fixed, Pathogen Extinct, then Host fixed*)
temp5=Select[temp4,#[[1,2]]>#[[3,1]]&];
AppendTo[out,{Length[temp1]-Length[temp2](*Host Polymorphic*),Length[temp2]-Length[temp3](*Host Fixed, Pathogen Polymorphic*),Length[temp3]-Length[temp4](*Host Fixed First*),Length[temp4]-Length[temp5](*Pathogen Fixed then Host Fixed*),Length[temp5](*Pathogen Fixed then Pathogen Extinct, then Host Fixed*)}];
];
(out/nRuns)
]

fixBarChart[file_,rep_,τ_]:=BarChart[fixOrder[file,rep][[Ordering[HetCompMtrx2[file,rep,τ][[;;,3]]]]],ChartLayout->"Stacked",ChartStyle->{Gray,Blue,Blue,Red,Green},Frame->True,FrameTicks->{{True,False},{True,False}},ImageSize->Medium,FrameTicksStyle->Directive[Black,14]]

Gray= Host and pathogen polymorphic at generation 250
Blue= Host fixed before pathogen
Red= Pathogen  fixes then host fixes
Green= Pathogen fixes, pathogen goes extinct, then host fixes
Runs (bars) are ordered in terms of increasing heterozygosity

(*GraphicsRow[{fixBarChart[fileE,2,250],fixBarChart[fileE,4,250]}]*)
fixBarChart[fileE,2,250]
Export[NotebookDirectory[]<>"Plots/Epi_DirSel_Rand.pdf",%,ImageResolution-> 300]

Model Comparison

constModel=Import["/Users/ailenemacpherson/Desktop/EcoEpiDrift/sims_11_1/const_11_20_2/MeanHet.csv"];
ecoModel=Import["/Users/ailenemacpherson/Desktop/EcoEpiDrift/sims_11_1/eco_11_6/MeanHet.csv"];
epiModel=Import["/Users/ailenemacpherson/Desktop/EcoEpiDrift/sims_11_1/epi_11_6/MeanHet.csv"];

constModel2[i_]:=constModel[[;;,{2i-1,2i}]]
ecoModel2[i_]:=ecoModel[[;;,{2i-1,2i}]]
epiModel2[i_]:=epiModel[[;;,{2i-1,2i}]]

Show[ListLinePlot[Table[constModel2[i],{i,2,100}],PlotStyle->Directive[Pink,Opacity[0.1]]],
ListLinePlot[constModel2[1],PlotStyle->Red],ListLinePlot[Table[ecoModel2[i],{i,2,100}],PlotStyle->Directive[Green,Opacity[0.1]]],
ListLinePlot[ecoModel2[1],PlotStyle->Darker[Green]],ListLinePlot[Table[epiModel2[i],{i,2,100}],PlotStyle->Directive[Lighter[Blue],Opacity[0.1]]],
ListLinePlot[epiModel2[1],PlotStyle->Blue],
Plot[x,{x,0,0.5},PlotStyle->{Black,Dashed}],AspectRatio->1,Frame->True,FrameTicks->{{True,False},{True,False}},FrameTicksStyle->Directive[Black,14],Axes->False]
Export[NotebookDirectory[]<>"Plots/ModelComparison.pdf",%,ImageResolution->300];