The Evolution of Sex Ratio Adjustment in the Presence of Sexually Antagonistic Selection

Gwylim S. Blackburn, Arianne Y. K. Albert, & Sarah P. Otto

Local stability analyses

Autosomal trait & modifier

Intro

Natural selection (s ) on female display trait genotypes:
(({{  TT  ,   Tt  ,   tt  }, {1, 1-hf*s, 1-s}}))
s=0 corresponds to no sexual anatagonism at the display trait locus

Natural and sexual selection on male display trait genotypes  (z):
(({{  TT  ,   Tt  ,   tt  }, {1-z, 1-hm*z, 1}}))

[In the text, s is defined as (Subscript[s, f]),while z is defined as (Subscript[s, m]).]

SRA by females when mating with TT, Tt, or tt males :
(({{      , TT, Tt      tt}, {g1
g2
g3
g4
g5
g6
g7
g8
g9
g10, m11
m21
m31
m41
m51
m61
m71
m81
m91
m101, m12      m13
m22      m23
m32      m33
m42      m43
m52      m53
m62     m63
m72      m73
m82      m83
m92      m93
m102   m103}}))

This matrix is written in a general form that can be used for both mate-SRA and self-SRA.  Appropriate choices of the (Subscript[m, ij]) allow conversion between this matrix and the (Subscript[D, ij]) used in the text.

Life cycle:
Natural selection, sexual selection, recombination, production of haploid gametes, pairing of gametes to create the new adult genotype frequencies

Female genotypes
g1 = TM/TM
g2 = TM/tM
g3 = tM/tM
g4 = TM/Tm
g5 = TM/tm
g6 = Tm/tM
g7 = tM/tm
g8 = Tm/Tm
g9 = Tm/tm
g10 = tm/tm

Although we keep the math general, in the text we assume that g5 = g6, so that there are no cis/trans differences in SRA.

Male genotypes
b1=TM/TM
b2=TM/tM
b3=tM/tM
b4=TM/Tm
b5=TM/tm
b6=Tm/tM
b7=tM/tm
b8=Tm/Tm
b9=Tm/tm
b10=tm/tm

Recursions [ENTER]

Off[General::Spell];
Off[General::spell1];

Natural selection in females

Mean fitness of females:

MFF=g1+g2(1-hf s)+g3(1-s)+g4+g5(1-hf s)+g6(1-hf s)+g7(1-s)+g8+g9(1-hf s)+g10(1-s);

Natural selection on different female genotypes:

gs1=g1/MFF;
gs2=g2(1-hf s)/MFF;
gs3=g3(1-s)/MFF;
gs4=g4/MFF;
gs5=g5(1-hf s)/MFF;
gs6=g6(1-hf s)/MFF;
gs7=g7(1-s)/MFF;
gs8=g8/MFF;
gs9=g9(1-hf s)/MFF;
gs10=g10(1-s)/MFF;

Factor[gs1+gs2+gs3+gs4+gs5+gs6+gs7+gs8+gs9+gs10]

Natural selection in males

Mean fitness of males:

MMF=b1(1-z)+b2(1-hm z)+b3+b4(1-z)+b5(1-hm z)+b6(1-hm z)+b7+b8(1-z)+b9(1-hm z)+b10;

Natural selection on different male genotypes:

bs1=b1(1-z)/MMF;
bs2=b2(1-hm z)/MMF;
bs3=b3/MMF;
bs4=b4(1-z)/MMF;
bs5=b5(1-hm z)/MMF;
bs6=b6(1-hm z)/MMF;
bs7=b7/MMF;
bs8=b8(1-z)/MMF;
bs9=b9(1-hm z)/MMF;
bs10=b10/MMF;

Factor[bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9+bs10]

Mating, recombination and syngamy

genotypes of daughters:

gr1=1/2*(1-r)gs5*1/2*bs2*m52+1/2*(1-r)gs5*1/2*r*bs6*m52+1/2*(1-r)gs5*1/2*(1-r)*bs5*m52+1/2*gs2*1/2*(1-r)*bs5*m22+1/2*gs2*1/2*r*bs6*m22+1/2*r*gs6*1/2*r*bs6*m62+1/2*r*gs6*1/2*bs2*m62+1/2*r*gs6*1/2*(1-r)*bs5*m62+gs1*1/2*r*bs6*m12+1/2*(1-r)gs5*1/2*bs4*m51+1/2*gs4*1/2*(1-r)*bs5*m42+1/2*gs4*1/2*r*bs6*m42+1/2*r*gs6*1/2*bs4*m61+1/2*gs4*1/2*bs4*m41+1/2*gs2*1/2*bs4*m21+1/2*gs4*1/2*bs2*m42+1/2*(1-r)gs5*bs1*m51+1/2*r*gs6*bs1*m61+gs1*1/2*(1-r)*bs5*m12+1/2*gs4*bs1*m41+gs1*1/2*bs4*m11+1/2*gs2*1/2*bs2*m22+1/2*gs2*bs1*m21+gs1*1/2*bs2*m12+gs1*bs1*m11;

gr2=1/2*r*gs5*1/2*bs2*m52+1/2*(1-r)gs5*1/2*bs2*m52+1/2*(1-r)gs5*1/2*(1-r)*bs6*m52+1/2*r*gs5*1/2*r*bs6*m52+1/2*(1-r)gs5*1/2*bs7*m53+1/2*(1-r)gs5*1/2*r*bs5*m52+1/2*r*gs5*1/2*(1-r)*bs5*m52+1/2*(1-r)gs5*bs3*m53+1/2*gs2*1/2*r*bs5*m22+1/2*gs2*1/2*(1-r)*bs5*m22+1/2*gs2*1/2*bs7*m23+1/2*gs2*1/2*(1-r)*bs6*m22+1/2*gs2*1/2*r*bs6*m22+1/2*gs7*1/2*(1-r)*bs5*m72+1/2*gs7*1/2*bs2*m72+1/2*gs7*1/2*r*bs6*m72+1/2*r*gs6*1/2*(1-r)*bs6*m62+1/2*(1-r)gs6*1/2*r*bs6*m62+1/2*r*gs6*1/2*bs2*m62+1/2*(1-r)gs6*1/2*bs2*m62+1/2*r*gs6*1/2*bs7*m63+1/2*r*gs6*1/2*r*bs5*m62+1/2*(1-r)gs6*1/2*(1-r)*bs5*m62+1/2*r*gs6*bs3*m63+gs3*1/2*(1-r)*bs5*m32+gs3*1/2*r*bs6*m32+gs1*1/2*(1-r)*bs6*m12+1/2*r*gs5*1/2*bs4*m51+1/2*gs4*1/2*r*bs5*m42+1/2*gs4*1/2*bs7*m43+1/2*gs4*1/2*(1-r)*bs6*m42+1/2*gs7*1/2*bs4*m71+1/2*(1-r)gs6*1/2*bs4*m61+1/2*gs2*1/2*bs4*m21+1/2*gs4*1/2*bs2*m42+1/2*gs4*bs3*m43+gs3*1/2*bs4*m31+1/2*r*gs5*bs1*m51+1/2*gs7*bs1*m71+1/2*(1-r)gs6*bs1*m61+gs1*1/2*r*bs5*m12+gs1*1/2*bs7*m13+1/2*gs2*bs2*m22+1/2*gs2*bs3*m23+gs3*1/2*bs2*m32+1/2*gs2*bs1*m21+gs1*1/2*bs2*m12+gs1*bs3*m13+gs3*bs1*m31;

gr3=1/2*r*gs5*1/2*bs2*m52+1/2*r*gs5*1/2*(1-r)*bs6*m52+1/2*r*gs5*1/2*bs7*m53+1/2*r*gs5*1/2*r*bs5*m52+1/2*r*gs5*bs3*m53+1/2*gs2*1/2*r*bs5*m22+1/2*gs2*1/2*bs7*m23+1/2*gs2*1/2*(1-r)*bs6*m22+1/2*gs7*1/2*r*bs5*m72+1/2*gs7*1/2*bs2*m72+1/2*gs7*1/2*(1-r)*bs6*m72+1/2*gs7*bs3*m73+1/2*gs7*1/2*bs7*m73+1/2*(1-r)gs6*1/2*(1-r)*bs6*m62+1/2*(1-r)gs6*1/2*bs2*m62+1/2*(1-r)gs6*1/2*bs7*m63+1/2*(1-r)gs6*1/2*r*bs5*m62+1/2*(1-r)gs6*bs3*m63+gs3*1/2*r*bs5*m32+gs3*1/2*bs7*m33+gs3*1/2*(1-r)*bs6*m32+1/2*gs2*1/2*bs2*m22+1/2*gs2*bs3*m23+gs3*1/2*bs2*m32+gs3*bs3*m33;

gr4=1/2*r*gs5*1/2*bs2*m52+1/2*r*gs5*1/2*r*bs6*m52+1/2*(1-r)gs5*1/2*(1-r)*bs6*m52+1/2*r*gs5*1/2*(1-r)*bs5*m52+1/2*(1-r)gs5*1/2*r*bs5*m52+1/2*r*gs6*1/2*bs9*m62+1/2*gs2*1/2*r*bs5*m22+1/2*gs2*1/2*bs9*m22+1/2*gs2*1/2*(1-r)*bs6*m22+1/2*gs9*1/2*r*bs6*m92+1/2*gs9*1/2*bs2*m92+1/2*gs9*1/2*(1-r)*bs5*m92+1/2*(1-r)gs5*1/2*bs9*m52+1/2*(1-r)gs6*1/2*r*bs6*m62+1/2*r*gs6*1/2*(1-r)*bs6*m62+1/2*(1-r)gs6*1/2*bs2*m62+1/2*(1-r)gs6*1/2*(1-r)*bs5*m62+1/2*r*gs6*1/2*r*bs5*m62+gs1*1/2*(1-r)*bs6*m12+1/2*r*gs5*1/2*bs4*m51+1/2*(1-r)gs5*1/2*bs4*m51+1/2*(1-r)gs5*bs8*m51+1/2*gs4*1/2*r*bs5*m42+1/2*gs4*1/2*(1-r)*bs5*m42+1/2*gs4*1/2*bs9*m42+1/2*gs4*1/2*(1-r)*bs6*m42+1/2*gs4*1/2*r*bs6*m42+1/2*gs9*1/2*bs4*m91+1/2*(1-r)gs6*1/2*bs4*m61+1/2*r*gs6*1/2*bs4*m61+1/2*r*gs6*bs8*m61+gs8*1/2*(1-r)*bs5*m82+gs8*1/2*r*bs6*m82+1/2*gs4*bs4*m41+1/2*gs4*bs8*m41+gs8*1/2*bs4*m81+1/2*gs2*1/2*bs4*m21+1/2*gs2*bs8*m21+1/2*gs4*1/2*bs2*m42+gs8*1/2*bs2*m82+1/2*r*gs5*bs1*m51+1/2*gs9*bs1*m91+1/2*(1-r)gs6*bs1*m61+gs1*1/2*r*bs5*m12+gs1*1/2*bs9*m12+1/2*gs4*bs1*m41+gs1*1/2*bs4*m11+gs1*bs8*m11+gs8*bs1*m81;

gr5=1/2*(1-r)gs5*1/2*bs2*m52+1/2*(1-r)gs5*r*bs6*m52+1/2*(1-r)gs5*1/2*bs7*m53+1/2*(1-r)gs5(1-r)*bs5*m52+1/2*r*gs6*1/2*bs9*m62+1/2*r*gs6*bs10*m63+1/2*gs2*1/2*(1-r)*bs5*m22+1/2*gs2*1/2*bs7*m23+1/2*gs2*1/2*bs9*m22+1/2*gs2*1/2*r*bs6*m22+1/2*gs2*bs10*m23+1/2*gs7*1/2*(1-r)*bs5*m72+1/2*gs7*1/2*bs2*m72+1/2*gs7*1/2*r*bs6*m72+1/2*gs9*1/2*r*bs6*m92+1/2*gs9*1/2*bs2*m92+1/2*gs9*1/2*(1-r)*bs5*m92+1/2*(1-r)gs5*1/2*bs9*m52+1/2*(1-r)gs5*bs10*m53+1/2*r*gs6*r*bs6*m62+1/2*r*gs6*1/2*bs2*m62+1/2*r*gs6*1/2*bs7*m63+1/2*r*gs6(1-r)*bs5*m62+gs10*1/2*r*bs6*m102+gs10*1/2*bs2*m102+gs10*1/2*(1-r)*bs5*m102+gs1*1/2*r*bs6*m12+1/2*(1-r)gs5*1/2*bs4*m51+1/2*gs4*1/2(1-r)*bs5*m42+1/2gs4*1/2*bs7*m43+1/2gs4*1/2*bs9*m42+1/2gs4*1/2r*bs6*m42+1/2gs4*bs10*m43+1/2gs7*1/2*bs4*m71+1/2gs9*1/2*bs4*m91+1/2r*gs6*1/2*bs4*m61+gs10*1/2*bs4*m101+1/2(1-r)gs5*bs1*m51+1/2gs7*bs1*m71+1/2gs9*bs1*m91+1/2r*gs6*bs1*m61+gs10*bs1*m101+gs1*1/2(1-r)*bs5*m12+gs1*1/2*bs7*m13+gs1*1/2*bs9*m12+gs1*bs10*m13;

gr6=1/2r*gs5*1/2*bs2*m52+1/2r*gs5(1-r)*bs6*m52+1/2r*gs5*1/2*bs7*m53+1/2r*gs5*r*bs5*m52+1/2r*gs5*bs3*m53+1/2(1-r)gs6*1/2*bs9*m62+1/2gs2*1/2r*bs5*m22+1/2gs2*1/2*bs9*m22+1/2gs2*1/2(1-r)*bs6*m22+1/2gs7*1/2r*bs5*m72+1/2gs7*1/2*bs9*m72+1/2gs7*1/2(1-r)*bs6*m72+1/2gs9*1/2(1-r)*bs6*m92+1/2gs9*1/2*bs2*m92+1/2gs9*1/2*bs7*m93+1/2gs9*1/2r*bs5*m92+1/2gs9*bs3*m93+1/2r*gs5*1/2*bs9*m52+1/2(1-r)gs6(1-r)*bs6*m62+1/2(1-r)gs6*1/2*bs2*m62+1/2(1-r)gs6*1/2*bs7*m63+1/2(1-r)gs6*r*bs5*m62+1/2(1-r)gs6*bs3*m63+gs3*1/2r*bs5*m32+gs3*1/2*bs9*m32+gs3*1/2(1-r)*bs6*m32+1/2r*gs5*1/2*bs4*m51+1/2r*gs5*bs8*m51+1/2gs4*1/2r*bs5*m42+1/2gs4*1/2*bs7*m43+1/2gs4*1/2(1-r)*bs6*m42+1/2gs7*1/2*bs4*m71+1/2gs7*bs8*m71+1/2(1-r)gs6*1/2*bs4*m61+1/2(1-r)gs6*bs8*m61+gs8*1/2r*bs5*m82+gs8*1/2*bs7*m83+gs8*1/2(1-r)*bs6*m82+1/2gs2*1/2*bs4*m21+1/2gs2*bs8*m21+1/2gs4*1/2*bs2*m42+1/2gs4*bs3*m43+gs3*1/2*bs4*m31+gs3*bs8*m31+gs8*1/2*bs2*m82+gs8*bs3*m83;

gr7=1/2*(1-r)gs5*1/2*bs2*m52+1/2(1-r)gs5*1/2(1-r)*bs6*m52+1/2r*gs5*1/2r*bs6*m52+1/2(1-r)gs5*1/2*bs7*m53+1/2r*gs5*1/2*bs7*m53+1/2(1-r)gs5*1/2r*bs5*m52+1/2r*gs5*1/2(1-r)*bs5*m52+1/2(1-r)gs5*bs3*m53+1/2(1-r)gs6*1/2*bs9*m62+1/2(1-r)gs6*bs10*m63+1/2gs2*1/2(1-r)*bs5*m22+1/2gs2*1/2*bs7*m23+1/2gs2*1/2*bs9*m22+1/2gs2*1/2r*bs6*m22+1/2gs2*bs10*m23+1/2gs7*1/2(1-r)*bs5*m72+1/2gs7*1/2r*bs5*m72+1/2gs7*1/2*bs2*m72+1/2gs7*1/2*bs9*m72+1/2gs7*1/2r*bs6*m72+1/2gs7*1/2(1-r)*bs6*m72+1/2gs7*bs10*m73+1/2gs7*bs3*m73+1/2gs7*bs7*m73+1/2gs9*1/2(1-r)*bs6*m92+1/2gs9*1/2*bs2*m92+1/2gs9*1/2*bs7*m93+1/2gs9*1/2r*bs5*m92+1/2gs9*bs3*m93+1/2r*gs5*1/2*bs9*m52+1/2r*gs5*bs10*m53+1/2r*gs6*1/2(1-r)*bs6*m62+1/2(1-r)gs6*1/2r*bs6*m62+1/2r*gs6*1/2*bs2*m62+1/2r*gs6*1/2*bs7*m63+1/2(1-r)gs6*1/2*bs7*m63+1/2r*gs6*1/2r*bs5*m62+1/2(1-r)gs6*1/2(1-r)*bs5*m62+1/2r*gs6*bs3*m63+gs10*1/2(1-r)*bs6*m102+gs10*1/2*bs2*m102+gs10*1/2*bs7*m103+gs10*1/2r*bs5*m102+gs10*bs3*m103+gs3*1/2(1-r)*bs5*m32+gs3*1/2*bs7*m33+gs3*1/2*bs9*m32+gs3*1/2r*bs6*m32+gs3*bs10*m33;

gr8=1/2(1-r)gs6*1/2(1-r)*bs6*m62+1/2(1-r)gs6*1/2*bs4*m61+1/2(1-r)gs6*1/2*bs9*m62+1/2(1-r)gs6*1/2r*bs5*m62+1/2(1-r)gs6*bs8*m61+1/2gs4*1/2(1-r)*bs6*m42+1/2gs4*1/2*bs4*m41+1/2gs4*1/2*bs9*m42+1/2gs4*1/2r*bs5*m42+1/2gs4*bs8*m41+1/2gs9*1/2(1-r)*bs6*m92+1/2gs9*1/2*bs4*m91+1/2gs9*1/2*bs9*m92+1/2gs9*1/2r*bs5*m92+1/2gs9*bs8*m91+1/2r*gs5*1/2(1-r)*bs6*m52+1/2r*gs5*1/2*bs4*m51+1/2r*gs5*1/2*bs9*m52+1/2r*gs5*1/2r*bs5*m52+1/2r*gs5*bs8*m51+gs8*1/2(1-r)*bs6*m82+gs8*1/2*bs4*m81+gs8*1/2*bs9*m82+gs8*1/2r*bs5*m82+gs8*bs8*m81;

gr9=1/2(1-r)gs5*1/2(1-r)*bs6*m52+1/2r*gs5*1/2r*bs6*m52+1/2(1-r)gs5*1/2*bs4*m51+1/2(1-r)gs5*1/2*bs9*m52+1/2r*gs5*1/2*bs9*m52+1/2(1-r)gs5*1/2r*bs5*m52+1/2r*gs5*1/2(1-r)*bs5*m52+1/2(1-r)gs5*bs8*m51+1/2(1-r)gs6*1/2*bs7*m63+1/2(1-r)gs6*bs10*m63+1/2gs4*1/2(1-r)*bs5*m42+1/2gs4*1/2*bs7*m43+1/2gs4*1/2*bs9*m42+1/2gs4*1/2r*bs6*m42+1/2gs4*bs10*m43+1/2gs7*1/2(1-r)*bs6*m72+1/2gs7*1/2*bs4*m71+1/2gs7*1/2*bs9*m72+1/2gs7*1/2r*bs5*m72+1/2gs7*bs8*m71+1/2gs9*1/2(1-r)*bs5*m92+1/2gs9*1/2r*bs5*m92+1/2gs9*1/2*bs4*m91+1/2gs9*1/2*bs7*m93+1/2gs9*1/2r*bs6*m92+1/2gs9*1/2(1-r)*bs6*m92+1/2gs9*bs10*m93+1/2gs9*bs8*m91+1/2gs9*bs9*m92+1/2r*gs5*1/2*bs7*m53+1/2r*gs5*bs10*m53+1/2r*gs6*1/2(1-r)*bs6*m62+1/2(1-r)gs6*1/2r*bs6*m62+1/2r*gs6*1/2*bs4*m61+1/2r*gs6*1/2*bs9*m62+1/2(1-r)gs6*1/2*bs9*m62+1/2r*gs6*1/2r*bs5*m62+1/2(1-r)gs6*1/2(1-r)*bs5*m62+1/2r*gs6*bs8*m61+gs10*1/2(1-r)*bs6*m102+gs10*1/2*bs4*m101+gs10*1/2*bs9*m102+gs10*1/2r*bs5*m102+gs10*bs8*m101+gs8*1/2(1-r)*bs5*m82+gs8*1/2*bs7*m83+gs8*1/2*bs9*m82+gs8*1/2r*bs6*m82+gs8*bs10*m83;

gr10=1/2(1-r)gs5*1/2(1-r)*bs5*m52+1/2(1-r)gs5*1/2*bs7*m53+1/2(1-r)gs5*1/2*bs9*m52+1/2(1-r)gs5*1/2r*bs6*m52+1/2(1-r)gs5*bs10*m53+1/2gs7*1/2(1-r)*bs5*m72+1/2gs7*1/2*bs7*m73+1/2gs7*1/2*bs9*m72+1/2gs7*1/2r*bs6*m72+1/2gs7*bs10*m73+1/2gs9*1/2(1-r)*bs5*m92+1/2gs9*1/2*bs7*m93+1/2gs9*1/2*bs9*m92+1/2gs9*1/2r*bs6*m92+1/2gs9*bs10*m93+1/2r*gs6*1/2(1-r)*bs5*m62+1/2r*gs6*1/2*bs7*m63+1/2r*gs6*1/2*bs9*m62+1/2r*gs6*1/2r*bs6*m62+1/2r*gs6*bs10*m63+gs10*1/2(1-r)*bs5*m102+gs10*1/2*bs7*m103+gs10*1/2*bs9*m102+gs10*1/2r*bs6*m102+gs10*bs10*m103;

