After reading this post, I wrote some code to simulate inbreeding.
We have a population of $N$ creatures. Each creature has two genes, which come in two forms: recessive (a) and dominant (A). The creatures' initial genomes are random. At each stage, the following happens, in order:
- Each creature with the recessive phenotype (aa) dies with probability $p_R$.
- Each creature with the dominant phenotype (AA, aA, or Aa) dies with probability $p_D<p_R$.
- Each creature then decides to have children with probability $q$. Each creature that decides to have children chooses a mate (see below) and has three children with that mate.
I wrote two different kinds of creature, incest-averse and incest-seking. Each time I ran the simulation, the population consisted of either all incest-averse or all invest-seeking creatures.
- An incest-averse creature chooses a mate randomly and uniformly from the set of its non-siblings.
- An incest-seeking creature chooses a mate randomly and uniformly from the set of its siblings. If it has no siblings, it just randomly picks anybody.
In practice, since each union produces 3 children, the incest-seeking creatures almost never fail to have siblings after a few generations.
What I found was that whichever type of creature I used, the results were the same: the frequency of the aa genotype plummeted to almost zero. I used the parameters $N=100$, $p_R=0.1$ $p_D=0.05$ $q=0.022$.
Presumably my model is in some way too simple. What am I not taking into account?