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.
20 """Test a Hookian chain with Bell model unfolding rate.
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.
26 Analytically, with a spring constant
28 .. math:: k = df/dx
30 and a pulling velocity
32 .. math:: v = dx/dt
36 .. math:: df/dt = df/dx dx/dt = kv \\;,
38 so
40 .. math:: f = kvt  + f_0 \\;.
42 Assuming :math:f_0 = 0 and an unfolding rate constant :math:K, the
43 population follows
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)
51 Where :math:\\rho \\equiv k_B T/\\Delta x and
52 :math:\\alpha \\equiv \\rho log(kv/K_0\\rho).
54 Events with such an exponentially increasing "hazard function" follow
55 the Gumbel (minimum) probability distribution
57   P(f) = K_0 p(f) = 1/\\rho exp((f-\\alpha)/\\rho - exp((f-\\alpha)/\\rho))
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.
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)
70 Now use reasonable protein parameters.
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)
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 """
85 from numpy import arange, exp, log, pi, sqrt
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
95 def probability_distribution(x, params):
96     """Gumbel (minimum) probability distribution.
98     .. math:: 1/\\rho exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
99     """
100     p = params  # convenient alias
101     p = abs(p)  # cannot normalize negative rho.
102     xs = (x - p) / p
103     return (1.0/p) * exp(xs - exp(xs))
106 class GumbelModelFitter (HistogramModelFitter):
107     """Gumbel (minimum) model fitter.
108     """
109     def model(self, params):
110         """A Gumbel (minimum) decay model.
112         Notes
113         -----
114         .. math:: y \\propto exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
115         """
116         self._model_data.counts = (
117             self.info['binwidth']*self.info['N']*probability_distribution(
118                 self._model_data.bin_centers, params))
119         return self._model_data
121     def guess_initial_params(self):
122         rho = sqrt(6) * self._data.std_dev / pi
123         alpha = self._data.mean + gamma_e * rho
124         return (alpha, rho)
126     def guess_scale(self, params):
127         return params
129     def model_string(self):
130         return 'p(x) ~ exp((x-alpha)/rho - exp((x-alpha)/rho))'
132     def param_string(self, params):
133         pstrings = []
134         for name,param in zip(['alpha', 'rho'], params):
135             pstrings.append('%s=%g' % (name, param))
136         return ', '.join(pstrings)
140 def bell_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
141               unfolding_distance=1, temperature=1/kB, spring_constant=1,
142               velocity=1, N=100):
144     rho = kB * temperature / unfolding_distance
146     w = 0.1 * rho # calculate bin width (in force)
147     force_mean = alpha - gamma_e * rho
148     theory = Histogram()
149     theory.bin_edges = arange(start=0, stop=3*force_mean, step=w)
150     theory.bin_centers = theory.bin_edges[:-1] + w/2
151     theory.counts = w*num_domains*N*probability_distribution(
152         theory.bin_centers, [alpha, rho])
153     theory.analyze()
155     max_force_step = w/10.0
157     param_string = (
158         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
159         '-s cantilever,hooke,%(spring_constant)g -N1 '
160         '-s folded,null -N%(num_domains)d -s unfolded,null '
161         '-k "folded,unfolded,bell,{%(unfolding_rate)g,%(unfolding_distance)g}" '
162         '-T %(temperature)g -q folded '
163         ) % locals()
165     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
166                            bin_edges=theory.bin_edges)
168     e = GumbelModelFitter(sim)
169     params = e.fit()
170     sim_alpha = params
171     sim_rho = abs(params)
172     for s,t,n in [(sim_alpha, alpha, 'alpha'), (sim_rho, rho, 'rho')]:
173         assert (s - t)/t < 0.1, 'simulation %s = %g != %g = %s' % (n,s,t,n)
174     return sim.residual(theory)
177 def test(threshold=0.2, **kwargs):
178     m = get_manager()()
179     sr = SawsimRunner(manager=m)
180     try:
181         residual = bell_rate(sawsim_runner=sr, **kwargs)
182         assert residual < threshold, residual
183     finally:
184         m.teardown()