PatWright / BatSimuPOP

Simulation of population changes over time to assess the use of temporal sampling
0 stars 0 forks source link

Migration. #3

Open BoPeng opened 6 years ago

BoPeng commented 6 years ago

With single directional migration from 1 -> 0, and with different migration rates for males and females, with existing size of subpopulations, the number of males and females in two populations are like:

before migration: m/f 0: 0/68 m/f 1: 37 155
after migration: m/f 0: 5/68 m/f 1: 32 155
before migration: m/f 0: 5/68 m/f 1: 32 155
after migration: m/f 0: 8/68 m/f 1: 29 155
before migration: m/f 0: 8/68 m/f 1: 29 155
after migration: m/f 0: 10/68 m/f 1: 27 155
before migration: m/f 0: 10/68 m/f 1: 27 155
after migration: m/f 0: 11/68 m/f 1: 26 155
before migration: m/f 0: 11/68 m/f 1: 26 155
after migration: m/f 0: 12/68 m/f 1: 25 155

this is certainly not what we want so we need to decide what to do here.

PatWright commented 6 years ago

Ideally we would have to start with a 50/50 sex ratio. We could maybe have a bidirectional migration rate. With this, would it be possible to create n male genotypes, before the start of the simulation in order to even the sex ratios?

BoPeng commented 6 years ago

The script can do anything, such as moving n males from population 1, or make some females male. I do not know the consequences though.

PatWright commented 6 years ago

I think it might be better to start with the 68 females from the first population, then create an additional 62 males before the start of the simulation. This will come to a total number of 130 which is close to the effective population size I calculated for that population.

For the second population, you could create 110 females and 228 males before the start of the simulation as this would give an even sex ratio for the source population and a total of 550 individuals which is close to the whole effective population size.

If a bidirectional migration rate is maintained, then the sex ratios should remain around 50/50.

Would this be possible?

BoPeng commented 6 years ago

ok, if I simply resize the population to sizes 130 and 338 with appropriate sex change, with migration between populations, we get

before migration: m/f 0: 62/68 m/f 1: 110 228
after migration: m/f 0: 72/68 m/f 1: 100 228
before migration: m/f 0: 62/68 m/f 1: 177 161
after migration: m/f 0: 69/69 m/f 1: 170 160
before migration: m/f 0: 68/62 m/f 1: 177 161
after migration: m/f 0: 84/62 m/f 1: 161 161
before migration: m/f 0: 55/75 m/f 1: 185 153
after migration: m/f 0: 72/75 m/f 1: 168 153
before migration: m/f 0: 54/76 m/f 1: 183 155
after migration: m/f 0: 69/76 m/f 1: 168 155

if we keep constant population size (e.g. 130 and 338 bats in two populations). Note that

  1. random mating will lead to roughly equal sex, so the initial 110 female/228 male will turn to ~160 females and 160 males.
  2. although migration by probability, migration is not balanced due to different population size. That is to say, 10% of males from subpop 0 to 1, which is less than the 10% of males from subpop 1 to 0 because subpop1 has more males.

This can be see from the following pattern

before migration: m/f 0: 62/68 m/f 1: 110 228
after migration: m/f 0: 69/69 m/f 1: 103 227
before migration: m/f 0: 65/73 m/f 1: 167 163
after migration: m/f 0: 78/73 m/f 1: 154 163
before migration: m/f 0: 80/71 m/f 1: 146 171
after migration: m/f 0: 100/71 m/f 1: 126 171
before migration: m/f 0: 76/95 m/f 1: 149 148
after migration: m/f 0: 82/96 m/f 1: 143 147
before migration: m/f 0: 100/78 m/f 1: 141 149
after migration: m/f 0: 108/78 m/f 1: 133 149

if we do not keep subpopulation sizes constant.

PatWright commented 6 years ago

I'm trying to test different scenarios, although some work not all of them do. For example:

Scenario 1: 250 generations (changed at the end of the script) Initial pop size (130,550) first event (170,550) at gen 100 second event x0.6 decline at gen 150 (minimum size possible 30)

def demoModel(gen, pop):
    if gen < 100:           #never change
        sz = 130, 550       #never change
    elif gen < 150:     #add generation for 2nd event
        sz = 170, 550
    else:
        sz = min(30,170*0.6**(gen-150)), 550

Scenario 2: 250 generations Initial pop size (130,550) first event x0.6 decline at gen 100 (minimum size possible 30) second event (170,550) at gen 150

BoPeng commented 6 years ago

what do you mean by "not work"? For

        sz = min(30,170*0.6**(gen-150)), 550

should not it be

        sz = max(30,170*0.6**(gen-150)), 550

?

PatWright commented 6 years ago

There seems to be something wrong with the script but I can't figure what. Also, in some of the tests I tried that should last 250 generations, will only go on for 150.

from simuOpt import setOptions
setOptions(quiet=True)
import simuPOP as sim
from simuPOP.utils import importPopulation

pop = importPopulation('GENEPOP', 'Genpop_2subpops.txt')

import random

pop.addInfoFields(['age'])

# age_cat is not imported as they are defined by age 
# 0, 1, 2 as juvenile and 2-above as adult
with open('Fstat_dat(age,agecat,sex).txt') as fs:
    for idx, line in enumerate(fs):
        if idx < 10:
            # skip the first few lines
            continue
        if line.startswith('1 '):
            age, age_cat, sex = line.strip().split()[-3:]
            pop.individual(idx-18).setSex(sim.MALE if sex == '2' else sim.FEMALE )
            pop.individual(idx-18).age = int(age)
            #pop.individual(idx-18).age_cat = int(age_cat)
        elif line.startswith('2 '):
            age_cat, sex = line.strip().split()[-2:]
            pop.individual(idx-18).setSex(sim.MALE if sex == '2' else sim.FEMALE )
            # no age info
            pop.individual(idx-18).age = random.randint(0, 17)
            #pop.individual(idx-18).age_cat = int(age_cat)

pop.setVirtualSplitter(
    sim.CombinedSplitter([
        sim.SexSplitter(),
        sim.InfoSplitter(field='age', cutoff=[3, 17] ),
    ]))

# population 0, from 68 female to 68 female + 62 male
pop.resize([130, 338], propagate=True)
# set sex
for idx in range(68, 130):
    pop.individual(idx, 0).setSex(sim.MALE)
# add 110-37=73 male
for idx in range(192, 192+73):
    pop.individual(idx, 1).setSex(sim.MALE)
for idx in range(192+73, 338):
    pop.individual(idx, 1).setSex(sim.FEMALE) 

migr = sim.Migrator(
    rate=[
        [0, 0.1 ],
        [0, 0.001 ],
        [0.1, 0],
        [0.001, 0],
       ],
    mode = sim.BY_PROBABILITY,
    subPops=[(0, 'Male'), (0, 'Female'), (1, 'Male'), (1, 'Female')],
    toSubPops=[0, 1],
)

# let us test the migration
pop.addInfoFields('migrate_to')

 def demoModel(gen, pop):
    if gen < 100:           #never change
        sz = 130, 550       #never change
    elif gen < 150:     #add generation for 2nd event
        sz = 170, 550
    else:
        sz = max(30,170*0.6**(gen-150)), 550 
    print(f"{gen} young/adult sp0 : {pop.subPopSize([0, 'age < 3'])} / {pop.subPopSize([0, '3 <= age < 17'])}" +
          f" sp1 : {pop.subPopSize([1, 'age < 3'])} / {pop.subPopSize([1, '3 <= age < 17'])}")
    return sz

pop1 = pop.clone()
pop1.evolve(
    preOps=[
        migr,
        sim.InfoExec('age += 1'),
    ],
    matingScheme=sim.HeteroMating([
        # only adult individuals with age >=3 will mate and produce
        # offspring. The age of offspring will be zero.
        sim.RandomMating(ops=[
            sim.MendelianGenoTransmitter()],
            subPops=[(sim.ALL_AVAIL,'3 <= age < 17')],
            weight=-0.1),
        # individuals with age < 17 will be kept, but might be removed due to
        # population size decline
        sim.CloneMating(subPops=[(sim.ALL_AVAIL, 'age < 3'), (sim.ALL_AVAIL, '3 <= age < 17')]),
        ],
        subPopSize=demoModel),
    postOps=[
        sim.Stat(popSize=True),
        sim.PyEval(r'f"{gen} {subPopSize}\n"'),
        sim.utils.Exporter(format='GENEPOP', step=10, output='!f"{gen}.pop"', gui='batch')
    ],
    gen=250
)
BoPeng commented 6 years ago

Please post error messages if you can. It is best for you to understand exactly what the script is doing before you modify it for other scenarios. My notebook has the step by step approach to test each models separately before merging them to a complete script. Perhaps you can do the same as well.

PatWright commented 6 years ago

I managed to sort out the issue. It was mainly an issue with python and a little thing with my script. Thanks.