MDM: describe behaviour when a pulse occurs during a period of continuous migration

[grahamgower]

The migration section of the population genetics model has this to say (thanks @apragsdale!):

Continuous migration rates are defined as the probability that parents in the "destination"
population are chosen from the "source" population. Migration rates are thus
per generation and must be less than one. Furthermore, if more than one source
population have continuous migration into the same destination population, the
sum of those migration rates must also be less than one, as rates define
probabilities. The probability that parents come from the same population is
just one minus the sum of incoming migration rates. On the other hand, pulse
migration events specify the instantaneous replacement of a given fraction of
individuals in a destination population by individuals with parents from
a source population. The fraction must be between zero and one, and if more
than one pulse occurs at the same time, those replacement events are applied
sequentially in the order that they are specified in the model.

So it's plausible that a simulator could support a model in which the sum of continuous migrations into a deme in a single generation is x < 1, and the sum of pulse migrations is y < 1, but x + y > 1. In particular, we can support this scenario in a continuous-time simulator like msprime. But for many simulators the pulse and continuous migrations will both be implemented using a single function or data structure, and in a discrete-time setting it may not be possible to decouple pulses and continuous migrations such that one occurs after the other. In this latter case, we must also have x + y <= 1.

Thoughts?

[jeromekelleher]

Is it biologically plausible to have x + y > 1? If not, I think we should rule it out.

[apragsdale]

The way I handled it in fwdpy11 (and maybe @molpopgen has thoughts about this as well, working with discrete generation models), is that when continuous migration occurs over a time period that includes a pulse, the continuous migration probabilities define the probability of choosing parents from each deme conditional on individuals not arriving via the pulse.

E.g. if deme 0 has continuous migration from deme 1 and deme 2 (say, rate 0.01 each), and then deme 3 contributes 50% in a generation, the migration row for deme zero in the non-pulse generations is [0.98, 0.01, 0.01, 0], and in the generation of the pulse, that row is [0.49, 0.005, 0.005, 0.5]. After the pulse, the migration matrix row changes back to just have continuous migration.

This way, I didn't have to worry about the cases where x + y > 1, and it relieves the need to do a bunch of verification across migration types.

[grahamgower]

That's nice and simple, I like it. It's equivalent to applying the pulse immediately after (more recently than) the migration, right? So it ought to match the current behaviour of msprime.Demography.from_demes().

[apragsdale]

It's necessarily concurrent in that generation, since it's all handled by adjusting the migration matrix. For example, if each of the other three populations provide migrants at rate 0.01, in that same example the mig matrix row without and with the pulse is [0.97, 0.01, 0.01, 0.01] and [0.485, 0.005, 0.005, 0.505], resp.

[sgravel]

I have a related question about the specification for pulse migrations at non-integer time values. Consider the model

time_units: generations
demes:
  - name: X
    epochs:
      - end_time: 1000
        start_size: 2000
  - name: A
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: B
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: C
    ancestors: [X]
    epochs:
      - start_size: 2000

pulses:
  - source: A
    dest: B
    proportion: 0.3
    time: 10
  - source: C
    dest: B
    proportion: 0.3
    time: 10.1

The possibly ambiguous value is the time: 10.1 in the last line.

As I understand things, if a continuous model such as Kingman’s coalescent, a pulse from A at generation tA=10 and a pulse from C at generation tC =10.1 are straightforward to interpret – they occur sequentially. There is some ambiguity as to what happens if two pulses occur at the exact same time, and the earlier suggestion was to treat them in order in which they appear in the text file. This is fine, but I don't love that solution as otherwise the order of things does not matter as I understand things, and it would be slightly preferable I think to keep it this way.

If such a model is passed on to a Wright-Fisher model, however, things are not so simple and I did not see an unambiguous specification. It may be that you want to simply forbid continuous times for discrete simulators, but of course this is limiting.

For example, TRACTS implements a discrete time model where recent migration pulses from different source populations can co-occur in a given generation, and fractional generations matter.

If I chose selfishly, I would like a specification such that:

1-The Wright-Fisher and Kingman models are comparable if simulated from the same model according to specification (i.e., at the end of generation 10, both models have the same proportion of lineages tracing back to sources A and C).

2-Both models work and are somewhat continuous in the continuous time parameter, if they accept continuous times.

In TRACTS, I had handled this issue by interpreting a pulse at a non-integer time t as a pulse spread over two generation, ceil(t) and floor(t), with proportions at each generation computed such that proportion of lineages that can be traced back to the pulse are as required by the pulse proportion.

However, this would not give directly comparable results to Kingman’s in the model above. To achieve consistency, one would need to similarly split pulses under Kingman to achieve compatibility. This is, admittedly, quite ugly.

A (small) benefit is that you would get back continuity in the time parameter, so that the following two models would give similar results:

time_units: generations
demes:
  - name: X
    epochs:
      - end_time: 1000
        start_size: 2000
  - name: A
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: B
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: C
    ancestors: [X]
    epochs:
      - start_size: 2000

pulses:
  - source: A
    dest: B
    proportion: 0.3
    time: 10.09
  - source: C
    dest: B
    proportion: 0.3
    time: 10.1

would be similar to

time_units: generations
demes:
  - name: X
    epochs:
      - end_time: 1000
        start_size: 2000
  - name: A
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: B
    ancestors: [X]
    epochs:
      - start_size: 2000
  - name: C
    ancestors: [X]
    epochs:
      - start_size: 2000

pulses:
  - source: A
    dest: B
    proportion: 0.3
    time: 10.1
  - source: C
    dest: B
    proportion: 0.3
    time: 10.09

Any thoughts on this? I don't really expect that you will want to change the Kingman coalescent specs to accommodate this somewhat niche case. Maybe the easiest would be to make a suggestion that discrete generation simulations issue a warning and use the convention above?

[molpopgen]

@sgravel -- I think this should be a separate issue.

Right now, I believe that we expect discrete simulators to apply events in order after rounding the times. For a discrete simulator, 10 and 10.1 will be 10 (hopefully).

[sgravel]

Thanks @molpopgen. Apologies for the incorrect issue management. I have created a new issue and hidden this one.

[molpopgen]

Thanks @molpopgen. Apologies for the incorrect issue management. I have created a new issue and hidden this one.

No worries! At least for me, it helps to think of issues in the smallest chunks possible.

[grahamgower]

The way I handled it in fwdpy11 (and maybe @molpopgen has thoughts about this as well, working with discrete generation models), is that when continuous migration occurs over a time period that includes a pulse, the continuous migration probabilities define the probability of choosing parents from each deme conditional on individuals not arriving via the pulse.

This needs some prose in the MDM. I don't think we have any checks in the reference implementation (or demes-python) that prohibits the "sum of pulses and migrations > 1".