In [1]:
# import the normal modules
import numpy as np
import matplotlib.pyplot as plt

# nice plotting
import seaborn as sns
sns.set()
import pboc_utils as pboc

Set up the parameters for our simulation

In [2]:
# biological parameters
k_on = 5 # rate of monomer addition, per min
k_off = 6 # rate of monomer subtraction, per min

# simulation parameters
dt = 0.01 # min
n_steps = 3000
n_filaments = 1000

Unregulated filaments

let's run our Gillespie simulation

In [3]:
# 2D array for saving lengths
L = np.zeros([n_filaments, n_steps])
L[:,0] = 1

# first loop through the filaments
for n in range(n_filaments):
    
    # then through time:
    for t in range(1, n_steps):
        
        rand = np.random.uniform()
        
        # check for production
        if rand > 1 - k_on*dt:
            L[n,t] = L[n,t-1] + 1
            
        # check for degradation
        elif rand < k_off*dt:
            # check that length is at least 2
            if L[n,t-1] > 1:
                L[n,t] = L[n,t-1] - 1
            else:
                L[n,t] = L[n,t-1]
                
        # otherwise the length stays the same
        else: 
            L[n,t] = L[n,t-1]

let's plot the trajectories

In [4]:
for n in range(n_filaments):
    plt.plot(np.arange(n_steps)*dt, L[n,:], "k", alpha = 0.01)
plt.xlabel("time (mins)")
plt.ylabel("MT length (monomers)");

now let's plot the distribution at the last time to see the steady state behavior

In [5]:
plt.hist(L[:,-1], np.arange(50));
plt.xlabel("MT length (monomers)")
plt.ylabel("counts")
Out[5]:
Text(0,0.5,'counts')

Length-dependent $k_{off}$, via Kip3

In [6]:
# 2D array for saving lengths
L1 = np.zeros([n_filaments, n_steps])
L1[:,0] = 1

# introducting a new biological parameter
k_bind = 0.1 # rate of Kip3 binding, per min

# first loop through the filaments
for n in range(n_filaments):
    
    # then through time:
    for t in range(1, n_steps):
        
        rand = np.random.uniform()
        
        # calculate new k_off, given Kip3 binding
        new_k_off = k_bind * L1[n,t-1]
        
        # check for production
        if rand > 1 - k_on*dt:
            L1[n,t] = L1[n,t-1] + 1
            
        # check for degradation
        elif rand < new_k_off*dt:
            # check that length is at least 2
            if L1[n,t-1] > 1:
                L1[n,t] = L1[n,t-1] - 1
            else:
                L1[n,t] = L1[n,t-1]
                
        # otherwise the length stays the same
        else: 
            L1[n,t] = L1[n,t-1]
In [7]:
for n in range(n_filaments):
    plt.plot(np.arange(n_steps)*dt, L1[n,:], "k", alpha = 0.01)
plt.xlabel("time (mins)")
plt.ylabel("MT length (monomers)");
In [8]:
plt.hist(L1[:,-1], bins=np.arange(25,75));
plt.xlabel("MT length (monomers)")
plt.ylabel("counts")
Out[8]:
Text(0,0.5,'counts')

Length-dependent $k_{on}$, via limited monomer pool

In [9]:
# 2D array for saving lengths
L2 = np.zeros([n_filaments, n_steps])
L2[:,0] = 1

# new parameter to limnit the number of monomers
n_monomers = 50 # limiting the number of monomers 

# first loop through the filaments
for n in range(n_filaments):
    
    # then through time:
    for t in range(1, n_steps):
        
        rand = np.random.uniform()
        
        # calculate new k_on, based on limited monomer pool
        new_k_on = k_on * (n_monomers - L2[n,t-1]) / n_monomers
        
        # check for production
        if rand > 1 - new_k_on*dt:
            L2[n,t] = L2[n,t-1] + 1
            
        # check for degradation
        elif rand < k_off*dt:
            # check that length is at least 2
            if L2[n,t-1] > 1:
                L2[n,t] = L2[n,t-1] - 1
            else:
                L2[n,t] = L2[n,t-1]
                
        # otherwise the length stays the same
        else: 
            L2[n,t] = L2[n,t-1]
In [10]:
for n in range(n_filaments):
    plt.plot(np.arange(n_steps)*dt, L2[n,:], "k", alpha = 0.01)
plt.xlabel("time (mins)")
plt.ylabel("MT length (monomers)");
In [11]:
plt.hist(L2[:,-1], bins=np.arange(25));
plt.xlabel("MT length (monomers)")
plt.ylabel("counts")
Out[11]:
Text(0,0.5,'counts')

Both length-dependent $k_{off}$ and length-dependent $k_{on}$

In [12]:
# 2D array for saving lengths
L3 = np.zeros([n_filaments, n_steps])
L3[:,0] = 1

# new parameter to limnit the number of monomers
n_monomers = 50 # limiting the number of monomers 

# first loop through the filaments
for n in range(n_filaments):
    
    # then through time:
    for t in range(1, n_steps):
        
        rand = np.random.uniform()
        
        # calculate new k_off, givin Kip3 binding
        new_k_off = k_bind * L3[n,t-1]
        
        # calculate new k_on, based on limited monomer pool
        new_k_on = k_on * (n_monomers - L3[n,t-1]) / n_monomers
        
        # check for production
        if rand > 1 - new_k_on*dt:
            L3[n,t] = L3[n,t-1] + 1
            
        # check for degradation
        elif rand < new_k_off*dt:
            # check that length is at least 2
            if L3[n,t-1] > 1:
                L3[n,t] = L3[n,t-1] - 1
            else:
                L3[n,t] = L3[n,t-1]
                
        # otherwise the length stays the same
        else: 
            L3[n,t] = L3[n,t-1]
In [13]:
for n in range(n_filaments):
    plt.plot(np.arange(n_steps)*dt, L3[n,:], "k", alpha = 0.01)
plt.xlabel("time (mins)")
plt.ylabel("MT length (monomers)");
In [14]:
plt.hist(L3[:,-1], bins=np.arange(10,40));
plt.xlabel("MT length (monomers)")
plt.ylabel("counts")
Out[14]:
Text(0,0.5,'counts')

Compare the four distributions

In [15]:
sorted_lengths = np.sort(L[:,-1])
sorted_lengths1 = np.sort(L1[:,-1])
sorted_lengths2 = np.sort(L2[:,-1])
sorted_lengths3 = np.sort(L3[:,-1])
y_values = np.linspace(start=0, stop=1, num=n_filaments)
plt.plot(sorted_lengths, y_values, '.')
plt.plot(sorted_lengths1, y_values, '.')
plt.plot(sorted_lengths2, y_values, '.')
plt.plot(sorted_lengths3, y_values, '.')
plt.xlabel("MT length (monomers)")
plt.ylabel("cumulative distribution")
plt.legend(["unregulated","Kip3","limited monomers","Kip3 and limited monomers"])
Out[15]:
<matplotlib.legend.Legend at 0x1c19d02588>