The evolution of condition-dependent self-fertilisation
with a fixed inbreeding depression component

Thomas Lesaffre
John Pannell
Charles Mullon

Last updated May 15th, 2026.

Definitions and initialisation

Zm={αAAm,αAam,αaam}; (* Mutant haplotype. *)
Z={αAA,αAa,αaa}; (* Resident haplotype. *)

(*** Condition depending on genotype at condition locus. **)

wAA=1;
wAa=1-s h;
waa=1-s;

wbar=wAA qAA+wAa qAa+waa qaa; (* Mean condition. *)

(*** Resident pollen types ***)

(** Haplotype frequencies before mutation... **)

PAMbef=((wAA qAA+(wAa/2)qAa)/wbar);
PaMbef=((waa qaa+(wAa/2)qAa)/wbar);

(** ... after mutation. **)

PAM = PAMbef(1-μA)+μa PaMbef;
PaM = PaMbef(1-μa)+μA PAMbef;

(*** Resident outcross ovule types ***)

FAMbef=qAA F0 wAA  (1- αAA) + (qAa/2)F0 wAa (1-αAa);
FaMbef=qaa F0 waa  (1- αaa) + (qAa/2)F0 wAa (1- αAa);

FAM = FAMbef(1-μA) + μa FaMbef;
FaM = FaMbef(1-μa)+ μA FAMbef;

(*** Mutant contributions ***)

w1i =F0 w (2α (1-δ)gAM gAm + (1-α)gAm PAM)+FAM gAm (w/wbar);
w2i = F0 w(2α (1-δ)gAM gam+ (1-α)gam PAM) + FAM gam(w/wbar); 
w3i =F0 w(2α(1-δ) gAm gaM + (1-α) gAm PaM) + FaM  gAm(w/wbar);
w4i = F0 w(2α(1-δ) gaM gam+ (1-α)gam PaM)+FaM gam (w/wbar);
w5i =F0 w α(1-δ) gAm gAm;
w6i =F0 w α(1-δ) gAm gam;
w7i=F0 w α(1-δ) gam gAm;
w8i =F0 w α(1-δ) gam gam;

Ww = {w1i,w2i,w3i,w4i,w5i,w6i,w7i,w8i};

(*** Class 1: AM/Am focal mutant. **)

(* Haplotype frequencies in gametes *)

gAM1 = ((1-μA)/2);
gaM1 = ((μA)/2);
gAm1 = ((1-μA)/2);
gam1 = ((μA)/2);

(* Selfing rate and condition *)

α1 = ((αAA+αAAm)/2);
w1 = wAA;

rep1={w->w1, α-> α1, gAM->gAM1, gaM->gaM1, gAm->gAm1, gam->gam1}; (* Replacement rule *)

Ww1 = Ww/. rep1;(* Mutant offspring produced. *)

(*** Class 2: AM/am focal mutant. ***)

gAM2 = (((1-r)(1-μA)+ r μa)/2);
gaM2 = (((1-r)μA+r(1-μa))/2);
gAm2 = ((r(1-μA)+(1-r)μa)/2);
gam2 = ((r μA+(1-r)(1-μa))/2);

α2 = ((αAa+αAam)/2);
w2 = wAa;

rep2={w->w2, α-> α2, gAM->gAM2, gaM->gaM2, gAm->gAm2, gam->gam2};

Ww2 = Ww/. rep2;(* Mutant offspring produced. *)

(*** Class 3: Am/aM focal mutant. ***)

gAM3 = ((r(1-μA)+(1-r)μa)/2);
gaM3 = ((r μA+(1-r)(1-μa))/2);
gAm3 = (((1-r)(1-μA)+r μa)/2);
gam3 = (((1-r)μA+r(1-μa))/2);

α3 = ((αAa+αAam)/2);
w3 = wAa;

rep3={w->w3, α-> α3, gAM->gAM3, gaM->gaM3, gAm->gAm3, gam->gam3};

Ww3 = Ww/. rep3;(* Mutant offspring produced. *)

(*** Class 4: aM/am focal mutant. ***)

gAM4 = ((μa)/2);
gaM4 = ((1-μa)/2);
gAm4 = ((μa)/2);
gam4 = ((1-μa)/2);

α4 = ((αaa+αaam)/2);
w4 = waa;

rep4={w->w4, α-> α4, gAM->gAM4, gaM->gaM4, gAm->gAm4, gam->gam4};

Ww4 = Ww/. rep4;(* Mutant offspring produced. *)
	
(*** Class 5: Am/Am focal mutant. ***)

gAM5 = 0;
gaM5 = 0;
gAm5 = 1-μA;
gam5 = μA;

α5 = αAAm;
w5 = wAA;

rep5={w->w5, α-> α5, gAM->gAM5, gaM->gaM5, gAm->gAm5, gam->gam5};

Ww5 = Ww/. rep5;(* Mutant offspring produced. *)

(*** Class 6: Am/am focal mutant. ***)

gAM6 = 0;
gaM6 = 0;
gAm6 = ((μa+1-μA)/2);
gam6 = ((μA+1-μa)/2);

α6 = αAam;
w6 = wAa;

rep6={w->w6, α-> α6, gAM->gAM6, gaM->gaM6, gAm->gAm6, gam->gam6};

Ww6 = Ww/. rep6;(* Mutant offspring produced. *)

(*** Class 7: am/Am focal mutant. ***)

gAM7 = 0;
gaM7 = 0;
gAm7 = ((μa+1-μA)/2);
gam7 = ((μA+1-μa)/2);

α7 = αAam;
w7 = wAa;

rep7={w->w7, α-> α7, gAM->gAM7, gaM->gaM7, gAm->gAm7, gam->gam7};

Ww7 = Ww/. rep7;(* Mutant offspring produced. *)

(*** Class 8: am/am focal mutant. ***)

gAM8 = 0;
gaM8 = 0;
gAm8 = μa;
gam8 = 1-μa;

α8 = αaam;
w8 = waa;

rep8={w->w8, α-> α8, gAM->gAM8, gaM->gaM8, gAm->gAm8, gam->gam8};

Ww8 = Ww/. rep8;(* Mutant offspring produced. *)

(***** Replacement rules *****)

WeakSel={μA->μA (ϵ^2), s->s ϵ, μa->μa (ϵ^2)};
PFdec={p->p0+ϵ p1+((ϵ^2)/2)p2+ ((ϵ^3)/6)p3, Cc->C0+ϵ C1+((ϵ^2)/2)C2+((ϵ^3)/6)C3};
Gen = {fAA->(((1-p))^2)+Cc, fAa-> 2p(1-p)-2 Cc, faa->(p^2)+Cc,qAA->(((1-p))^2)+Cc,qAa-> 2p(1-p)-2 Cc,qaa->(p^2)+Cc};
GenDec={fAA->fAA0+ϵ fAA1+ ((ϵ^2)/2)fAA2 + ((ϵ^3)/6)fAA3,fAa->fAa0+ϵ fAa1+ ((ϵ^2)/2)fAa2 + ((ϵ^3)/6)fAa3,faa->faa0+ϵ faa1+ ((ϵ^2)/2)faa2 + ((ϵ^3)/6)faa3};
mono = {αAAm->αAA, αAam->αAa, αaam->αaa};
NoCondDep = {αAA->α, αAa->α, αaa->α};

(***** Invasion matrix *****)

W= (1/(F0 (qAA wAA(1-αAA δ)+qAa wAa(1-αAa δ)+qaa waa(1-αaa δ)))){Ww1,Ww2,Ww3,Ww4,Ww5,Ww6,Ww7,Ww8} // Transpose;

Wneut = W/.mono; (* Invasion matrix under neutrality *)

(** Matrices of class-specific effects **)

DAA=D[W,αAAm]/.mono;
DAa=D[W,αAam]/.mono;
Daa=D[W,αaam]/.mono;

(***** Genetic dynamics at condition locus in resident population. *****)

PAfBef=((fAA wAA+fAa (wAa/2))/(fAA wAA+fAa wAa + faa waa));
PafBef=((faa waa+fAa (wAa/2))/(fAA wAA+fAa wAa + faa waa));

PAf = PAfBef(1-μA)+PafBef μa;
Paf = PafBef(1-μa)+PAfBef μA;

Wfbar= fAA wAA(1-αAA δ)+fAa wAa(1-αAa δ)+faa waa(1-αaa δ);

pAhet = ((1-μA+μa)/2);
pahet=((1-μa+μA)/2);