Factor[gr1+gr2+gr3+gr4+gr5+gr6+gr7+gr8+gr9+gr10/.{m11->1/2,m21->1/2,m31->1/2,m41->1/2,m51->1/2,m61->1/2,m71->1/2,m81->1/2,m91->1/2,m101->1/2,m12->1/2,m22->1/2,m32->1/2,m42->1/2,m52->1/2,m62->1/2,m72->1/2,m82->1/2,m92->1/2,m102->1/2,m13->1/2,m23->1/2,m33->1/2,m43->1/2,m53->1/2,m63->1/2,m73->1/2,m83->1/2,m93->1/2,m103->1/2}]

genotypes of sons:

br1=1/2(1-r)gs5*1/2*bs2(1-m52)+1/2(1-r)gs5*1/2r*bs6(1-m52)+1/2(1-r)gs5*1/2(1-r)*bs5(1-m52)+1/2gs2*1/2(1-r)*bs5(1-m22)+1/2gs2*1/2r*bs6(1-m22)+1/2r*gs6*1/2r*bs6(1-m62)+1/2r*gs6*1/2*bs2(1-m62)+1/2r*gs6*1/2(1-r)*bs5(1-m62)+gs1*1/2r*bs6(1-m12)+1/2(1-r)gs5*1/2*bs4(1-m51)+1/2gs4*1/2(1-r)*bs5(1-m42)+1/2gs4*1/2r*bs6(1-m42)+1/2r*gs6*1/2*bs4(1-m61)+1/2gs4*1/2*bs4(1-m41)+1/2gs2*1/2*bs4(1-m21)+1/2gs4*1/2*bs2(1-m42)+1/2(1-r)gs5*bs1(1-m51)+1/2r*gs6*bs1(1-m61)+gs1*1/2(1-r)*bs5(1-m12)+1/2gs4*bs1(1-m41)+gs1*1/2*bs4(1-m11)+1/2gs2*1/2*bs2(1-m22)+1/2gs2*bs1(1-m21)+gs1*1/2*bs2(1-m12)+gs1*bs1(1-m11);

br2=1/2r*gs5*1/2*bs2(1-m52)+1/2(1-r)gs5*1/2*bs2(1-m52)+1/2(1-r)gs5*1/2(1-r)*bs6(1-m52)+1/2r*gs5*1/2r*bs6(1-m52)+1/2(1-r)gs5*1/2*bs7(1-m53)+1/2(1-r)gs5*1/2r*bs5(1-m52)+1/2r*gs5*1/2(1-r)*bs5(1-m52)+1/2(1-r)gs5*bs3(1-m53)+1/2gs2*1/2r*bs5(1-m22)+1/2gs2*1/2(1-r)*bs5(1-m22)+1/2gs2*1/2*bs7(1-m23)+1/2gs2*1/2(1-r)*bs6(1-m22)+1/2gs2*1/2r*bs6(1-m22)+1/2gs7*1/2(1-r)*bs5(1-m72)+1/2gs7*1/2*bs2(1-m72)+1/2gs7*1/2r*bs6(1-m72)+1/2r*gs6*1/2(1-r)*bs6(1-m62)+1/2(1-r)gs6*1/2r*bs6(1-m62)+1/2r*gs6*1/2*bs2(1-m62)+1/2(1-r)gs6*1/2*bs2(1-m62)+1/2r*gs6*1/2*bs7(1-m63)+1/2r*gs6*1/2r*bs5(1-m62)+1/2(1-r)gs6*1/2(1-r)*bs5(1-m62)+1/2r*gs6*bs3(1-m63)+gs3*1/2(1-r)*bs5(1-m32)+gs3*1/2r*bs6(1-m32)+gs1*1/2(1-r)*bs6(1-m12)+1/2r*gs5*1/2*bs4(1-m51)+1/2gs4*1/2r*bs5(1-m42)+1/2gs4*1/2*bs7(1-m43)+1/2gs4*1/2(1-r)*bs6(1-m42)+1/2gs7*1/2*bs4(1-m71)+1/2(1-r)gs6*1/2*bs4(1-m61)+1/2gs2*1/2*bs4(1-m21)+1/2gs4*1/2*bs2(1-m42)+1/2gs4*bs3(1-m43)+gs3*1/2*bs4(1-m31)+1/2r*gs5*bs1(1-m51)+1/2gs7*bs1(1-m71)+1/2(1-r)gs6*bs1(1-m61)+gs1*1/2r*bs5(1-m12)+gs1*1/2*bs7(1-m13)+1/2gs2*bs2(1-m22)+1/2gs2*bs3(1-m23)+gs3*1/2*bs2(1-m32)+1/2gs2*bs1(1-m21)+gs1*1/2*bs2(1-m12)+gs1*bs3(1-m13)+gs3*bs1(1-m31);

br3=1/2*r*gs5*1/2*bs2(1-m52)+1/2*r*gs5*1/2*(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*bs7(1-m53)+1/2*r*gs5*1/2*r*bs5(1-m52)+1/2*r*gs5*bs3(1-m53)+1/2*gs2*1/2*r*bs5(1-m22)+1/2*gs2*1/2*bs7(1-m23)+1/2*gs2*1/2*(1-r)*bs6(1-m22)+1/2*gs7*1/2*r*bs5(1-m72)+1/2*gs7*1/2*bs2(1-m72)+1/2*gs7*1/2*(1-r)*bs6(1-m72)+1/2*gs7*bs3(1-m73)+1/2*gs7*1/2*bs7(1-m73)+1/2*(1-r)gs6*1/2*(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*bs2(1-m62)+1/2*(1-r)gs6*1/2*bs7(1-m63)+1/2*(1-r)gs6*1/2*r*bs5(1-m62)+1/2*(1-r)gs6*bs3(1-m63)+gs3*1/2*r*bs5(1-m32)+gs3*1/2*bs7(1-m33)+gs3*1/2*(1-r)*bs6(1-m32)+1/2*gs2*1/2*bs2(1-m22)+1/2*gs2*bs3(1-m23)+gs3*1/2*bs2(1-m32)+gs3*bs3(1-m33);

br4=1/2*r*gs5*1/2*bs2(1-m52)+1/2*r*gs5*1/2*r*bs6(1-m52)+1/2*(1-r)gs5*1/2*(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*(1-r)*bs5(1-m52)+1/2*(1-r)gs5*1/2*r*bs5(1-m52)+1/2*r*gs6*1/2*bs9(1-m62)+1/2*gs2*1/2*r*bs5(1-m22)+1/2*gs2*1/2*bs9(1-m22)+1/2*gs2*1/2*(1-r)*bs6(1-m22)+1/2*gs9*1/2*r*bs6(1-m92)+1/2*gs9*1/2*bs2(1-m92)+1/2*gs9*1/2*(1-r)*bs5(1-m92)+1/2*(1-r)gs5*1/2*bs9(1-m52)+1/2*(1-r)gs6*1/2*r*bs6(1-m62)+1/2*r*gs6*1/2*(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*bs2(1-m62)+1/2*(1-r)gs6*1/2*(1-r)*bs5(1-m62)+1/2*r*gs6*1/2*r*bs5(1-m62)+gs1*1/2*(1-r)*bs6(1-m12)+1/2*r*gs5*1/2*bs4(1-m51)+1/2*(1-r)gs5*1/2*bs4(1-m51)+1/2*(1-r)gs5*bs8(1-m51)+1/2*gs4*1/2*r*bs5(1-m42)+1/2*gs4*1/2*(1-r)*bs5(1-m42)+1/2*gs4*1/2*bs9(1-m42)+1/2*gs4*1/2*(1-r)*bs6(1-m42)+1/2*gs4*1/2*r*bs6(1-m42)+1/2*gs9*1/2*bs4(1-m91)+1/2*(1-r)gs6*1/2*bs4(1-m61)+1/2*r*gs6*1/2*bs4(1-m61)+1/2*r*gs6*bs8(1-m61)+gs8*1/2*(1-r)*bs5(1-m82)+gs8*1/2*r*bs6(1-m82)+1/2*gs4*bs4(1-m41)+1/2*gs4*bs8(1-m41)+gs8*1/2*bs4(1-m81)+1/2*gs2*1/2*bs4(1-m21)+1/2*gs2*bs8(1-m21)+1/2*gs4*1/2*bs2(1-m42)+gs8*1/2*bs2(1-m82)+1/2*r*gs5*bs1(1-m51)+1/2*gs9*bs1(1-m91)+1/2*(1-r)gs6*bs1(1-m61)+gs1*1/2*r*bs5(1-m12)+gs1*1/2*bs9(1-m12)+1/2*gs4*bs1(1-m41)+gs1*1/2*bs4(1-m11)+gs1*bs8(1-m11)+gs8*bs1(1-m81);

br5=1/2*(1-r)gs5*1/2*bs2(1-m52)+1/2*(1-r)gs5*r*bs6(1-m52)+1/2*(1-r)gs5*1/2*bs7(1-m53)+1/2*(1-r)gs5(1-r)*bs5(1-m52)+1/2*r*gs6*1/2*bs9(1-m62)+1/2*r*gs6*bs10(1-m63)+1/2*gs2*1/2*(1-r)*bs5(1-m22)+1/2*gs2*1/2*bs7(1-m23)+1/2*gs2*1/2*bs9(1-m22)+1/2*gs2*1/2*r*bs6(1-m22)+1/2*gs2*bs10(1-m23)+1/2*gs7*1/2*(1-r)*bs5(1-m72)+1/2*gs7*1/2*bs2(1-m72)+1/2*gs7*1/2*r*bs6(1-m72)+1/2*gs9*1/2*r*bs6(1-m92)+1/2*gs9*1/2*bs2(1-m92)+1/2*gs9*1/2*(1-r)*bs5(1-m92)+1/2*(1-r)gs5*1/2*bs9(1-m52)+1/2*(1-r)gs5*bs10(1-m53)+1/2*r*gs6*r*bs6(1-m62)+1/2*r*gs6*1/2*bs2(1-m62)+1/2*r*gs6*1/2*bs7(1-m63)+1/2*r*gs6(1-r)*bs5(1-m62)+gs10*1/2*r*bs6(1-m102)+gs10*1/2*bs2(1-m102)+gs10*1/2*(1-r)*bs5(1-m102)+gs1*1/2*r*bs6(1-m12)+1/2*(1-r)gs5*1/2*bs4(1-m51)+1/2*gs4*1/2*(1-r)*bs5(1-m42)+1/2*gs4*1/2*bs7(1-m43)+1/2*gs4*1/2*bs9(1-m42)+1/2*gs4*1/2*r*bs6(1-m42)+1/2*gs4*bs10(1-m43)+1/2*gs7*1/2*bs4(1-m71)+1/2*gs9*1/2*bs4(1-m91)+1/2*r*gs6*1/2*bs4(1-m61)+gs10*1/2*bs4(1-m101)+1/2*(1-r)gs5*bs1(1-m51)+1/2*gs7*bs1(1-m71)+1/2*gs9*bs1(1-m91)+1/2*r*gs6*bs1(1-m61)+gs10*bs1(1-m101)+gs1*1/2*(1-r)*bs5(1-m12)+gs1*1/2*bs7(1-m13)+gs1*1/2*bs9(1-m12)+gs1*bs10(1-m13);

br6=1/2*r*gs5*1/2*bs2(1-m52)+1/2*r*gs5(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*bs7(1-m53)+1/2*r*gs5*r*bs5(1-m52)+1/2*r*gs5*bs3(1-m53)+1/2*(1-r)gs6*1/2*bs9(1-m62)+1/2*gs2*1/2*r*bs5(1-m22)+1/2*gs2*1/2*bs9(1-m22)+1/2*gs2*1/2*(1-r)*bs6(1-m22)+1/2*gs7*1/2*r*bs5(1-m72)+1/2*gs7*1/2*bs9(1-m72)+1/2*gs7*1/2*(1-r)*bs6(1-m72)+1/2*gs9*1/2*(1-r)*bs6(1-m92)+1/2*gs9*1/2*bs2(1-m92)+1/2*gs9*1/2*bs7(1-m93)+1/2*gs9*1/2*r*bs5(1-m92)+1/2*gs9*bs3(1-m93)+1/2*r*gs5*1/2*bs9(1-m52)+1/2*(1-r)gs6(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*bs2(1-m62)+1/2*(1-r)gs6*1/2*bs7(1-m63)+1/2*(1-r)gs6*r*bs5(1-m62)+1/2*(1-r)gs6*bs3(1-m63)+gs3*1/2*r*bs5(1-m32)+gs3*1/2*bs9(1-m32)+gs3*1/2*(1-r)*bs6(1-m32)+1/2*r*gs5*1/2*bs4(1-m51)+1/2*r*gs5*bs8(1-m51)+1/2*gs4*1/2*r*bs5(1-m42)+1/2*gs4*1/2*bs7(1-m43)+1/2*gs4*1/2*(1-r)*bs6(1-m42)+1/2*gs7*1/2*bs4(1-m71)+1/2*gs7*bs8(1-m71)+1/2*(1-r)gs6*1/2*bs4(1-m61)+1/2*(1-r)gs6*bs8(1-m61)+gs8*1/2*r*bs5(1-m82)+gs8*1/2*bs7(1-m83)+gs8*1/2*(1-r)*bs6(1-m82)+1/2*gs2*1/2*bs4(1-m21)+1/2*gs2*bs8(1-m21)+1/2*gs4*1/2*bs2(1-m42)+1/2*gs4*bs3(1-m43)+gs3*1/2*bs4(1-m31)+gs3*bs8(1-m31)+gs8*1/2*bs2(1-m82)+gs8*bs3(1-m83);

br7=1/2*(1-r)gs5*1/2*bs2(1-m52)+1/2*(1-r)gs5*1/2*(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*r*bs6(1-m52)+1/2*(1-r)gs5*1/2*bs7(1-m53)+1/2*r*gs5*1/2*bs7(1-m53)+1/2*(1-r)gs5*1/2*r*bs5(1-m52)+1/2*r*gs5*1/2*(1-r)*bs5(1-m52)+1/2*(1-r)gs5*bs3(1-m53)+1/2*(1-r)gs6*1/2*bs9(1-m62)+1/2*(1-r)gs6*bs10(1-m63)+1/2*gs2*1/2*(1-r)*bs5(1-m22)+1/2*gs2*1/2*bs7(1-m23)+1/2*gs2*1/2*bs9(1-m22)+1/2*gs2*1/2*r*bs6(1-m22)+1/2*gs2*bs10(1-m23)+1/2*gs7*1/2*(1-r)*bs5(1-m72)+1/2*gs7*1/2*r*bs5(1-m72)+1/2*gs7*1/2*bs2(1-m72)+1/2*gs7*1/2*bs9(1-m72)+1/2*gs7*1/2*r*bs6(1-m72)+1/2gs7*1/2*(1-r)*bs6(1-m72)+1/2*gs7*bs10(1-m73)+1/2*gs7*bs3(1-m73)+1/2*gs7*bs7(1-m73)+1/2*gs9*1/2*(1-r)*bs6(1-m92)+1/2*gs9*1/2*bs2(1-m92)+1/2*gs9*1/2*bs7(1-m93)+1/2*gs9*1/2*r*bs5(1-m92)+1/2*gs9*bs3(1-m93)+1/2*r*gs5*1/2*bs9(1-m52)+1/2*r*gs5*bs10(1-m53)+1/2*r*gs6*1/2*(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*r*bs6(1-m62)+1/2*r*gs6*1/2*bs2(1-m62)+1/2*r*gs6*1/2*bs7(1-m63)+1/2*(1-r)gs6*1/2*bs7(1-m63)+1/2*r*gs6*1/2*r*bs5(1-m62)+1/2*(1-r)gs6*1/2*(1-r)*bs5(1-m62)+1/2*r*gs6*bs3(1-m63)+gs10*1/2*(1-r)*bs6(1-m102)+gs10*1/2*bs2(1-m102)+gs10*1/2*bs7(1-m103)+gs10*1/2*r*bs5(1-m102)+gs10*bs3(1-m103)+gs3*1/2*(1-r)*bs5(1-m32)+gs3*1/2*bs7(1-m33)+gs3*1/2*bs9(1-m32)+gs3*1/2*r*bs6(1-m32)+gs3*bs10(1-m33);

br8=1/2*(1-r)gs6*1/2*(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*bs4(1-m61)+1/2*(1-r)gs6*1/2*bs9(1-m62)+1/2*(1-r)gs6*1/2*r*bs5(1-m62)+1/2*(1-r)gs6*bs8(1-m61)+1/2*gs4*1/2*(1-r)*bs6(1-m42)+1/2*gs4*1/2*bs4(1-m41)+1/2*gs4*1/2*bs9(1-m42)+1/2*gs4*1/2*r*bs5(1-m42)+1/2*gs4*bs8(1-m41)+1/2*gs9*1/2*(1-r)*bs6(1-m92)+1/2*gs9*1/2*bs4(1-m91)+1/2*gs9*1/2*bs9(1-m92)+1/2*gs9*1/2*r*bs5(1-m92)+1/2*gs9*bs8(1-m91)+1/2*r*gs5*1/2*(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*bs4(1-m51)+1/2*r*gs5*1/2*bs9(1-m52)+1/2*r*gs5*1/2*r*bs5(1-m52)+1/2*r*gs5*bs8(1-m51)+gs8*1/2*(1-r)*bs6(1-m82)+gs8*1/2*bs4(1-m81)+gs8*1/2*bs9(1-m82)+gs8*1/2*r*bs5(1-m82)+gs8*bs8(1-m81);

br9=1/2*(1-r)gs5*1/2*(1-r)*bs6(1-m52)+1/2*r*gs5*1/2*r*bs6(1-m52)+1/2*(1-r)gs5*1/2*bs4(1-m51)+1/2*(1-r)gs5*1/2*bs9(1-m52)+1/2*r*gs5*1/2*bs9(1-m52)+1/2*(1-r)gs5*1/2*r*bs5(1-m52)+1/2*r*gs5*1/2*(1-r)*bs5(1-m52)+1/2*(1-r)gs5*bs8(1-m51)+1/2*(1-r)gs6*1/2*bs7(1-m63)+1/2*(1-r)gs6*bs10(1-m63)+1/2*gs4*1/2*(1-r)*bs5(1-m42)+1/2*gs4*1/2*bs7(1-m43)+1/2*gs4*1/2*bs9(1-m42)+1/2*gs4*1/2*r*bs6(1-m42)+1/2*gs4*bs10(1-m43)+1/2*gs7*1/2*(1-r)*bs6(1-m72)+1/2*gs7*1/2*bs4(1-m71)+1/2*gs7*1/2*bs9(1-m72)+1/2*gs7*1/2*r*bs5(1-m72)+1/2*gs7*bs8(1-m71)+1/2*gs9*1/2*(1-r)*bs5(1-m92)+1/2*gs9*1/2*r*bs5(1-m92)+1/2*gs9*1/2*bs4(1-m91)+1/2*gs9*1/2*bs7(1-m93)+1/2*gs9*1/2*r*bs6(1-m92)+1/2gs9*1/2*(1-r)*bs6(1-m92)+1/2*gs9*bs10(1-m93)+1/2*gs9*bs8(1-m91)+1/2*gs9*bs9(1-m92)+1/2*r*gs5*1/2*bs7(1-m53)+1/2*r*gs5*bs10(1-m53)+1/2*r*gs6*1/2*(1-r)*bs6(1-m62)+1/2*(1-r)gs6*1/2*r*bs6(1-m62)+1/2*r*gs6*1/2*bs4(1-m61)+1/2*r*gs6*1/2*bs9(1-m62)+1/2*(1-r)gs6*1/2*bs9(1-m62)+1/2*r*gs6*1/2*r*bs5(1-m62)+1/2*(1-r)gs6*1/2*(1-r)*bs5(1-m62)+1/2*r*gs6*bs8(1-m61)+gs10*1/2*(1-r)*bs6(1-m102)+gs10*1/2*bs4(1-m101)+gs10*1/2*bs9(1-m102)+gs10*1/2*r*bs5(1-m102)+gs10*bs8(1-m101)+gs8*1/2*(1-r)*bs5(1-m82)+gs8*1/2*bs7(1-m83)+gs8*1/2*bs9(1-m82)+gs8*1/2*r*bs6(1-m82)+gs8*bs10(1-m83);

br10=1/2*(1-r)gs5*1/2*(1-r)*bs5(1-m52)+1/2*(1-r)gs5*1/2*bs7(1-m53)+1/2*(1-r)gs5*1/2*bs9(1-m52)+1/2*(1-r)gs5*1/2*r*bs6(1-m52)+1/2*(1-r)gs5*bs10(1-m53)+1/2*gs7*1/2*(1-r)*bs5(1-m72)+1/2*gs7*1/2*bs7(1-m73)+1/2*gs7*1/2*bs9(1-m72)+1/2*gs7*1/2*r*bs6(1-m72)+1/2*gs7*bs10(1-m73)+1/2*gs9*1/2*(1-r)*bs5(1-m92)+1/2*gs9*1/2*bs7(1-m93)+1/2*gs9*1/2*bs9(1-m92)+1/2*gs9*1/2*r*bs6(1-m92)+1/2*gs9*bs10(1-m93)+1/2*r*gs6*1/2*(1-r)*bs5(1-m62)+1/2*r*gs6*1/2*bs7(1-m63)+1/2*r*gs6*1/2*bs9(1-m62)+1/2*r*gs6*1/2*r*bs6(1-m62)+1/2*r*gs6*bs10(1-m63)+gs10*1/2*(1-r)*bs5(1-m102)+gs10*1/2*bs7(1-m103)+gs10*1/2*bs9(1-m102)+gs10*1/2*r*bs6(1-m102)+gs10*bs10(1-m103);

Factor[br1+br2+br3+br4+br5+br6+br7+br8+br9+br10/.{m11->1/2,m21->1/2,m31->1/2,m41->1/2,m51->1/2,m61->1/2,m71->1/2,m81->1/2,m91->1/2,m101->1/2,m12->1/2,m22->1/2,m32->1/2,m42->1/2,m52->1/2,m62->1/2,m72->1/2,m82->1/2,m92->1/2,m102->1/2,m13->1/2,m23->1/2,m33->1/2,m43->1/2,m53->1/2,m63->1/2,m73->1/2,m83->1/2,m93->1/2,m103->1/2}]

standardize within daughters and sons so that the genotype frequencies = 1

fsum=gr1+gr2+gr3+gr4+gr5+gr6+gr7+gr8+gr9+gr10;

eqng1=gr1/fsum;
eqng2=gr2/fsum;
eqng3=gr3/fsum;
eqng4=gr4/fsum;
eqng5=gr5/fsum;
eqng6=gr6/fsum;
eqng7=gr7/fsum;
eqng8=gr8/fsum;
eqng9=gr9/fsum;
eqng10=gr10/fsum;

Factor[eqng1+eqng2+eqng3+eqng4+eqng5+eqng6+eqng7+eqng8+eqng9+eqng10]

msum=br1+br2+br3+br4+br5+br6+br7+br8+br9+br10;

eqnb1=br1/msum;
eqnb2=br2/msum;
eqnb3=br3/msum;
eqnb4=br4/msum;
eqnb5=br5/msum;
eqnb6=br6/msum;
eqnb7=br7/msum;
eqnb8=br8/msum;
eqnb9=br9/msum;
eqnb10=br10/msum;

Factor[eqnb1+eqnb2+eqnb3+eqnb4+eqnb5+eqnb6+eqnb7+eqnb8+eqnb9+eqnb10]

Mate-SRA

Summary of results

SRA should evolve to produce daughters in response to TT males and sons in response to tt males, as long as the population sex ratio remains near 1/2:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +2 (m42-m12) (1-pave) pave+ (m43-m13) (((1-pave))^2)}+
 (1-pave) pave {(OverBar[z])}{ (m41-m11) pave+ (m42-m12)  (1-2pave)-  (m43-m13) (1-pave)}

Polymorphism requires 2 (1-m) m (s(1-2hf)+z(1-2hm))+(1-2 H) (1-2 m) ΔD > 0 [low dominance coefficients] and the following:

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ<z<((s-hf s)/hm)-((H )/(hm ))θ
	
where θ=((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) )).  If the polymorphism is initially stable in the absence of SRA (ΔD=0), it will continue to be so as more and more extreme SRA evolves as long as the average population sex ratio, (OverBar[D]), stays near 1/2.

Definitions [ENTER]

The following sex ratio strategies are equivalent under mate SRA:

subNoRef={m21->m11,m22->m12,m23->m13,
m31->m11,m32->m12,m33->m13,
m51->m41,m52->m42,m53->m43,
m61->m41,m62->m42,m63->m43,
m71->m41,m72->m42,m73->m43,
m91->m81,m92->m82,m93->m83,
m101->m81,m102->m82,m103->m83};

In the absence of SRA, we would also have:  [Useful for some checks]

subNoSRA={m11->m,m12->m,m13->m};

Frequencies while the modifier is fixed on allele M:

subNoM={g3->1-g1-g2,g4->0,g5->0,g6->0,g7->0,g8->0,g9->0,g10->0,b3->1-b1-b2,b4->0,b5->0,b6->0,b7->0,b8->0,b9->0,b10->0};

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m12->m+H*ΔD,m13->m};

The population sex ratio among the offspring (M fixed):

fracFemale=Factor[bs1 m11+bs2 m12+bs3 m13/.subNoRef/.subNoM]

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficients (Ff, Fm) to simplify solutions:

subpF=Solve[{pf==g1+g2/2,Ff==1-(g2/(2 (1-pf) pf)),pm==b1+b2/2,Fm==1-(b2/(2 (1-pm) pm))},{g1,g2,b1,b2}]//Simplify//Flatten

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= ((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD));

pdif1=((-s+H s+hf s-H z+hm z) ΔD (2 m s-2 hf m s-2 (m^2) s+2 hf (m^2) s-2 hm m z+2 hm (m^2) z+H ΔD-2 H m ΔD) (-2 hf m s+2 hf (m^2) s+2 m z-2 hm m z-2 (m^2) z+2 hm (m^2) z-ΔD+H ΔD+2 m ΔD-2 H m ΔD))/(((2 m s-4 hf m s-2 (m^2) s+4 hf (m^2) s+2 m z-4 hm m z-2 (m^2) z+4 hm (m^2) z-ΔD+2 H ΔD+2 m ΔD-4 H m ΔD))^3);

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under mate-SRA:

eqng1/.subNoRef/.subNoM//Factor

eqng2/.subNoRef/.subNoM//Factor

eqnb1/.subNoRef/.subNoM//Factor

eqnb2/.subNoRef/.subNoM//Factor

recu={%%%%,%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]+1/2recu[[4]]-pm,1-((recu[[4]])/(2 (1-(recu[[3]]+1/2recu[[4]])) (recu[[3]]+1/2recu[[4]])))-Fm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm/.Fm->Fm*sm,{sm,0,1}]]]]/.sm->1

The second and fourth terms will be zero only if Ff and Fm are zero to this order.

solvemeAPPROX/.Ff->0/.Fm->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Finally, we can solve for the average frequency of T:

Solve[%[[3]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

pave0=pave/.%%[[3]]//Factor

which can be written as:

pave0= ((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD));

Factor[%%-pave0]

Equation 1:  The above is rewritten in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript:

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave^2+(m+H ΔD)(2pave(1-pave))+(m)(1-pave)^2),d==1-H},{m}]]

{m+ΔD,m+H ΔD,m}/.%/.(OverBar[D])->(m11*pave^2+m12(2pave(1-pave))+m13(1-pave)^2)/.ΔD->m11-m13/.H->(m12-m13)/(m11-m13)//Factor

Factor[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,0}]]/.sm->1]

