We take an individual which is heterozygous for 1000 loci. We represent this as a vector of 1000 loci with allele frequency 0.5

nloci=1000
sample=rep(0.5,nloci)
head(sample)
## [1] 0.5 0.5 0.5 0.5 0.5 0.5

If we inbred this individual, assuming loci are unlinked, this is the same as bernouli draws from an allele frequency from the previous generation

gen={}
gen[[1]]=rbinom(nloci,2,sample)/2
head(gen[[1]])
## [1] 0.0 0.0 0.5 0.5 0.5 0.5

We see that \(\sim 50\%\) of the loci remain heterozygous, as expected.

sum(gen[[1]]==0.5)
## [1] 502

For future generations, we use the allele frequencies of the previous generation’s parent.

gen[[2]]=rbinom(nloci,2,gen[[1]])/2

And again, we lose half.

sum(gen[[2]]==0.5)
## [1] 250

If we assume all 1000 of these loci are weakly deleterious and completely recessive, then we can just sum up the number of loci not homyzygous as an estimate of fitness. So our intial sample has a fitness of 1000.

fit=sum(sample!=1)
fit
## [1] 1000

And fitness of course declines by half every generation, approaching the expectation of 50%.

gen[[1]]=rbinom(nloci,2,sample)/2
gen[[2]]=rbinom(nloci,2,gen[[1]])/2
gen[[3]]=rbinom(nloci,2,gen[[2]])/2
gen[[4]]=rbinom(nloci,2,gen[[3]])/2
gen[[5]]=rbinom(nloci,2,gen[[4]])/2
gen[[6]]=rbinom(nloci,2,gen[[5]])/2
gen[[7]]=rbinom(nloci,2,gen[[6]])/2
fit=sapply(1:7, function(g) sum(gen[[g]]!=1))

plot(fit~c(1:7),xlab="generation",ylab="fitness",pch=19)

Now let’s do the same with a tetraploid. Same deal, but with sampling 4 alleles, so each generation after the initial heteorzygous sample can have allele freqeuncies of 0,0.25,0.5,0.75, and 1. We see that in the first few generations, inbreeding depression is much worse in the diploid, but by generation 7 we expect to see less of a difference.

gen4={}
gen4[[1]]=rbinom(nloci,4,sample)/4
head(gen4[[1]])
## [1] 0.75 0.50 0.75 0.50 0.50 0.50
gen4[[2]]=rbinom(nloci,4,gen4[[1]])/4
gen4[[3]]=rbinom(nloci,4,gen4[[2]])/4
gen4[[4]]=rbinom(nloci,4,gen4[[3]])/4
gen4[[5]]=rbinom(nloci,4,gen4[[4]])/4
gen4[[6]]=rbinom(nloci,4,gen4[[5]])/4
gen4[[7]]=rbinom(nloci,4,gen4[[6]])/4
fit4=sapply(1:7, function(g) sum(gen4[[g]]!=1))

plot(fit~c(1:7),xlab="generation",ylab="fitness",pch=19,ylim=c(450,1000))
points(fit4~c(1:7),col="chartreuse4",pch=19)
legend("topright",fill=c("black","chartreuse4"),legend=c("diploid","tetraploid"))

Now let’s try with partially recessive alleles. Mendel doesn’t change, but the impact on fitness does. We assume alleles come from some distribution of values between \(h=0\) (fully recessive) and \(h=0.5\) (additive). We’ll pick something that gives us a lot of strongly recessive alleles to start with.

curve(from=0,to=0.5,dbeta(x,shape1=0.25,shape2=1)/2,ylab="density",xlab="dominance (h)")

Now, in diploids, with a fitness model where at each locus an individual’s fitnes of the 0 (homozygous ancestral), 1 (heterozygous), and 2 (homozygous for the deleterious allele) genotype is written as 1, \(1-hs\), and, \(1-s\) where \(s\) is the selection coefficient, we can calculate new fitnesses.

h=rbeta(1000,0.25,1)/2
pgen={}
pgen[[1]]=rbinom(nloci,2,sample)/2
pgen[[2]]=rbinom(nloci,2,pgen[[1]])/2
pgen[[3]]=rbinom(nloci,2,pgen[[2]])/2
pgen[[4]]=rbinom(nloci,2,pgen[[3]])/2
pgen[[5]]=rbinom(nloci,2,pgen[[4]])/2
pgen[[6]]=rbinom(nloci,2,pgen[[5]])/2
pgen[[7]]=rbinom(nloci,2,pgen[[6]])/2
pfit=sapply(1:7, function(g) 1000-sum(pgen[[g]]==1)-sum(pgen[[g]][pgen[[g]]==0.5]*h[pgen[[g]]==0.5]))

plot(pfit~c(1:7),xlab="generation",ylab="fitness",pch=19,ylim=c(450,850),col="gray")
points(fit~c(1:7),pch=19)
legend("topright",fill=c("black","gray"),legend=c("complete recessive","partial recessive"))

Note that if we set \(h=0\) all alleles are completely recessive and we see no real difference from the previous case, as expected.

