Analysis of DMM data from bacteria

Import necessary libraries



In [1]:

    
import numpy as np
from scipy.optimize import leastsq, fmin_slsqp
from matplotlib.pylab import *
from matplotlib.colors import LogNorm
%matplotlib inline

Read the tables of $\mathcal{D}(\Delta t, q)$ and $\Delta t$ from files. Set the pixel size and thus the scale of wavenumbers $q$. Set the longest time index fit.



In [2]:

    
DDM = np.load('DDM_bacteria.npy')
DDM /= 4*DDM.shape[-1]**2
dts = np.load('dt_bacteria.npy')
pixelSize = 6.450/10. #in micrometre
qs = 2*np.pi/(2*DDM.shape[-1]*pixelSize) * np.arange(DDM.shape[-1])
tmax = len(DDM)

Scaling

Show $\mathcal{D}$ function of $q$ for various $\Delta t$ and select a safe range of $q$ for display.



In [3]:

    
for i in range(0,len(dts[:tmax]),5):
    plot(qs, DDM[i], 'o', color=cm.autumn(i/float(len(dts[:tmax]))))
xscale('log')
yscale('log')
ylabel(r'$\mathcal{D}$')
xlabel(r'$q\,(\mu m^{-1})$')

iqmin = 1+np.argmax(DDM[:tmax,1:].ptp(0))
iqmax = iqmin+np.where(DDM[:tmax,iqmin:].ptp(0)<2*DDM[:tmax,iqmin:].min(0))[0][0]
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
print iqmin, iqmax

Display $\mathcal{D}$ as function of $\Delta t$, $q\Delta t$ (ballistic scaling) and $q^2\Delta t$ (diffusive scaling).



In [4]:

    
iqs = np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)

figure(figsize=(18,6))
ax1 = subplot(1,3,1)
for i,iq in enumerate(iqs):
    plot(dts[:tmax], DDM[:tmax,iq], 'o', color=cm.autumn(i/10.))
xscale('log')
yscale('log')
ylabel(r'$\mathcal{D}(q,\Delta t)$')
xlabel(r'$\Delta t\,(s)$')

ax2 = subplot(1,3,2, sharey=ax1)
for i,iq in enumerate(iqs):
    plot(qs[iq]*dts[:tmax], DDM[:tmax,iq], 'o', color=cm.autumn(i/10.))
xscale('log')
xlabel(r'$q\Delta t\,(s/\mu m)$')
setp(ax2.get_yticklabels(), visible=False)
axvline(0.28)

ax2 = subplot(1,3,3, sharey=ax1)
for i,iq in enumerate(iqs):
    plot(qs[iq]**2*dts[:tmax], DDM[:tmax,iq], 'o', color=cm.autumn(i/10.))
xscale('log')
xlabel(r'$q^2\Delta t\,(s/\mu m^2)$')
setp(ax2.get_yticklabels(), visible=False)
axvline(21)









    Out[4]:





<matplotlib.lines.Line2D at 0x7f52f605d210>

Fitting

Define all functions necessary for the fit



In [6]:

    
def Pv(dts, q, tr, Z):
    a = dts / tr / (Z+1.)
    return 1./(a*Z) * np.sin(Z * np.arctan(a)) * (1 + a**2)**(-Z/2.)

ISF = lambda dts, q, td, tr, alpha, Z: np.exp(-dts/td) * ((1-alpha) + alpha * Pv(dts, q, tr, Z))
#bound A>0, B>0, tau_d>0, tau_r>0 and 0<alpha<1
LogDDM = lambda p, q, dts: np.log(np.abs(p[0]) * (1 - ISF(dts, q, np.abs(p[2]), np.abs(p[3]), alpha = 0.5*(1 + np.tanh(p[-2])), Z=p[-1])) + np.abs(p[1]));
simpleLogDDM = lambda p, dts: np.log(p[0] * (1-np.exp(-dts/p[2])) + p[1])

Choose initial parameters



