Alternation of Generations Model - trouble with subpopulations

[elissasoroj]

Hi again,

I am using an alternation of generations model in SLiM, which passes individuals between two subpopulations (p0 and p1 in my case). In generation 1 diploid sporophytes live in p1 and produce haploid gametophytes directly into p0, then the sporophytes die. In generation 2 the gametophytes then reproduce and place new diploid sporophytes back into p1 then the gametophytes die.

This seems to cause problems at the recapitation step in pyslim, because half the nodes are assigned to Population 0 and half are assigned to Population 1, but they really should all be treated as the same population by msprime.
I was not able to alter the population data in the tables, as I I got this error: raise FrozenInstanceError().
I also had no luck setting a migration matrix of [[0, 1], [1, 0] or other migration rates of < 1 for recapitation.

Any suggestions for how to proceed would be greatly appreciated.

Thanks!
~Elissa

[petrelharp]

Good question! Let's see - I imagine you're doing something like is described here:
https://tskit.dev/pyslim/docs/latest/tutorial.html#recapitation-with-more-than-one-population
... is that right? I'd like to see what code you're running to get those errors. And, it sounds like what you'd actually like to do is to merge the two populations into one for the recapitation step, which we could do with an msprime population_split - I'd be happy to suggest how to do this, but it'd help to see what you're doing, first?

[petrelharp]

Very nice model! I agree, definitely no-one is surviving more than one time step, and after thinking through the reproduction callbacks I don't see anything wrong. Now I wonder if you're seeing a selective sweep? I agree that the tree you sent this time looks a bit odd, but when I run this SLiM script and look at the resulting trees I don't see anything like what you're seeing. I'd need a reproducible example (like, exactly the same code that you ran to produce the figure you're seeing).

A small note - the second half of the s2 block (with maternal fitness effects) never runs, if I'm reading it correctly. Tangentially, instead of returning F from the alternating survival callbacks, you could set their fitnessScalings to 0 - e.g., line 67 could have

		p0.fitnessScaling = (hK / p0.individualCount);
		p1.fitnessScaling = 0;

[elissasoroj]

Thanks! I’ve attached a SLiM script with a proper chromosome (sorry I didn’t provide earlier, it’s just so unwieldy). It generated this ts (full and then sample of 20 ends):

I guess it could be a selective sweep but it seems too consistent between runs.

The scripts are kind of weird because they are generated in Python. If you’re wondering what the heck is going on with the genomic elements - I’ve written it so every third base in a coding region is empty, so we are simulating a proxy of synonymous and non-synonymous sites. And you’re right, although maternal fitness effects and random death are turned off in this particular model (because the constants are set to 0), they are still in there because they are options in the Python program that generated the script. Additionally, this is one of several different default models we have that simulate plant life cycles. The Python program that generates the SLiM script can use any of the default models, so although it would definitely be more efficient to use fitnessScaling, it is separated out in a survival() callback because we need to toggle it on and off for different models.

link to SLiM script

*sorry for the double post

[petrelharp]

Ah-ha: notice that all the uncoalesced lineages are even, and the ones that coalesce before recapitation are odd. Loading in the output from your slim script shows the same thing except the other way around (even ones coalesce):

ts = pyslim.load("bryo_branch4.trees")
t = ts.first()
for i in range(10):
    for n in ts.individual(i).nodes:
        print(n, t.parent(n))

gives

0 2492
1 -1
2 2184
3 -1
4 2184
5 -1
6 2511
7 -1
8 2813
9 -1
# etcetera

Since adjacent nodes are the two chromosomes of a single individual, it looks to me like your reproduction callbacks are making it so that the second chromosomes of each individual never coalesce. And, looking into things - it looks like this has something to do with addRecombinant for sporophytes adding in "empty" genomes. At first I thought that those weird branches are just the "empty" genomes, and so they'd not appear if I wrote out the tree sequence at an odd time instead of an even time, but that did not fix the problem.

So, I'm not quite sure where the problem is. If the problem is that you just need to ignore the second chromosomes of each individual, then simplify( ) will take care of that for you (but my quick try at that didn't fix the problem). Happily you don't need your big complicated script to debug this - you can make a very small example with like 5 individuals that runs for 4 generations (or something); if it were me I would (a) ask SLiM to please print out all parent-child relationships; and (b) remember everyone in every generation, so then I could compare what's in the tree sequence to what SLiM prints out to understand what's going on.

Let me know if the problem isn't clear.

[petrelharp]

Ah-ha, again: setting the population sizes to 10 and replacing the last blocks in your slim script with:

6 late() { // executes after selection occurs
	sim.treeSeqRememberIndividuals(sim.subpopulations.individuals)
		;
	sim.treeSeqOutput('test6.trees');
}
7 late() { // executes after selection occurs
	sim.treeSeqRememberIndividuals(sim.subpopulations.individuals)
		;
	sim.treeSeqOutput('test7.trees');
}

and then doing:

ts7 = pyslim.load("test7.trees")
ts6 = pyslim.load("test6.trees")
for ts in (ts6, ts7):
  for i in ts.individuals_alive_at(0):
    ind = ts.individual(i)
    print(ind.id, ind.metadata['subpopulation'])
  print("======")

gets me

0 0
1 0
2 0
3 0
4 0
=====
0 1
1 1
=====

which is telling me that in even time steps, you have only "haploid" gametes alive; so your weird branches are from those "empty" genomes; while in odd time steps you have diploids. The easy solution is to just output the tree sequence at an odd time step.

[elissasoroj]

Ah, of course! I had forgotten about the empty genomes in the haploids. Thank you very much for all your help - with the original question and all my follow up questions!