+# Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""Test a Hookian chain with Bell model unfolding rate.
+
+With the constant velocity experiment and a Hookian domain, the
+unfolding force is proportional to time, so we expect a peaked
+histogram.
+
+Analytically, with a spring constant
+
+.. math:: k = df/dx
+
+and a pulling velocity
+
+.. math:: v = dx/dt
+
+we have a loading rate
+
+.. math:: df/dt = df/dx dx/dt = kv \\;,
+
+so
+
+.. math:: f = kvt + f_0 \\;.
+
+Assuming :math:`f_0 = 0` and an unfolding rate constant :math:`K`, the
+population follows
+
+.. math::
+ dp/dt = -K_0 p = -p K_0 exp(f \\Delta x/k_B T)
+ dp/df = dp/dt dt/df = -p K_0/kv exp(f \\Delta x/k_B T)
+ = -p/\\rho exp(-\\alpha/\\rho) exp(f/\\rho)
+ = -p/\\rho exp((f-\\alpha)/\\rho)
+
+Where :math:`\\rho \\equiv k_B T/\\Delta x` and
+:math:`\\alpha \\equiv \\rho log(kv/K_0\\rho)`.
+
+Events with such an exponentially increasing "hazard function" follow
+the Gumbel (minimum) probability distribution
+
+ P(f) = K_0 p(f) = 1/\\rho exp((f-\\alpha)/\\rho - exp((f-\\alpha)/\\rho))
+
+which has a mean :math;`\\langle f \\rangle = \\alpha - \\gamma_e \\rho`
+and a variance :math:`\\sigma^2 = \pi^2 \\rho^2 / 6`, where
+:math:`\\gamma_e = 0.577\\ldots` is the Euler-Mascheroni constant.
+
+>>> test()
+>>> test(num_domains=5)
+>>> test(unfolding_rate=5)
+>>> test(unfolding_distance=5)
+>>> test(spring_constant=5)
+>>> test(velocity=5)
+
+Now use reasonable protein parameters.
+
+>>> test(num_domains=1, unfolding_rate=1e-3, unfolding_distance=1e-9,
+... temperature=300, spring_constant=0.05, velocity=1e-6)
+>>> test(num_domains=50, unfolding_rate=1e-3, unfolding_distance=1e-9,
+... temperature=300, spring_constant=0.05, velocity=1e-6)
+
+Problems with 50_1e-6_0.05_1e-3_1e-9_300's
+ z = K0/kv: 17106.1 (expected 20000.0)
+Strange banded structure too. Banding most pronounced for smaller
+forces. The banding is due to the double-unfolding problem. Reducing
+P from 1e-3 to 1e-5 (which takes a lot longer to run), gave
+50_1e-6_0.05_1e-3_1e-9_299, which looks much nicer.
+"""
+
+from numpy import arange, exp, log, pi, sqrt
+
+from ..constants import gamma_e, kB
+from ..histogram import Histogram
+from ..fit import HistogramModelFitter
+from ..manager import get_manager
+from ..sawsim import SawsimRunner
+from ..sawsim_histogram import sawsim_histogram
+
+
+def probability_distribution(x, params):
+ """Gumbel (minimum) probability distribution.
+
+ .. math:: 1/\\rho exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
+ """
+ p = params # convenient alias
+ p[1] = abs(p[1]) # cannot normalize negative rho.
+ xs = (x - p[0]) / p[1]
+ return (1.0/p[1]) * exp(xs - exp(xs))
+
+
+class GumbelModelFitter (HistogramModelFitter):
+ """Gumbel (minimum) model fitter.
+ """
+ def model(self, params):
+ """A Gumbel (minimum) decay model.
+
+ Notes
+ -----
+ .. math:: y \\propto exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
+ """
+ self._model_data.counts = (
+ self.info['binwidth']*self.info['N']*probability_distribution(
+ self._model_data.bin_centers, params))
+ return self._model_data
+
+ def guess_initial_params(self):
+ rho = sqrt(6) * self._data.std_dev / pi
+ alpha = self._data.mean + gamma_e * rho
+ return (alpha, rho)
+
+ def guess_scale(self, params):
+ return params
+
+ def model_string(self):
+ return 'p(x) ~ exp((x-alpha)/rho - exp((x-alpha)/rho))'
+
+ def param_string(self, params):
+ pstrings = []
+ for name,param in zip(['alpha', 'rho'], params):
+ pstrings.append('%s=%g' % (name, param))
+ return ', '.join(pstrings)
+
+
+
+def bell_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
+ unfolding_distance=1, temperature=1/kB, spring_constant=1,
+ velocity=1, N=100):
+ loading_rate = float(spring_constant * velocity)
+ rho = kB * temperature / unfolding_distance
+ alpha = rho * log(loading_rate / (unfolding_rate * rho))
+ w = 0.1 * rho # calculate bin width (in force)
+ force_mean = alpha - gamma_e * rho
+ theory = Histogram()
+ theory.bin_edges = arange(start=0, stop=3*force_mean, step=w)
+ theory.bin_centers = theory.bin_edges[:-1] + w/2
+ theory.counts = w*num_domains*N*probability_distribution(
+ theory.bin_centers, [alpha, rho])
+ theory.analyze()
+
+ max_force_step = w/10.0
+ max_time_step = max_force_step / loading_rate
+ param_string = (
+ '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
+ '-s cantilever,hooke,%(spring_constant)g -N1 '
+ '-s folded,null -N%(num_domains)d -s unfolded,null '
+ '-k "folded,unfolded,bell,{%(unfolding_rate)g,%(unfolding_distance)g}" '
+ '-T %(temperature)g -q folded '
+ ) % locals()
+
+ sim = sawsim_histogram(sawsim_runner, param_string, N=N,
+ bin_edges=theory.bin_edges)
+
+ e = GumbelModelFitter(sim)
+ params = e.fit()
+ sim_alpha = params[0]
+ sim_rho = abs(params[1])
+ for s,t,n in [(sim_alpha, alpha, 'alpha'), (sim_rho, rho, 'rho')]:
+ assert (s - t)/t < 0.1, 'simulation %s = %g != %g = %s' % (n,s,t,n)
+ return sim.residual(theory)
+
+
+def test(threshold=0.2, **kwargs):
+ m = get_manager()()
+ sr = SawsimRunner(manager=m)
+ try:
+ residual = bell_rate(sawsim_runner=sr, **kwargs)
+ assert residual < threshold, residual
+ finally:
+ m.teardown()
projects / sawsim.git / commitdiff