Cytoskeleton -- Solutions to problems

import numpy as np
import matplotlib.pyplot as plt
plt.ion()

## Set a bunch of parameters
h = 10.0   # barrier height (actual barrier heigth will be h/e)
dt = 0.01  # time step for integration of the Langevin equation
sqrt_dt = np.sqrt(2.0*dt)  # precompute the prefactor of the diffusive term

x_left = -10.0  # left boundary (reflecting)
x_right = 1.0   # right boundary (absorbing)
x = 0.5*(x_left+x_right) # initial position

x_distribution = []      # empty list to collect points for the spatial distribution
t_distribution = []      # empty list to collect the barrier crossing times

n_runs = 200              # number of runs
t_total = 1000000        # max time before breaking
batch = 10000            # batch size of time points to speed up random number generation


def force(x):
    # potential 10e^{x - e^{x}}
    #f = -10*np.exp(x-np.exp(x))*(1-np.exp(x))
    y = np.exp(x) # precalculate
    return -h*y*np.exp(-y)*(1-y)


def potential(x):
    return h*np.exp(x - np.exp(x))

## report n_runs times
for n in range(n_runs):
    escaped = False
    # loop over batches
    for j in range(t_total//batch):
        eta = np.random.normal(size=batch) # make batch gaussian random numbers
        for i in range(batch): # loop over the these
            # Langevin equation for the increment in position
            dx = dt*force(x) + sqrt_dt*eta[i]
            if x+dx < x_left: # reflecting boundary
                x = x_left - (x + dx - x_left)
            elif x+dx>x_right: # absorbing boundary
                escaped = True
                break          # -> break out of loop
            else:
                x += dx

            x_distribution.append(x)

        if escaped: # if we have hit the right boundary
            t = (j*batch+i)*dt
            t_distribution.append(t) # record time
            x = 0.5*(x_left+x_right)
            print("escape event, t=",t)
            break # start with  new run

# make histogram and plot spatial distribution
x_hist, bins = np.histogram(x_distribution, bins=np.linspace(x_left, x_right, 101))
bc = 0.5*(bins[1:]+bins[:-1])
bw = bins[1]-bins[0]

plt.figure()
dis = x_hist/bw/len(x_distribution) # normalized distribution
plt.plot(bc, dis, label='observed distribution')
# add boltzmann distribution scaled to overlap with the empirical distribution
plt.plot(bc, np.exp(-potential(bc))*dis[0], label='Boltzmann distribution')
plt.yscale('log')
plt.legend()

# plot temporal distribution
t_hist, t_bins = np.histogram(t_distribution, bins=30)
bc = 0.5*(t_bins[1:]+t_bins[:-1])
bw = np.diff(t_bins)
t_dis = t_hist/bw/len(t_distribution)
plt.figure()
plt.plot(bc, t_dis, label='escape time distribution')
plt.ylabel('time')
plt.yscale('log')
plt.legend()

print("Barrier height:", h)
print("Mean time:", np.mean(t_distribution))