from scipy.stats import poisson, nbinom
from scipy import __version__ as scipy_ver
print "scipy:", scipy_ver
np.set_printoptions(linewidth=80)
np.set_printoptions(precision=3)
from scipy import stats
def plot_freq(x, color='k', alpha=1):
bar(range(G), x, align="center", color=color, alpha=alpha)
xlabel("# mutations")
ylabel("frequency")
U = 0.003
s = 0.001
N = 1000
G = 25
def mutation_matrix(U):
mutation_rv = stats.poisson(U)
M = zeros((G,G))
for i in range(G):
for j in range(i+1):
M[i,j] = mutation_rv.pmf(i-j)
M[G-1,:] += 1-M.sum(axis=0)
assert (M.sum(axis=0)==1).all()
return M
M = mutation_matrix(U)
print M[:5,:5]
[[ 9.970e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00] [ 2.991e-03 9.970e-01 0.000e+00 0.000e+00 0.000e+00] [ 4.487e-06 2.991e-03 9.970e-01 0.000e+00 0.000e+00] [ 4.487e-09 4.487e-06 2.991e-03 9.970e-01 0.000e+00] [ 3.365e-12 4.487e-09 4.487e-06 2.991e-03 9.970e-01]]
gjfhjdgf
def selection_matrix(s):
S = diag((1-s)**arange(G))
return S
S = selection_matrix(s)
print S[:5,:5]
[[ 1. 0. 0. 0. 0. ] [ 0. 0.999 0. 0. 0. ] [ 0. 0. 0.998 0. 0. ] [ 0. 0. 0. 0.997 0. ] [ 0. 0. 0. 0. 0.996]]
E = M.dot(S)
print E[:5,:5]
[[ 9.970e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00] [ 2.991e-03 9.960e-01 0.000e+00 0.000e+00 0.000e+00] [ 4.487e-06 2.988e-03 9.950e-01 0.000e+00 0.000e+00] [ 4.487e-09 4.482e-06 2.985e-03 9.940e-01 0.000e+00] [ 3.365e-12 4.482e-09 4.478e-06 2.982e-03 9.930e-01]]
x_msb = stats.poisson(U/s).pmf(range(G))
x_free = zeros(G)
x_free[0] = 1
Evolve to MSB:
x = x_free.copy()
for _ in range(10000):
x = E.dot(x)
x /= x.sum()
plot_freq(x, color='w')
plot_freq(x_msb, color='b', alpha=0.5)
Evolve to first click of the ratchet:
x = x_free.copy()
t = 0
print x[:5]
while x[0] > 0:
t += 1
if t % N == 0: print t, x[:5]
x = E.dot(x)
x /= x.sum()
assert allclose(x.sum(),1)
n = multinomial(N, x)
assert len(n) == G
x = n/float(N)
print "CLICK"
print t, x[:5]
plot_freq(x, color='b', alpha=0.5)
[ 1. 0. 0. 0. 0.] CLICK 902 [ 0. 0.086 0.173 0.543 0.178]
def ratchet_click(s,U,N):
x = x_free.copy()
t = 0
M = mutation_matrix(U)
S = selection_matrix(s)
E = M.dot(S)
while x[0] > 0:
t += 1
x = E.dot(x)
x /= x.sum()
assert allclose(x.sum(),1)
n = multinomial(N, x)
assert len(n) == G
x = n/float(N)
return t
ratchet_click(s,U,N)
901
clicks = array([ratchet_click(s,U,N) for i in range(100)])
hist(clicks, normed=True)
xlabel("time to click")
ylabel("frequency");
Now using Diffusion equation [@Gordo2000]:
a = lambda x,x0,s: 0.6 * s * (1 - x/x0) * x
b = lambda x,N: x * (1-x) / N
G = lambda x,x0,s,N: exp(2 * N * 0.6 * s * x * (x/2 - x0) / x0)
from scipy.integrate import quad as integrate
intG = lambda z,w,x0,s,N: integrate(G, z, w, args=(x0,s,N))[0]
T0x0 = lambda x0: integrate( lambda x,x0,s,N: 2 * N / ( x * G(x,x0) ) * intG(0, x, x0, s, N), 0, x0, args=(x0, s, N))[0]
Tx01 = lambda x0: integrate( lambda x,x0,s,N: 2 * N / ( x * G(x,x0) ) * intG(0, x0, x0, s, N), x0, 1, args=(x0, s, N))[0]
T = lambda x0,s,N: T0x0(x0,s,N) + Tx01(x0,s,N)
print T(1)
clicks = array([[[ratchet_click(s,U,N) for i in range(100)] for s in [0.001, 0.01, 0.1]] for U in [0.0003, 0.003, 0.03]])
Ts = array([[[T(x_msb[0],s,N) for i in range(100)] for s in [0.001, 0.01, 0.1]] for U in [0.0003, 0.003, 0.03]])