Demographic model

The original paper used a very simple instant population growth model. Under the model assumption, a population with an initial population size \(N_{0}\) would evolve \(G_{0}\) generations, instantly expand its population size to \(N_{1}\) and evolve another \(G_{1}\) generations. Such a model can be easily implemented as follows:

def ins_expansion(gen):
    'An instant population growth model'
    if gen < G0:
        return N0
    else:
        return N1

Other demographic models could be implemented similarly. For example, an exponential population growth model that expand the population size from \(N_{0}\) to \(N_{1}\) in \(G_{1}\) generations could be defined as

def exp_expansion(gen):
    'An exponential population growth model'
    if gen < G0:
        return N0
    else:
        rate = (math.log(N1) - math.log(N0))/G1
        return int(N0 * math.exp((gen - G0) * rate))

That is to say, we first solve \(r\) from \(N_{1}=N_{0}\exp\left(rG_{1}\right)\) and then calculate \(N_{t}=N_{0}\exp\left(rG\right)\) for a given generation.

There is a problem here: the above definitions treat N0, G0, N1 and G1 as global variables. This is OK for small scripts but is certainly not a good idea for larger scripts especially when different parameters will be used. A better way is to wrap these functions by another function that accept N0, G0, N1 and G1 as parameters. That is demonstrated in Example reichDemo where a function demo_model is defined to return either an instant or an exponential population growth demographic function.

Example: A demographic function producer

>>> import simuPOP as sim
>>> import math
>>> def demo_model(model, N0=1000, N1=100000, G0=500, G1=500):
...     '''Return a demographic function
...     model: linear or exponential
...     N0:   Initial sim.population size.
...     N1:   Ending sim.population size.
...     G0:   Length of burn-in stage.
...     G1:   Length of sim.population expansion stage.
...     '''
...     def ins_expansion(gen):
...         if gen < G0:
...             return N0
...         else:
...             return N1
...     rate = (math.log(N1) - math.log(N0))/G1
...     def exp_expansion(gen):
...         if gen < G0:
...             return N0
...         else:
...             return int(N0 * math.exp((gen - G0) * rate))
...     if model == 'instant':
...         return ins_expansion
...     elif model == 'exponential':
...         return exp_expansion
...
>>> # when needed, create a demographic function as follows
>>> demo_func = demo_model('exponential', 1000, 100000, 500, 500)
>>> # sim.population size at generation 700
>>> print(demo_func(700))
6309

now exiting runScriptInteractively...

Download reichDemo.py

Note

The defined demographic functions return the total population size (a number) at each generation beacuse no subpopulation is considered. A list of subpopulation sizes should be returned if there are more than one subpopulations.