Collect[Numerator[%],ΔD,Factor]/Collect[Denominator[%],ΔD,Factor]

(((s(1-hf)-hm z)-H θ)/( (s(1-2 hf)+z(1-2 hm))+(1-2H) θ)) /.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD)),((2 (1-m) m ((1-hm) z-hf s)-(1-H)(1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD==0,z],ΔD,Factor]

{(((1-hf) s)/hm)-((H )/(hm ))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[2 (1-m) m ((1-hm) z-hf s)-(1-H)(1-2 m) ΔD==0,z],ΔD,Factor]

{((hf s)/((1-hm)))-((1-H )/(1-hm ))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that 2 (1 - m) m (s (1 - 2 hf) + z (1 - 2 hm)) - (1 - 2 H) (1 - 
    2 m) ΔD > 0 [low dominance coefficients] and

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ<z<((s-hf s)/hm)-((H )/(hm ))θ

* Negative numerators and denominators, which requires that 2 (1 - m) m (s (1 - 2 hf) + z (1 - 2 hm)) - (1 - 2 H) (1 - 
    2 m) ΔD < 0 [high dominance coefficients] and

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ>z>((s-hf s)/hm)-((H )/(hm ))θ

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqng1/.subNoM,g1],D[eqng1/.subNoM,g2],D[eqng1/.subNoM,b1],D[eqng1/.subNoM,b2]},{D[eqng2/.subNoM,g1],D[eqng2/.subNoM,g2],D[eqng2/.subNoM,b1],D[eqng2/.subNoM,b2]},{D[eqnb1/.subNoM,g1],D[eqnb1/.subNoM,g2],D[eqnb1/.subNoM,b1],D[eqnb1/.subNoM,b2]},{D[eqnb2/.subNoM,g1],D[eqnb2/.subNoM,g2],D[eqnb2/.subNoM,b1],D[eqnb2/.subNoM,b2]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[4]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Solve[%==0,constlam]//Simplify

The above has the sign of:

Factor[%*4 (1-m) m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable (λ < 1) requires that the denominator be positive, which requires that hf and hm be sufficiently small and the ΔD term not very negative.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat2={{D[eqng4,g4],D[eqng4,g5],D[eqng4,g6],D[eqng4,g7],D[eqng4,b4],D[eqng4,b5],D[eqng4,b6],D[eqng4,b7]},
{D[eqng5,g4],D[eqng5,g5],D[eqng5,g6],D[eqng5,g7],D[eqng5,b4],D[eqng5,b5],D[eqng5,b6],D[eqng5,b7]},
{D[eqng6,g4],D[eqng6,g5],D[eqng6,g6],D[eqng6,g7],D[eqng6,b4],D[eqng6,b5],D[eqng6,b6],D[eqng6,b7]},
{D[eqng7,g4],D[eqng7,g5],D[eqng7,g6],D[eqng7,g7],D[eqng7,b4],D[eqng7,b5],D[eqng7,b6],D[eqng7,b7]},
{D[eqnb4,g4],D[eqnb4,g5],D[eqnb4,g6],D[eqnb4,g7],D[eqnb4,b4],D[eqnb4,b5],D[eqnb4,b6],D[eqnb4,b7]},
{D[eqnb5,g4],D[eqnb5,g5],D[eqnb5,g6],D[eqnb5,g7],D[eqnb5,b4],D[eqnb5,b5],D[eqnb5,b6],D[eqnb5,b7]},
{D[eqnb6,g4],D[eqnb6,g5],D[eqnb6,g6],D[eqnb6,g7],D[eqnb6,b4],D[eqnb6,b5],D[eqnb6,b6],D[eqnb6,b7]},
{D[eqnb7,g4],D[eqnb7,g5],D[eqnb7,g6],D[eqnb7,g7],D[eqnb7,b4],D[eqnb7,b5],D[eqnb7,b6],D[eqnb7,b7]}
}/.subNoRef/.subNoM//Simplify;

The calculations are simplified by working in terms of the p's and F's:

mat2WeakSRA=mat2/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-mat2WeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm41,dm42,dm43},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

{((dm43 (1-2 m) (((1-pave))^2))/(4 (1-m) m))+((dm42 (1-2 m) (1-pave) pave)/(2 (1-m) m))+((dm41 (1-2 m) (pave^2))/(4 (1-m) m))};
%-stabcond//Factor

To leading order, (1-2 m) represents the difference between 1/2 and the population sex ratio [times 2]:

Factor[Normal[Series[2(1/2-fracFemale)/.subpF/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,0}]]/.sm->1]

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +2 (m42-m12) (1-pave) pave+ (m43-m13) (((1-pave))^2)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-mat2WeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm41,dm42,dm43,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm42->0/.dm43->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

newstabcond=Collect[stabcond-%%,{dm41,dm42,dm43},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +2 (m42-m12) (1-pave) pave+ (m43-m13) (((1-pave))^2)}+
(1/2) (1-pave) pave {((1-pave)(1-hf)+hf pave) s+(pave(1-hm)+hm (1-pave)) z}{ (m41-m11) pave+ (m42-m12)  (1-2pave)-  (m43-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(pave +(1-pave)(1-hf s))-(pave (1-hf s)+(1-pave)(1-s))]

that is, (OverBar[s]) is the average strength of selection acting against t in females:

(OverBar[s])=Factor[pave(1-(1-hf s))+(1-pave)((1-hf s)-(1-s))]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[(pave (1-hm z)+(1-pave))-(pave (1-z)+(1-pave)(1-hm z))]

that is, (OverBar[z]) is the average strength of selection acting against T in males:

(OverBar[z])=Factor[pave((1-hm z)-(1-z))+(1-pave)(1-(1-hm z))]

With these definitions, the total selection favoring T is (1/2) (OverBar[s])-(1/2)(OverBar[z]), which should be an equivalent means to derive pave (even when the sex ratio isn't near 1/2):

Factor[Normal[Series[(1/2) (OverBar[s])-(1/2)(OverBar[z])/.(OverBar[D])->fracFemale/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

pave/.Solve[%==0,pave]//Simplify

pave0

pave/.Solve[%==0,pave];
%-pave0//Factor

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +2 (m42-m12) (1-pave) pave+ (m43-m13) (((1-pave))^2)}+
 (1-pave) pave {(((OverBar[s])+(OverBar[z]))/2)}{ (m41-m11) pave+ (m42-m12)  (1-2pave)-  (m43-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because (OverBar[s]) and (OverBar[z]) must balance at equilibrium (when the sex ratio is near 1/2 so that θ is negligible), the above is equivalent to the following if we use the definition for pave at the equilibrium:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +2 (m42-m12) (1-pave) pave+ (m43-m13) (((1-pave))^2)}+
 (1-pave) pave {(OverBar[z])}{ (m41-m11) pave+ (m42-m12)  (1-2pave)-  (m43-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pave->pave0+sm*pave1/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form of the result used in the text.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when mated to TT males and if it increases investment in sons when mated to tt males.  Whether it is better to invest in daughters or sons when mated with Tt males depends on the current frequency of T.  If freq(T) is very low, then most matings are tt x tt, but a mother who has just mated with a Tt male is more likely to produce Tt children, and these will be very fit among the daughters (compared to the majority of tt females).

Self-SRA

Summary of results

SRA should evolve to produce daughters within TT females and sons within tt females, as long as the population sex ratio remains near 1/2:

	-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+((1-hm )pave+hm(1-pave))z pave (1-pave){(m41-m11) pave + (m51-m21)(r-pave) +
 (m61-m21) ((1-pave)-r)- (m71-m31 ) (1-pave) }

The recombination term drops out of the above if we assume no cis/trans effects, so that TM/tm and Tm/tM are equivalent (i.e., m51 = m61).

Polymorphism requires 2 (1 - m) m (s (1 - 2 hf) + z (1 - 2 hm)) - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and the following:

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ < z < ((s-hf s)/hm)-((H )/(hm ))θ	
		
where θ=((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) )).  If the polymorphism is initially stable in the absence of SRA (ΔD=0), it will continue to be so as more and more extreme SRA evolves as long as the average population sex ratio stays near 1/2.

Definitions [ENTER]

The following sex ratio strategies are equivalent under self-SRA:

subNoRef={m12->m11,m13->m11,
m22->m21,m23->m21,
m32->m31,m33->m31,
m42->m41,m43->m41,
m52->m51,m53->m51,
m62->m61,m63->m61,
m72->m71,m73->m71,
m82->m81,m83->m81,
m92->m91,m93->m91,
m102->m101,m103->m101};

In the absence of SRA, we would also have:  [Useful for some checks]

subNoSRA={m11->m,m21->m,m31->m};

Frequencies while the modifier is fixed on allele M:

subNoM={g3->1-g1-g2,g4->0,g5->0,g6->0,g7->0,g8->0,g9->0,g10->0,b3->1-b1-b2,b4->0,b5->0,b6->0,b7->0,b8->0,b9->0,b10->0};

The population sex ratio among the offspring (M fixed):

fracFemale=Factor[gs1 m11+gs2 m21+gs3 m31/.subNoRef/.subNoM]

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m21->m+H*ΔD,m31->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficients (Ff, Fm) to simplify solutions:

subpF=Solve[{pf==g1+g2/2,Ff==1-(g2/(2 (1-pf) pf)),pm==b1+b2/2,Fm==1-(b2/(2 (1-pm) pm))},{g1,g2,b1,b2}]//Simplify//Flatten

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= ((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD));

pdif1=((-s+H s+hf s-H z+hm z) ΔD (2 m s-2 hf m s-2 (m^2) s+2 hf (m^2) s-2 hm m z+2 hm (m^2) z+H ΔD-2 H m ΔD) (-2 hf m s+2 hf (m^2) s+2 m z-2 hm m z-2 (m^2) z+2 hm (m^2) z-ΔD+H ΔD+2 m ΔD-2 H m ΔD))/(((2 m s-4 hf m s-2 (m^2) s+4 hf (m^2) s+2 m z-4 hm m z-2 (m^2) z+4 hm (m^2) z-ΔD+2 H ΔD+2 m ΔD-4 H m ΔD))^3);

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under self-SRA:

eqng1/.subNoRef/.subNoM//Factor

eqng2/.subNoRef/.subNoM//Factor

eqnb1/.subNoRef/.subNoM//Factor

eqnb2/.subNoRef/.subNoM//Factor

recu={%%%%,%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]+1/2recu[[4]]-pm,1-((recu[[4]])/(2 (1-(recu[[3]]+1/2recu[[4]])) (recu[[3]]+1/2recu[[4]])))-Fm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm/.Fm->Fm*sm,{sm,0,1}]]]]/.sm->1

The second and forth terms will be zero only if Ff and Fm are zero to this order.

solvemeAPPROX/.Ff->0/.Fm->0//Factor

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[3]]==0,pave]//Simplify

We can then can solve for the difference in allele frequency between the sexes:

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave0=pave/.%%[[3]]//Factor

which can be written as:

((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD));

Factor[%-pave0]

Note that this is the same equilibrium as under mate-SRA.  As in the mate-SRA section, we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave^2+(m+H ΔD)(2pave(1-pave))+m(1-pave)^2)},{m}]]

{m+ΔD,m+H ΔD,m}/.%/.(OverBar[D])->(m11*pave^2+m12(2pave(1-pave))+m13(1-pave)^2)/.ΔD->m11-m13/.H->(m12-m13)/(m11-m13)//Factor

Factor[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

Collect[Numerator[%],ΔD,Factor]/Collect[Denominator[%],ΔD,Factor]

(((s(1-hf)-hm z)-H θ)/( (s(1-2 hf)+z(1-2 hm))+(1-2H) θ)) /.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{((2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD)),((2 (1-m) m ((1-hm) z-hf s)-(1-H) (1-2 m) ΔD)/(2 (1-m) m (s(1-2hf)+z(1-2hm))-(1-2 H) (1-2 m) ΔD))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[2 (1-m) m (s(1-hf)-hm z)+H (1-2 m) ΔD==0,z],ΔD,Factor]

{(((1-hf) s)/hm)-((H )/(hm ))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[2 (1-m) m ((1-hm) z-hf s)-(1-H) (1-2 m) ΔD==0,z],ΔD,Factor]

{((hf s)/((1-hm)))-((1-H )/(1-hm ))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires 2 (1 - m) m (s (1 - 2 hf) + z (1 - 2 hm)) - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ < z < ((s-hf s)/hm)-((H )/(hm ))θ

* Negative numerators and denominators, which requires that 2 (1 - m) m (s (1 - 2 hf) + z (1 - 2 hm)) - (1 - 2 H) (1 - 2 m) ΔD< 0 [high dominance coefficients]and

	((hf s)/(1-hm))-((1-H )/(1-hm ))θ > z > ((s-hf s)/hm)-((H )/(hm ))θ

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqng1/.subNoM,g1],D[eqng1/.subNoM,g2],D[eqng1/.subNoM,b1],D[eqng1/.subNoM,b2]},{D[eqng2/.subNoM,g1],D[eqng2/.subNoM,g2],D[eqng2/.subNoM,b1],D[eqng2/.subNoM,b2]},{D[eqnb1/.subNoM,g1],D[eqnb1/.subNoM,g2],D[eqnb1/.subNoM,b1],D[eqnb1/.subNoM,b2]},{D[eqnb2/.subNoM,g1],D[eqnb2/.subNoM,g2],D[eqnb2/.subNoM,b1],D[eqnb2/.subNoM,b2]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[4]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Solve[%==0,constlam]//Simplify

The above has the sign of:

Factor[%*4 (1-m) m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable (λ < 1) requires that the denominator be positive, which requires that hf and hm be sufficiently small and the ΔD term not be very negative.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat2={{D[eqng4,g4],D[eqng4,g5],D[eqng4,g6],D[eqng4,g7],D[eqng4,b4],D[eqng4,b5],D[eqng4,b6],D[eqng4,b7]},
{D[eqng5,g4],D[eqng5,g5],D[eqng5,g6],D[eqng5,g7],D[eqng5,b4],D[eqng5,b5],D[eqng5,b6],D[eqng5,b7]},
{D[eqng6,g4],D[eqng6,g5],D[eqng6,g6],D[eqng6,g7],D[eqng6,b4],D[eqng6,b5],D[eqng6,b6],D[eqng6,b7]},
{D[eqng7,g4],D[eqng7,g5],D[eqng7,g6],D[eqng7,g7],D[eqng7,b4],D[eqng7,b5],D[eqng7,b6],D[eqng7,b7]},
{D[eqnb4,g4],D[eqnb4,g5],D[eqnb4,g6],D[eqnb4,g7],D[eqnb4,b4],D[eqnb4,b5],D[eqnb4,b6],D[eqnb4,b7]},
{D[eqnb5,g4],D[eqnb5,g5],D[eqnb5,g6],D[eqnb5,g7],D[eqnb5,b4],D[eqnb5,b5],D[eqnb5,b6],D[eqnb5,b7]},
{D[eqnb6,g4],D[eqnb6,g5],D[eqnb6,g6],D[eqnb6,g7],D[eqnb6,b4],D[eqnb6,b5],D[eqnb6,b6],D[eqnb6,b7]},
{D[eqnb7,g4],D[eqnb7,g5],D[eqnb7,g6],D[eqnb7,g7],D[eqnb7,b4],D[eqnb7,b5],D[eqnb7,b6],D[eqnb7,b7]}
}/.subNoRef/.subNoM//Simplify;

The calculations are simplified by working in terms of the p's and F's:

mat2WeakSRA=mat2/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-mat2WeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm41,dm51,dm61,dm71},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

{((dm71 (1-2 m) (((1-pave))^2))/(4 (1-m) m))+(((dm51+dm61) (1-2 m) (1-pave) pave)/(4 (1-m) m))+((dm41 (1-2 m) (pave^2))/(4 (1-m) m))};
%%-%//Factor

To leading order, (1-2 m) represents the difference between 1/2 and the population sex ratio [times 2]:

Factor[Normal[Series[2(1/2-fracFemale)/.subpF/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,0}]]/.sm->1]

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-mat2WeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm41,dm51,dm61,dm71,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm51->0/.dm61->0/.dm71->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

newstabcond=Collect[stabcond-%%,{dm41,dm51,dm61,dm71},Factor]

newstabcond/.s->0/.z->0//Factor

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+(1/2) (1-pave) pave( (1-hf-pave+2 hf pave) s+(hm+pave-2 hm pave) z){(m41-m11) pave + (m51-m21)(r-pave) +
 (m61-m21) ((1-pave)-r)- (m71-m31 ) (1-pave) }/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(pave +(1-pave)(1-hf s))-(pave (1-hf s)+(1-pave)(1-s))]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[(pave (1-hm z)+(1-pave))-(pave (1-z)+(1-pave)(1-hm z))]

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+ (1-pave) pave {(((OverBar[s])+(OverBar[z]))/2)}{(m41-m11) pave + (m51-m21)(r-pave) +
 (m61-m21) ((1-pave)-r)- (m71-m31 ) (1-pave) }/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because (OverBar[s]) and (OverBar[z]) must balance at equilibrium (when the sex ratio is near 1/2 so that θ is negligible), the above is equivalent to the following if we use the definition for pave at the equilibrium:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+ (1-pave) pave {(OverBar[z])}{(m41-m11) pave + (m51-m21)(r-pave) +
 (m61-m21) ((1-pave)-r)- (m71-m31 ) (1-pave) }/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in the text.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when in TT females and if it increases investment in sons when in tt females.  In the text, we ignore cis/trans effects, so the SRA of TM/tm and Tm/tM are assumed equal (i.e., m51 = m61).  The recombination term in the above then drops out.  If we redefine the SRA strategies according to Table 2, the leading eigenvalue derived above for self-SRA becomes equivalent to the leading eigenvalue for mate-SRA.

