Script for the analysis of kinesin models (See Ref. [1])

References

[1] Fluctuating Currents in Stochastic Thermodynamics II. Energy Conversion and Nonequilibrium Response in Kinesin Models. Altaner, Wachtel and Vollmer (2015)

[2] Non-equilibrium fluctuations and mechanochemical couplings of a molecular motor. Lau, Lacoste, and Mallick, Phys. Rev. Lett., 99, 158102 (2007)

[3] Kinesin's Network of Chemomechanical Motor Cycles. Liepelt and Lipowsky. Phys. Rev. Lett. 98, 258102 (2007)

Set-up jupyter (ipython-notebook) and sympy libraries



In [ ]:

    
# inline plotting/interaction
%pylab inline 

from sympy import * # symbolic python
init_printing() # pretty printing

import numpy as np # numeric python

import time # timing, for performance monitoring

# activate latex text rendering
from matplotlib import rc
rc('text', usetex=True)

# use matplotlib for plotting
import matplotlib.pyplot as plt 
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

Initialize the kinesin model(s)

Overview of libraries:

models.py: contains the information about several kinesin models
cumulants.py: implements the algorithm described in [1]
lambdification.py: simplifies and transforms sympy expressions obtained for the models to lambda functions which are used to calculate numerical values



In [ ]:

    
# The lambda functions for the expressions discussed and analyzed in Ref. [1]
from lambdification import lau, ll



In [ ]:

    
# velocity, diffusion and normalized mechanical response for the kinesin model of Ref. [2]
velocityLa, diffusionLa, tmechLa = lau() #



In [ ]:

    
# velocity, hydrolysis rate, etc. for the 6-state [3] and 4-state [2] kinesin models
# NOTE: quick=False yields expressions that are numerically more stable, but take much more time and memory
vel6, hyd6, dif6, coupling6, invfano6, tmech6, \
vel4, hyd4, dif4, coupling4, invfano4, tmech4, \
vel_relerr, hyd_relerr, dif_relerr = ll(quick=True)

Set-up plot range and resolution

Full parameter space (Figs. 5 -- 9 in Ref. [1])



In [ ]:

    
[fmin, fmax, mumin, mumax] = [-30,30,-30,30] # boundaries of full phase diagram
resolution = 400 # plot resolution

plotarea = [fmin, fmax, mumin, mumax] # full plotarea

# prepare the plotting grid for kinesin 6 figures
xpts = linspace(fmin, fmax, resolution)
ypts = linspace(mumin, mumax, resolution)
X, Y  = meshgrid(xpts, ypts)

Zoomed parameter space (Fig. 10 in Ref [1])



In [ ]:

    
[fmin, fmax, mumin, mumax] = [0,20,5,30] # boundaries of zoomed phase diagram
resolution = 400 # plot resolution

zoomarea = [fmin, fmax, mumin, mumax] # zoomed plotarea

# prepare the plotting grid for kinesin 6 figures
xpts = linspace(fmin, fmax, resolution)
ypts = linspace(mumin, mumax, resolution)
XX, YY  = meshgrid(xpts, ypts)

Use the lambda functions to create data arrays

Lambda function evaluate the expressions at every point in the plotting grid



In [ ]:

    
# Full plot area
V6     = vel6(X,Y)
log_V6 = (np.log( np.abs( V6 ))/np.log(10))
H6     = hyd6(X,Y)
log_H6 = (np.log( np.abs( H6 ))/np.log(10))
D6     = dif6(X,Y)
log_D6 = (np.log( np.abs( D6 ))/np.log(10))
QTC6   = coupling6(X,Y)
F6     = invfano6(X,Y)
T6     = tmech6(X,Y)

V_relerr = vel_relerr(X,Y)
H_relerr = hyd_relerr(X,Y)
D_relerr = dif_relerr(X,Y)



In [ ]:

    
# Zoomed plot area (comparison plot Fig. 10)
VLa_comp = velocityLa(XX,YY)
TLa_comp = tmechLa(XX,YY)

V6_comp = vel6(XX,YY)
T6_comp = tmech6(XX,YY)

General Figure Settings



In [ ]:

    
## Setup for Figures (for instance for talks, slides, posters)

fig_size = (12,10) # in inch
fs = 25 # font size

font = {'family' : 'serif',
        'color'  : 'black',
        'weight' : 'normal',
        'size'   : fs,
        }

figdir = "plotdata/"

def pplot(g, t, logplot=True, highlight=0, crop=[], numerical=False):
    fig = figure(figsize=fig_size)
    
    if(numerical):
        G = g
        GG = G
        ccmm = cm.YlOrBr
    elif(logplot):
        G = np.log(np.abs(g(X,Y)))/np.log(10)
        GG = g(X,Y)
        ccmm = cm.gist_earth
    else:
        G = g(X,Y)
        GG = G
        ccmm = cm.coolwarm

    # the slicing parameter [::-1] reverses the y-axis before plotting
    
    im = imshow( G[::-1], cmap=ccmm, extent=plotarea ) # drawing the function
    if(len(crop)==2):
        im.set_clim(vmin=crop[0], vmax=crop[1])
        
    # adding the Contour lines with labels
    if( not numerical):
        cset1    = contour( X,Y, G, arange(-20,20,1),linewidths=1,linestyles="-",colors='black')
        stalling = contour( X,Y, GG,    [highlight], linewidths=3,linestyles="-",colors='white')
    
    # adding the colobar on the right
    cb = colorbar(im)
    
    # latex fashion title
    title(t, fontdict=font)
    xlabel('$f$', fontdict=font)
    ylabel(r'$\Delta \mu$', fontdict=font) # r'\newcommand' is a raw python string, meaning \n is not replaced by a newline, etc.
    
    #savefig(figdir+t+".png")
    return(fig)