In [16]:

    
figure(figsize=(8,6))
iq = 100
#test initial parameters
plot(dts[:tmax], DDM[:tmax,iq], 'o')
plot(dts[:tmax], DDM[:tmax,iq].ptp()*(1-ISF(dts[:tmax], qs[iq], 1/(0.28 * qs[iq]**2), 1./(21*qs[iq]), 0.65, 1)) + DDM[:tmax,iq].min(), 'k-')
xscale('log')
yscale('log')
ylabel(r'$\mathcal{D}$')
xlabel(r'$\Delta t\,(s)$')









    Out[16]:





<matplotlib.text.Text at 0x7f52fc03e510>

Perform the fit at each $q$. Generate the values of the fitted function



In [8]:

    
params = np.zeros((DDM.shape[-1], 6))
matrixfit = np.zeros(DDM[:tmax].T.shape)
for iq, ddm in enumerate(DDM[:tmax].T):
    params[iq] = leastsq(
        #function to minimize
        lambda p, q, dts, logd: LogDDM(p, q, dts) - logd,
        #initial parameters
        [ddm.ptp(), ddm.min(), 1/(0.28 * qs[iq]**2), 1./(21*qs[iq]), 0, 1],
        #data on which to perform minimization
        args=(qs[iq], dts[:tmax], np.log(ddm))
        )[0]
    matrixfit[iq] = np.exp(LogDDM(params[iq], qs[iq], dts[:tmax]))
params[:,:4] = np.abs(params[:,:4])
params[:,-2] = 0.5*(1 + np.tanh(params[:,-2]))









    



-c:8: RuntimeWarning: divide by zero encountered in double_scalars
-c:3: RuntimeWarning: divide by zero encountered in divide
-c:3: RuntimeWarning: invalid value encountered in multiply
/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:418: RuntimeWarning: Number of calls to function has reached maxfev = 1400.
  warnings.warn(errors[info][0], RuntimeWarning)
-c:7: RuntimeWarning: invalid value encountered in log

Show parameters $A(q)$ and $B(q)$ and deduce the range of $q$ where the fit is correct



In [9]:

    
figure(figsize=(8,6))
for p in params.T[:2]:
    plot(qs,p, 'o')
xscale('log')
yscale('log')
xlabel(r'$q\,(\mu m^{-1})$')
ylabel(r'$A(q),\, B(q)$')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
#ylim(1e-1,1e7)
print iqmin, iqmax

Show the results of the fit function of $q$ for various $\Delta t$ and function of $\Delta t$ for various $q$.



In [10]:

    
figure(figsize=(16,6))
ax1 = subplot(1,2,1)
for i in range(0,len(dts[:tmax]),5):
    plot(qs, DDM[i], 'o', color=cm.autumn(i/float(len(dts[:tmax]))))
    plot(qs, matrixfit[:,i], '-k')
xscale('log')
yscale('log')
ylabel(r'$\mathcal{D}$')
xlabel(r'$q\,(\mu m^{-1})$')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))

ax2 = subplot(1,2,2, sharey=ax1)
for i,iq in enumerate(np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)):
    plot(dts[:tmax], DDM[:tmax,iq], 'o', color=cm.autumn(i/10.))
    plot(dts[:tmax], matrixfit[iq], '-k')
xscale('log')
xlabel(r'$\Delta t\,(s)$')
setp(ax2.get_yticklabels(), visible=False)

ylim(1, 1e3)









    Out[10]:





(1, 1000.0)






    



/usr/lib/pymodules/python2.7/matplotlib/scale.py:90: RuntimeWarning: invalid value encountered in less_equal
  mask = a <= 0.0

Save (adapt to the folder you want to save to).



In [143]:

    
np.savetxt('bacteria.params', np.column_stack((qs, params)), header='q(nm)\tA\tB\tZ\ttr\ttd\talpha')
np.save('matrixfit_bacteria.npy', matrixfit)

Analyse fit results

Display the correlation functions as function of $\Delta t$, $q\Delta t$ (ballistic scaling) and $q^2\Delta t$ (diffusive scaling).



