-
Notifications
You must be signed in to change notification settings - Fork 2
Future directions
page containing a list of ideas (possibly sorted by relevance) for extension of the simulator in format:
Problem description: Consider current results for 100 loci on one chromosome. Suppose we wish to compare this to 100 loci under the same worst-genotype selection S, but spread across 10 chromosomes. Currently the only way to do this comparison is to simulate 10 loci on each chromosome. For some applications this is an 'unrealistically small' chromosome size.
Comment: There is no intrinsic cost to simulating large chromosome size - the computational bottleneck is the number of junctions, and these are produced as a function of Lambda irrespective of chromosome size. For further thoughts at this level see 'philosophy' below.
Potential solution: De-couple the number of selected loci from the chromosome size. Each chromosome i has some number of loci Li, and some number of selected loci N[S]i
Now we can, for the example situation above, simulate 10 chromosomes of 100 loci, each with 10 loci under selection. Now each chromosome has the same number of loci as in the original example (100) AND the selection acting on each selected locus is the same as in the original example - because S is spread over total 100 loci.
Minimum changes to current code for simplest case:
- Allow N[S] to be specified (in addition to L). [Assumption: all chromosomes are identical, N[S] will be distributed deterministically at regular intervals at positions P[S], starting eg for N[S]==2 one at either end of chromosome]
- Compute the fitness of individuals as a function of their state at positions P[S], rather than their whole genome state.
- Ensure the user knows P[S] (so they can interpret output).
Future proofing: Consider, eventually, the possibility of selection sj and P[S]j being user-specified across j from 1 to N[S]i for each chromosome.
Computational design issues for simplest case: "Compute the fitness of individuals as a function of their state at positions P[S], rather than their whole genome state." involves scanning state across chromosomes. Scanning across chromosomes is potentially computationally costly. If costly we should not do it twice if once is enough. Such scanning is already done when a gamete is produced by recombination at meiosis.
For the current implementation we could:
- During meiosis scanning, record the HI of each gamete produced.
- When calculating the fitness of an individual, use the stored HIs of the two gametes that make the individual [RestrictionA: This works as long as we do not want to implement a selection model that only depends on heterozygosity].
For the N[S],L decoupled implementation we could, by analogy with the above:
- During meiosis scanning record the state, at all P[S], of each gamete produced.
- When calculating the fitness of an individual, use the P[S] states of the two gametes that make the individual.
Comment: more simply, we could record the HI over all P[S], but this suffers from RestrictionA and also would be more difficult to extend to the possible future case of of selection sj being user-specified across j from 1 to N[S]i
Should development by in branch neutral, or branch 'decoupling L and N[S]' ?
Philosophy: As stated: there is no intrinsic cost to simulating large chromosome size - the computational bottleneck is the number of junctions, and these are produced as a function of Lambda irrespective of chromosome size. BUT, just because there is no cost does not mean it is useful to decouple L from N[S]. Why decouple? Consider the case where we wish to produce expectations for block size distributions using the simulator. Now we definitely want a large number of loci per chromosome - otherwise the block size distribution becomes highly discretised.
BUT, spreading S evenly over a large number of loci on a chromosome is not guaranteed to give a similar outcome as spreading S evenly over a small number of loci distributed across that same large chromosome.
Why?
If the global parameters lead to maintenance rather than collapse, and the majority of selective effects are due to association, then: because the 'neutral' loci between selected loci are each more closely linked with a selected locus than any selected locus is with another selected locus, it is possible than the selection 'neutral' loci experience through association will maintain a barrier even at the 'neutral' loci; at the same time the selected loci are 'stable' because the selection acting on them is greater than the nearly-neutral threshold [sims are stochastic]. In contrast, if S is spread across ALL loci of a large chromosome, then selection acting on them may be less than the nearly neutral threshold... now there are no 'stable' loci, and a quasi-stable genome wide barrier may slowly decay as random small blocks escape association. Returning to the 'neutral loci' case - of course we expect random (small) neutral blocks to escape as well, but their escape does not weaken the barrier (they are neutral). Once they escape they act neutrally and will drift to loss or fixation except, because sims are flanked by infinite source pops, they will never drift to fixation - there will instead be an escape-generating vs recombination&&drift-loss equilibrium. Returning once more to the case where all loci are (a little) selected: surely we expect this same escape-generating vs recombination&&drift-loss equilibrium? No, because each selected escape reduces the barrier: before an escapee is lost due to drift, it is possible the resulting reduced barrier will allow an increase of escapees... this is a positive feedback of escapes that could lead to long-term collapse.
will be developed in branch netural. Developmental version (with unit tests only) is available.
Implementation constrains:
- selected loci have to be: (LOCI - SELECTED_LOCI) % (SELECTED_LOCI - 2) == 0 to ensure correct distribution of selected loci over the chromosome.
- selected loci is a feature of every individual (maybe should be moved to every loci)
- method countB of Chromosome got overloaded, if number of selected loci is specified, it should computed number of selected B loci. It is not implemented yet, I have not figured out how to correctly iterated though blocks. This probably require remembering position on chromosome (which was not required before).
The output block output contain full information of the simulated world (without parameters and passed generations). It would be useful to resume simulation using block file, setting file that produced that output and generation of resume. The problem is that SEED has no way to be restored, the only possibility to reproduce resumed simulation would be to define new seed for resumed simulation.
Comment:If the problem is that we do not know the SEED state that would resume the SIM such that it places the random number generator in the same state as when it stopped, so that then we could not replicate what the SIM _would have done _if it had continued... well: it did not continue, so we will never need to replicate that.
What we might want to replicate is the full behavour of the interrupted and resumed sim. If we record the seed we resumed it with, and when we resumed, we can indeed replicate this sim.
developed in local branch load_blocks.
Requires parameterization of migration. It would be just extension of 0d case. User would define number of migrants emigrating from deme in given direction (i.e. 2 numbers for 1d, 4 numbers for 2d, etc.). Possibly with relative definitions of direction (towards center of HZ, left, right and opposite directions to the first direction). However, this extension would need a neat solution of different numbers of emigrants from each deme.
Allows to explore behavior of combined tension zone with capture of the zone by higher migration towards center.
Redefine a generation method to be able to discard demes filled purely by native individuals.
Allows travelling wave exploration.