1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
5 # Massimo Sandal <devicerandom@gmail.com>
6 # W. Trevor King <wking@drexel.edu>
8 # This file is part of Hooke.
10 # Hooke is free software: you can redistribute it and/or modify it
11 # under the terms of the GNU Lesser General Public License as
12 # published by the Free Software Foundation, either version 3 of the
13 # License, or (at your option) any later version.
15 # Hooke is distributed in the hope that it will be useful, but WITHOUT
16 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
18 # Public License for more details.
20 # You should have received a copy of the GNU Lesser General Public
21 # License along with Hooke. If not, see
22 # <http://www.gnu.org/licenses/>.
24 """The ``polymer_fit`` module proviews :class:`PolymerFitPlugin` and
25 serveral associated :class:`hooke.command.Command`\s for Fitting
26 velocity-clamp data to various polymer models (WLC, FJC, etc.).
30 from Queue import Queue
35 from scipy.optimize import newton
37 from ..command import Command, Argument, Success, Failure
38 from ..config import Setting
39 from ..curve import Data
40 from ..util.callback import is_iterable
41 from ..util.fit import PoorFit, ModelFitter
42 from ..util.peak import Peak
43 from ..util.si import join_data_label, split_data_label
44 from . import Plugin, argument_to_setting
45 from .curve import CurveArgument
46 from .vclamp import scale
49 kB = 1.3806503e-23 # Joules/Kelvin
53 """Hyperbolic cotangent.
57 >>> coth(1.19967874) # doctest: +ELLIPSIS
65 \coth(z) \equiv \frac{\exp(z) + \exp(-z)}{\exp(z) - \exp(-z)}
69 http://mathworld.wolfram.com/HyperbolicCotangent.html
71 return 1.0/numpy.tanh(z)
74 """Inverse hyperbolic cotangent.
78 >>> arccoth(1.19967874) # doctest: +ELLIPSIS
80 >>> arccoth(coth(numpy.pi)) # doctest: +ELLIPSIS
85 Inverting from the :func:`definition <coth>`.
88 z \equiv \atanh\left( \frac{1}{\coth(z)} \right)
90 return numpy.arctanh(1.0/z)
97 >>> Langevin(numpy.pi) # doctest: +ELLIPSIS
102 From `Wikipedia`_ or Hatfield [#]_
105 L(z) \equiv \coth(z) - \frac{1}{z}
108 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
110 .. [#]: J.W. Hatfield and S.R. Quake
111 "Dynamic Properties of an Extended Polymer in Solution."
112 Phys. Rev. Lett. 1999.
113 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
115 return coth(z) - 1.0/z
117 def inverse_Langevin(z, extreme=1.0 - 1e-8):
118 """Inverse Langevin function.
122 z : float or array_like
123 object whose inverse Langevin will be returned
125 value (close to one) for which we assume the inverse is +/-
126 infinity. This avoids problems with Newton-Raphson
127 convergence in regions with low slope.
131 >>> inverse_Langevin(Langevin(numpy.pi)) # doctest: +ELLIPSIS
133 >>> inverse_Langevin(Langevin(numpy.array([1,2,3]))) # doctest: +ELLIPSIS
138 We approximate the inverse Langevin via Newton's method
139 (:func:`scipy.optimize.newton`). For the starting point, we take
140 the first three terms from the Taylor series (from `Wikipedia`_).
143 L^{-1}(z) = 3z + \frac{9}{5}z^3 + \frac{297}{175}z^5 + \dots
146 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
149 ret = numpy.ndarray(shape=z.shape, dtype=z.dtype)
150 for i,z_i in enumerate(z):
151 ret[i] = inverse_Langevin(z_i)
157 # Catch stdout since sometimes the newton function prints exit
158 # messages to stdout, e.g. "Tolerance of %s reached".
160 tmp_stdout = StringIO.StringIO()
161 sys.stdout = tmp_stdout
162 ret = newton(func=lambda x: Langevin(x)-z,
163 x0=3*z + 9./5.*z**3 + 297./175.*z**5,)
165 output = tmp_stdout.getvalue()
168 def FJC_fn(x_data, T, L, a):
169 """The freely jointed chain model.
174 x values for which the WLC tension should be calculated.
176 temperature in Kelvin.
178 contour length in meters.
180 Kuhn length in meters.
189 >>> FJC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
190 ... T=300, L=15e-9, a=2.5e-10) # doctest: +ELLIPSIS
191 array([ 3.322...-12, 1.78...e-11, 4.889...e-11])
195 The freely jointed chain model is
198 F(x) = \frac{k_B T}{a} L^{-1}\left( \frac{x}{L} \right)
200 where :math:`L^{-1}` is the :func:`inverse_Langevin`, :math:`a` is
201 the Kuhn length, and :math:`L` is the contour length [#]_. For
202 the inverse formulation (:math:`x(F)`), see Ray and Akhremitchev [#]_.
205 .. [#]: J.W. Hatfield and S.R. Quake
206 "Dynamic Properties of an Extended Polymer in Solution."
207 Phys. Rev. Lett. 1999.
208 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
210 .. [#] C. Ray and B.B. Akhremitchev.
211 `"Fitting Force Curves by the Extended Freely Jointed Chain Model" <http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf>`
213 return kB*T / a * inverse_Langevin(x_data/L)
215 class FJC (ModelFitter):
216 """Fit the data to a freely jointed chain.
220 Generate some example data
223 >>> L = 35e-9 # meters
224 >>> a = 2.5e-10 # meters
225 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
226 >>> d_data = FJC_fn(x_data, T=T, L=L, a=a)
228 Fit the example data with a two-parameter fit (`L` and `a`).
231 ... 'temperature (K)': T,
232 ... 'x data (m)': x_data,
234 >>> model = FJC(d_data, info=info, rescale=True)
235 >>> outqueue = Queue()
236 >>> L,a = model.fit(outqueue=outqueue)
237 >>> fit_info = outqueue.get(block=False)
243 Fit the example data with a one-parameter fit (`L`). We introduce
244 some error in our fixed Kuhn length for fun.
246 >>> info['Kuhn length (m)'] = 2*a
247 >>> model = FJC(d_data, info=info, rescale=True)
248 >>> L = model.fit(outqueue=outqueue)
249 >>> fit_info = outqueue.get(block=False)
250 >>> print L # doctest: +ELLIPSIS
254 """To avoid invalid values of `L`, we reparametrize with `Lp`.
259 contour length in meters
264 rescaled version of `L`.
268 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
269 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
270 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
272 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
274 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
276 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
281 The rescaling is designed to limit `L` to values strictly
282 greater than the maximum `x` data value, so we use
285 L = [\exp(L_p))+1] * x_\text{max}
287 which has the inverse
290 L_p = \ln(L/x_\text{max}-1)
292 This will obviously effect the interpretation of the fit's covariance matrix.
294 x_max = self.info['x data (m)'].max()
295 return numpy.log(L/x_max - 1)
298 """Inverse of :meth:`Lp`.
303 rescaled version of `L`
308 contour length in meters
312 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
313 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
314 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
316 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
319 x_max = self.info['x data (m)'].max()
320 return (numpy.exp(Lp)+1)*x_max
322 def model(self, params):
323 """Extract the relavant arguments and use :func:`FJC_fn`.
327 params : list of floats
328 `[Lp, a]`, where the presence of `a` is optional.
330 # setup convenient aliases
331 T = self.info['temperature (K)']
332 x_data = self.info['x data (m)']
333 L = self.L(params[0])
337 a = self.info['Kuhn length (m)']
339 self._model_data[:] = FJC_fn(x_data, T, L, a)
340 return self._model_data
342 def fit(self, *args, **kwargs):
343 params = super(FJC, self).fit(*args, **kwargs)
344 if is_iterable(params):
345 params[0] = self.L(params[0]) # convert Lp -> L
346 params[1] = abs(params[1]) # take the absolute value of `a`
347 else: # params is a float
348 params = self.L(params) # convert Lp -> L
351 def guess_initial_params(self, outqueue=None):
352 """Guess initial fitting parameters.
357 A guess at the reparameterized contour length in meters.
359 A guess at the Kuhn length in meters. If `.info` has a
360 setting for `'Kuhn length (m)'`, `a` is not returned.
362 x_max = self.info['x data (m)'].max()
363 a_given = 'Kuhn length (m)' in self.info
365 a = self.info['Kuhn length (m)']
375 def inverse_FJC_PEG_fn(F_data, T=300, N=1, k=150., Lp=3.58e-10, Lh=2.8e-10, dG=3., a=7e-10):
376 """Inverse poly(ethylene-glycol) adjusted extended FJC model.
380 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
381 ... 'dG':3.0, 'a':7e-10}
382 >>> inverse_FJC_PEG_fn(F_data=200e-12, **kwargs) # doctest: +ELLIPSIS
387 The function is that proposed by F. Oesterhelt, et al. [#]_.
390 x(F) = N_\text{S} \cdot \left[
392 \frac{L_\text{planar}}{
393 \exp\left(-\Delta G/k_B T\right) + 1}
394 + \frac{L_\text{helical}}{
395 \exp\left(+\Delta G/k_B T\right) + 1}
397 \coth\left(\frac{Fa}{k_B T}\right)
400 + \frac{F}{K_\text{S}}
402 where :math:`N_\text{S}` is the number of segments,
403 :math:`K_\text{S}` is the segment elasticicty,
404 :math:`L_\text{planar}` is the ttt contour length,
405 :math:`L_\text{helical}` is the ttg contour length,
406 :math:`\Delta G` is the Gibbs free energy difference
407 :math:`G_\text{planar}-G_\text{helical}`,
408 :math:`a` is the Kuhn length,
409 and :math:`F` is the chain tension.
411 .. [#]: F. Oesterhelt, M. Rief, and H.E. Gaub.
412 "Single molecule force spectroscopy by AFM indicates helical
413 structure of poly(ethylene-glycol) in water."
414 New Journal of Physics. 1999.
415 doi: `10.1088/1367-2630/1/1/006 <http://dx.doi.org/10.1088/1367-2630/1/1/006>`
419 g = (dG - F_data*(Lp-Lh)) / kBT
421 return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
424 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
425 """Poly(ethylene-glycol) adjusted extended FJC model.
429 z : float or array_like
430 object whose inverse Langevin will be returned
432 value (close to one) for which we assume the inverse is +/-
433 infinity. This avoids problems with Newton-Raphson
434 convergence in regions with low slope.
438 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
439 ... 'dG':3.0, 'a':7e-10}
440 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(1e-11, **kwargs), **kwargs) # doctest: +ELLIPSIS
442 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs) # doctest: +ELLIPSIS
444 >>> FJC_PEG_fn(numpy.array([1e-10,2e-10,3e-10]), **kwargs) # doctest: +ELLIPSIS
445 array([ 5.207...e-12, 1.257...e-11, 3.636...e-11])
446 >>> kwargs['N'] = 123
447 >>> FJC_PEG_fn(numpy.array([1e-8,2e-8,3e-8]), **kwargs) # doctest: +ELLIPSIS
448 array([ 4.160...e-12, 9.302...e-12, 1.830...e-11])
452 We approximate the PEG adjusted eFJC via Newton's method
453 (:func:`scipy.optimize.newton`). For the starting point, we use
454 the standard FJC function with an averaged contour length.
456 kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
457 if is_iterable(x_data):
458 ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
459 for i,x in enumerate(x_data):
460 ret[i] = FJC_PEG_fn(x, **kwargs)
464 # Approximate with the simple FJC_fn()
467 while guess == numpy.inf:
468 guess = FJC_fn(x_data, T=T, L=L, a=a)
471 fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
472 ret = guess*abs(newton(func=fn, x0=1.0))
475 class FJC_PEG (ModelFitter):
476 """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
480 Generate some example data
482 >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
483 ... 'dG':3.0, 'a':7e-10}
484 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
485 >>> d_data = FJC_PEG_fn(x_data, **kwargs)
487 Fit the example data with a two-parameter fit (`N` and `a`).
490 ... 'temperature (K)': kwargs['T'],
491 ... 'x data (m)': x_data,
492 ... 'section elasticity (N/m)': kwargs['k'],
493 ... 'planar section length (m)': kwargs['Lp'],
494 ... 'helical section length (m)': kwargs['Lh'],
495 ... 'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
497 >>> model = FJC_PEG(d_data, info=info, rescale=True)
498 >>> outqueue = Queue()
499 >>> N,a = model.fit(outqueue=outqueue)
500 >>> fit_info = outqueue.get(block=False)
506 Fit the example data with a one-parameter fit (`N`). We introduce
507 some error in our fixed Kuhn length for fun.
509 >>> info['Kuhn length (m)'] = 2*kwargs['a']
510 >>> model = FJC_PEG(d_data, info=info, rescale=True)
511 >>> N = model.fit(outqueue=outqueue)
512 >>> fit_info = outqueue.get(block=False)
513 >>> print N # doctest: +ELLIPSIS
517 """To avoid invalid values of `L`, we reparametrize with `Lr`.
522 contour length in meters
527 rescaled version of `L`.
531 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
532 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
533 >>> model.Lr(20e-9) # doctest: +ELLIPSIS
535 >>> model.Lr(19e-9) # doctest: +ELLIPSIS
537 >>> model.Lr(21e-9) # doctest: +ELLIPSIS
539 >>> model.Lr(100e-9) # doctest: +ELLIPSIS
544 The rescaling is designed to limit `L` to values strictly
545 greater than zero, so we use
548 L = \exp(L_p) * x_\text{max}
550 which has the inverse
553 L_p = \ln(L/x_\text{max})
555 This will obviously effect the interpretation of the fit's covariance matrix.
557 x_max = self.info['x data (m)'].max()
558 return numpy.log(L/x_max)
561 """Inverse of :meth:`Lr`.
566 rescaled version of `L`
571 contour length in meters
575 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
576 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
577 >>> model.L(model.Lr(21e-9)) # doctest: +ELLIPSIS
579 >>> model.L(model.Lr(100e-9)) # doctest: +ELLIPSIS
582 x_max = self.info['x data (m)'].max()
583 return numpy.exp(Lr)*x_max
587 'T': self.info['temperature (K)'],
588 'k': self.info['section elasticity (N/m)'],
589 'Lp': self.info['planar section length (m)'],
590 'Lh': self.info['helical section length (m)'],
591 'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
594 def model(self, params):
595 """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
599 params : list of floats
600 `[N, a]`, where the presence of `a` is optional.
602 # setup convenient aliases
603 T = self.info['temperature (K)']
604 x_data = self.info['x data (m)']
605 N = self.L(params[0])
609 a = self.info['Kuhn length (m)']
610 kwargs = self._kwargs()
612 self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
613 return self._model_data
615 def fit(self, *args, **kwargs):
616 params = super(FJC_PEG, self).fit(*args, **kwargs)
617 if is_iterable(params):
618 params[0] = self.L(params[0]) # convert Nr -> N
619 params[1] = abs(params[1]) # take the absolute value of `a`
620 else: # params is a float
621 params = self.L(params) # convert Nr -> N
624 def guess_initial_params(self, outqueue=None):
625 """Guess initial fitting parameters.
630 A guess at the section count.
632 A guess at the Kuhn length in meters. If `.info` has a
633 setting for `'Kuhn length (m)'`, `a` is not returned.
635 x_max = self.info['x data (m)'].max()
636 a_given = 'Kuhn length (m)' in self.info
638 a = self.info['Kuhn length (m)']
641 f_max = self._data.max()
642 kwargs = self._kwargs()
644 x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
645 N = x_max / x_section;
652 def WLC_fn(x_data, T, L, p):
653 """The worm like chain interpolation model.
658 x values for which the WLC tension should be calculated.
660 temperature in Kelvin.
662 contour length in meters.
664 persistence length in meters.
673 >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
674 ... T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
675 array([ 1.717...e-12, 1.070...e-11, 4.418...e-11])
679 The function is the simple polynomial worm like chain
680 interpolation forumula proposed by C. Bustamante, et al. [#]_.
683 F(x) = \frac{k_B T}{p} \left[
684 \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
688 .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
689 "Entropic elasticity of lambda-phage DNA."
691 doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
694 scaled_data = x_data / L
695 return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
697 class WLC (ModelFitter):
698 """Fit the data to a worm like chain.
702 Generate some example data
705 >>> L = 35e-9 # meters
706 >>> p = 2.5e-10 # meters
707 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
708 >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
710 Fit the example data with a two-parameter fit (`L` and `p`).
713 ... 'temperature (K)': T,
714 ... 'x data (m)': x_data,
716 >>> model = WLC(d_data, info=info, rescale=True)
717 >>> outqueue = Queue()
718 >>> L,p = model.fit(outqueue=outqueue)
719 >>> fit_info = outqueue.get(block=False)
725 Fit the example data with a one-parameter fit (`L`). We introduce
726 some error in our fixed persistence length for fun.
728 >>> info['persistence length (m)'] = 2*p
729 >>> model = WLC(d_data, info=info, rescale=True)
730 >>> L = model.fit(outqueue=outqueue)
731 >>> fit_info = outqueue.get(block=False)
732 >>> print L # doctest: +ELLIPSIS
736 """To avoid invalid values of `L`, we reparametrize with `Lp`.
741 contour length in meters
746 rescaled version of `L`.
750 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
751 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
752 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
754 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
756 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
758 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
763 The rescaling is designed to limit `L` to values strictly
764 greater than the maximum `x` data value, so we use
767 L = [\exp(L_p))+1] * x_\text{max}
769 which has the inverse
772 L_p = \ln(L/x_\text{max}-1)
774 This will obviously effect the interpretation of the fit's covariance matrix.
776 x_max = self.info['x data (m)'].max()
777 return numpy.log(L/x_max - 1)
780 """Inverse of :meth:`Lp`.
785 rescaled version of `L`
790 contour length in meters
794 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
795 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
796 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
798 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
801 x_max = self.info['x data (m)'].max()
802 return (numpy.exp(Lp)+1)*x_max
804 def model(self, params):
805 """Extract the relavant arguments and use :func:`WLC_fn`.
809 params : list of floats
810 `[Lp, p]`, where the presence of `p` is optional.
812 # setup convenient aliases
813 T = self.info['temperature (K)']
814 x_data = self.info['x data (m)']
815 L = self.L(params[0])
819 p = self.info['persistence length (m)']
821 self._model_data[:] = WLC_fn(x_data, T, L, p)
822 return self._model_data
824 def fit(self, *args, **kwargs):
825 params = super(WLC, self).fit(*args, **kwargs)
826 if is_iterable(params):
827 params[0] = self.L(params[0]) # convert Lp -> L
828 params[1] = abs(params[1]) # take the absolute value of `p`
829 else: # params is a float
830 params = self.L(params) # convert Lp -> L
833 def guess_initial_params(self, outqueue=None):
834 """Guess initial fitting parameters.
839 A guess at the reparameterized contour length in meters.
841 A guess at the persistence length in meters. If `.info`
842 has a setting for `'persistence length (m)'`, `p` is not
845 x_max = self.info['x data (m)'].max()
846 p_given = 'persistence length (m)' in self.info
848 p = self.info['persistence length (m)']
858 class PolymerFitPlugin (Plugin):
859 """Polymer model (WLC, FJC, etc.) fitting.
862 super(PolymerFitPlugin, self).__init__(name='polymer_fit')
863 self._arguments = [ # For Command initialization
864 Argument('polymer model', default='WLC', help="""
865 Select the default polymer model for 'polymer fit'. See the
866 documentation for descriptions of available algorithms.
868 Argument('FJC Kuhn length', type='float', default=4e-10,
869 help='Kuhn length in meters'),
870 Argument('FJC_PEG Kuhn length', type='float', default=4e-10,
871 help='Kuhn length in meters'),
872 Argument('FJC_PEG elasticity', type='float', default=150.0,
873 help='Elasticity of a PEG segment in Newtons per meter.'),
874 Argument('FJC_PEG delta G', type='float', default=3.0, help="""
875 Gibbs free energy difference between trans-trans-trans (ttt) and
876 trans-trans-gauche (ttg) PEG states in units of kBT.
878 Argument('FJC_PEG L_helical', type='float', default=2.8e-10,
879 help='Contour length of PEG in the ttg state.'),
880 Argument('FJC_PEG L_planar', type='float', default=3.58e-10,
881 help='Contour length of PEG in the ttt state.'),
882 Argument('WLC persistence length', type='float', default=4e-10,
883 help='Persistence length in meters'),
886 Setting(section=self.setting_section, help=self.__doc__)]
887 for argument in self._arguments:
888 self._settings.append(argument_to_setting(
889 self.setting_section, argument))
890 argument.default = None # if argument isn't given, use the config.
891 self._arguments.extend([ # Non-configurable arguments
892 Argument(name='input distance column', type='string',
893 default='cantilever adjusted extension (m)',
895 Name of the column to use as the surface position input.
897 Argument(name='input deflection column', type='string',
898 default='deflection (N)',
900 Name of the column to use as the deflection input.
904 PolymerFitCommand(self), PolymerFitPeaksCommand(self),
905 TranslateFlatPeaksCommand(self),
908 def dependencies(self):
911 def default_settings(self):
912 return self._settings
915 class PolymerFitCommand (Command):
916 """Polymer model (WLC, FJC, etc.) fitting.
918 Fits an entropic elasticity function to a given chunk of the
919 curve. Fit quality compares the residual with the thermal noise
920 (a good fit should be 1 or less).
922 Because the models are complicated and varied, you should
923 configure the command by setting configuration options instead of
924 using command arguments. TODO: argument_to_setting()
926 def __init__(self, plugin):
927 super(PolymerFitCommand, self).__init__(
931 Argument(name='block', aliases=['set'], type='int', default=0,
933 Data block for which the fit should be calculated. For an
934 approach/retract force curve, `0` selects the approaching curve and
935 `1` selects the retracting curve.
937 Argument(name='bounds', type='point', optional=False, count=2,
939 Indicies of points bounding the selected data.
941 ] + plugin._arguments + [
942 Argument(name='output tension column', type='string',
943 default='polymer tension',
945 Name of the column (without units) to use as the polymer tension output.
947 Argument(name='fit parameters info name', type='string',
948 default='surface deflection offset',
950 Name (without units) for storing the fit parameters in the `.info` dictionary.
953 help=self.__doc__, plugin=plugin)
955 def _run(self, hooke, inqueue, outqueue, params):
956 for key,value in params.items():
957 if value == None: # Use configured default value.
958 params[key] = self.plugin.config[key]
959 scale(hooke, params['curve'], params['block']) # TODO: is autoscaling a good idea? (explicit is better than implicit)
960 data = params['curve'].data[params['block']]
961 # HACK? rely on params['curve'] being bound to the local hooke
962 # playlist (i.e. not a copy, as you would get by passing a
963 # curve through the queue). Ugh. Stupid queues. As an
964 # alternative, we could pass lookup information through the
966 model = params['polymer model']
967 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
968 new.info = copy.deepcopy(data.info)
970 new.info['columns'].append(
971 join_data_label(params['output tension column'], 'N'))
972 z_data = data[:,data.info['columns'].index(
973 params['input distance column'])]
974 d_data = data[:,data.info['columns'].index(
975 params['input deflection column'])]
976 start,stop = params['bounds']
977 tension_data,ps = self.fit_polymer_model(
978 params, z_data, d_data, start, stop, outqueue)
979 new.info[params['fit parameters info name']] = ps
980 new.info[params['fit parameters info name']]['model'] = model
981 new[:,-1] = tension_data
982 params['curve'].data[params['block']] = new
984 def fit_polymer_model(self, params, z_data, d_data, start, stop,
986 """Railyard for the `fit_*_model` family.
988 Uses the `polymer model` configuration setting to call the
989 appropriate backend algorithm.
991 fn = getattr(self, 'fit_%s_model'
992 % params['polymer model'].replace('-','_'))
993 return fn(params, z_data, d_data, start, stop, outqueue)
995 def fit_FJC_model(self, params, z_data, d_data, start, stop,
997 """Fit the data with :class:`FJC`.
1000 'temperature (K)': self.plugin.config['temperature'],
1001 'x data (m)': z_data[start:stop],
1003 if True: # TODO: optionally free persistence length
1004 info['Kuhn length (m)'] = (
1005 params['FJC Kuhn length'])
1006 model = FJC(d_data[start:stop], info=info, rescale=True)
1008 params = model.fit(outqueue=queue)
1009 if True: # TODO: if Kuhn length fixed
1011 a = info['Kuhn length (m)']
1015 T = info['temperature (K)']
1016 fit_info = queue.get(block=False)
1017 f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1018 f_data[start:stop] = FJC_fn(z_data[start:stop], T=T, L=L, a=a)
1019 return [f_data, fit_info]
1021 def fit_FJC_PEG_model(self, params, z_data, d_data, start, stop,
1023 """Fit the data with :class:`FJC_PEG`.
1026 'temperature (K)': self.plugin.config['temperature'],
1027 'x data (m)': z_data[start:stop],
1030 if True: # TODO: optionally free persistence length
1031 info['Kuhn length (m)'] = (
1032 params['FJC Kuhn length'])
1033 model = FJC_PEG(d_data[start:stop], info=info, rescale=True)
1035 params = model.fit(outqueue=queue)
1036 if True: # TODO: if Kuhn length fixed
1038 a = info['Kuhn length (m)']
1042 T = info['temperature (K)']
1043 fit_info = queue.get(block=False)
1044 f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1045 f_data[start:stop] = FJC_PEG_fn(z_data[start:stop], **kwargs)
1046 return [f_data, fit_info]
1048 def fit_WLC_model(self, params, z_data, d_data, start, stop,
1050 """Fit the data with :class:`WLC`.
1053 'temperature (K)': self.plugin.config['temperature'],
1054 'x data (m)': z_data[start:stop],
1056 if True: # TODO: optionally free persistence length
1057 info['persistence length (m)'] = (
1058 params['WLC persistence length'])
1059 model = WLC(d_data[start:stop], info=info, rescale=True)
1061 params = model.fit(outqueue=queue)
1062 if True: # TODO: if persistence length fixed
1064 p = info['persistence length (m)']
1068 T = info['temperature (K)']
1069 fit_info = queue.get(block=False)
1070 f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1071 f_data[start:stop] = WLC_fn(z_data[start:stop], T=T, L=L, p=p)
1072 return [f_data, fit_info]
1075 class PolymerFitPeaksCommand (Command):
1076 """Polymer model (WLC, FJC, etc.) fitting.
1078 Use :class:`PolymerFitCommand` to fit the each peak in a list of
1079 previously determined peaks.
1081 def __init__(self, plugin):
1082 super(PolymerFitPeaksCommand, self).__init__(
1083 name='polymer fit peaks',
1086 Argument(name='block', aliases=['set'], type='int', default=0,
1088 Data block for which the fit should be calculated. For an
1089 approach/retract force curve, `0` selects the approaching curve and
1090 `1` selects the retracting curve.
1092 Argument(name='peak info name', type='string',
1093 default='polymer peaks',
1095 Name for the list of peaks in the `.info` dictionary.
1097 Argument(name='peak index', type='int', count=-1, default=None,
1099 Index of the selected peak in the list of peaks. Use `None` to fit all peaks.
1101 ] + plugin._arguments,
1102 help=self.__doc__, plugin=plugin)
1104 def _run(self, hooke, inqueue, outqueue, params):
1105 data = params['curve'].data[params['block']]
1106 fit_command = [c for c in hooke.commands
1107 if c.name=='polymer fit'][0]
1110 p = copy.deepcopy(params)
1111 p['curve'] = params['curve']
1112 del(p['peak info name'])
1113 del(p['peak index'])
1114 for i,peak in enumerate(data.info[params['peak info name']]):
1115 if params['peak index'] == None or i in params['peak index']:
1116 p['bounds'] = [peak.index, peak.post_index()]
1117 p['output tension column'] = peak.name
1118 p['fit parameters info name'] = peak.name
1119 fit_command.run(hooke, inq, outq, **p)
1121 if not isinstance(ret, Success):
1125 class TranslateFlatPeaksCommand (Command):
1126 """Translate flat filter peaks into polymer peaks for fitting.
1128 Use :class:`~hooke.plugin.flatfilt.FlatPeaksCommand` creates a
1129 list of peaks for regions with large derivatives. For velocity
1130 clamp measurements, these regions are usually the rebound phase
1131 after a protein domain unfolds, the cantilever detaches, etc.
1132 Because these features occur after the polymer loading phase, we
1133 need to shift the selected regions back to align them with the
1134 polymer loading regions.
1136 def __init__(self, plugin):
1137 plugin_arguments = [a for a in plugin._arguments
1138 if a.name in ['input distance column',
1139 'input deflection column']]
1142 Argument(name='block', aliases=['set'], type='int', default=0,
1144 Data block for which the fit should be calculated. For an
1145 approach/retract force curve, `0` selects the approaching curve and
1146 `1` selects the retracting curve.
1148 ] + plugin_arguments + [
1149 Argument(name='input peak info name', type='string',
1150 default='flat filter peaks',
1152 Name for the input peaks in the `.info` dictionary.
1154 Argument(name='output peak info name', type='string',
1155 default='polymer peaks',
1157 Name for the ouptput peaks in the `.info` dictionary.
1159 Argument(name='end offset', type='int', default=-1,
1161 Number of points between the end of a new peak and the start of the old.
1163 Argument(name='start fraction', type='float', default=0.2,
1165 Place the start of the new peak at `start_fraction` from the end of
1166 the previous old peak to the end of the new peak. Because the first
1167 new peak will have no previous old peak, it uses a distance of zero
1171 super(TranslateFlatPeaksCommand, self).__init__(
1172 name='flat peaks to polymer peaks',
1173 arguments=arguments,
1174 help=self.__doc__, plugin=plugin)
1176 def _run(self, hooke, inqueue, outqueue, params):
1177 data = params['curve'].data[params['block']]
1178 z_data = data[:,data.info['columns'].index(
1179 params['input distance column'])]
1180 d_data = data[:,data.info['columns'].index(
1181 params['input deflection column'])]
1182 previous_old_stop = numpy.absolute(z_data).argmin()
1184 for i,peak in enumerate(data.info[params['input peak info name']]):
1185 next_old_start = peak.index
1186 stop = next_old_start + params['end offset']
1187 z_start = z_data[previous_old_stop] + params['start fraction']*(
1188 z_data[stop] - z_data[previous_old_stop])
1189 start = numpy.absolute(z_data - z_start).argmin()
1190 p = Peak('polymer peak %d' % i,
1192 values=d_data[start:stop])
1194 previous_old_stop = peak.post_index()
1195 data.info[params['output peak info name']] = new
1199 # def dist2fit(self):
1200 # '''Calculates the average distance from data to fit, scaled by
1201 # the standard deviation of the free cantilever area (thermal