Gene drive model—writing to tagL0, sample weighting
Jack Sissons
initialize() {
initializeSLiMModelType("nonWF");
initializeSLiMOptions(dimensionality="xy");
defineConstant("K", 1000); // scaling a pop of 100,000 down to 1000 (factor = 100)
initializeMutationRate(1.3e-6); // scaled
initializeMutationType("m1", 0.5, "f", 0.0); // neutral
initializeMutationType("m2", 0.5, "f", -0.1); // Gene drive
initializeMutationType("m3", 0.5, "f", 0.0); //S-locus mutation
m1.convertToSubstitution = T;
m2.convertToSubstitution = F;
m3.convertToSubstitution = T;
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElementType("g2", m3, 1.0); // S-locus mutations
initializeGenomicElement(g1, 0, 20000);
initializeGenomicElement(g2, 20001, 21000); // S-locus
initializeGenomicElement(g1, 21001, 99999);
initializeRecombinationRate(5e-6); //scaled
// spatial mate choice
initializeInteractionType(1, "xy", reciprocal=T, maxDistance=300);
i1.setInteractionFunction("n", 1.0, 10.0);
}
1 first() {
sim.addSubpop("p1", 50);
}
1: early() {
sim.killIndividuals(p1.subsetIndividuals(minAge=9));
p1.fitnessScaling = K / p1.individualCount;
}
//set up the map
1 late() {
p1.setSpatialBounds(c(0.0, 0.0, 701.0, 241.0));
mapImage = Image("Nelson_marina_greyscale.png");
map = p1.defineSpatialMap("Nelson marina", "xy", mapImage.floatK, valueRange=c(0.0, 1.0), colors=c("darkturquoise", "black"));
defineConstant("marina", map);
// set intial positions
for (ind in p1.individuals) {
ind.y = rnorm(1, 120.0, 50.0);
ind.x = rnorm(1, 350.0, 100.0);
}
}
2: first() {
i1.evaluate(p1);
}
// label gene drive carriers in tagL0 *****
2: late() {
ind = p1.individuals;
ind.tagL0 = ind.haploidGenome1.containsMarkerMutation(m2, 1000);
p1.individuals.color = ifelse(p1.individuals.tagL0, "red", "green");
}
// if over 2 months old, produce 200 offspring with random nearby mates. If no nearby mates, self-fertilise.
reproduction() {
if (individual.age >= 2)
{
neighbours = i1.nearestNeighbors(individual, count=1000);
closeness = 300 - i1.distance(individual, neighbours);
if (size(neighbours) >= 1)
{
for (mate in sample(neighbours, 200, replace=T, weights=closeness)); // does it know how to assign the correct closeness value?
subpop.addCrossed(individual, mate);
}
else subpop.addSelfed(individual, count=100);
}
}
// S-locus incompatibility - section 11.3 of manual.
modifyChild() {
spermSMuts = child.haploidGenome2.mutationsOfType(m3);
eggSMuts1 = parent1.haploidGenome1.mutationsOfType(m3);
eggSMuts2 = parent1.haploidGenome2.mutationsOfType(m3);
tick = sim.cycle;
rejectChance = 0.99 - (1.01 ^ -tick);
// compares s-locus of the sperm & egg and if they are identical, prevents fertilisation in 99% of cases (ramping).
if (identical(spermSMuts, eggSMuts1))
if (runif(1) < rejectChance)
return F;
if (identical(spermSMuts, eggSMuts2))
if (runif(1) < rejectChance)
return F;
// keep children in bounds and on structures
do {
do pos = parent1.spatialPosition + rnorm(2, 0, 50);
while (!p1.pointInBounds(pos));
}
while (marina.mapValue(pos) == 0.0);
child.setSpatialPosition(pos);
return T;
}
late() {
if (sim.cycle % 100 == 0)
{
inds = p1.individuals;
cat(sim.cycle + ": " + size(inds));
catn(" (" + max(inds.age) + ", " + mean(inds.age) + ")");
}
}
// gene drive release
1000 late()
{
released = sample(p1.haplosomes, 25);
released.addNewDrawnMutation(m2, 1000);
}
// convert heterozygous individuals to homozygous
modifyChild() {
mut = sim.mutationsOfType(m2);
if (size(mut) == 1) {
hasMutOnChromosome1 = child.haploidGenome1.containsMutations(mut);
hasMutOnChromosome2 = child.haploidGenome2.containsMutations(mut);
if (hasMutOnChromosome1 & !hasMutOnChromosome2)
child.haploidGenome2.addMutations(mut);
else if (hasMutOnChromosome2 & !hasMutOnChromosome1)
child.haploidGenome1.addMutations(mut);
}
return T;
}
//// gene drive carriers cannot produce eggs *****
1001: modifyChild() {
if (parent1.tagL0 == T) {
return F;
}
return T;
}
// Rainfall mortality: ~1 event every 36 months
1000: early() {
if (runif(1) < 0.027) {
severity = rnorm(1, 0.5, 0.1);
inds = p1.individuals;
popSize = p1.individualCount;
affected = p1.sampleIndividuals(asInteger(severity * popSize));
affected.fitnessScaling = 0.0;
catn(sim.cycle + ": 🌧️ mortality = " + severity * 100 + "%");
}
}
5000 late() {
sim.simulationFinished();
}
Ben Haller
I'm very glad you found the workshop useful. There's a bit of a SLiM community in New Zealand, by the way, and I think there was a somewhat recent thread on slim-discuss from them, if I recall correctly; you might search for that, if you want to connect with others. Actually – inspired by this, I just set up a new area in the slim-community organization online where people can set up their own regional groups of SLiM users, so maybe you NZ folks can ask for a discussion category there if that interests you. :-> See: https://github.com/orgs/slim-community/discussions/7 .
OK, onward to your questions.
1. Well, this indicates some kind of a logic error in your script. When it says that "tagL0 is not yet defined", that means what it sounds like it means! :-> You need to set the tagL0 property, on any given individual, before you access it. That error is there precisely to help you find logic errors like this in your script. Typically you want to set the value for each new individual when you create a new subpopulation, and for each new offspring. Exactly where you do that is up to you, as long as it gets set before you try to access its value. If you're puzzled because you think you're already doing this, you might debug the problem by adding a catn() call in each place where you set the tagL0 values and each place where you read them. You could turn on pedigree tracking with keepPedigrees=T, in initializeSLiMOptions(), and print the pedigreeID of the individual and the tag value you're setting/accessing. The output from your script will then show definitively that you're not setting the value before you try to use it – or, perhaps, will show that SLiM has a bug in this area, but I'd be fairly surprised since this facility does get tested by unit tests, as I recall. :->
2. If I understand correctly, you're using a call to sample() to choose a parent2 for a given parent1, out of a pool of candidate parent2 individuals. You would then pass all the parent2 individuals to sample(), and pass a vector of weights with the same length as that pool of candidates, providing the weight for each candidate. The sample() method would then use those weights. I think if you passed a vector of weights that didn't match the length of the candidate vector, sample() would throw an error. Does that answer your question? See the Eidos manual for the reference documentation on sample(); the reference doc is always the place to go to understand exactly what a given function/method does. But more broadly, note that choosing nearby neighbors of a focal individual based upon distances is exactly what the drawByStrength() method of InteractionType() is designed to do, so you might look at simply using that. Section 17.17 has an example of its use.
I hope this helps! I didn't look at your model, since it seemed like I didn't need to to answer your questions, but please follow up if you have more questions or I've misunderstood you. Happy modeling!
Cheers,
-B.
Benjamin C. Haller
Messer Lab
Cornell University
Jack Sissons wrote on 9/10/25 10:03 PM:
--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/d8dcf3d3-f517-4ccd-91c6-c78d5a767fd6n%40googlegroups.com.