from ipywidgets import interactive, IntSlider, FloatSlider, FloatLogSlider, fixed
import numpy as np
from matplotlib import pyplot as plt
function for Moran model (two alleles, mutations, with selection)
def moran(N = 100, mu12 = 0.0, mu21 = 0.0, relw = 1.0, fstart = 0.5, numsteps = 1000, numtrajs = 10):
w1 = relw * 1.0
w2 = 1.0
fs = np.zeros((numtrajs,numsteps))
fig, ax = plt.subplots(nrows = 1, ncols = 3, figsize = (20, 4))
for traj in range(numtrajs):
i = int(np.round(fstart*N))
for step in range(numsteps):
f = i/N
allelebirth = np.random.choice([1, 2], 1, p = [w1*f/(w1*f+w2*(1-f)), w2*(1-f)/(w1*f+w2*(1-f))])[0]
alleledeath = np.random.choice([1, 2], 1, p = [i/N, (N-i)/N])[0]
if (allelebirth == 1 and np.random.rand() >= mu12) or (allelebirth == 2 and np.random.rand() < mu21):
i = i + 1
if alleledeath == 1:
i = i - 1
fs[traj,step] = i/N
ax[0].plot(np.arange(numsteps) / N, fs[traj,:])
ax[0].set_ylim([0.0, 1.0])
ax[0].set_xlabel('time (generations)')
ax[0].set_ylabel('frequency of allele 1')
ax[1].plot(np.arange(numsteps) / N, np.mean(fs, axis = 0), color = 'black')
ax[1].set_ylim([0.0, 1.0])
ax[1].set_xlabel('time (generations)')
ax[1].set_ylabel('mean frequency of allele 1')
ax[2].plot(np.arange(numsteps) / N, np.var(fs, axis = 0), color = 'black')
ax[2].set_ylim(ymin = 0)
ax[2].set_xlabel('time (generations)')
ax[2].set_ylabel('variance of frequency of allele 1')
does it work in principle?
moran()
interactive(moran,
N = IntSlider(value = 100, min = 2, max = 200, continuous_update = False),
mu12 = fixed(0.0), mu21 = fixed(0.0), relw = fixed(1.0),
fstart = FloatSlider(value = 0.5, min = 0.0, max = 1.0, continuous_update = False),
numsteps = IntSlider(value = 1000, min = 10, max = 10000, continuous_update = False),
numtrajs = IntSlider(valud = 1, min = 1, max = 20, continuous_update = False))
interactive(moran,
N = IntSlider(value = 100, min = 2, max = 200, continuous_update = False),
mu12 = FloatSlider(value = 0.0, min = 0.0, max = 1.0, continuous_update = False),
mu21 = FloatSlider(value = 0.0, min = 0.0, max = 1.0, continuous_update = False),
relw = fixed(1.0),
fstart = FloatSlider(value = 0.5, min = 0.0, max = 1.0, continuous_update = False),
numsteps = IntSlider(value = 1000, min = 10, max = 10000, continuous_update = False),
numtrajs = IntSlider(valud = 1, min = 1, max = 20, continuous_update = False))
interactive(moran,
N = IntSlider(value = 100, min = 2, max = 200, continuous_update = False),
mu12 = fixed(0.0), mu21 = fixed(0.0),
relw = FloatLogSlider(value = 1.0, base = 10, min = -1, max = 1, continuous_update = False),
fstart = FloatSlider(value = 0.5, min = 0.0, max = 1.0, continuous_update = False),
numsteps = IntSlider(value = 1000, min = 10, max = 10000, continuous_update = False),
numtrajs = IntSlider(valud = 1, min = 1, max = 20, continuous_update = False))