Multilocus evolution

Here I develop exact recursions for a simple two-locus model with random mating, and show that recusions developed focussing on changes in genotype frequencies are equivalent to recursions using the Kirkpatrick et al. (2002) approach of focussing on changes in genetic associations.  The test case involves two biallelic loci i and j, each of which may be in a "1" state or "0" state.

First, I set up the population and calculate the initial genetic associations within it (section 2).  Then I work through selection (section 3), and show that associations after selection (D') are the same with both methods (section 3.2.3 and 3.2.4), and that the change in allele frequencies are the same (section 3.3).  After selection, I've recalculated the associations to the new allele frequencies, to give (D^(′* )) values (these don't exist in the Kirkpatrick paper), and shown these to be the same with both approaches (section 4.2).  These recalculated values are then transmitted to the next generation to create D''' values, which are the same with both approaches (section 5.2).  In section 5.2.2, I transmit associations immediately after selection, and then recalculate to the new allele frequencies.

This file contains a bit of a mix of things from quite general (and that should be easy to adapt when generalising this further), and things that are not general at all (in particular, things that depend on global variables or particular things about the system (locus names, in particular)).  This is the first large Mathematica file I've written, so it's a fair mix of styles.

. Some useful functions:

ExpandGrid: clone of R's expand.grid()
Rep: Simplify repeating an expression (repeats x n times).
Combinations: The "correct" order for cycling through n positions (I don't have a good definition for this).

ExpandGrid = Flatten[Outer[List, ##],Length[List[##]]-1] &;
Rep = Function[{x,n},Table[x,{n}]];
Combinations = Function[n,ExpandGrid @@ Rep[{1,0},n]];

. Genotypes and associations

.. Introduction & general set up

For a diploid organism with two loci there are four "positions" a gene may be in (two loci × two sexes of origin, and therefore four haplotypes.  These combine to give 16 possible diploid genotypes.  Under random mating, the diploid genotypes are constructed based on the product of the frequencies of their haplotypes, and the frequency of each haplotype is a function of the allele frequencies ((Subscript[p, i]) and (Subscript[p, j])) and the linkage disequilibrium between them ((Subscript[D, ij])).

haplotypes holds the frequencies of haplotypes, and genotypes holds the frequencies of genotypes, in terms of allele frequencies and the gametic phase linkage disequilibrium.  loci holds the names of the loci, for convenience.

loci = {i,j};
haplotypes = {(Subscript[p, i])(Subscript[p, j])+(Subscript[D, ij]), (Subscript[p, i])(1-(Subscript[p, j]))-(Subscript[D, ij]),(1-(Subscript[p, i]))(Subscript[p, j])-(Subscript[D, ij]),(1-(Subscript[p, i]))(1-(Subscript[p, j])) + (Subscript[D, ij])};
genotypes = Outer[Times,haplotypes,haplotypes] // Flatten;
Total[haplotypes] == Total[genotypes] == 1 // Factor

Let (Subscript[z, xy]) denote the position (or context) of a gene; z refers to the locus (i or j, above), x denotes the sex-of-carrier, and is typically omitted here, as the organisms are hermaphroditic, and y denotes the sex-of-origin; whether the gene was transmitted via eggs from the female parent (f) or through sperm from the male parent (m).

SOO = {f,m};
positions = Outer[Subscript[#2,#1] &, SOO, loci] // Flatten

Let (Subscript[X, 𝕚]) be an indicator variable that is 1 where the allele at position 𝕚 is 1, and 0 otherwise; the names of these variables are stored in Xs.

Xs = Thread[Subscript[X,positions]]

Construct a matrix of indicator variables, XX, with four columns correponding to positions, and 16 rows corresponding to genotypes.  To simplify things later, construct a set of four substitution rules, XXSubs, that each correspond to these 16 genotypes

XX = Combinations[4];
XXSubs = Thread[Xs -> Transpose[XX]];

.. Calculating associations

There are 16 different genetic associations that describe the population; these correspond to members of the power set of the set of positions.  The names of the different positions are stored in the vector Ds.  (Even though there are 16 of these, the four one-position associations are zero, and the zero-position association is one, by definition).

Ds=Thread[(Subscript[D, Subsets[positions]])]

The genetic association for a set of positions 𝔸 is
	(Subscript[D, 𝔸])=(Subscript[E, X])[(Subscript[ζ, 𝔸])]=(Subscript[E, X])[(∏_(𝕚⊆𝔸))((Subscript[X, 𝕚])-(Subscript[p, 𝕚]))]
so it is useful to define the (Subscript[ζ, 𝕚]) variables (zetas).

zetas = Xs - (Xs /. (Subscript[X, (Subscript[x_, _])])->(Subscript[p, x]));

The wrinkle is that (Subscript[p, i]) etc may occur in both the genotype frequencies and in relation to the indicator variables, but these need not be the same, so a function that generates the associations has to be willing to accept genotypes and ζs separately.  Associations is a function that calculates genetic assocaitions given a vector of genotype frequencies (constructed in the above way) and ζ values, by taking the sum of the products of genotype frequency and appropriate ζs.  

Because this uses XXSubs, Associations is not completely general; it should either accept XXSubs as an argument or automatically generate it.

Associations = Function[{genotypes,zetas},
(Total[genotypes * (# /. XXSubs)] & /@ (Times @@@ Subsets[zetas] ))//Factor];

Use this to calculate the initial genetic associations (D0):

D0 = Thread[Ds->Associations[genotypes,zetas]];
D0 // TableForm

.. Calculating genotypes

The inverse operation to calculating associations from genotypes is to calculate genotypes from associations.  Equation (4) in Kirkpatrick et al. (2002) gives the equation for doing this:
	f(X)=(∏_(𝕚∈𝔾))((Subscript[X, 𝕚])(Subscript[p, 𝕚])+(1-(Subscript[X, 𝕚]))(1-(Subscript[p, 𝕚])))+(∑_((Overscript[𝕌, *])⊆𝔾))((((Subscript[D, 𝕌])(-1))^((|𝕌|-(∑_(𝕚∈𝕌))(Subscript[X, 𝕚]))))(∏_(𝕚∈𝔾\𝕌))((Subscript[X, 𝕚])(Subscript[p, 𝕚])+(1-(Subscript[X, 𝕚]))(1-(Subscript[p, 𝕚]))))
(note that the starred term on the right does not include 𝕌=∅, and the term corresponding to this is the first term on the right ((((Subscript[D, ∅])(-1))^0)=1)).  

The function GenotpyesExpr takes a list of positions, and generates an expansion of the above equation in terms of X's and D's (genotypesExpr, for the diploid, two locus case.).

GenotypeExpr=Function[positions,
Total[Function[𝕌,
(Subscript[D, 𝕌])(-1)^(Length[𝕌]-Sum[(Subscript[X, 𝕚]),{𝕚,𝕌}])Product[(Subscript[X, 𝕚])(Subscript[p, 𝕚[[1]]])+(1-(Subscript[X, 𝕚]))(1-(Subscript[p, 𝕚[[1]]])),{𝕚,Complement[positions,𝕌]}]]/@Subsets[positions]]];

Substituting the current associations (D0) and the matrix of indicator variables (XXSubs) into this expression recovers the original genotype frequencies, in terms of p's and (Subscript[D, ij]).  The second expression shows that the recovered genotypes are correct.

genotypeExpr = GenotypeExpr[positions]
(genotypeExpr/.D0 /. XXSubs) == genotypes//Simplify

This same approach can be used to generate the haplotype frequencies, too, but using just two of the positions.

GenotypeExpr[positions[[1;;2]]] /. D0 /.  Thread[Xs[[1;;2]] -> Transpose[Combinations[2]]]
% == haplotypes

. Selection

.. Specifying selection regimes

With selection, I'm attempting to be as general as possible, but also have a reasonable test case.  The general case allows for completely arbitrary selection regimes, and the specific case uses the epistatic selection model in Appendix 1 of Kirkpatrick et al. 2002.

... Genotypes-base selection

Wpoly is a polynomial in terms of indicator variables (the Xs from above), and some selection coefficients (c).

Wpoly = Total[Thread[Subscript[c,Subsets[positions]]] *
(Times @@@ Subsets[Xs])]

Selection regimes may be specified in terms of c parameters: epiSubs is an epistatic selection scheme, and neuSubs is a neutral scheme.

The epistatic selection regime is a set of c parameter, epiSubs, and the neutral selection regime is

epiSubs = Join[{(Subscript[c, {}])->1},
Thread[ {(Subscript[c, {(Subscript[i, f])}]), (Subscript[c, {(Subscript[i, m])}])} ->(Subscript[s, i])],
Thread[ {(Subscript[c, {(Subscript[j, f])}]), (Subscript[c, {(Subscript[j, m])}])} ->(Subscript[s, j])],
Thread[ {(Subscript[c, {(Subscript[i, f]),(Subscript[j, f])}]), (Subscript[c, {(Subscript[i, f]),(Subscript[j, m])}]), (Subscript[c, {(Subscript[j, f]),(Subscript[i, m])}]), (Subscript[c, {(Subscript[i, m]),(Subscript[j, m])}])} ->(Subscript[e, ij])],
{(Subscript[c, {(Subscript[i, f]),(Subscript[i, m])}])->0,(Subscript[c, {(Subscript[j, f]),(Subscript[j, m])}])->0},
(Subscript[c, #])->0 &/@ Subsets[positions,{3,4}]];
neuSubs = Prepend[(Subscript[c, #])->0 &/@ Subsets[positions,{1,4}],(Subscript[c, {}])->1];

These can be substituted into the fitness polynomial to give simplified fitness equations:

WpolyEpi = Wpoly /. epiSubs // Simplify
WpolyNeu = Wpoly /. neuSubs // Simplify

These can be converted into traditional-looking w values (illustrated below for the epistatic case).  The row names indicate the state of the i and j loci for the maternally-inhereted positions, and the column names for the paternally inhereted positions.

TableForm[Partition[WpolyEpi /.XXSubs,4],TableHeadings->Rep[Combinations[2],2]]

... Converting to Kirkpatrick et al. coefficients:

The "c" selection coefficients above cannot be directly used; they have to be converted into "a" selection coefficients.  This is the conversion described in Appendix 1 in the paper.  First construct a fitness polynomial in Xs and D values (Kpoly), using equation (7)
	W(X) = (OverBar[w])(1 + (∑_(𝕌⊆𝕎))(Subscript[a, 𝕌])((Subscript[ζ, 𝕌]) - (Subscript[D, 𝕌])))
Unfortunately, and for reasons I cannot explain, I have to set (Subscript[D, ∅]) to zero, when it should really be 1.  I'm not sure if this is something wrong with the paper or my interpretation though.  For the rest of the workbook, (Subscript[D, ∅])=1.  Note  that I don't use the "1+" in the equation above either (these problems are definately related)

First, write out the polynomial in terms of a coefficients and associations:

coefA = Thread[(Subscript[a, Subsets[positions]])];
Kpoly = (OverBar[w]) Total[coefA*(Times @@@ Subsets[zetas] - Ds /. (Subscript[D, {}])->0)]

Set these equations to be equal to one another, and collect different combinations of X's and sove this system of equations "a" coefficients in terms of "c" coefficients. The selection coefficients can be solved by manually folding through the list, or by passing it directly to Solve.  These give the same answer, but for larger numbers of loci one may be preferable to the other; here I'm just using Solve, but see the previous version for a manual fold.  coefACSubs contains substitutions for "a" coefficients in terms of "c" coeffients.

coefACSubs = Solve[Wpoly - Kpoly == 0 /. XXSubs,coefA]//Flatten;

The messy (Subscript[a, ∅]) coefficient simplifies to 1, for any selection regime:

(Subscript[a, {}])/.coefACSubs /. (OverBar[w])->Total[genotypes*(Wpoly/.XXSubs)] /. D0 // Factor

(this uses the expression for (OverBar[w]) that is generated in the next section.) Given this, it makes sense to replace (Subscript[a, ∅]) with 1 in the substitutions.

coefACSubs = ReplacePart[coefACSubs,{1,2}->1];

The substitutions above can be used to find the Kirkpatrick coefficients for the epistatic selection model  (kEpiSubs).

kEpiSubs = coefACSubs /. epiSubs // Factor;
kEpiSubs // TableForm
coefACSubs /. neuSubs // Factor

.. Applying selection

There are two steps here; apply selection to generate new genotype frequencies and calculate the associations among these genotypes.  Using the traditional (genotype-based) approach, first calculate the new genotype frequencies, then calculate the associations.  Using the second, the new associations are calculated, then the genotypes are calculated from that.

Below, I'm presenting the easy step for each first (changing the genotypes for the genotype-based approach and changing the associations for the associations-based  approach), then the opposite step and demonstrating equivalency between the approaches.

... Genotypes approach, changing genotypes

The change in genotypes here is very easy to calculate (genotypesP is a vector of post association genotype frequencies [P for prime]):

genotypesP = genotypes * (Wpoly /. XXSubs) / (OverBar[w]);
wbar = Total[genotypes * (Wpoly /. XXSubs)] // Factor;

genotypesPEpi is a vector of the 16 genotype frequencies for the epistatic selection model:

genotypesPEpi = genotypesP/.(OverBar[w])->wbar/. epiSubs

... Associations-based, changing associations

This is less straightforward, and there's a fair bit of information on this in the previous version of this file.

The function ReduceD removes duplicated positions, recursively applying Kirkpatrick's equation (5):
	(Subscript[D, 𝕌𝕚𝕚])=(Subscript[p, 𝕚])(Subscript[q, 𝕚])(Subscript[D, 𝔸]) + (1-2(Subscript[p, 𝕚]))(Subscript[D, 𝔸𝕚])

ReSortD = Function[x,x/.(Subscript[D, 𝔸_]):>(Subscript[D, SortBy[𝔸,#[[2]]&]])];
ReduceD = Function[x,
x//.(Subscript[D, {𝕛___, (Subscript[𝕚_, c_]),(Subscript[𝕚_, c_]),𝕜___}]):>(Subscript[D, {𝕛,𝕜}])(Subscript[p, 𝕚])(1-(Subscript[p, 𝕚]))+(1-2(Subscript[p, 𝕚]))(Subscript[D, {𝕛,(Subscript[𝕚, c]),𝕜}])];

Next, the function assocSelection applies selection to associations, for one set A.  Mapping this over all sets generates the associations after selection, in terms of the original allele frequencies (Dprime).

assocSelection = Function[𝔸,
(Subscript[D, 𝔸])+Total[Function[𝕌,(Subscript[a, 𝕌])((Subscript[D, Join[𝔸,𝕌]])-(Subscript[D, 𝔸]) (Subscript[D, 𝕌]))] /@ Subsets[positions]]];
Dprime =assocSelection /@ Subsets[positions] //ReSortD// ReduceD;
DprimeSubs = Thread[Ds->Dprime];

This is a monsterous list of associations, but the associations for the epistatsis model (DprimeEpi) are easier to display:

DprimeEpi = Dprime /. kEpiSubs  /. (OverBar[w])->(wbar /. epiSubs)/. D0 // Simplify;
DprimeEpiSubs = Thread[Ds -> DprimeEpi];
DprimeEpiSubs // TableForm

... Genotypes-based, changing associations

The associations calculated above (DprimeEpi) are the same as the associations among loci from the genotypes-based approach:

DprimeEpi2 = Associations[genotypesPEpi,zetas] // Simplify;
And @@ Thread[DprimeEpi ==DprimeEpi2] // Simplify

... Associations-based, determining genotypes

The approach in section 2.3 can be used to calculate the genotype frequencies from the genetic associations:

(genotypeExpr/.DprimeEpiSubs/.XXSubs)==(genotypesP/.(OverBar[w])->wbar/. epiSubs)//Simplify

.. Comparison of changes in genotype frequencies

After selection, the "1" allele frequency at the ith locus has changed from (Subscript[p, i]) to
	(Subscript[p, i]^′)=(Subscript[p, i])+Δ (Subscript[p, i])=(Subscript[p, i])+(Subscript[D, i]^′)
Since transmission via meiosis will not change allele frequencies further, it is worth confirming that the methods agree on allele frequencies:  pPrimeGenotype is a vector containing the new allele frequencies calculated from the genotype-based approach, and pPrimeAssoc contains the allele frequencies calculated from the associations-based approach (using just the epistatic model).  The last line shows that these are equivalent.

pPrimeGenotype={Total[((Subscript[X, (Subscript[i, f])])+(Subscript[X, (Subscript[i, m])])/.XXSubs)*genotypesPEpi/2]//Factor,
Total[((Subscript[X, (Subscript[j, f])])+(Subscript[X, (Subscript[j, m])])/.XXSubs)*genotypesPEpi/2]//Factor}
pPrimeAssoc={(Subscript[p, i])+(Subscript[D, {(Subscript[i, f])}]),(Subscript[p, j])+(Subscript[D, {(Subscript[j, f])}])}/.DprimeEpiSubs
pPrimeGenotype == pPrimeAssoc // Simplify

pPrimeSubs is a set of substitutions for these new allele frequencies

pPrimeSubs = Thread[{(Subscript[p, i]),(Subscript[p, j])}->pPrimeAssoc]

. Recalculating associations (pre-transmission)

One strategy is to recalculate genetic associations immediately after selection, since meiosis will not change the allele frequencies, which are our reference values.  This converts the D' associations from above to (D^(′* )) associations.  First, calculate the new value of the associations.  Then compare these against manually calculated associations among the post selection genotype frequencies

.. Computing recalculated associations:

The assocations after recalculation are:
	(Subscript[D, 𝔸]^(′* ))=(∑_(𝕌⊆𝔸))((Subscript[D, 𝔸\𝕌]^′) (∏_(𝕚∈𝕌))((Subscript[p, 𝕚])-(Subscript[p, 𝕚]^′)))
	=(∑_(𝕌⊆𝔸))((Subscript[D, 𝔸\𝕌]^′) (∏_(𝕚∈𝕌))(-(Subscript[D, 𝕚]^′)))
The vector DpsDprime holds these recalculated post-selection associations ((D^(′* ))), in terms of D' associations.

RecalcD=Function[𝔸,Total[(Subscript[D, Complement[𝔸,#]])Product[-(Subscript[D, {𝕚}]),{𝕚,#}] & /@Subsets[𝔸]]//ReSortD];
DpsDprime=RecalcD /@ Subsets[positions]

Dps holds recalculated associations in terms of associations and allele frequencies before selection began (and is very large), wth DpsEpi being the epistatic case. DpsSubs, DpsEpiSubs are substitution versions of Dps and DpsEpi.  To keep things managable, I have substituted in the values of the associations before selection was applied (D0).

(* Dps =DpsDprime /. DprimeSubs /. D0; *)
(* DpsSubs = Thread[Ds->Dps] *)
DpsEpi = DpsDprime /. DprimeEpiSubs;
DpsEpiSubs = Thread[Ds->DpsEpi];

.. Comparing associations

Note first the zero- and one-position associations now have their expected values that they should take by definition: (Subscript[D, ∅]) is one since the product over the empty set is one, and the one-position associations, (Subscript[D, 𝕚]), are zero, since E[(Subscript[ζ, 𝕚])]=0 (i.e. the expected deviation of the average allele frequency from the average allele frequency is zero).

DpsEpiSubs[[1;;5]]

The other associations do not have such a simple form, but can be shown to be the same as the associations among post-selection genotypes, with the reference values set to the new allele frequencies.

DpsEpi2=Associations[genotypesPEpi,zetas /.pPrimeSubs ];
DpsEpi==DpsEpi2//Simplify

. Transmission

Following selection, meiosis occurs, which changes the asociations again.  With a random mating, associations involving two sexes of origin will be zero (except for (Subscript[D, (Subscript[i, f])(Subscript[j, f])(Subscript[i, m])(Subscript[j, m])])=(Subscript[D, (Subscript[i, f])(Subscript[j, f])])(Subscript[D, (Subscript[i, f])(Subscript[j, f])])=(Subscript[D, ij]^2)), associations with one position are zero by definition, leaving only (Subscript[D, (Subscript[i, f])(Subscript[j, f])])and (Subscript[D, (Subscript[i, f])(Subscript[j, f])]) with interesting values.  Here, I'm working through things to show that these are the same, and that (Subscript[D, (Subscript[i, f])(Subscript[j, f])])and (Subscript[D, (Subscript[i, f])(Subscript[j, f])]) are the same as (Subscript[D, ij]) in the genotypes-based approach.

.. Genotypes-based

xm is a manually constructed a transition matrix that gives the probability of generating one of the four haplotypes (i,j) ->(1,1),(1,0),(0,1),(0.0) from one of the 16 genotypes.  These transition probabilities involve the recombination rate, r.

xm = {
{1, 0, 0, 0}, {1/2 ,1/2, 0, 0}, {1/2, 0, 1/2, 0}, {(1-r)/2, r/2, r/2, (1-r)/2},
{1/2, 1/2, 0,0}, {0, 1, 0, 0}, {r/2, (1-r)/2, (1-r)/2, r/2}, {0, 1/2, 0, 1/2},
{1/2,0,1/2,0}, {r/2, (1-r)/2, (1-r)/2, r/2}, {0, 0, 1, 0}, {0, 0, 1/2, 1/2},
{(1-r)/2, r/2, r/2, (1-r)/2}, {0, 1/2, 0, 1/2}, {0, 0, 1/2, 1/2}, {0, 0, 0, 1}};
xm // MatrixForm;

haplotypesPP is a vector of the frequencies of the haplotypes in the gamete pool leading to the next generation:

haplotypesPP = genotypesP . xm;

The three key variables needed to describe genotypes ((Subscript[p, i]), (Subscript[p, j]) and (Subscript[D, ij])) can be calculated from these haplotype frequencies.  recGenotype is a set of substitutions containing these three variables after one entire generation:

recGenotype={(Subscript[p, i])->(Subscript[x, 1])+(Subscript[x, 2]),
(Subscript[p, j])->(Subscript[x, 1])+(Subscript[x, 3]),
(Subscript[D, ij])->(Subscript[x, 1])(Subscript[x, 4])-(Subscript[x, 2])(Subscript[x, 3])}/.Thread[Thread[(Subscript[x, Range[4]])]->haplotypesPP];

The above expression is very long, and it is simper just to consider the epistasis case:

recGenotypeEpi =recGenotype/.(OverBar[w])->wbar/. epiSubs//Simplify

Note that transmission does not change the allele frequencies:

(Thread[(Subscript[p, loci])]/.recGenotypeEpi)==pPrimeGenotype // Simplify

Use the recursion substitutions above to transform the population through one complete generation and generate the 16 dipoid genotypes:

genotypesPP= genotypes/.recGenotype;
genotypesPPEpi= genotypes/.recGenotypeEpi;

However, it might be more natural to express this using substitutions based on the previous allele frequencies and disequilbria, plus some deviation ((Subscript[Δp, i]), (Subscript[Δp, j]), or (Subscript[ΔD, ij])):

nextGenSubs = {(Subscript[p, i])->(Subscript[p, i])+(Subscript[Δp, i]),(Subscript[p, j])->(Subscript[p, j])+(Subscript[Δp, j]), (Subscript[D, ij])->(Subscript[D, ij])+(Subscript[ΔD, ij])};
genotypesPP2 = genotypes/.nextGenSubs;

This makes the genetic associations easier to interpret.  Where ζ variables still include the previous allele frequencies, the associations after transmission but before recalculation (D'') are:

Dpp=Associations[genotypesPP2, zetas];
Thread[Ds->Dpp]//TableForm

To calculate the associations after transmission and recalculation (D'''), use the new allele frequencies ((Subscript[p, i])+(Subscript[Δp, i])) in the zetas:

Dppp=Associations[genotypesPP2, zetas /. nextGenSubs];
Thread[Ds->Dppp]//TableForm

.. Associations-based

There are two effects of transmission to work through; those that have been recalculated and those that have not.  Some of these are straightforward (particulary those that have only one sex of origin), but some of these are fairly messy and involved at the moment.  Much of the approach is shared between post and pre recalculation.

Transmission is described by equation (12):
	(Subscript[D, 𝔸]^(′′))=(∑_(𝕌: U=A))(Subscript[t, 𝔸←𝕌])(Subscript[D, 𝕌]^′)
The first step is determining the sets of positions 𝕌 from which 𝔸 might be transmitted.  The function TransPos takes a set of positions 𝔸 and returns a list of sets of positions 𝕌, that includes sex-of-carrier information.

TransPos=Function[𝔸,
If[Length[𝔸]==0,{{}},
MapThread[(Subscript[#1[[1]], {#1[[2]],#2}])&,{𝔸,#}]&/@ExpandGrid@@Rep[{f,m},Length[𝔸]]]];

To calculate the probability that a set of positions pos are transmitted together, consider positions with the same sex of carrier independently.  For neighbouring loci i and j, the recombination rate is (Subscript[r, ij]) (in theory this can depend on the sex of carrier too).  Therefore, the probability of transmission of positions at loci i and j with the same sex of origin is (1-(Subscript[r, ij]))/2, and with different sexes of origin with (Subscript[r, ij])/2.  Once the transmission probabilities within each sex of carrier are computed, the product over all sexes of carrier is taken to give the final transmission coefficient.

PosSplitSOC takes a set of positions pos and returns a nested list with posititions divided into sex of carrier (e.g. {(Subscript[i, {f,m}]),(Subscript[j, {m,m}])} becomes {{(Subscript[i, {f,m}])},{(Subscript[j, {m,m}])}}).
Trans applies the algorithm described above, over all sexes of carrier, using the helper function Trans1SOC to recursively account for all positions.  There is a possibility that the way this is set up will require some odd manipulation of the produced recombination rates (this should not be an issue for the two loci case).

PosSplitSOC=Function[pos,
Cases[pos,(Subscript[_, {#,_}])]&/@(#[[2,1]]&/@pos//Union)];
Trans=Function[pos,
Times@@Trans1SOC/@PosSplitSOC[pos]];
Trans1SOC = Function[x,
If[Length[x]<2,1/2,
(If[x[[1,2,2]]===x[[2,2,2]],1-r,r]/.r->(Subscript[r, {x[[1,1]],x[[2,1]]}]))*Trans1SOC[Rest[x]]]];

For example, consider linkage disequilibrium ((Subscript[D, {(Subscript[i, f]),(Subscript[j, f])}])) vs trans-linkage disequilibrium ((Subscript[D, {(Subscript[i, f]),(Subscript[j, m])}])).  The transmission coefficients are:

(Subscript[t, {{(Subscript[i, f]),(Subscript[j, f])}←#}])->Trans[#]&/@TransPos[{(Subscript[i, f]),(Subscript[j, f])}]//TableForm
(Subscript[t, {{(Subscript[i, f]),(Subscript[j, m])}←#}])->Trans[#]&/@TransPos[{(Subscript[i, f]),(Subscript[j, m])}]//TableForm

The second thing required to calculate are the associations (Subscript[D, 𝕌]).  If 𝕌 has one sex of carrier, this is just (Subscript[D, 𝕌]) (e.g. (Subscript[D, {(Subscript[i, ff])(Subscript[j, ff])}])=(Subscript[D, {(Subscript[i, mf]),(Subscript[j, mf])}])), but if 𝕌 has two sexes of carrier, this is not correct; this association is an association among different individuals.  Two sex of carrier associations are computed in a similar way to the one sex of carrier associations: summing over the product of genotype frequencies and ζ variables.  However, the genotype frequencies for these still require calculation.

With random mating, the genotype frequencies are independent among individuals; i.e. the genotype frequency across two sex of carrers is the product of the genotype frequencies in each sex of carrier.  The genotype frequency within each sex of carrier is calculated using GenotypeExpr, as above.  The function Association2SOC calculates these associations, using several smaller helper functions:
	Association is identical to Assocaitions (section 2.2 above), but it accepts an XXSubs argument, since the indicator variables will involve two sexes of carrier.
	Assoc2SOCGenotypes calculates the "genotype frequency" for a genotype present across two individuals (at positions pos, which contain both sex of carrier and sex of origin in their context), in terms of allele frequencies and associations with one sex of carrier.
	Assoc2SOCZeta generates ζ expressions from which to calculate the associations.
	Assoc2SOCSubs generates all possible combinations of X variables, for substituting.
	
Finally, the associations created by Association2SOC would contain sex-of-carrier information, but no association contains multiple sexes of carrier, so it can be removed.  The function OneSOCD can be applied to any expression to strip the sex of carrier information; this is the last step before returning from Association2SOC.

Association2SOC will generate the correct associations even with one sex of carrier, but to avoid needless work, it short-cuts and just creates (Subscript[D, 𝕌]).

Association2SOC=Function[pos,
If[Length[PosSplitSOC[pos]]<2,(Subscript[D, pos]),
Association[Assoc2SOCGenotypes[pos],Assoc2SOCZeta[pos],Assoc2SOCSubs[pos]]]//OneSOCD//ReSortD];
Association=Function[{genotypes,zetas,XXSubs},
Total[genotypes*(zetas/.XXSubs)]//Factor];
Assoc2SOCGenotypes=Function[pos,Times @@ GenotypeExpr /@ PosSplitSOC[pos]/.Assoc2SOCSubs[pos]];
Assoc2SOCZeta=Function[pos,Times@@pos/. (Subscript[loc_, con_]):>(Subscript[X, (Subscript[loc, con])])-(Subscript[p, loc])];
Assoc2SOCSubs=Function[pos,Thread[Thread[(Subscript[X, pos])]->Transpose[Combinations[Length[pos]]]]];
OneSOCD=Function[x,x//.(Subscript[D, {𝕚___,(Subscript[j_, {SOC_,SOO_}]),𝕜___}]):>(Subscript[D, {𝕚,(Subscript[j, SOO]),𝕜}])];

Putting these together gives the transmission function, Transmission:

Transmission=Function[pos,Total[(Trans[#1] Association2SOC[#1]&)/@TransPos[pos]]];

Using this, it is straightforward to generate DppAssoc; the associations after transmission (note that since there is only one recombination rate, I'm substituting r back in for (Subscript[r, ij]))

DppAssoc=Transmission/@Subsets[positions]/.(Subscript[r, {i,j}])->r//Simplify;
Thread[Ds->DppAssoc]//TableForm

... Transmission of recalculated associations

If recalculation is done after selection, but before transmission, then the D values that appear in DppAssoc are (D^(′* )) associations: substitute in these values, for the epistasis model (DpsEpiSubs):

DpppSubs=Thread[Ds->(DppAssoc/.DpsEpiSubs)];
DpppSubs//TableForm

This is consistent with the genotype-based approach:
	Zero position association (Subscript[D, ∅]) stays 1
	One-position associations (Subscript[D, 𝕚]) stay zero
	Two-position and three position associations with two sexes of carrier are zero.
This leaves just the two one sex of carrier associations (linkage disequilibrium) and the four locus assocaition.  The linkage disequilibrium measures ((Subscript[D, {(Subscript[i, f]),(Subscript[j, f])}])and (Subscript[D, {(Subscript[i, m]),(Subscript[j, m])}])) are the same as those calculated with the genotype approach:

((Subscript[D, {(Subscript[i, f]),(Subscript[j, f])}])/.DpppSubs)==((Subscript[D, {(Subscript[i, m]),(Subscript[j, m])}])/.DpppSubs)==((Subscript[D, ij])/.recGenotypeEpi)//Simplify

As is the four-position association; this is equal to ((((Subscript[D, ij])+(Subscript[ΔD, ij])))^2) (see Dppp, section 5.1):

(Subscript[D, positions])==((((Subscript[D, {(Subscript[i, f]),(Subscript[j, f])}])))^2)/.DpppSubs//Simplify

This can be seen more clearly in the general expressions for D'', above:

(Subscript[D, positions])==((((Subscript[D, {(Subscript[i, f]),(Subscript[j, f])}])))^2)/.Thread[Ds->DppAssoc]//Simplify

... Transmission of non-recalculated associations

These do not immediately yield the final associations expected, so these must be compared against Dpp, calculated above (associations among transmitted positions, in terms of the old genotype associations).  To keep things simple above, the post transmission associations D'' (Dpp) were expressed in terms of changes in allele frequencies and changes in the linkage disequilibrium:

Thread[Ds->Dpp]//TableForm

Substitute in values for (Subscript[Δp, i]), (Subscript[Δp, j]) and (Subscript[ΔD, ij]) (by manipulation of recGenotypeEpi):

DppEpi=Dpp/.Thread[{(Subscript[Δp, i]),(Subscript[Δp, j]),(Subscript[ΔD, ij])}->(#[[2]]-#[[1]]&/@recGenotypeEpi)];

Substitute in assocaitions after selection (D', DprimeEpiSubs) into DppAssoc:

DppEpiAssoc=DppAssoc/.DprimeEpiSubs;

And compare these, showing that they are the same:

DppEpi==DppEpiAssoc//Simplify

Finally, recalculate the post-selection, post-transmission associations, and show that these are the the same as the associations that were recalculated before transmission, and the associations calculated from the genotype-based approach:

(RecalcD/@Subsets[positions])/.Thread[Ds->Dpp]//Simplify
%==Dppp