If we allowed cis/trans effects, then whether it is better to invest in daughters or sons when in Tt females depends on the current frequency of T and the recombination rate.  If r is small, the m modifier is more likely to spread if TM/tm (x5) mothers produce sons (m linked to t) and Tm/tM (x6) mothers produce daughters (m linked to T).  If r is large, both dm51 and dm61 terms approach the sign of (1/2-pave), then when freq(T) is very low, most matings are tt x tt, but a Tt mother is more likely to produce Tt children, and these will be very fit among the daughters (compared to the majority of tt females).  So heterozygous mothers should evolve to invest in daughters when freq(T)< 1/2 and should invest in sons if freq(T)>1/2.

X-linked trait & modifier

Intro

Traits:
A male trait with alleles T (less viable or attractive in males), t (more viable or attractive in males with dominance h [range 0 to 1])
Female-expressed sex ratio adjustment (SRA) trait with alelles M (wild), m (mutant with different levels of expression depending on display trait genotype). 

The T locus could control any kind of trait that has a sex-specific fitness effect (e.g. hip width, body size, mating rate) and can also include traits that are more tradionally thought of as 'sexually selected' (e.g. horn/antler size, tail length, etc). Anything that causes differential effects on the probability of leaving descendents in the two sexes.

Natural selection (s ) on female display trait genotypes:
(({{  TT  ,   Tt  ,   tt  }, {1, 1-h*s, 1-s}}))
s=0 corresponds to no sexual anatagonism at the display trait locus

Natural and sexual selection on male display trait genotypes  (z ):
(({{  T  ,   t  }, {1-z, 1}}))

[In the text, s is defined as (Subscript[s, f]),while z is defined as (Subscript[s, m]).]

SRA by females when mating with t or T males :
(({{      ,   T  ,   t  }, {x1
x2
x3
x4
x5
x6
x7
x8
x9
x10, m11
m21
m31
m41
m51
m61
m71
m81
m91
m101, m12
m22
m32
m42
m52
m62
m72
m82
m92
m102}}))

This matrix is written in a general form that can be used for both mate-SRA and self-SRA.  Appropriate choices of the (Subscript[m, ij]) allow conversion between this matrix and the (Subscript[D, ij]) used in the text.

Life cycle:
Natural selection, sexual selection, recombination, production of haploid gametes, pairing of gametes to create the new adult genotype frequencies

Female genotypes
x1=TM/TM
x2=TM/tM
x3=tM/tM
x4=TM/Tm
x5=TM/tm
x6=Tm/tM
x7=tM/tm
x8=Tm/Tm
x9=Tm/tm
x10=tm/tm

Male genotypes
y1=TM
y2=tM
y3=Tm
y4=tm

Recursions [ENTER]

Off[General::Spell];
Off[General::spell1];

Natural selection in females

Mean fitness of females:

MFF=x1+x2(1-h s)+x3(1-s)+x4+x5(1-h s)+x6(1-h s)+x7(1-s)+x8+x9(1-h s)+x10(1-s);

Natural selection on different female genotypes:

xs1=x1/MFF;
xs2=x2(1-h s)/MFF;
xs3=x3(1-s)/MFF;
xs4=x4/MFF;
xs5=x5(1-h s)/MFF;
xs6=x6(1-h s)/MFF;
xs7=x7(1-s)/MFF;
xs8=x8/MFF;
xs9=x9(1-h s)/MFF;
xs10=x10(1-s)/MFF;

Factor[xs1+xs2+xs3+xs4+xs5+xs6+xs7+xs8+xs9+xs10]

Natural selection in males

Mean fitness of males:

MFM=y1(1-z)+y2+y3(1-z)+y4;

Natural selection on different male genotypes:

ys1=y1(1-z)/MFM;
ys2=y2/MFM;
ys3=y3(1-z)/MFM;
ys4=y4/MFM;

Factor[ys1+ys2+ys3+ys4]

Mating, recombination and syngamy

genotypes of daughters:

xr1=xs1 ys1 m11+1/2xs2 ys1 m21+1/2xs4 ys1 m41+1/2 (xs5 (1-r)ys1) m51+(1/2(xs6 r)ys1)m61;

xr2=xs3 ys1 m31+xs1 ys2 m12+1/2xs2 ys1 m21+1/2xs7 ys1 m71+1/2xs2 ys2 m22+1/2xs4 ys2 m42+1/2(xs6(1-r)ys1)m61+(1/2(xs5 r)ys1)m51+1/2(xs5(1-r)ys2)m52+(1/2(xs6 r)ys2)m62;

xr3=xs3 ys2 m32+1/2xs2 ys2 m22+1/2xs7 ys2 m72+1/2(xs6(1-r)ys2)m62+(1/2(xs5 r)ys2)m52;

xr4=xs8 ys1 m81+xs1 ys3 m11+1/2xs4 ys1 m41+1/2xs9 ys1 m91+1/2xs2 ys3 m21+1/2xs4 ys3 m41+1/2(xs5(1-r)ys3) m51+(1/2(xs6 r)ys3)m61+1/2(xs6(1-r)ys1)m61+(1/2(xs5 r)ys1)m51;

xr5=xs10 ys1 m101+xs1 ys4 m12+1/2xs7 ys1 m71+1/2xs9 ys1 m91+1/2xs2 ys4 m22+1/2xs4 ys4 m42+1/2(xs5(1-r)ys4)m52+(1/2(xs6 r)ys4)m62+1/2(xs5(1-r)ys1)m51+(1/2(xs6 r)ys1)m61;

xr6=xs8 ys2 m82+xs3 ys3 m31+1/2xs4 ys2 m42+1/2xs9 ys2 m92+1/2xs2 ys3 m21+1/2xs7 ys3 m71+1/2(xs6(1-r)ys2)m62+(1/2(xs5 r)ys2)m52+1/2(xs6(1-r)ys3)m61+(1/2(xs5 r)ys3)m51;

xr7=xs10 ys2 m102+xs3 ys4 m32+1/2xs7 ys2 m72+1/2xs9 ys2 m92+1/2xs2 ys4 m22+1/2xs7 ys4 m72+1/2(xs5(1-r)ys2)m52+(1/2(xs6 r)ys2)m62+1/2(xs6(1-r)ys4)m62+(1/2(xs5 r)ys4)m52;

xr8=xs8 ys3 m81+1/2xs4 ys3 m41+1/2xs9 ys3 m91+1/2(xs6(1-r)ys3)m61+(1/2(xs5 r)ys3)m51;

xr9=xs10 ys3 m101+xs8 ys4 m82+1/2xs7 ys3 m71+1/2xs9 ys3 m91+1/2xs4 ys4 m42+1/2xs9 ys4 m92+1/2(xs5(1-r)ys3)m51+(1/2(xs6 r)ys3)m61+1/2(xs6 (1-r)ys4)m62+(1/2(xs5 r)ys4)m52;

xr10=xs10 ys4 m102+1/2xs7 ys4 m72+1/2xs9 ys4 m92+1/2(xs5(1-r)ys4)m52+(1/2(xs6 r)ys4)m62;

Factor[xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10/.{m11->1/2,m21->1/2,m31->1/2,m41->1/2,m51->1/2,m61->1/2,m71->1/2,m81->1/2,m91->1/2,m101->1/2,m12->1/2,m22->1/2,m32->1/2,m42->1/2,m52->1/2,m62->1/2,m72->1/2,m82->1/2,m92->1/2,m102->1/2}]

genotypes of sons:

yr1=xs1(ys1(1-m11)+ys2(1-m12)+ys3(1-m11)+ys4(1-m12))+1/2 xs2 (ys1(1-m21)+ys2(1-m22)+ys3(1-m21)+ys4(1-m22))+1/2xs4(ys1(1-m41)+ys2(1-m42)+ys3(1-m41)+ys4(1-m42))+1/2xs5(1-r)(ys1(1-m51)+ys2(1-m52)+ys3(1-m51)+ys4(1-m52))+1/2xs6 r(ys1(1-m61)+ys2(1-m62)+ys3(1-m61)+ys4(1-m62));

yr2=1/2xs2(ys1(1-m21)+ys2(1-m22)+ys3(1-m21)+ys4(1-m22))+xs3(ys1(1-m31)+ys2(1-m32)+ys3(1-m31)+ys4(1-m32))+1/2xs6(1-r)(ys1(1-m61)+ys2(1-m62)+ys3(1-m61)+ys4(1-m62))+1/2xs5 r(ys1(1-m51)+ys2(1-m52)+ys3(1-m51)+ys4(1-m52))+1/2xs7(ys1(1-m71)+ys2(1-m72)+ys3(1-m71)+ys4(1-m72));

yr3=1/2xs4(ys1(1-m41)+ys2(1-m42)+ys3(1-m41)+ys4(1-m42))+1/2xs6(1-r)(ys1(1-m61)+ys2(1-m62)+ys3(1-m61)+ys4(1-m62))+1/2xs5 r(ys1(1-m51)+ys2(1-m52)+ys3(1-m51)+ys4(1-m52))+xs8(ys1(1-m81)+ys2(1-m82)+ys3(1-m81)+ys4(1-m82))+1/2xs9(ys1(1-m91)+ys2(1-m92)+ys3(1-m91)+ys4(1-m92));

yr4= 1/2xs5(1-r)(ys1(1-m51)+ys2(1-m52)+ys3(1-m51)+ys4(1-m52))+1/2xs6 r(ys1(1-m61)+ys2(1-m62)+ys3(1-m61)+ys4(1-m62))+1/2xs7(ys1(1-m71)+ys2(1-m72)+ys3(1-m71)+ys4(1-m72))+1/2xs9(ys1(1-m91)+ys2(1-m92)+ys3(1-m91)+ys4(1-m92))+xs10(ys1(1-m101)+ys2(1-m102)+ys3(1-m101)+ys4(1-m102));

Factor[yr1+yr2+yr3+yr4/.{m11->1/2,m21->1/2,m31->1/2,m41->1/2,m51->1/2,m61->1/2,m71->1/2,m81->1/2,m91->1/2,m101->1/2,m12->1/2,m22->1/2,m32->1/2,m42->1/2,m52->1/2,m62->1/2,m72->1/2,m82->1/2,m92->1/2,m102->1/2}]

Factor[xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10+yr1+yr2+yr3+yr4]

standardize within daughters and sons so that the genotype frequencies = 1

fsum=xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10;

eqnx1=xr1/fsum;
eqnx2=xr2/fsum;
eqnx3=xr3/fsum;
eqnx4=xr4/fsum;
eqnx5=xr5/fsum;
eqnx6=xr6/fsum;
eqnx7=xr7/fsum;
eqnx8=xr8/fsum;
eqnx9=xr9/fsum;
eqnx10=xr10/fsum;

Factor[eqnx1+eqnx2+eqnx3+eqnx4+eqnx5+eqnx6+eqnx7+eqnx8+eqnx9+eqnx10]

msum=yr1+yr2+yr3+yr4;

eqny1=yr1/msum;
eqny2=yr2/msum;
eqny3=yr3/msum;
eqny4=yr4/msum;

Factor[eqny1+eqny2+eqny3+eqny4]

Mate-SRA

Summary of results

SRA should evolve to produce daughters in response to TT males and sons in response to tt males, as long as the population sex ratio remains near 1/2:

	-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
((((OverBar[z])-((ΔD )/((OverBar[D]) )))/3))(1-pave) pave ((m41-m11)-(m42-m12))


Note that (z-((ΔD )/((OverBar[D]) ))) must remain positive for the system to remain polymorphic. 

Polymorphism requires  (1-2 h) > 0 [low dominance coefficients] and the following:

	2 h s+((ΔD )/((OverBar[D]) )) < z <  2 (1-h) s+((ΔD )/((OverBar[D]) ))

SRA is expected to drive ΔD to a more and more positive value, until z hits the lower boundary and the polymorphism is lost.

Definitions [ENTER]

The following sex ratio strategies are equivalent under mate-SRA:

subNoRef={m21->m11,m31->m11,m22->m12,m32->m12,m51->m41,m61->m41,m71->m41,m52->m42,m62->m42,m72->m42,m91->m81,m101->m81,m92->m82,m102->m82};

Frequencies while the modifier is fixed on allele M:

subNoM={x3->1-x1-x2,x4->0,x5->0,x6->0,x7->0,x8->0,x9->0,x10->0,y2->1-y1,y3->0,y4->0};

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m12->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in females (Ff) to simplify solutions:

subpF=Solve[{pf==x1+x2/2,Ff==1-(x2/(2 (1-pf) pf)),pm==y1},{x1,x2,y1}]//Simplify//Flatten

fracFemale=ys1*m11+ys2*m12/.subNoM;

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (( m (2 s(1-h)-z)+ ΔD)/(2 m s(1-2 h)));

pdif1=(((2 h m s-m z+ΔD) (m z-ΔD) (-2 m s+2 h m s+m z-ΔD))/(8 (((-1+2 h))^2) (m^3) (s^2)));

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under mate-SRA:

Factor[eqnx1/.subNoM/.subNoRef]

Factor[eqnx2/.subNoM/.subNoRef]

Factor[eqny1/.subNoM/.subNoRef]

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]-pm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm,{sm,0,1}]]]]/.sm->1

The second term will be zero only if Ff is zero to this order (implying it is of order (sm^2)).

solvemeAPPROX/.Ff->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Finally, we can solve for the average frequency of T:

Solve[%[[3]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

pave/.%%[[3]]//Factor

Collect[{Numerator[%],Denominator[%]},ΔD,Factor]

which can be written as:

pave0= (( m (2 s(1-h)-z)+ ΔD)/(2 m s(1-2 h)));
Factor[%%%-pave0]

The above is rewritten in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript:

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave+m(1-pave))},{m}]]

{m+ΔD,m}/.%/.(OverBar[D])->(m11*pave+m12(1-pave))/.ΔD->m11-m12//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

(((s(1-h)-(1/2)z)+(1/2) ((ΔD)/(OverBar[D])))/(s(1-2 h)));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{((m (2 (1-h) s-z)+ΔD)/(2 (1-2 h) m s)),((m (z-2 h s)-ΔD)/(2 (1-2 h) m s))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[m (2 (1-h) s-z)+ΔD==0,z],ΔD,Factor]

{2(1-h) s+((ΔD )/((OverBar[D]) ))}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[m (z-2 h s)-ΔD==0,z],ΔD,Factor]

{2h s+((ΔD )/((OverBar[D]) ))}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that (1-2 h) > 0 [low dominance coefficients] and

	2 h s+((ΔD )/((OverBar[D]) )) < z <  2 (1-h) s+((ΔD )/((OverBar[D]) ))

* Negative numerators and denominators, which requires that (1-2 h) < 0 [high dominance coefficients] and

	2 h s+((ΔD )/((OverBar[D]) )) > z >  2 (1-h) s+((ΔD )/((OverBar[D]) ))

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnx1/.subNoM,x1],D[eqnx1/.subNoM,x2],D[eqnx1/.subNoM,y1]},{D[eqnx2/.subNoM,x1],D[eqnx2/.subNoM,x2],D[eqnx2/.subNoM,y1]},{D[eqny1/.subNoM,x1],D[eqny1/.subNoM,x2],D[eqny1/.subNoM,y1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

For the general model, there is a natural constraint (x3/(1-y1)+x1/(y1)==1, but we aren't bothering to use this in this weak selection/weak SRA approximation):

Factor[charpoly/.lambda->0]

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Flatten[Solve[%==0,constlam]]//Factor

The above has the sign of:

Factor[%*3  m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable (λ < 1) requires that the denominator be positive, which requires that h be sufficiently small.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat={{D[eqnx4,x4],D[eqnx4,x5],D[eqnx4,x6],D[eqnx4,x7],D[eqnx4,y3],D[eqnx4,y4]},{D[eqnx5,x4],D[eqnx5,x5],D[eqnx5,x6],D[eqnx5,x7],D[eqnx5,y3],D[eqnx5,y4]},{D[eqnx6,x4],D[eqnx6,x5],D[eqnx6,x6],D[eqnx6,x7],D[eqnx6,y3],D[eqnx6,y4]},{D[eqnx7,x4],D[eqnx7,x5],D[eqnx7,x6],D[eqnx7,x7],D[eqnx7,y3],D[eqnx7,y4]},{D[eqny3,x4],D[eqny3,x5],D[eqny3,x6],D[eqny3,x7],D[eqny3,y3],D[eqny3,y4]},{D[eqny4,x4],D[eqny4,x5],D[eqny4,x6],D[eqny4,x7],D[eqny4,y3],D[eqny4,y4]}}/.subNoRef/.subNoM//Simplify;

The calculations are simplified by working in terms of the p's and F's:

matWeakSRA=mat/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm42},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

(((dm41 (1-2 m) pave)/(3 (1-m) m))+((dm42(1-2 m) (1-pave))/(3 (1-m) m)));
%%-%//Factor

To leading order, (1-2 m) represents the difference between 1/2 and the population sex ratio [times 2]:

Factor[Normal[Series[2(1/2-fracFemale)/.subpF/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,0}]]/.sm->1]

so that stabcond can be written as:

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm42,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm42->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

The leading eigenvalue can then be written as (to leading order):

newstabcond=Collect[stabcond-%%,{dm41,dm42},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
{(1/3) (4 (1-h-pave+2 h pave) s-z+2 ΔD)} (1-pave) pave ((m41-m11)-(m42-m12))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(pave +(1-pave)(1-h s))-(pave (1-h s)+(1-pave)(1-s))]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[1-(1-z)]

and κ as the difference between the transmission of the T and the t allele in males:

κ=((ΔD)/(OverBar[D]))/.(OverBar[D])->fracFemale/.directional;
Factor[Normal[Series[%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

With these definitions, the total selection favoring T is (2/3) (OverBar[s])+((κ )/3)-(1/3)(OverBar[z]), which should be an equivalent means to derive pave (even when the sex ratio isn't near 1/2):

Factor[Normal[Series[(2/3) (OverBar[s])+((κ )/3)-(1/3)(OverBar[z])/.(OverBar[D])->fracFemale/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42,{sm,0,1}]]/.sm->1]

pave/.Solve[%==0,pave];
%-pave0//Factor

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
(((4 (OverBar[s])+κ-(OverBar[z]) )/3))(1-pave) pave ((m41-m11)-(m42-m12))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because the total selection favoring T is (2/3) (OverBar[s])+((κ )/3)-(((OverBar[z]) )/3) , which must equal zero at equilibrium, we can use the substitution pave -> pave0 to show that the above equals:

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
((((OverBar[z])-((ΔD)/(OverBar[D])) )/3))(1-pave) pave ((m41-m11)-(m42-m12))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pave->pave0/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

We know that for the equilibrium to be valid and stable, we must have (1-2h) > 0 and 

	2 h s<z-((ΔD)/(OverBar[D]))<2 (1-h) s
	
which for (OverBar[D]) ~ 1/2 implies that z-((ΔD)/(OverBar[D])) is positive and bounded by:

	2 h s<z-((ΔD)/(OverBar[D]))<2 (1-h) s

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when mated with T males and if it increases investment in sons when mated with t males.

Self-SRA

Summary of results

SRA should evolve to produce daughters within TT females and sons within tt females, as long as the population sex ratio remains near 1/2:

	-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+
(4/3)(z/2) (1-pave) pave { ( (m41-m11) pave +
 (m51-m21)  (-pave +r )+ (m61-m21)  ((1-pave)-r )-(m71-m31) (1-pave))}

The recombination term drops out of the above if we assume no cis/trans effects, so that TM/tm and Tm/tM are equivalent (i.e., m51 = m61).

Polymorphism requires  2 (1 - m) m s (1 - 2 h) - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and the following:

	2 h s-(1-H) θ <  z  <  2 (1-h) s-H θ

where θ=((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) )).  As long as the average population sex ratio stays near 1/2, SRA should be able to increase indefinitely.

Definitions [ENTER]

The following sex ratio strategies are equivalent under self-SRA:

subNoRef={m12->m11,m22->m21,m32->m31,m42->m41,m52->m51,m62->m61,m72->m71,m82->m81,m92->m91,m102->m101};

Frequencies while the modifier is fixed on allele M:

subNoM={x3->1-x1-x2,x4->0,x5->0,x6->0,x7->0,x8->0,x9->0,x10->0,y2->1-y1,y3->0,y4->0};

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m21->m+H*ΔD,m31->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in females (Ff) to simplify solutions:

subpF=Solve[{pf==x1+x2/2,Ff==1-(x2/(2 (1-pf) pf)),pm==y1},{x1,x2,y1}]//Simplify//Flatten

fracFemale=xs1*m11+xs2*m21+xs3*m31/.subNoM;

We will show that the equilibrium trait frequency is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD));

pdif1=((2 m s-2 h m s-2 (m^2) s+2 h (m^2) s-m z+(m^2) z+H ΔD-2 H m ΔD) (-2 h m s+2 h (m^2) s+m z-(m^2) z-ΔD+H ΔD+2 m ΔD-2 H m ΔD) (m s z-2 h m s z-(m^2) s z+2 h (m^2) s z-s ΔD+h s ΔD+H s ΔD+m z ΔD-2 H m z ΔD))/(((2 m s-4 h m s-2 (m^2) s+4 h (m^2) s-ΔD+2 H ΔD+2 m ΔD-4 H m ΔD))^3);

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under self-SRA:

Factor[eqnx1/.subNoM/.subNoRef]

Factor[eqnx2/.subNoM/.subNoRef]

Factor[eqny1/.subNoM/.subNoRef]

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]-pm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm,{sm,0,1}]]]]/.sm->1

The second term will be zero only if Ff is zero to this order.

solvemeAPPROX/.Ff->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[3]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave/.%%[[3]]//Factor

Collect[{Numerator[%],Denominator[%]},ΔD,Factor]

which can be written as:

pave0= (((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD));

Factor[%%%-pave0]

Note that this is the same equilibrium as under mate-SRA.  As in the mate-SRA section, we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave^2+(m+H ΔD)(2pave(1-pave))+m(1-pave)^2)},{m}]]

