That’s really cool! I’d be really curious to see the lesson plans for this. Perhaps we could publish them at PS.
The software from the Nature 2004 paper is not available. I’ve written some software of my own but its not in any shape to share at this time.
My suggestion is to use SLIM.
They have an excellent forward simulator and an active community. While it doesn’t usually track genealogical ancestry, it might be flexible enough configure for that purpose. I’d check with their discussion group.
The other option is to hold off from a full fledged simulator, and simulate focused questions. For example,
In a well mixed population of size N, how many generations before a randomly chosen individual becomes a universal ancestor?
That’s pretty easy to code up in python. Here it is in pseudocode:
import numpy as np
N = 10000
generation = np.zeros(N)
generation[0] = 1
for x in xrange(num_generations):
men = generation[N/2:]
father = (np.rand(N) * N/2).floor()
father = np.take(men, father)
women = generation[:N/2]
mother = (np.rand(N) * N/2).floor()
mother = np.take(women, mother)
generation = mother + father
if np.alltrue(generation > 0 ):
print("Generations to universal ancestor: ", x+1)
break
if np.alltrue(generation == 0):
print("Generations till extinction: ", x +1)
break
I’m calling this pseudocode because I haven’t tested it yet, but something very close to this will work.
You can have them compare that with the theoretical expectation (it is very close). They can also measure the variance there is between multiple runs (it is exceedingly low). You can also ask what chance there is of extinction of a lineage (very close to the theoretical value). You can also ask how this adjusts if you change the mating algorithm (usually not much affect at all). You can ask what the distribution of times is for multiple runs (it will look like a truncated and skewed gaussian).
One hint is that computing the most recent universal ancestor of a given population is much harder. That requires you to track every signal person in the whole simulation. I do NOT recommend this as it isn’t the most interesting quantity and it is technically challenging to compute. Instead, the average time to universal ancestry is a better place to focus.
Perhaps @evograd, @davecarlson or @Joe_Felsenstein has some suggestions I don’t know about.