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 constant unfolding rate.
22 With the constant velocity experiment and a Hookian domain, the
23 unfolding force is proportional to time, so we expect exponential
24 population decay, even with multiple unfolding domains.
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 With an unfolding rate constant :math:K, the population follows
44 .. math::
45   dp/dt = Kp
46   p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \\;.
48 Therefore, a histogram of unfolding vs. force :math:p(f) normalized
49 to :math:p(0)=1 should follow
51 .. math::
52   -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
54 where :math:w is the binwidth and :math:N is the number of tested
55 domains.
57 >>> test()
58 """
60 from numpy import arange, exp, log
62 from ..histogram import Histogram
63 from ..fit import HistogramModelFitter
64 from ..manager import get_manager
65 from ..sawsim import SawsimRunner
66 from ..sawsim_histogram import sawsim_histogram
69 def probability_distribution(x, params):
70     """Exponential decay decay probability distribution.
72     .. math:: 1/\\tau \\cdot exp(-x/\\tau)
73     """
74     p = params  # convenient alias
75     p = abs(p)  # cannot normalize negative tau.
76     return (1.0/p) * exp(-x/p)
79 class ExponentialModelFitter (HistogramModelFitter):
80     """Exponential decay model fitter.
81     """
82     def model(self, params):
83         """An exponential decay model.
85         Notes
86         -----
87         .. math:: y \\propto \\cdot e^{-t/\\tau}
88         """
89         self._model_data.counts = (
90             self.info['binwidth']*self.info['N']*probability_distribution(
91                 self._model_data.bin_centers, params))
92         return self._model_data
94     def guess_initial_params(self):
95         return [self._data.mean]
97     def guess_scale(self, params):
98         return [self._data.mean]
100     def model_string(self):
101         return 'p(x) = A exp(-t/tau)'
103     def param_string(self, params):
104         pstrings = []
105         for name,param in zip(['tau'], params):
106             pstrings.append('%s=%g' % (name, param))
107         return ', '.join(pstrings)
111 def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
112                   spring_constant=1, velocity=1, N=100):
115     w = 0.1 * tau  # calculate bin width (in force)
116     A = w*num_domains*N / tau
117     theory = Histogram()
118     # A exp(-x/tau) = 0.001
119     #        -x/tau = log(0.001/A)
120     #             x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
121     theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
122     theory.bin_centers = theory.bin_edges[:-1] + w/2
123     theory.counts = w*num_domains*N*probability_distribution(
124         theory.bin_centers, [tau])
125     theory.analyze()
127     max_time_step = tau/10.0
129     param_string = (
130         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
131         '-s cantilever,hooke,%(spring_constant)g -N1 '
132         '-s folded,null -N%(num_domains)d -s unfolded,null '
133         '-k folded,unfolded,const,%(unfolding_rate)g -q folded '
134         ) % locals()
135     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
136                            bin_edges=theory.bin_edges)
138     e = ExponentialModelFitter(sim)
139     params = e.fit()
140     sim_tau = abs(params)
141     for s,t,n in [(sim_tau, tau, 'tau')]:
142         assert (s - t)/t < 0.1, 'simulation %s = %g != %g = %s' % (n,s,t,n)
143     return sim.residual(theory)
146 def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
147          velocity=[1,2], threshold=0.2):
148     m = get_manager()()
149     sr = SawsimRunner(manager=m)
150     try:
151         for nd in num_domains:
152             for ur in unfolding_rate:
153                 for sc in spring_constant:
154                     for v in velocity:
155                         residual = constant_rate(
156                             sawsim_runner=sr, num_domains=nd,
157                             unfolding_rate=ur, spring_constant=sc, velocity=v)
158                         assert residual < threshold, residual
159     finally:
160         m.teardown()