Create some actual plots

Six-state model of Ref [2]



In [ ]:

    
# Some exemplary figures

pplot(vel6,"Velocity 6-state (Fig. 5b)", crop=[-20,2])
pplot(dif6,"Diffusion 6-state (Fig. 7a)", crop=[-20,2])
pplot(coupling6,"Current ratio 6-state (Fig. 5c)",logplot=False, crop=[-2,2], highlight=1)
pplot(tmech6,"Mechanical Response 6-state (Fig. 8)",logplot=False, crop=[-2,2],highlight=1)

# Efficiency of energy conversion, see Fig. 6
P_mech = X*vel6(X,Y) # mechanical power
P_chem = Y*hyd6(X,Y) # chemical power

inversion = np.sign(P_mech)*np.sign(P_chem)
eff = P_mech/P_chem
EFF = np.minimum((eff)**inversion,(eff)**(-inversion))

pplot(EFF,"Efficiency",crop=[0,1],numerical=True)
print("Done")

Corresponding plots in four-state model of Ref [1] (Figures are not in the paper, only relative errors are shown in Fig. 9)



In [ ]:

    
pplot(vel4,"Velocity 4-state", crop=[-20,2])
pplot(dif4,"Diffusion 4-state", crop=[-20,2])
pplot(coupling4,"Current ratio 4-state",logplot=False, crop=[-2,2], highlight=1)
pplot(tmech4,"Normalized mechanical response 4-state",logplot=False, crop=[-2,2],highlight=1)
print("Done")

Some plots for the model of Ref [3] (Note the deviation to the models above in physically unrealistic parameter regions)



In [ ]:

    
pplot(velocityLa,"Velocity Lau et. al. (2007)", crop=[-20,2])
pplot(diffusionLa,"Diffusion Lau et. al. (2007)", crop=[-20,2])
pplot(tmechLa,"Mechanical Response Lau et. al. (2007)",logplot=False, crop=[-2,2],highlight=1)
print("Done")

Code that generates Fig. 10 in Ref. [1]



In [ ]:

    
# Colors
tango_chameleon3 = "#4e9a06"
tango_sky3 = "#204a87"

#Fonts
overlay = {'family' : 'serif',
        'color'  : 'white',
        'weight' : 'normal',
        'size'   : 55,
        }

a=1.5 # scaling factor: figure size/font size

f = figure(figsize=(3.5*a,1.8*a))
gs = GridSpec(100,100)


axLL = f.add_subplot(gs[:,0:43])
im = axLL.imshow( T6_comp[::-1], cmap=cm.coolwarm, extent=zoomarea, aspect=1 ) # drawing the function
# adding the Contour lines with labels
stalling = contour( XX,YY, V6_comp, [0], linewidths=1.4,linestyles="-",colors='white')

FDR = contour( XX,YY, T6_comp , [1], linewidths=1,linestyles="--",colors=tango_chameleon3)
for c in FDR.collections:
    c.set_dashes([(0, (6, 3))])
NDR = contour( XX,YY, T6_comp , [0], linewidths=1,linestyles="--",colors=tango_sky3)
for c in NDR.collections:
    c.set_dashes([(0, (1.0, 1.0))])

#axLL.set_title(r'Liepelt, Lipowsky,' '\n' r'Phys. Rev. Lett. \textbf{98} (2007)', fontsize=5)
text(.5, 27, r"Liepelt \& Lipowsky" '\n' r"Phys.~Rev.~Lett.~\textbf{98} (2007)", fontdict=overlay,fontsize=6,color='black')  # u'Unícôdè ∫triŋs'


xlabel('$f$', labelpad = 2)
ylabel(r'$\Delta\mu$', labelpad=2) #labelpad: move label closer to axis

im.set_clim(vmin=-2, vmax=2)


axLa = f.add_subplot(gs[:,32:100])
im = axLa.imshow( TLa_comp [::-1], cmap=cm.coolwarm, extent=zoomarea, aspect=1 ) # drawing the function
# adding the Contour lines with labels
stalling = contour( XX,YY, VLa_comp, [0], linewidths=1.4,linestyles="-",colors='white')
FDR = contour( XX,YY, TLa_comp , [1], linewidths=1,linestyles="--",colors=tango_chameleon3)
for c in FDR.collections:
    c.set_dashes([(0, (6, 3))])
NDR = contour( XX,YY, TLa_comp , [0], linewidths=1,linestyles="--",colors=tango_sky3)
for c in NDR.collections:
    c.set_dashes([(0, (1.0, 1.0))])
text(.5, 27, r"Lau, Lacoste \& Mallick" '\n' r"Phys.~Rev.~Lett.~\textbf{99} (2007)", fontdict=overlay,fontsize=6,color='black')  # u'Unícôdè ∫triŋs'

## Modify labels
setp( axLa.get_yticklabels(), visible=False) # hide y labels on La plot
#axLa.set_xticks(arange(5,25,5)) # don't display zero
axLL.set_xticks(arange(0,20,5)) # don't display "20"

xlabel('$f$', labelpad = 2)
im.set_clim(vmin=-2, vmax=2)


cb=colorbar(im,ticks=arange(-2,2.1,.5))



In [ ]: