In [1]:
# import modules
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats 

# nice plotting
import seaborn as sns
import pboc_utils as pboc

define our parameters

In [33]:
# parameters
r = 0.5 # min^-1 
gamma = 0.1 # min^-1

# array parameters
M = 20 # maximum number of mRNAs we keep track of
n_steps = 500

# array to store mRNA probabilities
mRNA = np.zeros([M, n_steps])
mRNA[0,0] = 1

# time step
dt = 0.1 # min
In [34]:
# loop through time
for t in range(1, n_steps): 
    
    # loop through mRNAS of 1 to M-1
    for m in range(1, M-1):
        mRNA[m,t] = mRNA[m,t-1] + r*dt*mRNA[m-1,t-1] + gamma*(m+1)*dt*mRNA[m+1,t-1] - \
                    r*dt*mRNA[m,t-1] - gamma*m*dt*mRNA[m,t-1]
        
    # special case for 0 mRNAS
    mRNA[0,t] = mRNA[0,t-1] + gamma*dt*mRNA[1,t-1] - r*dt*mRNA[0,t-1]
    
    # special case for maximum number of mRNAS, M
    mRNA[-1,t] = mRNA[-1,t-1] + r*dt*mRNA[-2,t-1] - gamma*(M-1)*dt*mRNA[-1,t-1]
In [35]:
# plot the results
pboc.bar_plot(mRNA, n_slices=6, dy=dt, x_label="mRNA count", y_label="time (mins)" )
Out[35]:
(<Figure size 720x576 with 1 Axes>,
 <matplotlib.axes._subplots.Axes3DSubplot at 0x1c18171748>)
In [36]:
mRNA[-1,-1]
Out[36]:
9.647730285443068e-07

Let's see if this distribution is Poisson

In [37]:
plt.plot(mRNA[:,-1],'.')
plt.xlabel("mRNA count")
plt.ylabel("frequency")

# plot poisson distrubtion to compare
n_mRNAs = np.arange(M)
lam = r/gamma
poisson_dist = stats.poisson.pmf(n_mRNAs,lam)

plt.plot(poisson_dist,'.')
plt.legend(["simulation","poisson"])
Out[37]:
<matplotlib.legend.Legend at 0x1c18064f98>