In [11]:

    
iqs = np.logspace(np.log10(iqmin),np.log10(iqmax),10, base=10).astype(int)

figure(figsize=(18,6))
ax1 = subplot(1,3,1)
for i,iq in enumerate(iqs):
    plot(dts[:tmax], 1- (DDM[:tmax,iq]-params[iq,1])/params[iq,0], 'o', color=cm.autumn(i/10.))
    plot(dts[:tmax], ISF(dts[:tmax], qs[iq], params[iq,2], params[iq,3], params[iq,4], params[iq,5]), '-k')
xscale('log')
ylabel(r'$f(q,\Delta t)$')
xlabel(r'$\Delta t\,(s)$')

ax2 = subplot(1,3,2, sharey=ax1)
for i,iq in enumerate(iqs):
    plot(qs[iq]*dts[:tmax], 1- (DDM[:tmax,iq]-params[iq,1])/params[iq,0], 'o', color=cm.autumn(i/10.))
    plot(qs[iq]*dts[:tmax], ISF(dts[:tmax], qs[iq], params[iq,2], params[iq,3], params[iq,4], params[iq,5]), '-k')
xscale('log')
xlabel(r'$q\Delta t\,(s/\mu m)$')
setp(ax2.get_yticklabels(), visible=False)

ax2 = subplot(1,3,3, sharey=ax1)
for i,iq in enumerate(iqs):
    plot(qs[iq]**2*dts[:tmax], 1- (DDM[:tmax,iq]-params[iq,1])/params[iq,0], 'o', color=cm.autumn(i/10.))
    plot(qs[iq]**2*dts[:tmax], ISF(dts[:tmax], qs[iq], params[iq,2], params[iq,3], params[iq,4], params[iq,5]), '-k')
xscale('log')
xlabel(r'$q^2\Delta t\,(s/\mu m^2)$')
setp(ax2.get_yticklabels(), visible=False)

ylim(-0.1,1.1)









    Out[11]:





(-0.1, 1.1)

Display the fit parameters relative to the velocity distribution



In [13]:

    
figure(figsize=(8,12))
ax1 = subplot(2,1,1)
plot(qs, params[:,-2], 'o')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
ylabel(r'$\alpha$')
ylim(0,1)
xscale('log')

ax2 = subplot(2,1,2, sharex=ax1)
plot(qs, params[:,-1], 'o')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
ylabel(r'$Z$')
ylim(-5,5)
xlabel(r'$q\, (\mu m^{-1})$')
setp(ax2.get_xticklabels(), visible=False)









    Out[13]:





[None, None, None, None, None, None]

Perform and show the fit on $\tau_r$ and $\tau_d$, print the mean velocity $\bar{v}$ and the diffusion coefficient $D$.



In [25]:

    
v = np.exp(-leastsq(
    lambda p, q, tr: p[0] - np.log(q) - np.log(tr),
    [-1.],
    args=(qs[iqmin:iqmax], params[iqmin:iqmax,3])
    )[0][0])
D = np.exp(-leastsq(
    lambda p, q, td: p[0] - 2*np.log(q) - np.log(td),
    [-1.],
    args=(qs[iqmin:iqmax], params[iqmin:iqmax,2])
    )[0][0])

plot(qs, params[:,2], 'or', label=r'$\tau_d$')
plot(qs, 1/(D*qs**2), '-k')
plot(qs, params[:,3], 'ob', label=r'$\tau_r$')
plot(qs, 1/(v*qs), '-k')
axvspan(qs[iqmin], qs[iqmax], color=(0.9,0.9,0.9))
xscale('log')
yscale('log')
xlabel(r'$q\, (\mu m^{-1})$')
ylabel(r'$\tau_d,\,\tau_r\, (s)$')
ylim(1e-3,1e7)
legend(loc='upper center',bbox_to_anchor = (0.65, 1))
print u'D = {:.02f} µm²/s'.format(D)
print u'v = {:.02f} µm/s'.format(v)









    



D = 0.20 µm²/s
v = 20.12 µm/s



In [22]:

    
legend?



In [ ]: