|
| 1 | +// Keywords: monogamy, monogamous mating, nonWF, non-Wright-Fisher |
| 2 | + |
| 3 | +initialize() { |
| 4 | + defineConstant("K", 500); // carrying capacity |
| 5 | + defineConstant("R_AGE_M", 3); // minimum age of reproduction (male) |
| 6 | + defineConstant("R_AGE_F", 4); // minimum age of reproduction (female) |
| 7 | + defineConstant("FECUN", 0.2); // mean fecundity per female per tick |
| 8 | + |
| 9 | + initializeSLiMModelType("nonWF"); |
| 10 | + initializeSex(); |
| 11 | +} |
| 12 | +1 first() { |
| 13 | + sim.addSubpop("p1", K); |
| 14 | + p1.individuals.age = rdunif(K, min=0, max=15); // initial variation in age |
| 15 | + p1.individuals.tag = -1; // mark all individuals as unmated |
| 16 | +} |
| 17 | +first() { |
| 18 | + // find mated individuals whose mate has died, and mark them as unmated |
| 19 | + mated_individuals = p1.individuals; |
| 20 | + mated_individuals = mated_individuals[mated_individuals.tag >= 0]; |
| 21 | + |
| 22 | + if (size(mated_individuals) > 0) |
| 23 | + { |
| 24 | + tags = mated_individuals.tag; |
| 25 | + tag_counts = tabulate(tags, maxbin=size(mated_individuals)-1); |
| 26 | + tags_to_fix = which(tag_counts == 1); |
| 27 | + unmated_indices = match(tags_to_fix, tags); |
| 28 | + mated_individuals[unmated_indices].tag = -1; |
| 29 | + } |
| 30 | + |
| 31 | + // find the next tag value to use for new mating pairs |
| 32 | + next_tag = max(p1.individuals.tag) + 1; |
| 33 | + |
| 34 | + // find unmated individuals that are of reproductive age |
| 35 | + unmated_F = p1.subsetIndividuals(sex="F", tag=-1, minAge=R_AGE_F); |
| 36 | + unmated_M = p1.subsetIndividuals(sex="M", tag=-1, minAge=R_AGE_M); |
| 37 | + |
| 38 | + // pair individuals randomly; some individuals may be left unpaired |
| 39 | + pair_count = min(size(unmated_F), size(unmated_M)); |
| 40 | + unmated_F = sample(unmated_F, pair_count, replace=F); |
| 41 | + unmated_M = sample(unmated_M, pair_count, replace=F); |
| 42 | + |
| 43 | + for (f in unmated_F, m in unmated_M, tag in seqLen(pair_count) + next_tag) |
| 44 | + { |
| 45 | + f.tag = tag; |
| 46 | + m.tag = tag; |
| 47 | + } |
| 48 | +} |
| 49 | +reproduction() { |
| 50 | + // find the subset of individuals that have a mate |
| 51 | + mated_F = p1.subsetIndividuals(sex="F"); |
| 52 | + mated_F = mated_F[mated_F.tag >= 0]; |
| 53 | + |
| 54 | + mated_M = p1.subsetIndividuals(sex="M"); |
| 55 | + mated_M = mated_M[mated_M.tag >= 0]; |
| 56 | + |
| 57 | + // look up the male for each female, by tag |
| 58 | + male_indices = match(mated_F.tag, mated_M.tag); |
| 59 | + mated_M = mated_M[male_indices]; |
| 60 | + |
| 61 | + pair_count = size(mated_F); |
| 62 | + |
| 63 | + // produce offspring from each mated pair |
| 64 | + for (f in mated_F, |
| 65 | + m in mated_M, |
| 66 | + c in rpois(pair_count, FECUN), |
| 67 | + new_tag in seqLen(pair_count)) |
| 68 | + { |
| 69 | + // re-tag paired individuals to compact tags down |
| 70 | + f.tag = new_tag; |
| 71 | + m.tag = new_tag; |
| 72 | + |
| 73 | + offspring = p1.addCrossed(f, m, count=c); |
| 74 | + offspring.tag = -1; // mark offspring as unmated |
| 75 | + } |
| 76 | + |
| 77 | + self.active = 0; // deactivate for the rest of the tick ("big bang") |
| 78 | +} |
| 79 | +early() { |
| 80 | + // density-dependent population regulation |
| 81 | + p1.fitnessScaling = K / p1.individualCount; |
| 82 | +} |
| 83 | +10000 late() { } |
0 commit comments