h=rep(0,1000)
pgen={}
pgen[[1]]=rbinom(nloci,2,sample)/2
pgen[[2]]=rbinom(nloci,2,pgen[[1]])/2
pgen[[3]]=rbinom(nloci,2,pgen[[2]])/2
pgen[[4]]=rbinom(nloci,2,pgen[[3]])/2
pgen[[5]]=rbinom(nloci,2,pgen[[4]])/2
pgen[[6]]=rbinom(nloci,2,pgen[[5]])/2
pgen[[7]]=rbinom(nloci,2,pgen[[6]])/2
pfit=sapply(1:7, function(g) 1000-sum(pgen[[g]]==1)-sum(pgen[[g]][pgen[[g]]==0.5]*h[pgen[[g]]==0.5]))

plot(pfit~c(1:7),xlab="generation",ylab="fitness",pch=19,ylim=c(450,850),col="gray")
points(fit~c(1:7),pch=19)
legend("topright",fill=c("black","gray"),legend=c("complete recessive","partial recessive"))

Now we do the same for tetraploids. Only here, additive is \(h=0.25\) so our partially recessive distribution changes. Now we see no difference between tetraploids and diploids.

h=rbeta(1000,0.25,1)/4
pgen4={}
pgen4[[1]]=rbinom(nloci,4,sample)/4
head(gen4[[1]])
## [1] 0.75 0.50 0.75 0.50 0.50 0.50
pgen4[[2]]=rbinom(nloci,4,pgen4[[1]])/4
pgen4[[3]]=rbinom(nloci,4,pgen4[[2]])/4
pgen4[[4]]=rbinom(nloci,4,pgen4[[3]])/4
pgen4[[5]]=rbinom(nloci,4,pgen4[[4]])/4
pgen4[[6]]=rbinom(nloci,4,pgen4[[5]])/4
pgen4[[7]]=rbinom(nloci,4,pgen4[[6]])/4
pfit4=sapply(1:7, function(g) 1000-sum(pgen[[g]]==1)-sum(pgen[[g]][pgen[[g]]==0.25]*h[pgen[[g]]==0.25])-2*sum(pgen[[g]][pgen[[g]]==0.5]*h[pgen[[g]]==0.5])-3*sum(pgen[[g]][pgen[[g]]==0.75]*h[pgen[[g]]==0.75]))

plot(pfit4~jitter(c(1:7),.5),xlab="generation",ylab="fitness",pch=19,ylim=c(450,1000),col="chartreuse")
points(pfit~jitter(c(1:7),.5),pch=19,col="gray")
points(fit~jitter(c(1:7),.5),pch=19,col="black")
points(fit4~jitter(c(1:7),.5),pch=19,col="chartreuse4")

legend("topright",fill=c("chartreuse4","chartreuse","black","gray"),legend=c("4X recessive","4X partial recessive","2X recessive","2X partial recessive"))

And this is independent of the details of the distribution of \(h\). First a new distribution:

curve(from=0,to=0.5,dbeta(x,shape1=1,shape2=1)/2,ylab="density",xlab="dominance (h)")

And now repeat calcs.

h=rbeta(1000,1,1)/2

pgen={}
pgen[[1]]=rbinom(nloci,2,sample)/2
pgen[[2]]=rbinom(nloci,2,pgen[[1]])/2
pgen[[3]]=rbinom(nloci,2,pgen[[2]])/2
pgen[[4]]=rbinom(nloci,2,pgen[[3]])/2
pgen[[5]]=rbinom(nloci,2,pgen[[4]])/2
pgen[[6]]=rbinom(nloci,2,pgen[[5]])/2
pgen[[7]]=rbinom(nloci,2,pgen[[6]])/2
pfit=sapply(1:7, function(g) 1000-sum(pgen[[g]]==1)-sum(pgen[[g]][pgen[[g]]==0.5]*h[pgen[[g]]==0.5]))

h=rbeta(1000,1,1)/4
pgen4={}
pgen4[[1]]=rbinom(nloci,4,sample)/4
head(gen4[[1]])
## [1] 0.75 0.50 0.75 0.50 0.50 0.50
pgen4[[2]]=rbinom(nloci,4,pgen4[[1]])/4
pgen4[[3]]=rbinom(nloci,4,pgen4[[2]])/4
pgen4[[4]]=rbinom(nloci,4,pgen4[[3]])/4
pgen4[[5]]=rbinom(nloci,4,pgen4[[4]])/4
pgen4[[6]]=rbinom(nloci,4,pgen4[[5]])/4
pgen4[[7]]=rbinom(nloci,4,pgen4[[6]])/4
pfit4=sapply(1:7, function(g) 1000-sum(pgen[[g]]==1)-sum(pgen[[g]][pgen[[g]]==0.25]*h[pgen[[g]]==0.25])-2*sum(pgen[[g]][pgen[[g]]==0.5]*h[pgen[[g]]==0.5])-3*sum(pgen[[g]][pgen[[g]]==0.75]*h[pgen[[g]]==0.75]))

plot(pfit4~jitter(c(1:7),.5),xlab="generation",ylab="fitness",pch=19,ylim=c(450,1000),col="chartreuse")
points(pfit~jitter(c(1:7),.5),pch=19,col="gray")
points(fit~jitter(c(1:7),.5),pch=19,col="black")
points(fit4~jitter(c(1:7),.5),pch=19,col="chartreuse4")

legend("topright",fill=c("chartreuse4","chartreuse","black","gray"),legend=c("4X recessive","4X partial recessive","2X recessive","2X partial recessive"))