fAANext= (1/Wfbar)(fAA  wAA (αAA(1-δ)(((1-μA))^2)+ (1-αAA)(1-μA) PAf) + fAa wAa(αAa (1-δ) pAhet pAhet + pAhet (1-αAa) PAf) + faa waa (αaa(1-δ) ((μa)^2)+ (1-αaa) μa PAf));
fAaNext =(1/Wfbar)(fAA  wAA(2 μA(1-μA)αAA(1-δ) + (1-αAA)((1-μA)Paf + μA PAf))+ fAa wAa (2 αAa (1-δ)pAhet  pahet + (1-αAa) (pAhet Paf+pahet PAf))+ faa waa (2 μa(1-μa)αaa(1-δ)+(1-αaa)((1-μa)PAf+μa Paf)));
faaNext=(1/Wfbar)(fAA wAA (αAA(1-δ) ((μA)^2)+(1-αAA) μA Paf)+fAa wAa(αAa (1-δ)pahet pahet+pahet(1-αAa)Paf) + faa waa (αaa(1-δ)(((1-μa))^2)+(1-αaa)(1-μa) Paf));

ΔfAA=fAANext-fAA;
ΔfAa=fAaNext-fAa;
Δfaa=faaNext - faa;
Δp=Δfaa+((ΔfAa)/2)/.Gen;
ΔC=faaNext-(((faaNext+(fAaNext/2)))^2)- Cc /.Gen;

(*** Numerical sanity check ***)

num={r->0.5, μA->0.001,μa->0.000, s->0.05, h->.25, αAA->1, αAa->1, αaa->0, F0->1, δ->0.4695};

freq={qAA->fAA, qAa->fAa, qaa->faa}/.FindRoot[ΔfAA==0&& ΔfAa==0&& Δfaa ==0/.num, {{fAA,0.998}, {fAa,0.001},{faa,0.001}}];

Wn = Wneut/.num/.freq;
Eigenvectors[Wn][[1]];

Qn =((%)/(Total[%])) // Chop
Eigenvectors[Wn // Transpose][[1]] //Chop;
Vn=((%)/(%.Qn));
S={Vn.DAA.Qn, Vn.DAa.Qn, Vn.Daa.Qn}/.freq/.num

Numerical investigations

Gradient iteration approach

Aa ={0.5,0.5,0.5}; (** Initial trait values. **)
pAa={-100,-100,-100}; (** Values in the previous time-step used for stopping criterion and so initialised to something large. **)
vA={Aa};

tol=(10^(-12)); (** Tolerance for stopping. **)
num={r->0.5, μA->0.001,μa->0.001, s->0.05, h->0.1,  F0->1, δ->.47}; (** Parameter values. **)

While[Total[(((Aa-pAa))^2)] > tol (* While the selfing rates keep changing ... *),
rA=(#[[1]]->#[[2]]&)/@Transpose@{{αAA,αAa,αaa}, Aa}; (* Replacement rule for selfing rates. *)

(*** RESIDENT EQUILIBRIUM ***)

Δf = {ΔfAA, ΔfAa, Δfaa}/.num /.rA; (* Change in genotypic frequencies at condition locus in resident pop. *)
ff={0.25,0.5,0.25}; (* Initial genotypic freqs. *)
pf={-100,-100,-100}; (* Genotypic frequencies in previous gen. for stopping criterion. *)

While[Total[(((ff-pf))^2)] > tol (* While genotypic frequencies keep changing ... *),
rf=(#[[1]]->#[[2]]&)/@Transpose@{{fAA,fAa,faa}, ff};(* Replacement rule for frequencies. *)
pf=ff; (* Update freqs. in previous gen. *)
ff+=Δf/.rf  (* Genotypic frequencies in the next generation. *)];

rq=(#[[1]]->#[[2]]&)/@Transpose@{{qAA,qAa,qaa},ff};  (* Replacement rule for genotypic frequencies at equilibrium. *)

(*** SELECTION GRADIENT ***)

Wn = Wneut/.num/.rq /.rA;
Qn = Eigenvectors[Wn][[1]];
Qn /=Total[Qn] // Chop;
Vn = Eigenvectors[Wn // Transpose][[1]] //Chop;
Vn /= Vn.Qn;
Sn = {Vn.DAA.Qn, Vn.DAa.Qn, Vn.Daa.Qn}/.rq/.num /.rA;

(*** CHANGE IN TRAIT VALUES ***)

pAa=Aa; (* Update trait values in previous timestep. *)
Aa+= dA{((Sn[[1]])/qAA),(( Sn[[2]])/qAa), ((Sn[[3]])/qaa)}/.rq/.num /.rA/.dA->0.05; (* Trait values move in the direction of the gradients. *)
(* Here, gradients are divided by frequencies of each genotype to speed up the process, because selection on rare genotypes is very weak. This does not change the sign of gradients, so this does not affect the equilibrium reached. *)
Aa = Clip[#,{0,1}] &/@ Aa; (* Trait values are clipped to be kept within bounds. *)
vA =Append[vA,Aa];
];