{m+ΔD,m+H ΔD,m}/.%/.(OverBar[D])->(m11*pave^2+m12(2pave(1-pave))+m13(1-pave)^2)/.ΔD->m11-m13/.H->(m12-m13)/(m11-m13)//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

((((1-h) s-1/2z)-H ΔD((( (OverBar[D])-1/2))/((1-(OverBar[D])) (OverBar[D]))))/((1-2 h) s+(1-2 H) ΔD ((( (OverBar[D])-1/2))/((1-(OverBar[D])) (OverBar[D])))));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{(((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD)),(((1-m) m (z-2 h s)-(1-H) (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[(1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD==0,z],ΔD,Factor]

{2 (1-h) s-H θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[(1-m) m (z-2 h s)-(1-H) (1-2 m) ΔD==0,z],ΔD,Factor]

{2 h s-(1-H) θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that 2 (1 - m) m s (1 - 2 h) - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and

	2 h s-(1-H) θ <  z  <  2 (1-h) s-H θ

* Negative numerators and denominators, which requires that 2 (1 - m) m s (1 - 2 h) - (1 - 2 H) (1 - 2 m) ΔD < 0 [high dominance coefficients] and

	2 h s-(1-H) θ >  z  >  2 (1-h) s-H θ

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnx1/.subNoM,x1],D[eqnx1/.subNoM,x2],D[eqnx1/.subNoM,y1]},{D[eqnx2/.subNoM,x1],D[eqnx2/.subNoM,x2],D[eqnx2/.subNoM,y1]},{D[eqny1/.subNoM,x1],D[eqny1/.subNoM,x2],D[eqny1/.subNoM,y1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Flatten[Solve[%==0,constlam]]//Factor

The above has the sign of:

Factor[%*3 (1-m) m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable requires that the denominator be positive, which requires that h be sufficiently small and the ΔD term not be very negative.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat={{D[eqnx4,x4],D[eqnx4,x5],D[eqnx4,x6],D[eqnx4,x7],D[eqnx4,y3],D[eqnx4,y4]},{D[eqnx5,x4],D[eqnx5,x5],D[eqnx5,x6],D[eqnx5,x7],D[eqnx5,y3],D[eqnx5,y4]},{D[eqnx6,x4],D[eqnx6,x5],D[eqnx6,x6],D[eqnx6,x7],D[eqnx6,y3],D[eqnx6,y4]},{D[eqnx7,x4],D[eqnx7,x5],D[eqnx7,x6],D[eqnx7,x7],D[eqnx7,y3],D[eqnx7,y4]},{D[eqny3,x4],D[eqny3,x5],D[eqny3,x6],D[eqny3,x7],D[eqny3,y3],D[eqny3,y4]},{D[eqny4,x4],D[eqny4,x5],D[eqny4,x6],D[eqny4,x7],D[eqny4,y3],D[eqny4,y4]}}/.subNoRef/.subNoM//Simplify;

The calculations are simplified by working in terms of the p's and F's:

matWeakSRA=mat/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm51,dm61,dm71},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

(((dm71 (1-2 m) (((1-pave))^2))/(3 (1-m) m))+(((dm51+dm61) (1-2 m) (1-pave) pave)/(3 (1-m) m))+((dm41 (1-2 m) (pave^2))/(3 (1-m) m)));
%%-%//Factor

so that stabcond can be written as:

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm51,dm61,dm71,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm51->0/.dm61->0/.dm71->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

Leaving us with:

newstabcond=Collect[stabcond-%%,{dm41,dm51,dm61,dm71},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-(4/3)((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+
(4/3)(z/2) (1-pave) pave { ( (m41-m11) pave +
 (m51-m21)  (-pave +r )+ (m61-m21)  ((1-pave)-r )-(m71-m31) (1-pave))}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when in TT females and if it increases investment in sons when in tt females.    In the text, we ignore cis/trans effects, so the SRA of TM/tm and Tm/tM are assumed equal (i.e., m51 = m61).  The recombination term in the above then drops out.  

If we allowed cis/trans effects, then whether it is better to invest in daughters or sons when in Tt females depends on the current frequency of T and the recombination rate.  If r is small, the m modifier is more likely to spread if TM/tm (x5) mothers produce sons (m linked to t) and Tm/tM (x6) mothers produce daughters (m linked to T).  If r is large, both dm51 and dm61 terms approach the sign of (1/2-pave), then when freq(T) is very low, most matings are tt x tt, but a Tt mother is more likely to produce Tt children, and these will be very fit among the daughters (compared to the majority of tt females).  So heterozygous mothers should evolve to invest in daughters when freq(T)< 1/2 and should invest in sons if freq(T)>1/2.

Note that at the polymorphic equilibrium we must have 2/3 (OverBar[s]) equal in strength to 1/3 (OverBar[z]), so that ((2 z)/3) is the same as ((2 (OverBar[s])+(OverBar[z]))/3).

X-linked trait & autosomal modifier

Intro

Traits:
A male trait with alleles T (less viable or attractive in males), t (more viable or attractive in males).
Female-expressed sex ratio adjustment (SRA) trait with alelles M (wild), m (mutant with different levels of expression depending on display trait genotype). 

The T locus could control any kind of trait that has a sex-specific fitness effect (e.g. hip width, body size, mating rate) and can also include traits that are more tradionally thought of as 'sexually selected' (e.g. horn/antler size, tail length, etc). Anything that causes differential effects on the probability of leaving descendents in the two sexes.

Natural selection (s ) on female display trait genotypes:
(({{  TT  ,   Tt  ,   tt  }, {1, 1-h*s, 1-s}}))
s=0 corresponds to no sexual anatagonism at the display trait locus

Natural and sexual selection on male display trait genotypes  (z ):
(({{  T  ,   t  }, {1-z, 1}}))

[In the text, s is defined as (Subscript[s, f]),while z is defined as (Subscript[s, m]).]

SRA by females when mating with t or T males :
(({{      ,   T  ,   t  }, {x1
x2
x3
x4
x5
x6
x7
x8
x9
x10, m11
m21
m31
m41
m51
m61
m71
m81
m91
m101, m12
m22
m32
m42
m52
m62
m72
m82
m92
m102}}))

This matrix is written in a general form that can be used for both mate-SRA and self-SRA.  Appropriate choices of the (Subscript[m, ij]) allow conversion between this matrix and the (Subscript[D, ij]) used in the text.

Life cycle:
Natural selection, sexual selection, recombination, production of haploid gametes, pairing of gametes to create the new adult genotype frequencies

Female genotypes
x1 = TM/TM
x2 = TM/tM
x3 = tM/tM
x4 = TM/Tm
x5 = TM/tm
x6 = Tm/tM
x7 = tM/tm
x8 = Tm/Tm
x9 = Tm/tm
x10 = tm/tm

Male genotypes
y1 = TMM
y2 = tMM
y3 = TMm
y4 = tMm
y5 = Tmm
y6 = tmm

Recursions [ENTER]

Off[General::Spell];
Off[General::spell1];

Natural selection in females

Mean fitness of females:

MFF=x1+x2(1-h s)+x3(1-s)+x4+x5(1-h s)+x6(1-h s)+x7(1-s)+x8+x9(1-h s)+x10(1-s);

Natural selection on different female genotypes:

xs1=x1/MFF;
xs2=x2(1-h s)/MFF;
xs3=x3(1-s)/MFF;
xs4=x4/MFF;
xs5=x5(1-h s)/MFF;
xs6=x6(1-h s)/MFF;
xs7=x7(1-s)/MFF;
xs8=x8/MFF;
xs9=x9(1-h s)/MFF;
xs10=x10(1-s)/MFF;

Factor[xs1+xs2+xs3+xs4+xs5+xs6+xs7+xs8+xs9+xs10]

Natural selection in males

Mean fitness of males:

MFM=y1(1-z)+y2+y3(1-z)+y4+y5(1-z)+y6;

Natural selection on different male genotypes:

ys1=y1(1-z)/MFM;
ys2=y2/MFM;
ys3=y3(1-z)/MFM;
ys4=y4/MFM;
ys5=y5(1-z)/MFM;
ys6=y6/MFM;

Factor[ys1+ys2+ys3+ys4+ys5+ys6]

Mating, recombination and syngamy

genotypes of daughters (recombination is included in the female genotypes, but will be set to 0.5 in the analysis):

xr1=m11 xs1 ys1+1/2 m11 xs1 ys3+1/2 m21 xs2 ys1+1/4 m21 xs2 ys3+1/2 m41 xs4 ys1+1/4 m41 xs4 ys3+ 1/2(m51 xs5(1-r)+m61 xs6 r)ys1+1/4(m51 xs5(1-r)+m61 xs6 r)ys3;

xr2=m12 xs1 ys2+1/2 m12 xs1 ys4+1/2 m22 xs2 ys2+1/4 m22 xs2 ys4+1/2 m21 xs2 ys1+1/4 m21 xs2 ys3+m31 xs3 ys1+1/2 m31 xs3 ys3+1/2 m42 xs4 ys2+1/4 m42 xs4 ys4+1/2(m52 xs5(1-r)+m62 xs6 r)ys2+1/4(m52 xs5(1-r)+m62 xs6 r)ys4+1/2(m61 xs6(1-r)+m51 xs5 r)ys1+1/4(m61 xs6(1-r)+m51 xs5 r)ys3+1/2 m71 xs7 ys1+1/4 m71 xs7 ys3;

xr3=1/2 m22 xs2 ys2 +1/4 m22 xs2 ys4+m32 xs3 ys2+1/2 m32 xs3 ys4+1/2(m62 xs6(1-r)+m52 xs5 r)ys2+1/4(m62 xs6(1-r)+m52 xs5 r)ys4+1/2 m72 xs7 ys2+1/4 m72 xs7 ys4;

xr4= 1/2 m11 xs1 ys3+m11 xs1 ys5+1/4 m21 xs2 ys3+1/2 m21 xs2 ys5+1/4 m41 xs4 ys3+1/2 m41 xs4 ys5+1/2 m41 xs4 ys1+1/4 m41 xs4 ys3+1/4(m51 xs5(1-r)+m61 xs6 r)ys3+1/2(m51 xs5(1-r)+m61 xs6 r)ys5+1/2(m61 xs6(1-r)+m51 xs5 r)ys1+1/4(m61 xs6(1-r)+m51 xs5 r)ys3+m81 xs8 ys1+1/2 m81 xs8 ys3+1/2 m91 xs9 ys1+1/4 m91 xs9 ys3;

xr5=1/2 m12 xs1 ys4+m12 xs1 ys6+1/4 m22 xs2 ys4+1/2 m22 xs2 ys6 +1/4 m42 xs4 ys4 +1/2 m42 xs4 ys6+1/4(m52 xs5(1-r)+m62 xs6 r)ys4+1/2(m52 xs5(1-r)+m62 xs6 r)ys6+1/2(m51 xs5(1-r)+m61 xs6 r)ys1+1/4(m51 xs5(1-r)+m61 xs6 r)ys3+1/2 m71 xs7 ys1+1/4 m71 xs7 ys3+1/2 m91 xs9 ys1+1/4 m91 xs9 ys3+m101 xs10 ys1+1/2 m101 xs10 ys3;

xr6=1/2 m42 xs4 ys2+1/4 m42 xs4 ys4+1/2(m62 xs6(1-r)+m52 xs5 r) ys2+1/4(m62 xs6(1-r)+m52 xs5 r)ys4+m82 xs8 ys2+1/2 m82 xs8 ys4+1/2 m92 xs9 ys2+1/4 m92 xs9 ys4+1/4 m21 xs2 ys3+1/2 m21  xs2 ys5+1/2 m31 xs3 ys3+m31 xs3 ys5+1/4(m61 xs6(1-r)+m51 xs5 r)ys3+1/2(m61 xs6(1-r)+m51 xs5 r)ys5+1/4 m71 xs7 ys3+1/2 m71 xs7 ys5;

xr7=1/4 m22 xs2 ys4+1/2 m22 xs2 ys6+1/2 m32 xs3 ys4+m32 xs3 ys6+1/2(m52 xs5(1-r)+m62 xs6 r)ys2+1/4(m52 xs5(1-r)+m62 xs6 r)ys4+1/4(m62 xs6(1-r)+m52 xs5 r)ys4+1/2(m62 xs6(1-r)+m52 xs5 r)ys6+1/4 m72 xs7 ys4+1/2 m72 xs7 ys6+1/2 m72  xs7 ys2+1/4 m72 xs7 ys4+1/2 m92 xs9 ys2+1/4 m92 xs9 ys4+m102 xs10 ys2+1/2 m102 xs10 ys4;

xr8=1/4 m41 xs4 ys3+1/2 m41  xs4 ys5+1/4(m61 xs6(1-r)+m51 xs5 r)ys3+1/2(m61 xs6(1-r)+m51 xs5 r)ys5+1/2 m81 xs8 ys3+m81 xs8 ys5+1/4 m91 xs9 ys3+1/2 m91 xs9 ys5;

xr9=1/4 m42 xs4 ys4+1/2 m42  xs4 ys6+1/4(m51 xs5(1-r)+m61 xs6 r)ys3+1/2(m51 xs5(1-r)+m61 xs6 r)ys5+1/4(m62 xs6(1-r)+m52 xs5 r)ys4+1/2(m62 xs6(1-r)+m52 xs5 r)ys6+1/4 m71 xs7 ys3+1/2 m71  xs7 ys5+1/2 m82 xs8 ys4+m82 xs8 ys6+1/4 m92 xs9 ys4+1/2 m92 xs9 ys6+1/4 m91 xs9 ys3+1/2 m91  xs9 ys5+1/2 m101 xs10 ys3+m101 xs10 ys5;

xr10=1/4(m52 xs5(1-r)+m62 xs6 r)ys4+1/2(m52 xs5(1-r)+m62 xs6 r)ys6+1/4 m72 xs7 ys4+1/2 m72  xs7 ys6+1/4 m92 xs9 ys4+1/2 m92 xs9 ys6+1/2 m102 xs10 ys4+m102 xs10 ys6;

Factor[xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10/.{m11->m,m21->m,m31->m,m41->m,m51->m,m61->m,m71->m,m81->m,m91->m,m101->m,m12->m,m22->m,m32->m,m42->m,m52->m,m62->m,m72->m,m82->m,m92->m,m102->m}]

genotypes of sons:

yr1=xs1 (1-m11) ys1+xs1 (1-m12) ys2+1/2 xs1 (1-m11) ys3+1/2 xs1 (1-m12) ys4+1/2 xs2 (1-m21) ys1+1/2 xs2 (1-m22) ys2+1/4 xs2 (1-m21) ys3+1/4 xs2 (1-m22) ys4+1/2 xs4 (1-m41) ys1+1/2 xs4 (1-m42) ys2+1/4 xs4 (1-m41) ys3+1/4 xs4 (1-m42) ys4+1/2(xs5 (1-m51) (1-r)+xs6 (1-m61) r)ys1+1/2(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys2+1/4(xs5 (1-m51) (1-r)+xs6 (1-m61) r)ys3+1/4(xs5 (1-m52) (1-r)+xs6 (1-m62) r)ys4;

yr2=1/2 xs2 (1-m21) ys1+1/2 xs2 (1-m22) ys2+1/4 xs2 (1-m21) ys3+1/4 xs2 (1-m22) ys4+xs3 (1-m31) ys1+xs3 (1-m32) ys2+1/2 xs3 (1-m31) ys3+1/2 xs3 (1-m32) ys4+1/2(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys1+1/2(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys2+1/4(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys3+1/4(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys4+1/2 xs7 (1-m71) ys1+1/2 xs7 (1-m72) ys2+1/4 xs7 (1-m71) ys3+1/4 xs7 (1-m72) ys4;

yr3=1/2 xs1 (1-m11) ys3+1/2 xs1 (1-m12) ys4+xs1 (1-m11) ys5+xs1 (1-m12) ys6+1/4 xs2 (1-m21) ys3+1/4 xs2 (1-m22) ys4+1/2 xs2 (1-m21) ys5+1/2 xs2 (1-m22) ys6+1/4 xs4 (1-m41) ys3+1/4 xs4 (1-m42) ys4+1/2 xs4 (1-m41) ys5+1/2 xs4 (1-m42) ys6+1/2 xs4 (1-m41) ys1+1/2 xs4 (1-m42) ys2+1/4 xs4 (1-m41) ys3+1/4 xs4 (1-m42) ys4+1/4(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys3+1/4(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys4+1/2(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys5+1/2(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys6+1/2(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys1+1/2(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys2+1/4(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys3+1/4(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys4+xs8 (1-m81) ys1+xs8 (1-m82) ys2+1/2 xs8 (1-m81) ys3+1/2 xs8 (1-m82) ys4+1/2 xs9 (1-m91) ys1+1/2 xs9 (1-m92) ys2+1/4 xs9 (1-m91) ys3+1/4 xs9 (1-m92) ys4;

yr4=1/4 xs2 (1-m21) ys3+1/4 xs2 (1-m22) ys4+1/2 xs2 (1-m21) ys5+1/2 xs2 (1-m22) ys6+1/2 xs3 (1-m31) ys3+1/2 xs3 (1-m32) ys4+xs3 (1-m31) ys5+xs3 (1-m32) ys6+1/2(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys1+1/2(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys2+1/4(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys3+1/4(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys4+1/4(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys3+1/4(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys4+1/2(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys5+1/2(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys6+1/2 xs7 (1-m71) ys1+1/2 xs7 (1-m72) ys2+1/4 xs7 (1-m71) ys3+1/4 xs7 (1-m72) ys4+1/4 xs7 (1-m71) ys3+1/4 xs7 (1-m72) ys4+1/2 xs7 (1-m71) ys5+1/2 xs7 (1-m72) ys6+1/2 xs9 (1-m91) ys1+1/2 xs9 (1-m92) ys2+1/4 xs9 (1-m91) ys3+1/4 xs9 (1-m92) ys4+xs10 (1-m101) ys1+xs10 (1-m102) ys2+1/2 xs10 (1-m101) ys3+1/2 xs10 (1-m102) ys4;

yr5=1/4 xs4 (1-m41) ys3+1/4 xs4 (1-m42) ys4+1/2 xs4 (1-m41) ys5+1/2 xs4 (1-m42) ys6+1/4(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys3+1/4(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys4+1/2(xs6 (1-m61)(1-r)+xs5 (1-m51) r)ys5+1/2(xs6 (1-m62)(1-r)+xs5 (1-m52) r)ys6+1/2 xs8 (1-m81) ys3+1/2 xs8 (1-m82) ys4+xs8 (1-m81) ys5+xs8 (1-m82) ys6+1/4 xs9 (1-m91) ys3+1/4 xs9 (1-m92) ys4+1/2 xs9 (1-m91) ys5+1/2 xs9 (1-m92) ys6;

yr6=1/4(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys3+1/4(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys4+1/2(xs5 (1-m51 )(1-r)+xs6 (1-m61) r)ys5+1/2(xs5 (1-m52 )(1-r)+xs6 (1-m62) r)ys6+1/4 xs7 (1-m71) ys3+1/4 xs7 (1-m72) ys4+1/2 xs7 (1-m71) ys5+1/2xs7 (1-m72) ys6+1/4 xs9 (1-m91) ys3+1/4 xs9 (1-m92) ys4+1/2 xs9 (1-m91) ys5+1/2 xs9 (1-m92) ys6+1/2 xs10 (1-m101) ys3+1/2 xs10 (1-m102) ys4+xs10 (1-m101) ys5+xs10 (1-m102) ys6;

Factor[yr1+yr2+yr3+yr4+yr5+yr6/.{m11->m,m21->m,m31->m,m41->m,m51->m,m61->m,m71->m,m81->m,m91->m,m101->m,m12->m,m22->m,m32->m,m42->m,m52->m,m62->m,m72->m,m82->m,m92->m,m102->m}]

Factor[xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10+yr1+yr2+yr3+yr4+yr5+yr6]

standardize within daughters and sons so that the genotype frequencies = 1

fsum=xr1+xr2+xr3+xr4+xr5+xr6+xr7+xr8+xr9+xr10;

eqnx1=xr1/fsum;
eqnx2=xr2/fsum;
eqnx3=xr3/fsum;
eqnx4=xr4/fsum;
eqnx5=xr5/fsum;
eqnx6=xr6/fsum;
eqnx7=xr7/fsum;
eqnx8=xr8/fsum;
eqnx9=xr9/fsum;
eqnx10=xr10/fsum;

Factor[eqnx1+eqnx2+eqnx3+eqnx4+eqnx5+eqnx6+eqnx7+eqnx8+eqnx9+eqnx10]

msum=yr1+yr2+yr3+yr4+yr5+yr6;

eqny1=yr1/msum;
eqny2=yr2/msum;
eqny3=yr3/msum;
eqny4=yr4/msum;
eqny5=yr5/msum;
eqny6=yr6/msum;

Factor[eqny1+eqny2+eqny3+eqny4+eqny5+eqny6]

Mate-SRA

Summary of results

The equilibrium and internal stability are the same as with an X-linked modifier (as expected).

SRA should evolve to produce daughters in response to TT males and sons in response to tt males, as long as the population sex ratio remains near 1/2, but only up to a point:

	- (1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))*((m41-m11)  pave+(m42-m12)  (1-pave))+
(1/5)( (1-pave) pave  (z-2 ((ΔD)/(OverBar[D])))((m41-m11)  -(m42-m12) ))

when (z-2 ((ΔD)/(OverBar[D])))becomes zero, then the system will remain polymorphic but SRA will stop evolving (an intermediate ESS).

->  SRA should reach an ESS rather than increase, even while the system remains polymorphic.

Polymorphism requires  (1-2 h) > 0 [low dominance coefficients] and the following:

		2 h s+((ΔD)/(OverBar[D])) < z <  2 (1-h) s+((ΔD)/(OverBar[D]))
		
(SRA is expected to drive ΔD to a more and more positive value, but it will hit the ESS first).

Definitions [ENTER]

The following sex ratio strategies are equivalent under mate-SRA:

subNoRef={m21->m11,m31->m11,m22->m12,m32->m12,m51->m41,m61->m41,m71->m41,m52->m42,m62->m42,m72->m42,m91->m81,m101->m81,m92->m82,m102->m82};

Frequencies while the modifier is fixed on allele M:

subNoM={x3->1-x1-x2,x4->0,x5->0,x6->0,x7->0,x8->0,x9->0,x10->0,y2->1-y1,y3->0,y4->0,y5->0,y6->0};

directional={m11->m+ΔD,m12->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in females (Ff) to simplify solutions:

subpF=Solve[{pf==x1+x2/2,Ff==1-(x2/(2 (1-pf) pf)),pm==y1},{x1,x2,y1}]//Simplify//Flatten

fracFemale=ys1*m11+ys2*m12/.subNoM;

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (( m (2 s(1-h)-z)+ ΔD)/(2 m s(1-2 h)));

pdif1=(((m z-ΔD) (-2 m s+2 h m s+m z-ΔD) (2 h m s-m z+ΔD))/(8 (((-1+2 h))^2) (m^3) (s^2)));

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under mate-SRA:

Factor[eqnx1/.subNoM/.subNoRef]

Factor[eqnx2/.subNoM/.subNoRef]

Factor[eqny1/.subNoM/.subNoRef]

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]-pm}/.subpF]]

Assuming weak selection:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm,{sm,0,1}]]]]/.sm->1

The second term will be zero only if Ff is zero to this order.

solvemeAPPROX/.Ff->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[3]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave/.%%[[3]]//Factor

Collect[{Numerator[%],Denominator[%]},ΔD,Factor]

which can be written as:

pave0= (( m (2 s(1-h)-z)+ ΔD)/(2 m s(1-2 h)));

Factor[%%%-pave0]

Note that we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave+m(1-pave))},{m}]]

{m+ΔD,m}/.%/.(OverBar[D])->(m11*pave+m12(1-pave))/.ΔD->m11-m12//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

(((2 s(1-h)-z)+(( ΔD)/(OverBar[D])))/(2 s(1-2 h)));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{((m (2 (1-h) s-z)+ΔD)/(2 (1-2 h) m s)),((m (z-2 h s)-ΔD)/(2 (1-2 h) m s))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[m (2 (1-h) s-z)+ΔD==0,z],ΔD,Factor]

{2 (1-h) s+((ΔD)/(OverBar[D]))}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[m (z-2 h s)-ΔD==0,z],ΔD,Factor]

{2 h s+((ΔD)/(OverBar[D]))}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that (1-2 h) > 0 [low dominance coefficients] and

	2 h s+((ΔD)/(OverBar[D]))< z <  2 (1-h) s+((ΔD)/(OverBar[D]))

* Negative numerators and denominators, which requires that (1-2 h)  < 0 [high dominance coefficients] and

	2 h s+((ΔD)/(OverBar[D]))> z >  2 (1-h) s+((ΔD)/(OverBar[D]))

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnx1/.subNoM,x1],D[eqnx1/.subNoM,x2],D[eqnx1/.subNoM,y1]},{D[eqnx2/.subNoM,x1],D[eqnx2/.subNoM,x2],D[eqnx2/.subNoM,y1]},{D[eqny1/.subNoM,x1],D[eqny1/.subNoM,x2],D[eqny1/.subNoM,y1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

For the general model, there is a natural constraint (x3/(1-y1)+x1/(y1)==1, but we aren't bothering to use this in this weak selection/weak SRA approximation):

Factor[charpoly/.lambda->0]

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Flatten[Solve[%==0,constlam]]//Factor

The above has the sign of:

Factor[%*3 m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable (λ < 1) requires that the denominator be positive, which requires that h be sufficiently small.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat={{D[eqnx4,x4],D[eqnx4,x5],D[eqnx4,x6],D[eqnx4,x7],D[eqnx4,y3],D[eqnx4,y4]},{D[eqnx5,x4],D[eqnx5,x5],D[eqnx5,x6],D[eqnx5,x7],D[eqnx5,y3],D[eqnx5,y4]},{D[eqnx6,x4],D[eqnx6,x5],D[eqnx6,x6],D[eqnx6,x7],D[eqnx6,y3],D[eqnx6,y4]},{D[eqnx7,x4],D[eqnx7,x5],D[eqnx7,x6],D[eqnx7,x7],D[eqnx7,y3],D[eqnx7,y4]},{D[eqny3,x4],D[eqny3,x5],D[eqny3,x6],D[eqny3,x7],D[eqny3,y3],D[eqny3,y4]},{D[eqny4,x4],D[eqny4,x5],D[eqny4,x6],D[eqny4,x7],D[eqny4,y3],D[eqny4,y4]}}/.subNoRef/.subNoM//Simplify;

Because the modifier and trait are on different loci, r is set to 1/2.

r=1/2;

The calculations are simplified by working in terms of the p's and F's:

matWeakSRA=mat/.subpF//Factor;

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm42},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

(((dm41 (1-2 m) pave)/(4 (1-m) m))+((dm42(1-2 m) (1-pave))/(4 (1-m) m)));
%%-%//Factor

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.m43->m13+sm*dm43/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term, which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm42,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm42->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

newstabcond=Collect[stabcond-%%,{dm41,dm42},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
(1/5) (1-pave) pave {(4 s-4 h s-4 pave s+8 h pave s-z)}{ (m41-m11)- (m42-m12)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(pave +(1-pave)(1-h s))-(pave (1-h s)+(1-pave)(1-s))]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[1-(1-z)]

and κ as the difference between the transmission of the T and the t allele in males:

κ=((ΔD)/(OverBar[D]))/.(OverBar[D])->fracFemale/.directional;
Factor[Normal[Series[%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
(((4 (OverBar[s])-(OverBar[z]) ))/5) (1-pave) pave { (m41-m11)- (m42-m12)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because the total selection favoring T is (2/3) (OverBar[s])+((κ )/3)-(((OverBar[z]) )/3) , which must equal zero at equilibrium, we can use the substitution pave -> pave0 to show that the above equals:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) pave +(m42-m12) (1-pave)}+
((((OverBar[z]) -2((ΔD)/(OverBar[D]))))/5) (1-pave) pave { (m41-m11)- (m42-m12)}/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m42->m12+dm42/.m43->m13+dm43;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pave->pave0/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid and SRA is absent, the modifier will invade if it increases investment in daughters when mated with T males and if it increases investment in sons when mated with t males.  In this case, however, SRA is limited and only evolves to the point that z-2 ((ΔD)/(OverBar[D]))=0 (assuming an equal sex ratio).

We know that for the equilibrium to be valid and stable, we must have (1-2h)> 0 and 

	2 h s<z-((ΔD)/(OverBar[D]))<2 (1-h) s
	
which for (OverBar[D]) ~ 1/2 implies that z - ((ΔD)/(OverBar[D])) is positive and bounded by:

	2 h s<z-((ΔD)/(OverBar[D]))<2 (1-h) s
	
But this doesn't require that z-2((ΔD)/(OverBar[D])) be positive.  The following is a counter-example:

numpar={z->1/100,s->1/100,h->1/10,ΔD->1/300,dm41->0.001,dm42->0};

We must first solve for the equilibrium and the value of m that allows the sex ratio to be ~1/2:

recu/.directional/.numpar//Simplify

fracFemale/.directional/.numpar//Simplify

numx=NSolve[{%%[[1]]==x1,%%[[2]]==x2,%%[[3]]==y1,%==1/2},{m,x1,x2,y1},25]

The solution of interest:

numx=numx[[3]]

{pave0,1-pave0,fracFemale,z-4ΔD}/.directional/.numpar/.numx//N

While we would expect mothers mated to T to produce more daughters (dm41 > 0), such alleles don't invade

Collect[N[- (1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))*((m41-m11)  pave+(m42-m12)  (1-pave))+
(1/5)( (1-pave) pave  (z-2 ((ΔD)/(OverBar[D])))((m41-m11)  -(m42-m12) ))/.(OverBar[D])->fracFemale/.
m41->m11+dm41/.m42->m12+dm42/.pave->pave0/.directional/.numpar/.numx,25],dm41]

The numerically exact leading eigenvalue has the same order of magnitude:

N[Eigenvalues[mat/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.numpar/.numx/.sm->1]-1]

It is still the case at this equilibrium that T is favored in females and t is favored in males by natural selection:

(xs1-x1)+(xs2-x2)/2/.subNoM/.numx/.numpar

(ys1-y1)/.subNoM/.numx/.numpar

Overall the freq of T goes down during natural selection:

2/3*%%+1/3*%

Here we look at the change in allele frequency during reproduction:

(eqnx1-x1)+(eqnx2-x2)/2/.subNoRef/.directional/.subNoM/.numx/.{z->0,s->0,h->0}/.numpar

(eqny1-y1)/.subNoRef/.directional/.subNoM/.numx/.{z->0,s->0,h->0}/.numpar

Overall the freq of T goes up up during reproduction because males bearing T alleles are more likely to produce daughters:

2/3*%%+1/3*%

The eigenvector associated with the leading eigenvalue tells us the direction in which the new m modifier evolves:

x4 = TM/Tm
x5 = TM/tm
x6 = Tm/tM
x7 = tM/tm
y3 = TMm
y4 = tMm

N[Eigenvectors[mat/.m41->m11+sm*dm41/.m42->m12+sm*dm42/.directional/.numpar/.numx/.sm->1]][[1]]

The above points in a direction associated with an excess of females:

eigenSR=%.{1,1,1,1,0,0}/%.{1,1,1,1,1,1}

Calculating the genotypic frequencies among daughters and sons along this eigenvector:

{{%%[[1]],%%[[2]]+%%[[3]],%%[[4]]}/Drop[%%,-2].{1,1,1,1},Drop[%%,4]/Drop[%%,4].{1,1}}

{{x1,x2,1-x1-x2},{y1,1-y1}}/.numx

The modifier also becomes associated with TT daughters and T sons.

%%-%

The above associations with TT daughters and T sons were repeated in other similar cases (e.g., dm41 -> 0, dm42 -> -0.001), but the effects of the modifier on the sex ratio depended on the action of the modifier.

With more intermediate values of h, however, the polymorphism tended to be lost before the ESS was reached:

numpar={z->1/100,s->1/100,h->1/3,ΔD->1/500};

We must first solve for the equilibrium and the value of m that allows the sex ratio to be ~1/2:

recu/.directional/.numpar//Simplify

fracFemale/.directional/.numpar//Simplify

numx=NSolve[{%%[[1]]==x1,%%[[2]]==x2,%%[[3]]==y1,%==1/2},{m,x1,x2,y1},25]

None of the above solutions allow a polymorphic equilibrium.

Self-SRA

Summary of results

The equilibrium and internal stability are the same as with an X-linked modifier (as expected).

SRA should evolve to produce daughters within TT females and sons within tt females, as long as the population sex ratio remains near 1/2:

	- (1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))*((m41-m11) (pave^2)+((m51-m21)+(m61-m21))pave(1-pave)+(m71-m31) (((1-pave))^2))+(3/5)(1-pave) pave z((m41-m11) pave+(((m51-m21)+(m61-m21))/2) (1-2 pave)-(m71-m31) (1-pave))
	

Polymorphism requires  2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD > 0 [low dominance coefficients] and the following:

	2 h s-H θ <  z  <  2 (1-h) s-H θ,
	
where θ=((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) )).  As long as the average population sex ratio stays near 1/2, SRA should be able to increase indefinitely.

Definitions [ENTER]

The following sex ratio strategies are equivalent under self-SRA:

subNoRef={m12->m11,m22->m21,m32->m31,m42->m41,m52->m51,m62->m61,m72->m71,m82->m81,m92->m91,m102->m101};

Frequencies while the modifier is fixed on allele M:

subNoM={x3->1-x1-x2,x4->0,x5->0,x6->0,x7->0,x8->0,x9->0,x10->0,y2->1-y1,y3->0,y4->0,y5->0,y6->0};

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m21->m+H*ΔD,m31->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in females (Ff) to simplify solutions:

subpF=Solve[{pf==x1+x2/2,Ff==1-(x2/(2 (1-pf) pf)),pm==y1},{x1,x2,y1}]//Simplify//Flatten

fracFemale=xs1*m11+xs2*m21+xs3*m31/.subNoM;

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD));

pdif1=((2 m s-2 h m s-2 (m^2) s+2 h (m^2) s-m z+(m^2) z+H ΔD-2 H m ΔD) (-2 h m s+2 h (m^2) s+m z-(m^2) z-ΔD+H ΔD+2 m ΔD-2 H m ΔD) (m s z-2 h m s z-(m^2) s z+2 h (m^2) s z-s ΔD+h s ΔD+H s ΔD+m z ΔD-2 H m z ΔD))/(((2 m s-4 h m s-2 (m^2) s+4 h (m^2) s-ΔD+2 H ΔD+2 m ΔD-4 H m ΔD))^3);

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under self-SRA:

Factor[eqnx1/.subNoM/.subNoRef]

Factor[eqnx2/.subNoM/.subNoRef]

Factor[eqny1/.subNoM/.subNoRef]

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]+1/2recu[[2]]-pf,1-((recu[[2]])/(2 (1-(recu[[1]]+1/2recu[[2]])) (recu[[1]]+1/2recu[[2]])))-Ff,recu[[3]]-pm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Ff->Ff*sm,{sm,0,1}]]]]/.sm->1

The second term will be zero only if Ff is zero to this order.

solvemeAPPROX/.Ff->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[1]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[3]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave/.%%[[3]]//Factor

Collect[{Numerator[%],Denominator[%]},ΔD,Factor]

which can be written as:

pave0= (((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD));

Factor[%%%-pave0]

Note that this is the same equilibrium as under mate-SRA.  As in the mate-SRA section, we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave^2+(m+H ΔD)(2pave(1-pave))+m(1-pave)^2)},{m}]]

{m+ΔD,m+H ΔD,m}/.%/.(OverBar[D])->(m11*pave^2+m12(2pave(1-pave))+m13(1-pave)^2)/.ΔD->m11-m13/.H->(m12-m13)/(m11-m13)//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

((((1-h) s-1/2z)-H θ)/((1-2 h) s+(1-2 H)θ))/. θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{(((1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD)),(((1-m) m (z-2 h s)-(1-H) (1-2 m) ΔD)/(2 (1-m) m s(1-2 h)-(1-2 H) (1-2 m) ΔD))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[(1-m) m (2 s(1-h)-z)+H (1-2 m) ΔD==0,z],ΔD,Factor]

{2(1-h) s-H θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[(1-m) m (z-2 h s)-(1-H) (1-2 m) ΔD==0,z],ΔD,Factor]

{2h s-H θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that 2 (1 - m) m s (1 - 2 h) - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and

	2 h s-H θ <  z  <  2 (1-h) s-H θ

* Negative numerators and denominators, which requires that 2 (1 - m) m s (1 - 2 h) - (1 - 2 H) (1 - 2 m) ΔD < 0 [high dominance coefficients] and

	2 h s-H θ >  z  >  2 (1-h) s-H θ

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnx1/.subNoM,x1],D[eqnx1/.subNoM,x2],D[eqnx1/.subNoM,y1]},{D[eqnx2/.subNoM,x1],D[eqnx2/.subNoM,x2],D[eqnx2/.subNoM,y1]},{D[eqny1/.subNoM,x1],D[eqny1/.subNoM,x2],D[eqny1/.subNoM,y1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Flatten[Solve[%==0,constlam]]//Factor

The above has the sign of:

Factor[%*3 (1-m) m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable (λ < 1) requires that the denominator be positive, which requires that hf and hm be sufficiently small and the ΔD term not be very negative.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat={{D[eqnx4,x4],D[eqnx4,x5],D[eqnx4,x6],D[eqnx4,x7],D[eqnx4,y3],D[eqnx4,y4],D[eqnx4,y5],D[eqnx4,y6]},{D[eqnx5,x4],D[eqnx5,x5],D[eqnx5,x6],D[eqnx5,x7],D[eqnx5,y3],D[eqnx5,y4],D[eqnx5,y5],D[eqnx5,y6]},{D[eqnx6,x4],D[eqnx6,x5],D[eqnx6,x6],D[eqnx6,x7],D[eqnx6,y3],D[eqnx6,y4],D[eqnx6,y5],D[eqnx6,y6]},{D[eqnx7,x4],D[eqnx7,x5],D[eqnx7,x6],D[eqnx7,x7],D[eqnx7,y3],D[eqnx7,y4],D[eqnx7,y5],D[eqnx7,y6]},{D[eqny3,x4],D[eqny3,x5],D[eqny3,x6],D[eqny3,x7],D[eqny3,y3],D[eqny3,y4],D[eqny3,y5],D[eqny3,y6]},{D[eqny4,x4],D[eqny4,x5],D[eqny4,x6],D[eqny4,x7],D[eqny4,y3],D[eqny4,y4],D[eqny4,y5],D[eqny4,y6]},
{D[eqny5,x4],D[eqny5,x5],D[eqny5,x6],D[eqny5,x7],D[eqny5,y3],D[eqny5,y4],D[eqny5,y5],D[eqny5,y6]},{D[eqny6,x4],D[eqny6,x5],D[eqny6,x6],D[eqny6,x7],D[eqny6,y3],D[eqny6,y4],D[eqny6,y5],D[eqny6,y6]}}/.subNoRef/.subNoM//Simplify;

r=1/2;

The calculations are simplified by working in terms of the p's and F's:

matWeakSRA=mat/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm51,dm61,dm71},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

(((dm71 (1-2 m) (((1-pave))^2))/(4 (1-m) m))+(((dm51+dm61) (1-2 m) (1-pave) pave)/(4 (1-m) m))+((dm41 (1-2 m) (pave^2))/(4 (1-m) m)));
%%-%//Factor

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[8]*lambda-matWeakSRA]/.m41->m11+sm*dm41/.m51->m21+sm*dm51/.m61->m21+sm*dm61/.m71->m31+sm*dm71/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm41,dm51,dm61,dm71},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm41->0/.dm51->0/.dm61->0/.dm71->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

newstabcond=Collect[stabcond-%%,{dm41,dm51,dm61,dm71},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+
((2(s-h s-pave s+2 h pave s+z))/5) (1-pave) pave ((m41-m11) pave+(((m51-m21)+(m61-m21))/2) (1-2 pave)-(m71-m31) (1-pave))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(pave +(1-pave)(1-h s))-(pave (1-h s)+(1-pave)(1-s))]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[(1)-(1-z)]

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+
((2((OverBar[s])+(OverBar[z])))/5) (1-pave) pave ((m41-m11) pave+(((m51-m21)+(m61-m21))/2) (1-2 pave)-(m71-m31) (1-pave))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because the total selection favoring T is (2/3) (OverBar[s])-(((OverBar[z]) )/3), which must equal zero at equilibrium (assuming a nearly equal sex ratio), we can use the substitution pave -> pave0 to show that the above equals:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m41-m11) (pave^2) +((((m51-m21))/2)+(((m61-m21))/2)) 2(1-pave) pave+ (m71-m31 )(((1-pave))^2)}+
((3(OverBar[z]))/5) (1-pave) pave ((m41-m11) pave+(((m51-m21)+(m61-m21))/2) (1-2 pave)-(m71-m31) (1-pave))/.(OverBar[D])->fracFemale/.m41->m11+dm41/.m51->m21+dm51/.m61->m21+dm61/.m71->m31+dm71;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pave->pave0+sm*pave1/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm51->sm*dm51/.dm61->sm*dm61/.dm71->sm*dm71/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of - ((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when in TT females and if it increases investment in sons when in tt females.  Whether it is better to invest in daughters or sons when in Tt females depends on the current frequency of T.  When freq(T) is very low, most matings are tt x t, but a Tt mother is more likely to produce Tt children, and these will be very fit among the daughters (compared to the majority of tt females).  So heterozygous mothers should evolve to invest in daughters when freq(T)< 1/2 and should invest in sons if freq(T)>1/2.

Z-linked trait & modifier

Because sex ratio control is assumed to be governed by the maternal Z, sex ratios are strongly driven toward producing all sons, which inherit their mother's Z chromosome (Hamilton 1967).   SRA has no effect to leading order;  SRA modification is proportional to the strength of selection on the trait, but sex chromosome drive is not and is stronger.

Z-linked trait & autosomal modifier

Intro

Natural selection (s ) on female display trait genotypes:
(({{  T  ,   t  }, {1, 1-s}}))

Natural and sexual selection on male display trait genotypes  (z ):
(({{  TT  ,   Tt  ,   tt  }, {1-z, 1-h*z, 1}}))

[In the text, s is defined as (Subscript[s, f]),while z is defined as (Subscript[s, m]).]

SRA by females when mating with t or T males :
(({{      ,   TT  ,  T t  }, {w1
w2
w3
w4
w5
w6, m11
m21
m31
m41
m51
m61, m12
m22
m32
m42
m52
m62}})({{  tt  }, {m13
m23
m33
m43
m53
m63}}))

This matrix is written in a general form that can be used for both mate-SRA and self-SRA.  Appropriate choices of the (Subscript[m, ij]) allow conversion between this matrix and the (Subscript[D, ij]) used in the text.

Life cycle:
Natural selection, sexual selection, recombination, production of haploid gametes, pairing of gametes to create the new adult genotype frequencies

Female genotypes
w1=TMM
w2=tMM
w3=TMm
w4=tMm
w5=Tmm
w6=tmm

Male genotypes
z1=TM/TM
z2=TM/tM
z3=tM/tM
z4=TM/Tm
z5=TM/tm
z6=Tm/tM
z7=tM/tm
z8=Tm/Tm
z9=Tm/tm
z10=tm/tm

Recursions [ENTER]

Off[General::spell]
Off[General::spell1]

Natural selection in females

Mean fitness of females:

MFF=w1+w2(1-s)+w3+w4(1-s)+w5+w6(1-s);

Natural selection on different female genotypes:

ws1=w1/MFF;
ws2=w2(1-s)/MFF;
ws3=w3/MFF;
ws4=w4(1-s)/MFF;
ws5=w5/MFF;
ws6=w6(1-s)/MFF;

Factor[ws1+ws2+ws3+ws4+ws5+ws6]

Natural selection in males

Mean fitness of males:

MMF=z1(1-z)+z2(1-h z)+z3+z4(1-z)+z5(1-h z)+z6(1-h z)+z7+z8(1-z)+z9(1-h z)+z10;

Natural selection on different male genotypes:

zs1=z1(1-z)/MMF;
zs2=z2(1-h z)/MMF;
zs3=z3/MMF;
zs4=z4(1-z)/MMF;
zs5=z5(1-h z)/MMF;
zs6=z6(1-h z)/MMF;
zs7=z7/MMF;
zs8=z8(1-z)/MMF;
zs9=z9(1-h z)/MMF;
zs10=z10/MMF;

Factor[zs1+zs2+zs3+zs4+zs5+zs6+zs7+zs8+zs9+zs10]

Mating, recombination and syngamy

genotypes of daughters (recombination is included in the female genotypes, but will be set to 0.5 in the analysis):

wr1=zs1 m11 ws1+zs1 m21 ws2+1/2 zs1 m31 ws3+1/2 zs1 m41 ws4+1/2 zs2 m12 ws1+1/2 zs2 m22 ws2+1/4 zs2 m32 ws3+1/4 zs2 m42 ws4+1/2 zs4 m11 ws1+1/2 zs4 m21 ws2+1/4 zs4 m31 ws3+1/4 zs4 m41 ws4+1/2(zs5 m12(1-r)+zs6 m12 r)ws1+1/2(zs5 m22(1-r)+zs6 m22 r)ws2+1/4(zs5 m32(1-r)+zs6 m32 r)ws3+1/4(zs5 m42(1-r)+zs6 m42 r)ws4;

wr2=1/2 zs2 m12 ws1+1/2 zs2 m22 ws2+1/4 zs2 m32 ws3+1/4 zs2 m42 ws4+zs3 m13 ws1+zs3 m23 ws2+1/2 zs3 m33 ws3+1/2 zs3 m43 ws4+1/2(zs6 m12(1-r)+zs5 m12 r)ws1+1/2(zs6 m22(1-r)+zs5 m22 r)ws2+1/4(zs6 m32(1-r)+zs5 m32 r)ws3+1/4(zs6 m42(1-r)+zs5 m42 r)ws4+1/2 zs7 m13 ws1+1/2 zs7 m23 ws2+1/4 zs7 m33 ws3+1/4 zs7 m43 ws4;

wr3=1/2 zs1 m31 ws3+1/2 zs1 m41 ws4+zs1 m51 ws5+zs1 m61 ws6+1/4 zs2 m32 ws3+1/4 zs2 m42 ws4+1/2 zs2 m52 ws5+1/2 zs2 m62 ws6+1/4 zs4 m31 ws3+1/4 zs4 m41 ws4+1/2 zs4 m51 ws5+1/2 zs4 m61 ws6+1/2 zs4 m11 ws1+1/2 zs4 m21 ws2+1/4 zs4 m31 ws3+1/4 zs4 m41 ws4+1/4(zs5 m32(1-r)+zs6 m32 r)ws3+1/4(zs5 m42(1-r)+zs6 m42 r)ws4+1/2(zs5 m52(1-r)+zs6 m52 r)ws5+1/2(zs5 m62(1-r)+zs6 m62 r)ws6+1/2(zs6 m12(1-r)+zs5 m12 r)ws1+1/2(zs6 m22(1-r)+zs5 m22 r)ws2+1/4(zs6 m32(1-r)+zs5 m32 r)ws3+1/4(zs6 m42(1-r)+zs5 m42 r)ws4+zs8 m11 ws1+zs8 m21 ws2+1/2 zs8 m31 ws3+1/2 zs8 m41 ws4+1/2 zs9 m12 ws1+1/2 zs9 m22 ws2+1/4 zs9 m32 ws3+1/4 zs9 m42 ws4;

wr4=1/4 zs2 m32 ws3+1/4 zs2 m42 ws4+1/2 zs2 m52 ws5+1/2 zs2 m62 ws6+1/2 zs3 m33 ws3+1/2 zs3 m43 ws4+zs3 m53 ws5+zs3 m63 ws6+1/2(zs5 m12(1-r)+zs6 m12 r)ws1+1/2(zs5 m22(1-r)+zs6 m22 r)ws2+1/4(zs5 m32(1-r)+zs6 m32 r)ws3+1/4(zs5 m42(1-r)+zs6 m42 r)ws4+1/4(zs6 m32(1-r)+zs5 m32 r)ws3+1/4(zs6 m42(1-r)+zs5 m42 r)ws4+1/2(zs6 m52(1-r)+zs5 m52 r)ws5+1/2(zs6 m62(1-r)+zs5 m62 r)ws6+1/2 zs7 m13 ws1+1/2 zs7 m23 ws2+1/4 zs7 m33 ws3+1/4 zs7 m43 ws4+1/4 zs7 m33 ws3+1/4 zs7 m43 ws4+1/2 zs7 m53 ws5+1/2 zs7 m63 ws6+1/2 zs9 m12 ws1+1/2 zs9 m22 ws2+1/4 zs9 m32 ws3+1/4 zs9 m42 ws4+zs10 m13 ws1+zs10 m23 ws2+1/2 zs10 m33 ws3+1/2 zs10 m43 ws4;

wr5=1/4 zs4 m31 ws3+1/4 zs4 m41 ws4+1/2 zs4 m51 ws5+1/2 zs4 m61 ws6+1/4(zs6 m32(1-r)+zs5 m32 r)ws3+1/4(zs6 m42(1-r)+zs5 m42 r)ws4+1/2(zs6 m52(1-r)+zs5 m52 r)ws5+1/2(zs6 m62(1-r)+zs5 m62 r)ws6+1/2 zs8 m31 ws3+1/2 zs8 m41 ws4+zs8 m51 ws5+zs8 m61 ws6+1/4 zs9 m32 ws3+1/4 zs9 m42 ws4+1/2 zs9 m52 ws5+1/2 zs9 m62 ws6;

wr6=1/4(zs5 m32(1-r)+zs6 m32 r)ws3+1/4(zs5 m42(1-r)+zs6 m42 r)ws4+1/2(zs5 m52(1-r)+zs6 m52 r)ws5+1/2(zs5 m62(1-r)+zs6 m62 r)ws6+1/4 zs7 m33 ws3+1/4 zs7 m43 ws4+1/2 zs7 m53 ws5+1/2zs7 m63 ws6+1/4 zs9 m32 ws3+1/4 zs9 m42 ws4+1/2 zs9 m52 ws5+1/2 zs9 m62 ws6+1/2 zs10 m33 ws3+1/2 zs10 m43 ws4+zs10 m53 ws5+zs10 m63 ws6;

Factor[wr1+wr2+wr3+wr4+wr5+wr6/.{m11->m,m12->m,m13->m,m21->m,m22->m,m23->m,m31->m,m32->m,m33->m,m41->m,m42->m,m43->m,m51->m,m52->m,m53->m,m61->m,m62->m,m63->m}]

genotypes of sons:

zr1=zs1 (1-m11)ws1+1/2 zs1 (1-m31)ws3+1/2 zs2 (1-m12)ws1+1/4 zs2 (1-m32)ws3+1/2 zs4 (1-m11)ws1+1/4 zs4 (1-m31)ws3+ 1/2(zs5 (1-m12) (1-r)+zs6 (1-m12) r)ws1+1/4(zs5 (1-m32) (1-r)+zs6 (1-m32) r)ws3;

zr2=zs1 (1-m21)ws2+1/2 zs1 (1-m41)ws4+1/2 zs2 (1-m22)ws2+1/4 zs2 (1-m42)ws4+1/2 zs2 (1-m12)ws1+1/4 zs2 (1-m32)ws3+zs3 (1-m13)ws1+1/2 zs3 (1-m33)ws3+1/2 zs4 (1-m21)ws2+1/4 zs4 (1-m41)ws4+1/2(zs5 (1-m22) (1-r)+zs6 (1-m22) r)ws2+1/4(zs5 (1-m42) (1-r)+zs6 (1-m42) r)ws4+1/2(zs6 (1-m12) (1-r)+zs5 (1-m12) r)ws1+1/4(zs6 (1-m32) (1-r)+zs5 (1-m32) r)ws3 +1/2 zs7 (1-m13)ws1+1/4 zs7 (1-m33)ws3;

zr3=1/2 zs2 (1-m22)ws2 +1/4 zs2 (1-m42)ws4+zs3 (1-m23)ws2+1/2 zs3 (1-m43)ws4+1/2(zs6 (1-m22) (1-r)+zs5 (1-m22) r)ws2+1/4(zs6 (1-m42) (1-r)+zs5 (1-m42) r)ws4+1/2 zs7 (1-m23)ws2+1/4 zs7 (1-m43)ws4;

zr4= 1/2 zs1 (1-m31)ws3+zs1 (1-m51)ws5+1/4 zs2 (1-m32)ws3+1/2 zs2 (1-m52)ws5+1/4 zs4 (1-m31)ws3+1/2 zs4 (1-m51)ws5+1/2 zs4 (1-m11)ws1+1/4 zs4 (1-m31)ws3+1/4(zs5 (1-m32) (1-r)+zs6 (1-m32) r)ws3+1/2(zs5 (1-m52) (1-r)+zs6 (1-m52) r)ws5+1/2(zs6 (1-m12) (1-r)+zs5 (1-m12) r)ws1+1/4(zs6 (1-m32) (1-r)+zs5 (1-m32) r)ws3+zs8 (1-m11)ws1+1/2 zs8 (1-m31)ws3+1/2 zs9 (1-m12)ws1+1/4 zs9 (1-m32)ws3;

zr5=1/2 zs1 (1-m41)ws4+zs1 (1-m61)ws6+1/4 zs2 (1-m42)ws4+1/2 zs2 (1-m62)ws6 +1/4 zs4 (1-m41)ws4 +1/2 zs4 (1-m61)ws6+1/4(zs5 (1-m42) (1-r)+zs6 (1-m42) r)ws4+1/2(zs5 (1-m62) (1-r)+zs6 (1-m62) r)ws6+1/2(zs5 (1-m12) (1-r)+zs6 (1-m12) r)ws1+1/4(zs5 (1-m32) (1-r)+zs6 (1-m32) r)ws3+1/2 zs7 (1-m13)ws1+1/4 zs7 (1-m33)ws3+1/2 zs9 (1-m12)ws1+1/4 zs9 (1-m32)ws3+zs10 (1-m13)ws1+1/2 zs10 (1-m33)ws3;

zr6=1/2 zs4 (1-m21)ws2+1/4 zs4 (1-m41)ws4+1/2(zs6 (1-m22) (1-r)+zs5 (1-m22) r) ws2+1/4(zs6 (1-m42) (1-r)+zs5 (1-m42) r)ws4+zs8 (1-m21)ws2+1/2 zs8 (1-m41)ws4+1/2 zs9 (1-m22)ws2+1/4 zs9 (1-m42)ws4+1/4 zs2 (1-m32)ws3+1/2 zs2 (1-m52)ws5+1/2 zs3 (1-m33)ws3+zs3 (1-m53)ws5+1/4(zs6 (1-m32) (1-r)+zs5 (1-m32) r)ws3+1/2(zs6 (1-m52) (1-r)+zs5 (1-m52) r)ws5+1/4 zs7 (1-m33)ws3+1/2 zs7 (1-m53)ws5;

zr7=1/4 zs2 (1-m42)ws4+1/2 zs2 (1-m62)ws6+1/2 zs3 (1-m43)ws4+zs3 (1-m63)ws6+1/2(zs5 (1-m22) (1-r)+zs6 (1-m22) r)ws2+1/4(zs5 (1-m42) (1-r)+zs6 (1-m42) r)ws4+1/4(zs6 (1-m42) (1-r)+zs5 (1-m42) r)ws4+1/2(zs6 (1-m62) (1-r)+zs5 (1-m62) r)ws6+1/4 zs7 (1-m43)ws4+1/2 zs7 (1-m63)ws6+1/2 zs7 (1-m23)ws2+1/4 zs7 (1-m43)ws4+1/2 zs9 (1-m22)ws2+1/4 zs9 (1-m42)ws4+zs10 (1-m23)ws2+1/2 zs10 (1-m43)ws4;

zr8=1/4 zs4 (1-m31)ws3+1/2 zs4 (1-m51)ws5+1/4(zs6 (1-m32) (1-r)+zs5 (1-m32) r)ws3+1/2(zs6 (1-m52) (1-r)+zs5 (1-m52) r)ws5+1/2 zs8 (1-m31)ws3+zs8 (1-m51)ws5+1/4 zs9 (1-m32)ws3+1/2 zs9 (1-m52)ws5;

zr9=1/4 zs4 (1-m41)ws4+1/2 zs4 (1-m61)ws6+1/4(zs5 (1-m32) (1-r)+zs6 (1-m32) r)ws3+1/2(zs5 (1-m52) (1-r)+zs6 (1-m52) r)ws5+1/4(zs6 (1-m42) (1-r)+zs5 (1-m42) r)ws4+1/2(zs6 (1-m62) (1-r)+zs5 (1-m62) r)ws6+1/4 zs7 (1-m33)ws3+1/2 zs7 (1-m53)ws5+1/2 zs8 (1-m41)ws4+zs8 (1-m61)ws6+1/4 zs9 (1-m42)ws4+1/2 zs9 (1-m62)ws6+1/4 zs9 (1-m32)ws3+1/2 zs9 (1-m52)ws5 +1/2 zs10 (1-m33)ws3+zs10 (1-m53)ws5;

zr10=1/4(zs5 (1-m42) (1-r)+zs6 (1-m42) r)ws4+1/2(zs5 (1-m62) (1-r)+zs6 (1-m62) r)ws6+1/4 zs7 (1-m43)ws4+1/2 zs7 (1-m63)ws6+1/4 zs9 (1-m42)ws4+1/2 zs9 (1-m62)ws6+1/2 zs10 (1-m43)ws4+zs10 (1-m63)ws6;

Factor[zr1+zr2+zr3+zr4+zr5+zr6+zr7+zr8+zr9+zr10/.{m11->m,m12->m,m13->m,m21->m,m22->m,m23->m,m31->m,m32->m,m33->m,m41->m,m42->m,m43->m,m51->m,m52->m,m53->m,m61->m,m62->m,m63->m}]

Factor[wr1+wr2+wr3+wr4+wr5+wr6+zr1+zr2+zr3+zr4+zr5+zr6+zr7+zr8+zr9+zr10]

standardize within daughters and sons so that the genotype frequencies = 1

msum=zr1+zr2+zr3+zr4+zr5+zr6+zr7+zr8+zr9+zr10;

eqnz1=zr1/msum;
eqnz2=zr2/msum;
eqnz3=zr3/msum;
eqnz4=zr4/msum;
eqnz5=zr5/msum;
eqnz6=zr6/msum;
eqnz7=zr7/msum;
eqnz8=zr8/msum;
eqnz9=zr9/msum;
eqnz10=zr10/msum;

Factor[eqnz1+eqnz2+eqnz3+eqnz4+eqnz5+eqnz6+eqnz7+eqnz8+eqnz9+eqnz10]

fsum=wr1+wr2+wr3+wr4+wr5+wr6;

eqnw1=wr1/fsum;
eqnw2=wr2/fsum;
eqnw3=wr3/fsum;
eqnw4=wr4/fsum;
eqnw5=wr5/fsum;
eqnw6=wr6/fsum;

Factor[eqnw1+eqnw2+eqnw3+eqnw4+eqnw5+eqnw6]

Mate-SRA

Summary of results

The equilibrium and internal stability are the same as with an Z-linked modifier (as expected).

SRA should evolve to produce daughters in response to TT males and sons in response to tt males, as long as the population sex ratio remains near 1/2:
		- (1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))( (pave^2) (m31-m11)+2 (1-pave) pave (m32-m12)+(((1-pave))^2)(m33-m13))+
(3/5) (1-pave)pave s( pave (m31-m11)+(1-2 pave)(m32-m12)- (1-pave)(m33-m13))

Polymorphism requires  2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD > 0 [low dominance coefficients] and the following:

	(s/(2 (1-h)))-(((1-H) )/(2 (1-h)))θ < z < (s/(2 h))-((H )/(2 h ))θ

where θ=((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) )).   As long as the average population sex ratio stays near 1/2, SRA should be able to increase indefinitely.

Definitions [ENTER]

The following sex ratio strategies are equivalent under mate-SRA:

subNoRef={m21->m11,m22->m12,m23->m13,m41->m31,m42->m32,m43->m33,m61->m51,m62->m52,m63->m53};

Frequencies while the modifier is fixed on allele M:

subNoM={z3->1-z1-z2,z4->0,z5->0,z6->0,z7->0,z8->0,z9->0,z10->0,w2->1-w1,w3->0,w4->0,w5->0,w6->0};

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m12->m+H*ΔD,m13->m};

The population sex ratio among the offspring (M fixed):

fracFemale=Factor[zs1 m11+zs2 m12+zs3 m13/.subNoRef/.subNoM]

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in males (Fm) to simplify solutions:

subpF=Solve[{pf==w1,pm==z1+z2/2,Fm==1-(z2/(2 (1-pm) pm))},{w1,z1,z2}]//Simplify//Flatten

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (((1-m) m (s-2 h z)+H(1-2 m) ΔD)/(2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD));

pdif1=((m s-(m^2) s-2 h m z+2 h (m^2) z+H ΔD-2 H m ΔD) (-m s+(m^2) s+2 m z-2 h m z-2 (m^2) z+2 h (m^2) z-ΔD+H ΔD+2 m ΔD-2 H m ΔD) (m s z-2 h m s z-(m^2) s z+2 h (m^2) s z-s ΔD+2 H s ΔD+m s ΔD-2 H m s ΔD+h z ΔD-H z ΔD))/(((2 m z-4 h m z-2 (m^2) z+4 h (m^2) z-ΔD+2 H ΔD+2 m ΔD-4 H m ΔD))^3);

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under self-SRA:

eqnw1/.subNoRef/.subNoM//Factor

eqnz1/.subNoRef/.subNoM//Factor

eqnz2/.subNoRef/.subNoM//Factor

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]-pf,recu[[2]]+1/2recu[[3]]-pm,1-((recu[[3]])/(2 (1-(recu[[2]]+1/2recu[[3]])) (recu[[2]]+1/2recu[[3]])))-Fm}/.subpF]]

Assuming weak selection:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Fm->Fm*sm,{sm,0,1}]]]]/.sm->1

The third term will be zero only if Fm is zero to this order.

solvemeAPPROX/.Fm->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[2]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[1]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave/.%%[[3]]//Factor

Collect[{Numerator[%],Denominator[%]},ΔD,Factor]

which can be written as:

pave0= (((1-m) m (s-2 h z)+H(1-2 m) ΔD)/(2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD));

Factor[%%%-pave0]

Note that this is the same equilibrium as under mate-SRA.  As in the mate-SRA section, we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave^2+(m+H ΔD)(2pave(1-pave))+m(1-pave)^2)},{m}]]

{m+ΔD,m+H ΔD,m}/.%/.(OverBar[D])->(m11*pave^2+m12(2pave(1-pave))+m13(1-pave)^2)/.ΔD->m11-m13/.H->(m12-m13)/(m11-m13)//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

((( 1/2 s-h z)-H θ)/((1-2 h) z+(1-2 H) θ))/.θ->( (OverBar[D])-1/2)((ΔD)/((1-(OverBar[D])) (OverBar[D])));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{(((1-m) m (s-2 h z)+H(1-2 m) ΔD)/(2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD)),(((1-m) m (2 z(1-h)-s)-(1-H)(1-2 m) ΔD)/(2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[(1-m) m (s-2 h z)+H(1-2 m) ΔD==0,z],ΔD,Factor]

{(s/(2h))-((H )/(2 h ))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[(1-m) m (2 z(1-h)-s)-(1-H)(1-2 m) ΔD==0,z],ΔD,Factor]

{(s/(2(1-h)))-(((1-H) )/(2 (1-h)))θ}/.θ->((OverBar[D])-(1/2)) ((ΔD)/((1-(OverBar[D])) (OverBar[D]) ))/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that 2 (1 - 2 h) (1 - m) m z - (1 - 2 H) (1 - 2 m) ΔD > 0 [low dominance coefficients] and

	(s/(2 (1-h)))-(((1-H) )/(2 (1-h)))θ< z < (s/(2 h))-((H )/(2 h ))θ

* Negative numerators and denominators, which requires that 2 (1 - 2 h) (1 - m) m z - (1 - 2 H) (1 - 2 m) ΔD < 0 [high dominance coefficients] and

	(s/(2 (1-h)))-(((1-H) )/(2 (1-h)))θ> z > (s/(2 h))-((H )/(2 h ))θ

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnz1/.subNoM,z1],D[eqnz1/.subNoM,z2],D[eqnz1/.subNoM,w1]},{D[eqnz2/.subNoM,z1],D[eqnz2/.subNoM,z2],D[eqnz2/.subNoM,w1]},{D[eqnw1/.subNoM,z1],D[eqnw1/.subNoM,z2],D[eqnw1/.subNoM,w1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

whatami=constlam/.Flatten[Solve[%==0,constlam]]//Simplify

Factor[%*3 (1-m) m/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable requires that the denominator be positive, 2 (1-2 h) (1-m) m z-(1-2 H) (1-2 m) ΔD > 0, which requires that h be sufficiently small and the ΔD term not be very positive.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat2={{D[eqnz4,z4],D[eqnz4,z5],D[eqnz4,z6],D[eqnz4,z7],D[eqnz4,w3],D[eqnz4,w4]},{D[eqnz5,z4],D[eqnz5,z5],D[eqnz5,z6],D[eqnz5,z7],D[eqnz5,w3],D[eqnz5,w4]},{D[eqnz6,z4],D[eqnz6,z5],D[eqnz6,z6],D[eqnz6,z7],D[eqnz6,w3],D[eqnz6,w4]},{D[eqnz7,z4],D[eqnz7,z5],D[eqnz7,z6],D[eqnz7,z7],D[eqnz7,w3],D[eqnz7,w4]},{D[eqnw3,z4],D[eqnw3,z5],D[eqnw3,z6],D[eqnw3,z7],D[eqnw3,w3],D[eqnw3,w4]},{D[eqnw4,z4],D[eqnw4,z5],D[eqnw4,z6],D[eqnw4,z7],D[eqnw4,w3],D[eqnw4,w4]}}/.subNoRef/.subNoM//Simplify;

r=1/2;

The calculations are simplified by working in terms of the p's and F's:

mat2WeakSRA=mat2/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-mat2WeakSRA]/.m31->m11+sm*dm31/.m32->m12+sm*dm32/.m33->m13+sm*dm33/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm31,dm32,dm33},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

{((dm33 (1-2 m) (((1-pave))^2))/(4 (1-m) m))+((dm32 (1-2 m) (1-pave) pave)/(2 (1-m) m))+((dm31 (1-2 m) (pave^2))/(4 (1-m) m))};
Factor[%-%%]

To leading order, (1-2 m) represents the difference between 1/2 and the population sex ratio [times 2]:

Factor[Normal[Series[2(1/2-fracFemale)/.subpF/.m31->m11+sm*dm31/.m32->m12+sm*dm32/.m33->m13+sm*dm33/.directional/.lambda->1+constlam*sm/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2,{sm,0,0}]]/.sm->1]

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) (pave^2) +2 (m32-m12) (1-pave) pave+ (m33-m13) (((1-pave))^2)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m31->m11+sm*dm31/.m32->m12+sm*dm32/.m33->m13+sm*dm33/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-mat2WeakSRA]/.m31->m11+sm*dm31/.m32->m12+sm*dm32/.m33->m13+sm*dm33/.directional/.lambda->1+constlam*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm31,dm32,dm33,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm31->0/.dm32->0/.dm33->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

Leaving us with:

newstabcond=Collect[stabcond-%%,{dm31,dm32,dm33,ΔD},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) (pave^2) +2 (m32-m12) (1-pave) pave+ (m33-m13) (((1-pave))^2)}+
((2(s+h z+pave z-2 h pave z))/5) (1-pave) pave { (m31-m11) pave+ (m32-m12)  (1-2pave)-  (m33-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m32->m12+dm32/.m33->m13+dm33;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm32->dm32*sm/.dm33->dm33*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(1)-(1-s)]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[(pave (1-h z)+(1-pave))-(pave (1-z)+(1-pave)(1-h z))]

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) (pave^2) +2 (m32-m12) (1-pave) pave+ (m33-m13) (((1-pave))^2)}+
((2(s+(OverBar[z])))/5) (1-pave) pave { (m31-m11) pave+ (m32-m12)  (1-2pave)-  (m33-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m32->m12+dm32/.m33->m13+dm33;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm32->dm32*sm/.dm33->dm33*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because the total selection favoring T is (1/3) (OverBar[s])-((2 )/3)(OverBar[z]) , which must equal zero at equilibrium (assuming a 50:50 sex ratio), we can use the substitution pave -> pave0 to show that the above equals:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) (pave^2) +2 (m32-m12) (1-pave) pave+ (m33-m13) (((1-pave))^2)}+
((3s)/5) (1-pave) pave { (m31-m11) pave+ (m32-m12)  (1-2pave)-  (m33-m13) (1-pave)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m32->m12+dm32/.m33->m13+dm33;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+sm*pave1/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm32->dm32*sm/.dm33->dm33*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

The first component reflects selection for returning the sex ratio toward 1/2: These all have the sign of -((OverBar[D])-(1/2)), so if the fraction of females is less than 1/2, then the proportion of females will be selected to increase, and vice versa.

The second component reflects selection for SRA:  Assuming that the equilibrium is valid, the modifier will invade if it increases investment in daughters when mated to TT males and if it increases investment in sons when mated to tt males.  Whether it is better to invest in daughters or sons when mated with Tt males depends on the current frequency of T.  If freq(T) is very low, then most matings are tt x tt, but a mother who has just mated with a Tt male is more likely to produce Tt children, and these will be very fit among the daughters (compared to the majority of tt females).

Self-SRA

Summary of results

The equilibrium and internal stability are the same as with an Z-linked modifier (as expected).

SRA should evolve to produce daughters within TT females and sons within tt females, as long as the population sex ratio remains near 1/2, but only up to a point:

	-(1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))((m31-m11) pave+dm41 (1-pave))+
	((1/5) (s-4 ΔD)(1-pave) pave((m31-m11)-(m41-m21))  )

when (s - 4 ΔD) becomes zero, then the system will remain polymorphic but SRA will stop evolving (an intermediate ESS).

->  SRA should reach an ESS rather than increase, even while the system remains polymorphic.

Polymorphism requires  (1-2 h) > 0 [low dominance coefficients] and the following:

	(s/(2 (1-h)))-((1 )/(2 (1-h) ))((ΔD)/((1-(OverBar[D])) ))<z<(s/(2 h))-((1 )/(2 h ))((ΔD)/((1-(OverBar[D])) ))
	
SRA is expected to drive ΔD to a more and more positive value, but it will hit the ESS first.

Definitions [ENTER]

The following sex ratio strategies are equivalent under self-SRA:

subNoRef={m12->m11,m13->m11,m22->m21,m23->m21,m32->m31,m33->m31,m42->m41,m43->m41,m52->m51,m53->m51,m62->m61,m63->m61};

Frequencies while the modifier is fixed on allele M:

subNoM={z3->1-z1-z2,z4->0,z5->0,z6->0,z7->0,z8->0,z9->0,z10->0,w2->1-w1,w3->0,w4->0,w5->0,w6->0};

The population sex ratio among the offspring (M fixed):

fracFemale=Factor[ws1 m11+ws2 m21/.subNoRef/.subNoM]

Assuming small differences in the sex ratio strategies:

directional={m11->m+ΔD,m21->m};

We will transform the recursions into allele frequencies in the two sexes (pf, pm) and Wright's coefficient in males (Fm) to simplify solutions:

subpF=Solve[{pf==w1,pm==z1+z2/2,Fm==1-(z2/(2 (1-pm) pm))},{w1,z1,z2}]//Simplify//Flatten

We will show that the equilibrium trait frequency and the difference in freq(T) between the sexes is approximately the following when the population sex ratio is near 1/2 and selection is weak.

pave0= (((1-m) (s-2 h z)-ΔD)/(2 (1-2 h) (1-m) z));

pdif1=(((s-m s-2 h z+2 h m z-ΔD) (-s+m s+ΔD) (-s+m s+2 z-2 h z-2 m z+2 h m z+ΔD))/(8 (((-1+2 h))^2) (((-1+m))^3) (z^2)));

Polymorphic equilibrium (M fixed)

The recursions when fixed for allele M under self-SRA:

eqnw1/.subNoRef/.subNoM//Factor

eqnz1/.subNoRef/.subNoM//Factor

eqnz2/.subNoRef/.subNoM//Factor

recu={%%%,%%,%};

The following works in terms of the allele frequencies (pf, pm) and Wright's coefficients (Ff, Fm), and sets the values of these variables in the next generation to those in the previous generation to find the equilibrium:

solveme=Numerator[Factor[{recu[[1]]-pf,recu[[2]]+1/2recu[[3]]-pm,1-((recu[[3]])/(2 (1-(recu[[2]]+1/2recu[[3]])) (recu[[2]]+1/2recu[[3]])))-Fm}/.subpF]]

To simplify these equilibrium equations, we assume weak selection and a weak difference in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).
Under these assumptions, we know that the difference in allele frequencies between the sexes (pdif) and the departures from Hardy-Weinberg (Ff and Fm) will be small, but the average allele frequency (pave) need not be:

solvemeAPPROX=Numerator[Factor[Normal[Series[solveme/.s->s*sm/.z->z*sm/.directional/.ΔD->ΔD*sm/.pf->pave-pdif*sm/2/.pm->pave+pdif*sm/2/.Fm->Fm*sm,{sm,0,1}]]]]/.sm->1

The third term will be zero only if Fm is zero to this order.

solvemeAPPROX/.Fm->0//Factor

We can then can solve for the difference in allele frequency between the sexes:

Solve[%[[2]]==0,pdif]//Simplify//Flatten

%%/.%//Factor

Solve[%[[1]]==0,pave]//Simplify

pdif1=pdif/.%%%/.%[[3]]//Factor

Finally, we can solve for the average frequency of T:

pave/.%%[[3]]//Factor

which can be written as:

pave0= (((1-m) (s-2 h z)-ΔD)/(2 (1-2 h) (1-m) z));

Factor[%%-pave0]

Note that this is the same equilibrium as under mate-SRA.  As in the mate-SRA section, we can rewrite the average allele frequency in terms of the average sex ratio ((OverBar[D])), ΔD, and the dominance of the modifier, as defined in the manuscript (Equation 1):

Flatten[Solve[{(OverBar[D])==((m+ΔD)*pave+m(1-pave))},{m}]]

{m+ΔD,m}/.%/.(OverBar[D])->(m11*pave+m12(1-pave))/.ΔD->m11-m12//Factor

Simplify[Normal[Series[pave0/.%%/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,0}]]/.sm->1]

(((1/2 s-h z)- ΔD/(1-(OverBar[D]))/2)/(z(1-2 h)));
Factor[%%-%]

We next examine when the equilibrium allele frequency, pave, is a valid equilibrium (falling between 0 and 1):

{pave0,1-pave0} equal:

{(((1-m) (s-2 h z)-ΔD)/(2 (1-2 h) (1-m) z)),(((1-m) (2 z(1-h)-s)+ΔD)/(2 (1-2 h) (1-m) z))};
%-{pave0,1-pave0}//Factor

To be valid, the numerators must have the same sign (the denominator, being the sum, will then have that sign).  Thus, we solve for values of z that allow the numerators to be equal in sign:

Collect[z/.Solve[(1-m) (s-2 h z)-ΔD==0,z],ΔD,Factor]

{(s/(2 h))-((1 )/(2 h ))((ΔD)/((1-(OverBar[D])) ))}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,1}]]/.sm->1]

Collect[z/.Solve[(1-m) (2 z(1-h)-s)+ΔD==0,z],ΔD,Factor]

{(s/(2 (1-h)))-((1 )/(2 (1-h) ))((ΔD)/((1-(OverBar[D])) ))}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[%-%%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->dm41*sm/.dm42->dm42*sm/.dm43->dm43*sm,{sm,0,1}]]/.sm->1]

For both pf and 1-pf to be positive, we require:
* Positive numerators and denominators, which requires that (1-2 h) > 0 [low dominance coefficients] and

	(s/(2 (1-h)))-((1 )/(2 (1-h) ))((ΔD)/((1-(OverBar[D])) ))<z<(s/(2 h))-((1 )/(2 h ))((ΔD)/((1-(OverBar[D])) ))

* Negative numerators and denominators, which requires that (1-2 h) < 0 [high dominance coefficients]and

	(s/(2 (1-h)))-((1 )/(2 (1-h) ))((ΔD)/((1-(OverBar[D])) ))>z>(s/(2 h))-((1 )/(2 h ))((ΔD)/((1-(OverBar[D])) ))

Stability of polymorphic equilibrium (M fixed)

The local stability matrix with the M allele fixed:

mateq=Simplify[{{D[eqnz1/.subNoM,z1],D[eqnz1/.subNoM,z2],D[eqnz1/.subNoM,w1]},{D[eqnz2/.subNoM,z1],D[eqnz2/.subNoM,z2],D[eqnz2/.subNoM,w1]},{D[eqnw1/.subNoM,z1],D[eqnw1/.subNoM,z2],D[eqnw1/.subNoM,w1]}}/.subNoRef];

The characteristic polynomial:

charpoly=Det[mateq-lambda*IdentityMatrix[3]];

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[charpoly/.lambda->1+constlam*sm/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm,{sm,0,1}]]/.sm->1]

constlam/.Solve[%==0,constlam]//Simplify

The above has the sign of:

Factor[%*3(1-m)/(pave(1-pave))/.pave->pave0]

which is the opposite in sign to the denominator for pave0.

Factor[%/Denominator[pave0]]

Thus, for the equilibrium to be stable requires that the denominator be positive, (1-2h) > 0, which requires that h be sufficiently small.

Stability to the introduction of a new SRA allele (m rare)

The local stability matrix once the m allele has been introduced:

mat={{D[eqnz4,z4],D[eqnz4,z5],D[eqnz4,z6],D[eqnz4,z7],D[eqnz4,w3],D[eqnz4,w4]},{D[eqnz5,z4],D[eqnz5,z5],D[eqnz5,z6],D[eqnz5,z7],D[eqnz5,w3],D[eqnz5,w4]},{D[eqnz6,z4],D[eqnz6,z5],D[eqnz6,z6],D[eqnz6,z7],D[eqnz6,w3],D[eqnz6,w4]},{D[eqnz7,z4],D[eqnz7,z5],D[eqnz7,z6],D[eqnz7,z7],D[eqnz7,w3],D[eqnz7,w4]},{D[eqnw3,z4],D[eqnw3,z5],D[eqnw3,z6],D[eqnw3,z7],D[eqnw3,w3],D[eqnw3,w4]},{D[eqnw4,z4],D[eqnw4,z5],D[eqnw4,z6],D[eqnw4,z7],D[eqnw4,w3],D[eqnw4,w4]}}/.subNoRef/.subNoM//Simplify;

r=1/2;

The calculations are simplified by working in terms of the p's and F's:

mat2WeakSRA=mat/.subpF//Factor;

To simplify the leading eigenvalue, we assume weak selection and weak differences in SRA (all small terms are proportional to "sm", a small term, with respect to which we perform a Taylor Series).  
Note: As we let sm -> 0, we expect the leading eigenvalue to approach one (neutral case), so we represent the leading eigenvalue as a power series of the form λ = 1 + (some constant)*sm + (some other constant)*(sm^2)...

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-mat2WeakSRA]/.lambda->1+constlam*sm/.m31->m11+sm*dm31/.m41->m21+sm*dm41/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

stabcond=Collect[constlam/.Solve[%==0,constlam],{dm31,dm41},Factor]

Each of these terms is proportional to the effect of the modifier times (1-2 m):

{((dm41 (1-2 m) (1-pave))/(4 (1-m) m))+((dm31 (1-2 m) pave)/(4 (1-m) m))};
%%-%//Factor

so that stabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) pave +(m41-m21) (1-pave)}/.(OverBar[D])->fracFemale;

Factor[Normal[Series[stabcond-%/.subpF/.m31->m11+sm*dm31/.m41->m21+sm*dm41/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2,{sm,0,1}]]/.sm->1]

To leading order, (1-2 m) represents the difference between 1/2 and the population sex ratio [times 2]:

Factor[Normal[Series[2(1/2-fracFemale)/.subpF/.m31->m11+sm*dm31/.m41->m21+sm*dm41/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2,{sm,0,0}]]/.sm->1]

Thus, invasion will occur if modifiers increase allocation to the rarer sex (i.e., to daughters if the sex ratio is below 1/2).

Assuming that the sex ratio is near 1/2, we can get the next order term (proportional to (sm^2)), which clarifies the evolution of SRA:

Factor[Normal[Series[Det[IdentityMatrix[6]*lambda-mat2WeakSRA]/.lambda->1+constlam*sm^2/.m31->m11+sm*dm31/.m41->m21+sm*dm41/.directional/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Fm->Fm*sm^2/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

stabcond=Collect[constlam/.Flatten[Solve[%==0,constlam]],{dm31,dm41,ΔD},Factor]

The terms not obviously involving the effects of the modifier cancel out when we substitute for pave and pdif:

Factor[stabcond/.dm31->0/.dm41->0]

Factor[Normal[Series[%/.pave->pave0+pave1*sm/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.m->(m-1/2)*sm+1/2,{sm,0,2}]]/.sm->1]

Leaving us with:

newstabcond=Collect[stabcond-%%,{dm31,dm41, ΔD},Factor]

The terms not involving s or z are related to the Fisherian sex ratio selection identified above, so that newstabcond can be written as:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) pave +(m41-m21) (1-pave)}+
+(((4 h z+4 pave z-8 h pave z-s))/5) (1-pave) pave { (m31-m11)-  (m41-m21)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m41->m21+dm41;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm41->dm41*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

If we define (OverBar[s]) as the difference between the fitness of the more fit T and the less fit t allele in females:

(OverBar[s])=Factor[(1)-(1-s)]

and (OverBar[z]) as the difference between the fitness of the more fit t and the less fit T allele in males:

(OverBar[z]) =Factor[(pave (1-h z)+(1-pave))-(pave (1-z)+(1-pave)(1-h z))]

and κ as the difference between the transmission of the t and the T allele in females:

κ=((ΔD)/(1-(OverBar[D])))/.(OverBar[D])->fracFemale/.directional;
Factor[Normal[Series[%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42/.m->diff*sm+1/2,{sm,0,1}]]/.sm->1]

With these definitions, the total selection favoring T is (1/3) (OverBar[s])-((κ )/3)-(2/3)(OverBar[z]), which should be an equivalent means to derive pave (even when the sex ratio isn't near 1/2):

Factor[Normal[Series[(1/3) (OverBar[s])-((κ )/3)-(2/3)(OverBar[z])/.(OverBar[D])->fracFemale/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm41->sm*dm41/.dm42->sm*dm42,{sm,0,1}]]/.sm->1]

pave/.Solve[%==0,pave];
%-pave0//Factor

Then the modifier can be more simply written as in the text (note that the SRA matrix is written in a different form here to encompass both mate-SRA and self-SRA, see Intro):

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) pave +(m41-m21) (1-pave)}+
+(((4 (OverBar[z])-s))/5) (1-pave) pave { (m31-m11)-  (m41-m21)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m41->m21+dm41;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm41->dm41*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

Because the total selection favoring T is (1/3) (OverBar[s])-((κ )/3)-(2/3)(OverBar[z]) , which must equal zero at equilibrium, we can use the substitution pave -> pave0 to show that the above equals:

-((((OverBar[D])-(1/2)))/(2(OverBar[D])(1-(OverBar[D])))){(m31-m11) pave +(m41-m21) (1-pave)}+
+((( s-2((ΔD)/(1-(OverBar[D])))))/5) (1-pave) pave { (m31-m11)-  (m41-m21)}/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m41->m21+dm41;

Factor[Normal[Series[newstabcond-%/.subpF/.directional/.pf->pave-pdif/2/.pm->pave+pdif/2/.pave->pave0+sm*pave1/.pdif->pdif*sm/.Ff->Ff*sm^2/.Fm->Fm*sm^2/.pdif->pdif1+pdif2*sm^2/.s->s*sm/.z->z*sm/.ΔD->ΔD*sm/.dm31->dm31*sm/.dm41->dm41*sm/.m->diff*sm+1/2,{sm,0,2}]]/.sm->1]

This is the form used in Appendix C.

We know that for the equilibrium to be valid and stable, we must have (1-2h)> 0 and 

	(s/(2 (1-h)))-((1 )/(2 (1-h) ))((ΔD)/((1-(OverBar[D])) ))<z<(s/(2 h))-((1 )/(2 h ))((ΔD)/((1-(OverBar[D])) ))
		
which for (OverBar[D]) ~ 1/2 implies that s -2ΔD is positive and must lie above 2*z*h:

	(s/(2 (1-h)))-((ΔD)/((1-h)))<z<(s/(2 h))-((ΔD)/h)
	
But this doesn't require that s-4ΔD be positive.  The following is a counter-example:

numpar={z->1/100,s->1/100,h->1/10,ΔD->35/10000,dm31->0.001,dm41->0};

We must first solve for the equilibrium and the value of m that allows the sex ratio to be ~1/2:

recu/.directional/.numpar//Simplify

fracFemale/.directional/.numpar//Simplify

numx=NSolve[{%%[[1]]==w1,%%[[2]]==z1,%%[[3]]==z2,%==1/2},{m,w1,z1,z2},25]

The solution of interest:

numx=numx[[3]]

{pave0,1-pave0,fracFemale,s-4ΔD}/.directional/.numpar/.numx//N

While we would expect mothers mated to T to produce more daughters (dm31 > 0), such alleles don't invade

Collect[N[-(1/2)((((OverBar[D])-(1/2)))/((OverBar[D])(1-(OverBar[D]))))((m31-m11) pave+dm41 (1-pave))+((1/5) ((m31-m11)-(m41-m21)) (1-pave) pave (s-4 ΔD))/.(OverBar[D])->fracFemale/.m31->m11+dm31/.m41->m21+dm41/.pave->pave0/.directional/.numpar/.numx,25],dm41]

The numerically exact leading eigenvalue has the same order of magnitude:

N[Eigenvalues[mat/.m31->m11+dm31/.m41->m21+dm41/.directional/.numpar/.numx/.sm->1]-1]

It is still the case at this equilibrium that T is favored in females and t is favored in males by natural selection:

(ws1-w1)/.subNoM/.numx/.numpar

(zs1-z1)+(zs2-z2)/2/.subNoM/.numx/.numpar

Overall the freq of T goes up during natural selection:

1/3*%%+2/3*%

Here we look at the change in allele frequency during reproduction:

(eqnw1-w1)/.subNoRef/.directional/.subNoM/.numx/.{z->0,s->0,h->0}/.numpar

(eqnz1-z1)+(eqnz2-z2)/2/.subNoRef/.directional/.subNoM/.numx/.{z->0,s->0,h->0}/.numpar

Overall the freq of T goes down up during reproduction because males bearing T alleles are more likely to produce daughters:

1/3*%%+2/3*%

The eigenvector associated with the leading eigenvalue tells us the direction in which the new m modifier evolves:

z4=TM/Tm
z5=TM/tm
z6=Tm/tM
z7=tM/tm
w3=TMm
w4=tMm

N[Eigenvectors[mat/.m31->m11+dm31/.m41->m21+dm41/.directional/.numpar/.numx/.sm->1]][[1]]

The above points in a direction associated with an excess of females:

eigenSR=%.{0,0,0,0,1,1}/%.{1,1,1,1,1,1}

Calculating the genotypic frequencies among sons and daughters along this eigenvector:

{{%%[[1]],%%[[2]]+%%[[3]],%%[[4]]}/Drop[%%,-2].{1,1,1,1},Drop[%%,4]/Drop[%%,4].{1,1}}

{{z1,z2,1-z1-z2},{w1,1-w1}}/.numx

The modifier also becomes associated with tt sons and t daughters.

%%-%

The above associations with  tt sons and t daughters were repeated in other similar cases (e.g., dm31 -> 0, dm41 -> -0.001), but the effects of the modifier on the sex ratio depended on the action of the modifier.

With more intermediate values of h, however, the polymorphism tended to be lost before the ESS was reached:

numpar={z->1/100,s->1/100,h->4/10,ΔD->1/500};

We must first solve for the equilibrium and the value of m that allows the sex ratio to be ~1/2:

recu/.directional/.numpar//Simplify

fracFemale/.directional/.numpar//Simplify

numx=NSolve[{%%[[1]]==w1,%%[[2]]==z1,%%[[3]]==z2,%==1/2},{m,w1,z1,z2},25]

None of the above solutions allow a polymorphic equilibrium.