pysawsim/test/bell_rate.py

   1 # Copyright (C) 2008-2010  W. Trevor King <wking@drexel.edu>
   2 #
   3 # This program is free software: you can redistribute it and/or modify
   4 # it under the terms of the GNU General Public License as published by
   5 # the Free Software Foundation, either version 3 of the License, or
   6 # (at your option) any later version.
   7 #
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  15 #
  16 # The author may be contacted at <wking@drexel.edu> on the Internet, or
  17 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
  18 # Philadelphia PA 19104, USA.
  19
  20 """Test a Hookian chain with Bell model unfolding rate.
  21
  22 With the constant velocity experiment and a Hookian domain, the
  23 unfolding force is proportional to time, so we expect a peaked
  24 histogram.
  25
  26 Analytically, with a spring constant
  27
  28 .. math:: k = df/dx
  29
  30 and a pulling velocity
  31
  32 .. math:: v = dx/dt
  33
  34 we have a loading rate
  35
  36 .. math:: df/dt = df/dx dx/dt = kv \\;,
  37
  38 so
  39
  40 .. math:: f = kvt  + f_0 \\;.
  41
  42 Assuming :math:`f_0 = 0` and an unfolding rate constant :math:`K`, the
  43 population follows
  44
  45 .. math::
  46   dp/dt = -K_0 p = -p K_0 exp(f \\Delta x/k_B T)
  47   dp/df = dp/dt dt/df = -p K_0/kv exp(f \\Delta x/k_B T)
  48         = -p/\\rho exp(-\\alpha/\\rho) exp(f/\\rho)
  49         = -p/\\rho exp((f-\\alpha)/\\rho)
  50
  51 Where :math:`\\rho \\equiv k_B T/\\Delta x` and
  52 :math:`\\alpha \\equiv \\rho log(kv/K_0\\rho)`.
  53
  54 Events with such an exponentially increasing "hazard function" follow
  55 the Gumbel (minimum) probability distribution
  56
  57   P(f) = K_0 p(f) = 1/\\rho exp((f-\\alpha)/\\rho - exp((f-\\alpha)/\\rho))
  58
  59 which has a mean :math;`\\langle f \\rangle = \\alpha - \\gamma_e \\rho`
  60 and a variance :math:`\\sigma^2 = \pi^2 \\rho^2 / 6`, where
  61 :math:`\\gamma_e = 0.577\\ldots` is the Euler-Mascheroni constant.
  62
  63 >>> test()
  64 >>> test(num_domains=5)
  65 >>> test(unfolding_rate=5)
  66 >>> test(unfolding_distance=5)
  67 >>> test(spring_constant=5)
  68 >>> test(velocity=5)
  69
  70 Now use reasonable protein parameters.
  71
  72 >>> test(num_domains=1, unfolding_rate=1e-3, unfolding_distance=1e-9,
  73 ...      temperature=300, spring_constant=0.05, velocity=1e-6)
  74 >>> test(num_domains=50, unfolding_rate=1e-3, unfolding_distance=1e-9,
  75 ...      temperature=300, spring_constant=0.05, velocity=1e-6)
  76
  77 Problems with 50_1e-6_0.05_1e-3_1e-9_300's
  78   z = K0/kv:      17106.1 (expected 20000.0)
  79 Strange banded structure too.  Banding most pronounced for smaller
  80 forces.  The banding is due to the double-unfolding problem.  Reducing
  81 P from 1e-3 to 1e-5 (which takes a lot longer to run), gave
  82 50_1e-6_0.05_1e-3_1e-9_299, which looks much nicer.
  83 """
  84
  85 from numpy import arange, exp, log, pi, sqrt
  86
  87 from ..constants import gamma_e, kB
  88 from ..histogram import Histogram
  89 from ..fit import HistogramModelFitter
  90 from ..manager import get_manager
  91 from ..sawsim import SawsimRunner
  92 from ..sawsim_histogram import sawsim_histogram
  93
  94
  95 def probability_distribution(x, params):
  96     """Gumbel (minimum) probability distribution.
  97
  98     .. math:: 1/\\rho exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
  99
 100     Integral
 101
 102     .. math:: -exp(-exp((x-\\alpha)/\\rho))
 103
 104     So integrated over the range x = [0,\\infty]
 105
 106     .. math:: -exp(-\\infty) - (-exp(-exp(-\\alpha/\\rho)))
 107       = exp(-exp(-\\alpha/\\rho)))
 108     """
 109     p = params  # convenient alias
 110     p[1] = abs(p[1])  # cannot normalize negative rho.
 111     xs = (x - p[0]) / p[1]
 112     return (exp(exp(-p[0]/p[1]))/p[1]) * exp(xs - exp(xs))
 113
 114
 115 class GumbelModelFitter (HistogramModelFitter):
 116     """Gumbel (minimum) model fitter.
 117     """
 118     def model(self, params):
 119         """A Gumbel (minimum) decay model.
 120
 121         Notes
 122         -----
 123         .. math:: y \\propto exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
 124         """
 125         self._model_data.counts = (
 126             self.info['binwidth']*self.info['N']*probability_distribution(
 127                 self._model_data.bin_centers, params))
 128         return self._model_data
 129
 130     def guess_initial_params(self):
 131         rho = sqrt(6) * self._data.std_dev / pi
 132         alpha = self._data.mean + gamma_e * rho
 133         return (alpha, rho)
 134
 135     def guess_scale(self, params):
 136         return params
 137
 138     def model_string(self):
 139         return 'p(x) ~ exp((x-alpha)/rho - exp((x-alpha)/rho))'
 140
 141     def param_string(self, params):
 142         pstrings = []
 143         for name,param in zip(['alpha', 'rho'], params):
 144             pstrings.append('%s=%g' % (name, param))
 145         return ', '.join(pstrings)
 146
 147
 148
 149 def bell_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
 150               unfolding_distance=1, temperature=1/kB, spring_constant=1,
 151               velocity=1, N=100):
 152     loading_rate = float(spring_constant * velocity)
 153     rho = kB * temperature / unfolding_distance
 154     alpha = rho * log(loading_rate / (unfolding_rate * rho))
 155     w = 0.2 * rho # calculate bin width (in force)
 156     force_mean = alpha - gamma_e * rho
 157     theory = Histogram()
 158     theory.bin_edges = arange(start=0, stop=max(force_mean,0)+3*rho, step=w)
 159     theory.bin_centers = theory.bin_edges[:-1] + w/2
 160     theory.counts = w*num_domains*N*probability_distribution(
 161         theory.bin_centers, [alpha, rho])
 162     theory.analyze()
 163
 164     max_force_step = w/10.0
 165     max_time_step = max_force_step / loading_rate
 166     param_string = (
 167         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
 168         '-s cantilever,hooke,%(spring_constant)g -N1 '
 169         '-s folded,null -N%(num_domains)d -s unfolded,null '
 170         '-k "folded,unfolded,bell,{%(unfolding_rate)g,%(unfolding_distance)g}" '
 171         '-T %(temperature)g -q folded '
 172         ) % locals()
 173
 174     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
 175                            bin_edges=theory.bin_edges)
 176
 177     e = GumbelModelFitter(sim)
 178     params = e.fit()
 179     sim_alpha = params[0]
 180     sim_rho = abs(params[1])
 181     for s,t,n in [(sim_alpha, alpha, 'alpha'), (sim_rho, rho, 'rho')]:
 182         assert abs(s - t)/w < 3, (
 183             'simulation %s = %g != %g = %s (bin width = %g)' % (n,s,t,n,w))
 184     return sim.residual(theory)
 185
 186
 187 def test(threshold=0.2, **kwargs):
 188     m = get_manager()()
 189     sr = SawsimRunner(manager=m)
 190     try:
 191         residual = bell_rate(sawsim_runner=sr, **kwargs)
 192         assert residual < threshold, residual
 193     finally:
 194         m.teardown()