Load stuff

library(lattice)
library(cowplot)
library(pbapply)
library(tidyverse)
library(gridExtra)
library(viridis)

Model

Let \(x\) be the frequency of ab10 and \(1-x\) the frequency of n10 in a panmictic population of infinite size. In the absence of drive, the frequency of ab10 in the next generation is simply \(x'=x^2+\frac{1}{2}*2x(1-x)\), as \(\frac{1}{2}\) of the alleles produced by heterozygotes will be ab10.

We then introduce drive, such that heterozygotes for ab10 produce ab10:n10 in the ration \(d\):\((1-d)\) during female meioisis.

Note this means we have to keep track of frequencies in pollen (\(p\)) and ovules (\(o\)) separately.

Finally, we introduce selection in females (\(s\)) and dominance (\(h\)), such that n10 homozygotes have relative fitness of \(1\), heterozygotes have fitness \(1-hs\), and ab10 homozygotes have fitness of \(1-s\). On top of this viability, selection against male pollen competition acts similarly with strength \(t\) and dominance \(h\), such that population mean fitness of females (\(\bar{w}_{f}\)) and males (\(\bar{w}_{m}\)) are:

\(\bar{w}_{f}=x^2(1-s)+2x(1-x)(1-hs)+(1-x)^2\)

\(\bar{w}_{m}=x^2(1-t)+2x(1-x)(1-ht)+(1-x)^2\)

We assume a starting population with genotype frequencies \(X=(1-x)^2\), \(Y=2x(1-x)\), and \(Z=(1-x)^2\) of the n10 homozygote, the heterozygote, and the ab10 homozygote, respectively.

Pollen

Assuming no drive in males, and current genotpye frequencies \(X,Y,Z\) the frequency in pollen contributing to the next generation will be:

\(p = \frac{Z(1-t)}{\bar{w}_m} + \frac{Y(1-ht)}{2\bar{w}_m}\)

Ovules

With drive in females, and current frequency \(o\) of ab10 in ovules, the frequency in ovules contributing to the next generation will be:

\(o = \frac{Z(1-s)}{\bar{w}_f} + \frac{Yd(1-hs)}{\bar{w}_f}\)

Equilibrium

After a bit of algebra, and solving for when the frequency in the next generation \(x'=x\) I get:

\(\frac{x(1-t)+(1-x)(1-ht)}{2\Big(x^2(1-t)+(1-x)^2+2x(1-x)(1-ht)\Big)}+\frac{x(1-s)+2(1-x)d(1-sh)}{2\Big(x^2(1-s)+(1-x)^2+2x(1-x)(1-hs)\Big)}-1=0\)

Still doesn’t help much, as we’ve got five variables. So I’ll just evaluate this numerically. We’ll assume drive, so \(d>0.5\), and that observed frequency is \(0.05<x<0.25\).

d=600:800/1000
s=5:10/10
h=c(0:500/1000)
t=0.06
x=(10:25)/100
bob=expand.grid(d,s,h,t,x)
colnames(bob)=c("d","s","h","t","x")
bob<-mutate(bob,g=1/2*(2*(h*s - 1)*d*(x - 1) - (s - 1)*x)/(2*(h*s - 1)*(x - 1)*x - (s - 1)*x^2 + (x - 1)^2) + 1/2*((h*t - 1)*(x - 1) - (t - 1)*x)/(2*(h*t - 1)*(x - 1)*x - (t - 1)*x^2 + (x - 1)^2) - 1)
## Warning: package 'bindrcpp' was built under R version 3.3.2
sue=filter(bob,abs(bob$g)<0.0001)
ggplot(sue,aes(h,d,color=s))+geom_point()

Simulate

Next try: let’s simulate for 1000 generations and track final allele frequency.

Function

Simulate

Graph