Update to TeXLive 2011 needs protected caption subrefs

[thesis.git] / src / sawsim / discussion.tex

\section{Results and Discussion}
\label{sec:sawsim:results}

\subsection{Force curves generated by simulation}
\label{sec:sawsim:results:force-curves}

\Cref{fig:sawsim:sim-sawtooth} shows three simulated force curves from
pulling a polymer composed of eight identical protein molecules using
parameters from typical experimental settings.  The order of the peaks
in the force curves reflects the temporal sequence of the unfolding
events instead of the positions of the protein molecules in the
polymer\citep{li00}.  As observed experimentally
(\cref{fig:expt-sawtooth}), the forces at which identical protein
molecules unfold fluctuate, revealing the stochastic nature of protein
unfolding since no instrumental noise is included in the simulation.

After aquiring a series of experimental unfolding curves, we need to
fit the data to an explanatory model.  For velocity-clamp experiments
(\cref{sec:procedure,sec:sawsim:velocity-clamp}), we extract unfolding
forces from the sawtooth curves (\cref{sec:hooke}) and generate
histograms of unfolding forces.  Then we construct a parameterized
model of the experimental system (\cref{tab:sawsim:model}).  We can
then run \insilico\ experiments mimicking our \invitro\ experiments
(\cref{fig:sawsim:sim-hist}).  We extract the model parameters which
provide the best fit using a ``fit quality'' metric and a nonlinear
optimization routine (or a full parameter space sweep, for
low-dimensional parameter spaces).

\begin{table}[btp]
  \begin{center}
    \subfloat[][]{\label{tab:sawsim:domains}
      \begin{tabular}{l l l l}
        \toprule
        \multicolumn{4}{c}{Domain states} \\
        \midrule
        Domain name & Initial count & Tension model & Model parameters \\
        \midrule
        AFM cantilever & 1 & Hooke (\cref{eq:sawsim:hooke}) &
          $k_c=0.05\U{N/m}$ \\
        Folded I27 & 8 & WLC (\cref{eq:sawsim:wlc}) &
          $p=3.9\U{\AA}$, $L=5.1\U{nm}$ \\
        Unfolded I27 & 0 & WLC (\cref{eq:sawsim:wlc}) &
          $p=3.9\U{\AA}$, $L=33.8\U{nm}$ \\
        \bottomrule
      \end{tabular}} \\
    \subfloat[][]{\label{tab:sawsim:transitions}    
      \begin{tabular}{l l l l l}
        \toprule
        \multicolumn{5}{c}{Transition rates} \\
        \midrule
        Transition & Source & Target & Rate model & Model parameters \\
        \midrule
        Unfolding & Folded I27 & Unfolded I27 & Bell (\cref{eq:sawsim:bell}) &
          $k_{u0}=3.3\E{-4}\U{s$^{-1}$}$, $\Delta x=0.35\U{nm}$. \\
        \bottomrule
      \end{tabular}}
    \caption{\protect\subref{tab:sawsim:domains} Model for
      I27\textsubscript{8} domain states and
      \protect\subref{tab:sawsim:transitions} transitions between them
      (compare with \cref{fig:sawsim:domains}).  The models and
      parameters are those given by \citet{carrion-vazquez99b}.
      \citet{carrion-vazquez99b} don't list their cantilever spring
      constant (or, presumably, use it in their simulations), but we
      can estimate it from the rebound slope in their Figures~2.a and
      2.b, see
      \cref{fig:sawsim:kappa-sawteeth}.\label{tab:sawsim:model}}
  \end{center}
\end{table}

Because the unfolding behavious of an individual sawtooth curve is
stochastic (\cref{fig:sawsim:sim-sawtooth}), we cannot directly
compare single curves in our fit quality metric.  Instead, we gather
many experimental and simulated curves, and compare the aggregate
properties.  For velocity-clamp experiments, the usual aggregate
property used for comparison is a histogram of unfolding
forces\citep{carrion-vazquez99b} (\cref{fig:sawsim:sim-hist}).
Defining and extracting “unfolding force” is suprisingly complicated
(\cref{sec:hooke:unfolding-force}), but basically it is the highest
tension force achieved by the chain before an unfolding event (the
drops in the sawtooth).  The final drop is not an unfolding event, it
is the entire chain breaking away from the cantilever tip, severing
the connection between the substrate and the cantilever.

\begin{figure}
  \begin{center}
    \subfloat[][]{\asyinclude{figures/sim-sawtooth/sim-sawtooth}%
      \label{fig:sawsim:sim-sawtooth}%
    }
    \hspace{.25in}%
    \subfloat[][]{\asyinclude{figures/sim-hist/sim-hist}%
      \label{fig:sawsim:sim-hist}%
    }
    \caption{\protect\subref{fig:sawsim:sim-sawtooth} Three simulated force
      curves from pulling a polymer of eight identical protein
      molecules.  The simulation was carried out using the parameters:
      pulling speed $v=1\U{$\mu$m/s}$, cantilever spring constant
      $\kappa_c=50\U{pN/nm}$, temperature $T=300\U{K}$, persistence
      length of unfolded proteins $p_u=0.40\U{nm}$, $\Delta
      x_u=0.225\U{nm}$, and $k_{u0}=5\E{-5}\U{s$^{-1}$}$.  The contour
      length between the two linking points on a protein molecule is
      $L_{f1}=3.7\U{nm}$ in the folded form and $L_{u1}=28.1\U{nm}$ in
      the unfolded form.  These parameters are those of ubiquitin
      molecules connected through the N-C
      termini\citep{chyan04,carrion-vazquez03}.  Detachment from the
      tip or substrate is assumed to occur at a force of $400\U{pN}$.
      In experiments, detachments have been observed to occur at a
      variety of forces.  For clarity, the green and blue curves are
      offset by $200$ and $400\U{pN}$ respectively.
      \protect\subref{fig:sawsim:sim-hist} The distribution of the unfolding
      forces from $400$ simulated force curves ($3200$ data points)
      such as those shown in \protect\subref{fig:sawsim:sim-sawtooth}. The
      frequency is normalized by the total number of points, \ie, the
      height of each bin is equal to the number of data points in that
      bin divided by the total number of data
      points.\label{fig:sawsim:sim-all}}
  \end{center}
\end{figure}

\subsection{The supramolecular scaffold}
\label{sec:sawsim:results:scaffold}

Analysis of the mechanical unfolding data is complicated by the
dependence of the average unfolding force on the unfolding order due
to the serial linkage of the molecules.  Under an external stretching
force $F$, the probability of some domain unfolding in a polymer with
$N_f$ folded domains is $N_fP_1$ (\cref{eq:sawsim:prob-n}), which is
higher than the unfolding probability for a single molecule $P_1$.
Consequently, the average unfolding force is lower for the earlier
unfolding events when $N_f$ is larger, and the force should increase
as more and more molecules become unfolded.  However, there is a
competing factor that opposes this trend.  As the protein molecules
unfold, the chain becomes softer and the force loading rate becomes
lower when the pulling speed is constant.  This lower loading rate
leads to a decrease in the unfolding force (in the no-loading limit,
all unfolding events occur at a tension of $0\U{N}$).  The dependence
of the average unfolding force on the unfolding order is the result of
these two opposing effects.  \Cref{fig:sawsim:order-dep} shows the
dependence of the average unfolding force on the unfolding force peak
order (the temporal order of unfolding events) for four polymers with
$4$, $8$, $12$, and $16$ identical protein molecules.  The effect of
polymer chain softening dominates the initial unfolding events, and
the average unfolding force decreases as more molecules unfold.  After
several molecules have unfolded, the softening for each additional
unfolding event becomes less significant, the change in unfolding
probability becomes dominant, and the unfolding force increases upon
each subsequent unfolding event\citep{zinober02}.

We validate this explanation by calculating the unfolding force
probability distribution's dependence on the two competing factors.
The rate of unfolding events with respect to force is
\begin{align}
  r_{uF} &= -\deriv{F}{N_f}
    = -\frac{\dd N_f/\dd t}{\dd F/\dd t}
    = \frac{N_f k_u}{\kappa v} \\
    &= \frac{N_f k_{u0}}{\kappa v}\exp\p({\frac{F\Delta x_u}{k_B T}})
    = \frac{1}{\rho}\exp\p({\frac{F-\alpha}{\rho}}) \;,
\end{align}
where $N_f$ is the number of folded domain,
$\kappa=(1/\kappa_c+N_u/\kappa_\text{WLC})^{-1}$ is the spring
constant of the cantilever-polymer system ($\kappa_\text{WLC}$ is the
effective spring constant of one unfolded domain, assumed constant for
a particular polymer/cantilever combination), $\kappa v$ is the force
loading rate, and $k_u$ is the unfolding rate constant
(\cref{eq:sawsim:bell}).  In the last expression, $\rho\equiv
k_BT/\Delta x_u$, and $\alpha\equiv-\rho\ln(N_fk_{u0}\rho/\kappa v)$.
The event probability density for events with an exponentially
increasing likelihood function follows the Gumbel (minimum)
probability density\citep{NIST:gumbel} with $\rho$ and $\alpha$ being
the scale and location parameters, respectively\citep{hummer03}
\begin{equation}
  \mathcal{P}(F) = \frac{1}{\rho} \exp\p[{\frac{F-\alpha}{\rho}
                                           -\exp\p({\frac{F-\alpha}{\rho}})
                                           }] \;.  \label{eq:sawsim:gumbel}
\end{equation}
The distribution has a mean $\avg{F}=\alpha-\gamma_e\rho$ and a
variance $\sigma^2 = \pi^2\rho^2/6$, where $\gamma_e=0.577\ldots$ is
the Euler-Mascheroni constant.  Therefore, the unfolding force
distribution has a variance
\begin{equation}
  \sigma^2 = \frac{\p({\frac{\pi k_BT}{\Delta x_u}})^2}{6} \;,
  \label{eq:sawsim:variance}
\end{equation}
and and average\citep{brockwell02,hummer03}
% TODO: verify brockwell equivalence (p465)
\begin{equation}
  \avg{F(i)} = \frac{k_BT}{\Delta x_u}
               \p[{\ln\p({\frac{\kappa v\Delta x_u}{N_fk_{u0}k_BT}})
                   -\gamma_e}] \;,  \label{eq:sawsim:order-dep}
\end{equation}
where $N_f$ and $\kappa$ depend on the domain index $i=N_u$.  Curves
based on this formula fit the simulated data remarkably well
considering the effective WLC\index{WLC} stiffness $\kappa_\text{WLC}$
is the only fitted parameter, and that the actual WLC stiffness is not
constant, as we have assumed here, but a non-linear function of $F$.
\citet{dudko08} derived a formula for the loading rate for a WLC, but
as far as I know, nobody has found an analytical form for the
unfolding force histograms produced under such a variable loading
rate.
\nomenclature{$r_{uF}$}{Unfolding loading rate (newtons per second)}
\nomenclature{$\gamma_e$}{Euler-Macheroni constant, $\gamma_e=0.577\ldots$}

From \cref{fig:sawsim:order-dep}, we see that the proper way to
process data from mechanical unfolding experiments is to group the
curves according to the length of the polymer and to perform
statistical analysis separately for peaks with the same unfolding
order.  However, in most experiments, the tethering of the polymer to
the AFM tip is by nonspecific adsorption; as a result, the polymers
being stretched between the tip and the substrate have various
lengths\citep{1li00}.  In addition, the interactions between the tip
and the surface often cause irregular features in the beginning of the
force curve (\cref{fig:expt-sawtooth}), making the identification of
the first peak uncertain\citep{carrion-vazquez00}.  Furthermore, it is
often difficult to acquire a large amount of data in single molecule
experiments.  These difficulties make the aforementioned data analysis
approach unfeasible for many mechanical unfolding experiments.  As a
result, the values of all force peaks from polymers of different
lengths are often pooled together for statistical analysis.  To assess
the errors caused by such pooling, simulation data were analyzed using
different pooling methods and the results were compared.
\Cref{fig:sawsim:sim-hist} shows that, for a polymer with eight
protein molecules, the average unfolding force is $281\U{pN}$ with a
standard deviation of $25\U{pN}$ when all data is pooled.  If only the
first peaks in the force curves are analyzed, the average force is
$279\U{pN}$ with a standard deviation of $22\U{pN}$.  While for the
fourth and eighth peaks, the average force are $275\U{pN}$ and
$300\U{pN}$, respectively, and the standard deviations are $23\U{pN}$
and $25\U{pN}$, respectively.  As expected from the Gumbel
distribution, the width of the unfolding force distribution (insets in
\cref{fig:sawsim:order-dep}) is only weakly effected by unfolding
order, but the average unfolding force can be quite different for the
same protein because of the differences in unfolding order and polymer
length.

\begin{figure}
  \begin{center}
  \asyinclude{figures/order-dep/order-dep}
  \caption{The dependence of the unfolding force on the temporal
    unfolding order for four polymers with $4$, $8$, $12$, and $16$
    identical protein domains.  Each point in the figure is the
    average of $400$ data points.  The first point in each curve
    represents the average of only the first peak in each of the $400$
    simulated force curves, the second point represents the average of
    only the second peak, and so on.  The solid lines are fits of
    \cref{eq:sawsim:order-dep} to the simulated data, with best fit
    $\kappa_\text{WLC}=203$, $207$, $161$, and $157\U{pN/nm}$,
    respectively, for lengths $4$ through $16$.  The insets show the
    force distributions of the first, fourth, and eighth peaks, left
    to right, for the polymer with eight protein domains.  The
    parameters used for generating the data were the same as those
    used for \cref{fig:sawsim:sim-sawtooth}, except for the number of
    domains.  The histogram insets were normalized in the same way as
    in \cref{fig:sawsim:sim-hist}.\label{fig:sawsim:order-dep}}
  \end{center}
\end{figure}


\subsection{The effect of cantilever force constant}
\label{sec:sawsim:cantilever}

In mechanical unfolding experiments, the ability to observe the
unfolding of a single protein molecule depends on the tension drop
after an unfolding event such that another molecule does not unfold
immediately.  The magnitude of this drop is determined by many
factors, including the magnitude of the unfolding force, the contour
and persistence lengths of the protein polymer, the contour length
increase from unfolding, and the stiffness (force constant) of the
cantilever.  Among these, the effect of the cantilever force constant
is particularly interesting because cantilevers with a wide range of
force constants are available.  In addition, different single molecule
manipulation techniques, such as the AFM and laser tweezers, differ
mainly in the range of the spring constants of their force
transducers\citep{walton08}.  \Cref{fig:sawsim:kappa-sawteeth} shows
the simulated force curves from pulling an octamer of protein
molecules using cantilevers with different force constants, while
other parameters are identical.  For this model protein, the
appearance of the force curve does not change much until the force
constant of the cantilever reaches a certain value
($\kappa_c=50\U{pN/nm}$).  When $\kappa_c$ is lower than this value,
the individual unfolding events become less identifiable.  In order to
observe individual unfolding events, the cantilever needs to have a
force constant high enough so that the bending at the maximum force is
small in comparison with the contour length increment from the
unfolding of a single molecule.  \Cref{fig:sawsim:kappa-sawteeth} also
shows that the back side of the force peaks becomes more tilted as the
cantilever becomes softer.  This is due to the fact that the extension
(end-to-end distance) of the protein polymer has a large sudden
increase as the tension rebalances after an unfolding event.

It should also be mentioned that the contour length increment from
each unfolding event is not equal to the distance between adjacent
peaks in the force curve because the chain is never fully stretched.
This contour length increase can only be obtained by fitting the curve
to WLC\index{WLC} or other polymer models (\cref{fig:expt-sawtooth}).

\begin{figure}
\begin{center}
\asyinclude{figures/kappa-sawteeth/kappa-sawteeth}
\caption{Simulated force curves obtained from pulling a polymer with
  eight protein molecules using cantilevers with different force
  constants $\kappa_c$.  Parameters used in generating these curves
  are the same as those used in \cref{fig:sawsim:sim-all}, except the
  cantilever force constant.  Successive force curves are offset by
  $300\U{pN}$ for clarity.\label{fig:sawsim:kappa-sawteeth}}
\end{center}
\end{figure}

\subsection{Determination of $\Delta x_u$ and $k_{u0}$}
\label{sec:sawsim:results:fitting}

As mentioned in \cref{sec:sawsim:results:force-curves}, fitting
experimental unfolding force histograms to simulated histograms allows
you to extract best-fit parameters for your simulation model.  For
example, if you have Bell model unfolding
(\cref{sec:sawsim:rate:bell}), your two fitting parameters are the
zero-force unfolding rate $k_{u0}$ and the distance $\Delta x_u$ from
the native state to the transition state.  \cref{fig:sawsim:v-dep}
shows the dependence of the unfolding force on the pulling speed for
different values of $k_{u0}$ and $\Delta x_u$.  As expected, the
unfolding force increases linearly with the pulling speed in the
linear-log plot\citep{evans99}.  While the magnitude of the unfolding
forces is affected by both $k_{u0}$ and $\Delta x_u$, the slope of
speed dependence is primarily determined by $\Delta x_u$
(\cref{eq:sawsim:order-dep}).  \Cref{fig:sawsim:width-v-dep} shows
that the width of the unfolding force distribution is very sensitive
to $\Delta x_u$, as expected from the Gumbel distribution
(\cref{eq:sawsim:variance}).  To obtain the values of $k_{u0}$ and
$\Delta x_u$ for the protein, the pulling speed dependence and the
distribution of the unfolding forces from simulation, such as those
shown in \cref{fig:sawsim:v-dep} and the insets of
\cref{fig:sawsim:width-v-dep}, are compared with the experimentally
measured results.  The values of $k_{u0}$ and $\Delta x_u$ that
provide the best match are designated as the parameters describing the
protein under study.  Since $k_{u0}$ and $\Delta x_u$ affect the
unfolding forces differently, the values of both parameters can be
determined simultaneously.  The data used in plotting
\cref{fig:sawsim:all-v-dep} includes all force peaks from the
simulated force curves because most experimental data is analyzed that
way.  % TODO: all?  most data analyzed what way?

In most published literature, $k_{u0}$ and $\Delta x_u$ were fit by
carrying out simulations using a handful of possible unfolding
parameters and selected the best fit by eye%
%\citep{us,carrion-vazques99b,schlierf06}
.  This approach does not allow estimation of uncertainties in the
fitting parameters, as shown by \citet{best02}.  A more rigorous
approach involves quantifying the quality of fit between the
experimental and simulated force distributions, allowing the use of a
numerical minimization algorithm to pick the best fit parameters.  We
use the Jensen--Shannon divergence\citep{sims09,lin91}, a measure of
the similarity between two probability distributions.
\begin{equation}
  D_\text{JS}(p_e, p_s)
    = D_\text{KL}(p_e, p_m) + D_\text{KL}(p_s, p_m) \;,  \label{eq:sawsim:D_JS}
\end{equation}
where $p_e(i)$ and $p_s(i)$ are the the values of the $i^\text{th}$
bin in the experimental and simulated unfolding force histograms,
respectively.  $D_\text{KL}$ is the Kullback-Leibler divergence
\begin{equation}
  D_\text{KL}(p_p,p_q)
    = \sum_i p_p(i) \log_2\p({\frac{p_p(i)}{p_q(i)}}) \;,  \label{eq:sawsim:D_KL}
\end{equation}
where the sum is over all unfolding force histogram bins.  $p_m$ is
the symmetrized probability distribution
\begin{equation}
  p_m(i) \equiv [p_e(i)+p_s(i)]/2 \;.
\end{equation}
% DONE: Mention inter-histogram normalization? no.
%  For experiments carried out over a series of pulling velocities, we
%  simply sum residuals computed for each velocity, although it would
%  also be reasonable to weight this sum according to the number of
%  experimental unfolding events recorded for each velocity.

The major advantage of the Jensen--Shannon divergence is that
$D_\text{JS}$ is bounded ($0\le D_\text{JS}\le 1$) regardless of the
experimental and simulated histograms.  For comparison, Pearson's
$\chi^2$ test,
\begin{equation}
  D_{χ^2} = \sum_i \frac{(p_e(i)-p_s(i))^2}{p_s(i)}) \;,  \label{eq:sawsim:X2}
\end{equation}
is infinite if there is a bin for which $p_e(i)>0$ but $p_s(i)=0$.

\Cref{fig:sawsim:fit-space} shows the Jensen--Shannon divergence
calculated using \cref{eq:sawsim:D_JS} between an experimental data
set and simulation results obtained using a range of values of
$k_{u0}$ and $\Delta x_u$.  There is an order of magnitude range of
$k_{u0}$ that produce reasonable fits to experimental data
(\cref{fig:sawsim:fit-space}), which is consistent with the results
\citet{best02} obtained using a chi-square test.  The values of
$k_{u0}$ and $\Delta x_u$ can be determined to higher precision by
using both the pulling speed dependent data and the unfolding force
distribution, as well as any relevant information about the protein
from other sources.

\begin{figure}
  \begin{center}
  \subfloat[][]{\asyinclude{figures/v-dep/v-dep}%
    \label{fig:sawsim:v-dep}%
  } \\
  \subfloat[][]{\asyinclude{figures/v-dep/v-dep-sd}%
    \label{fig:sawsim:width-v-dep}%
  }
  \caption{\protect\subref{fig:sawsim:v-dep} The dependence of the
    unfolding forces on the pulling speed for three different model
    protein molecules characterized by the parameters $k_{u0}$ and
    $\Delta x_u$.  The polymer length is eight molecules, and each
    symbol is the average of $3200$ data points.
    \protect\subref{fig:sawsim:width-v-dep} The dependence of standard
    deviation of the unfolding force distribution on the pulling speed
    for the simulation data shown in
    \protect\subref{fig:sawsim:v-dep}, using the same symbols.  The
    insets show the force distribution histograms for the three
    proteins at the pulling speed of $1\U{$\mu$m/s}$.  The left,
    middle and right histograms are for the proteins represented by
    the top, middle, and bottom lines in
    \protect\subref{fig:sawsim:v-dep},
    respectively.\label{fig:sawsim:all-v-dep}}
  \end{center}
\end{figure}

\begin{figure}
  \begin{center}
  \asyinclude{figures/fit-space/fit-valley}
  \caption{Fit quality between an experimental data set and simulated
    data sets obtained using various values of unfolding rate
    parameters $k_{u0}$ and $\Delta x_u$.  The experimental data are
    from octameric ubiquitin pulled at $1\U{$\mu$m/s}$\citep{chyan04},
    and the other model parameters are the same as those in
    \cref{fig:sawsim:sim-all}.  The best fit parameters are $\Delta
    x_u=0.17\U{nm}$ and $k_{u0}=1.2\E{-2}\U{s$^{-1}$}$.  The
    simulation histograms were built from $400$ pulls at for each
    parameter pair.\label{fig:sawsim:fit-space}}
  \end{center}
\end{figure}

\subsection{Features}
\label{sec:sawsim:features}

\sawsim\ is a great improvement over existing work in this field.
\citet{best02} are the only authors to mention such automatic
simulation comparisons, and their $\chi^2$ fit only compares mean
unfolding forces over a range of speeds.  They calculate $\avg{F}$
through an iterative method, and assume a standard deviation of
$20\U{pN}$ on their simulated $\avg{F}$.  \sawsim, by comparison,
makes full use of your experimental histograms, which you specify in a
plain-text histogram file:
\begin{center}
\begin{spacing}{1}
\begin{Verbatim}
#HISTOGRAM: -v 6e-7
#Force (N)      Unfolding events
1.4e-10 1
1.5e-10 0
\ldots
3e-10   116
3.1e-10 18
3.2e-10 1

#HISTOGRAM: -v 8e-7
#Force (N)      Unfolding events
1.4e-10 0
1.5e-10 3
\ldots
3.2e-10 50
3.3e-10 13

#HISTOGRAM: -v 1e-6
#Force (N)      Unfolding events
1.5e-10 2
1.6e-10 3
\ldots
3.3e-10 24
3.4e-10 2
\end{Verbatim}
\end{spacing}
\end{center}

Each \sawsim\ run simulates a single sawtooth curve, so you need to
run many \sawsim\ instances to generate your simulated histograms.  To
automate this task, \sawsim\ comes with a \citet{Python} wrapping
library (\pysawsim), which provides convenient programmatic and
command line interfaces for generating and manipulating \sawsim\ runs.
For example, to compare the experimental histograms listed above with
simulated data over a 50-by-50 grid of $k_{u0}$ and $\Delta x$, you
would use something like
\begin{Verbatim}[samepage]
$ sawsim_hist_scan.py -f '-s cantilever,hooke,0.05 -N1 -s folded,null -N8
     -s "unfolded,wlc,{0.39e-9,28e-9}" -k "folded,unfolded,bell,{%g,x%g}"
     -q folded' -r '[1e-5,1e-3,50],[0.1e-9,1e-9,50]' --logx histograms.txt
\end{Verbatim}
That's a bit of a mouthful, so let's break it down.  Without the
\sawsim\ template (\Verb+-f ...+), we can focus on the comparison
options:
\begin{Verbatim}[samepage]
$ sawsim_hist_scan.py \ldots -r '[1e-5,1e-3,50],[0.1e-9,1e-9,50]' --logx histograms.txt
\end{Verbatim}
This sets up a two-parameter sweep, with the first parameter going
from $1\E{-5}$ to $1\E{-3}$ in 50 logarithmic steps, and the second
going from $0.1\E{-9}$ to $1\E{-9}$ in 50 linear steps.  The
\sawsim\ template defines the simulation model
(\cref{fig:sawsim:domains,tab:sawsim:model}), and \Verb+%g+ marks the
location where the swept parameters will be inserted.

Behind the scenes, \pysawsim\ is spawning several concurrent
\sawsim\ processes to take advantage of any parallel processing
facilities you may have access to (e.g.~multiple cores, MPI, PBS,
\ldots).  A 50-by-50 grid with 400 runs per pixel at about one second
per \sawsim\ pull would take arount 12 days of serial execution.
Moving the simulation to the departments' 16 core file server cuts
that execution time down to 18 hours, which will easily complete over
a quiet weekend.  Using MPI on the departments' 15 box, dual core
computer lab, the simulation would finish overnight.
\nomenclature{MPI}{Message passing interface, a parallel computing
  infrastructure}
\nomenclature{PBS}{Portable batch system, a parallel computing
  infrastructure.  You should be able to distinguish this from the
  other PBS (phosphate buffered saline) based on the context}

\sawsim\ also takes advantage of a number of optimizations for faster
execution.  One of the bottlenecks in the \sawsim\ code is the TODO:
interpolating tree.

\subsection{Testing}
\label{sec:sawsim:testing}

Once a body of code reaches a certain level of complication, it
becomes difficult to convince others (or yourself) that it's actually
working correctly.  In order to test \sawsim, I've developed a test
suite that compares simulated unfolding force histograms with
analytical histograms for a number of situations where solving for the
analytical histogram is possible.

\section{Review of current research}

There are two main approaches to modeling protein domain unfolding
under tension: Bell's and Kramers'\citep{schlierf06,hummer03,dudko06}.
Bell introduced his model in the context of cell
adhesion\citep{bell78}, but it has been widely used to model
mechanical unfolding in
proteins\citep{rief97a,carrion-vazquez99b,schlierf06} due to it's
simplicity and ease of use\citep{hummer03}.  Kramers introduced his
theory in the context of thermally activated barrier crossings, which
is how we use it here.


\subsection{Who's who}

The field of mechanical protein unfolding is developing along three main branches.
Some groups are predominantly theoretical,
\begin{itemize}
  \item Evans, University of British Columbia (Emeritus) \\
    \url{http://www.physics.ubc.ca/php/directory/research/fac-1p.phtml?entnum=55}
  \item Thirumalai, University of Maryland \\
    \url{http://www.marylandbiophysics.umd.edu/}
  \item Onuchic, University of California, San Diego \\
    \url{http://guara.ucsd.edu/}
  \item Hyeon, Chung-Ang University (Onuchic postdoc, Thirumalai postdoc?) \\
    \url{http://physics.chem.cau.ac.kr/} \\
  \item Dietz (Rief grad) \\
    \url{http://www.hd-web.de/}
  \item Hummer and Szabo, National Institute of Diabetes and Digestive and Kidney Diseases \\
    \url{http://intramural.niddk.nih.gov/research/faculty.asp?People_ID=1615}
    \url{http://intramural.niddk.nih.gov/research/faculty.asp?People_ID=1559}
\end{itemize}
and the experimentalists are usually either AFM based
\begin{itemize}
  \item Rief, Technischen Universität München \\
    \url{http://cell.e22.physik.tu-muenchen.de/gruppematthias/index.html}
  \item Fernandez, Columbia University \\
    \url{http://www.columbia.edu/cu/biology/faculty/fernandez/FernandezLabWebsite/}
  \item Oberhauser, University of Texas Medical Branch (Fernandez postdoc) \\
    \url{http://www.utmb.edu/ncb/Faculty/OberhauserAndres.html}
  \item Marszalek, Duke University (Fernandez postdoc) \\
    \url{http://smfs.pratt.duke.edu/homepage/lab.htm}
  \item Guoliang Yang, Drexel University \\
    \url{http://www.physics.drexel.edu/~gyang/}
  \item Wojcikiewicz, University of Miami \\
    \url{http://chroma.med.miami.edu/physiol/faculty-wojcikiewicz_e.htm}
\end{itemize}
or laser-tweezers based
\begin{itemize}
  \item Bustamante, University of California, Berkley \\
    \url{http://alice.berkeley.edu/}
  \item Forde, Simon Fraser University \\
    \url{http://www.sfu.ca/fordelab/index.html}
\end{itemize}

\subsection{Evolution of unfolding modeling}

Evans introduced the saddle-point Kramers' approximation in a protein unfolding context 1997 (\citet{evans97} Eqn.~3).
However, early work on mechanical unfolding focused on the simper Bell model\citep{rief97a}.%TODO
In the early `00's, the saddle-point/steepest-descent approximation to Kramer's model (\citet{hanggi90} Eqn.~4.56c) was introduced into our field\citep{dudko03,hyeon03}.%TODO
By the mid `00's, the full-blown double-integral form of Kramer's model (\citet{hanggi90} Eqn.~4.56b) was in use\citep{schlierf06}.%TODO

There have been some tangential attempts towards even fancier models.
\citet{dudko03} attempted to reduce the restrictions of the single-unfolding-path model.
\citet{hyeon03} attempted to measure the local roughness using temperature dependent unfolding.

\subsection{History of simulations}

Early molecular dynamics (MD) work on receptor-ligand breakage by Grubmuller 1996 and Izrailev 1997 (according to Evans 1997).
\citet{evans97} introduce a smart Monte Carlo (SMC) Kramers' simulation.

\subsection{History of experimental AFM unfolding experiments}

\begin{itemize}
  \item \citet{rief97a}: 
\end{itemize}

\subsection{History of experimental laser tweezer unfolding experiments}

\begin{itemize}
  \item \citet{izrailev97}: 
\end{itemize}

\section{Single-domain proteins under constant loading}

eq:sawsim:order-dep

Let $x$ be the end to end distance of the protein, $t$ be the time since loading began, $F$ be tension applied to the protein, $P$ be the surviving population of folded proteins.
Make the definitions
\begin{align}
  v &\equiv \deriv{t}{x} && \text{the pulling velocity} \\
  k &\equiv \deriv{x}{F} && \text{the loading spring constant} \\
  P_0 &\equiv P(t=0) && \text{the initial number of folded proteins} \\
  D &\equiv P_0 - P && \text{the number of dead (unfolded) proteins} \\
  \kappa &\equiv -\frac{1}{P} \deriv{t}{P} && \text{the unfolding rate}
\end{align}
\nomenclature{$\equiv$}{Defined as (\ie equivalent to)}
The proteins are under constant loading because
\begin{equation}
  \deriv{t}{F} = \deriv{x}{F}\deriv{t}{x} = kv\;,
\end{equation}
a constant, since both $k$ and $v$ are constant (\citet{evans97} in the text on the first page, \citet{dudko06} in the text just before Eqn.~4).

The instantaneous likelyhood of a protein unfolding is given by $\deriv{F}{D}$, and the unfolding histogram is merely this function discretized over a bin of width $W$(This is similar to \citet{dudko06} Eqn.~2, remembering that $\dot{F}=kv$, that their probability density is not a histogram ($W=1$), and that their pdf is normalized to $N=1$).
\begin{equation}
  h(F) \equiv \deriv{\text{bin}}{F}
    = \deriv{F}{D} \cdot \deriv{\text{bin}}{F}
    = W \deriv{F}{D}
    = -W \deriv{F}{P}
    = -W \deriv{t}{P} \deriv{F}{t}
    = \frac{W}{vk} P\kappa \label{eq:unfold:hist}
\end{equation}
Solving for theoretical histograms is merely a question of taking your chosen $\kappa$, solving for $P(f)$, and plugging into Eqn. \ref{eq:unfold:hist}.
We can also make a bit of progress solving for $P$ in terms of $\kappa$ as follows:
\begin{align}
  \kappa &\equiv -\frac{1}{P} \deriv{t}{P} \\
  -\kappa \dd t \cdot \deriv{t}{F} &= \frac{\dd P}{P} \\
  \frac{-1}{kv} \int \kappa \dd F &= \ln(P) + c \\
  P &= C\exp{\p({\frac{-1}{kv}\integral{}{}{F}{\kappa}})} \;, \label{eq:P}
\end{align}
where $c \equiv \ln(C)$ is a constant of integration scaling $P$.

\subsection{Constant unfolding rate}

In the extremely weak tension regime, the proteins' unfolding rate is independent of tension, we have
\begin{align}
  P &= C\exp{\p({\frac{-1}{kv}\integral{}{}{F}{\kappa}})}
     = C\exp{\p({\frac{-1}{kv}\kappa F})}
     = C\exp{\p({\frac{-\kappa F}{kv}})} \\
  P(0) &\equiv P_0 = C\exp(0) = C \\
  h(F) &= \frac{W}{vk} P \kappa
     = \frac{W\kappa P_0}{vk} \exp{\p({\frac{-\kappa F}{kv}})}
\end{align}
So, a constant unfolding-rate/hazard-function gives exponential decay.
Not the most earth shattering result, but it's a comforting first step, and it does show explicitly the dependence in terms of the various unfolding-specific parameters.

\subsection{Bell model}

Stepping up the intensity a bit, we come to Bell's model for unfolding
(\citet{hummer03} Eqn.~1 and the first paragraph of \citet{dudko06} and \citet{dudko07}).
\begin{equation}
  \kappa = \kappa_0 \cdot \exp\p({\frac{F \dd x}{k_B T}})
         = \kappa_0 \cdot \exp(a F) \;,
\end{equation}
where we've defined $a \equiv \dd x/k_B T$ to bundle some constants together.
The unfolding histogram is then given by
\begin{align}
  P &= C\exp\p({\frac{-1}{kv}\integral{}{}{F}{\kappa}})
     = C\exp\p[{\frac{-1}{kv} \frac{\kappa_0}{a} \exp(a F)}]
     = C\exp\p[{\frac{-\kappa_0}{akv}\exp(a F)}] \\
  P(0) &\equiv P_0 = C\exp\p({\frac{-\kappa_0}{akv}}) \\
  C &= P_0 \exp\p({\frac{\kappa_0}{akv}}) \\
  P &= P_0 \exp\p\{{\frac{\kappa_0}{akv}[1-\exp(a F)]}\} \\
  h(F) &= \frac{W}{vk} P \kappa
     = \frac{W}{vk} P_0 \exp\p\{{\frac{\kappa_0}{akv}[1-\exp(a F)]}\} \kappa_0 \exp(a F)
     = \frac{W\kappa_0 P_0}{vk} \exp\p\{{a F + \frac{\kappa_0}{akv}[1-\exp(a F)]}\} \label{eq:unfold:bell_pdf}\;.
\end{align}
The $F$ dependent behavior reduces to
\begin{equation}
  h(F) \propto \exp\p[{a F - b\exp(a F)}] \;,
\end{equation}
where $b \equiv \kappa_0/akv \equiv \kappa_0 k_B T / k v \dd x$ is
another constant rephrasing.

This looks similar to the Gompertz / Gumbel / Fisher-Tippett
distribution, where
\begin{align}
  p(x) &\propto z\exp(-z) \\
  z &\equiv \exp\p({-\frac{x-\mu}{\beta}}) \;,
\end{align}
but we have
\begin{equation}
  p(x) \propto z\exp(-bz) \;.
\end{equation}
Strangely, the Gumbel distribution is supposed to derive from an
exponentially increasing hazard function, which is where we started
for our derivation.  I haven't been able to find a good explaination
of this discrepancy yet, but I have found a source that echos my
result (\citet{wu04} Eqn.~1).  TODO: compare \citet{wu04} with
my successful derivation in \cref{sec:sawsim:results-scaffold}.

Oh wait, we can do this:
\begin{equation}
  p(x) \propto z\exp(-bz) = \frac{1}{b} z'\exp(-z')\propto z'\exp(-z') \;,
\end{equation}
with $z'\equiv bz$.  I feel silly...  From
\href{Wolfram}{http://mathworld.wolfram.com/GumbelDistribution.html},
the mean of the Gumbel probability density
\begin{equation}
  P(x) = \frac{1}{\beta} \exp\p[{\frac{x-\alpha}{\beta}
                                 -\exp\p({\frac{x-\alpha}{\beta}})
                                 }]
\end{equation}
is given by $\mu=\alpha-\gamma\beta$, and the variance is
$\sigma^2=\frac{1}{6}\pi^2\beta^2$, where $\gamma=0.57721566\ldots$ is
the Euler-Mascheroni constant.  Selecting $\beta=1/a=k_BT/\dd x$,
$\alpha=-\beta\ln(\kappa\beta/kv)$, and $F=x$ we have
\begin{align}
  P(F)
    &= \frac{1}{\beta} \exp\p[{\frac{F+\beta\ln(\kappa\beta/kv)}{\beta}
                               -\exp\p({\frac{F+\beta\ln(\kappa\beta/kv)}
                                             {\beta}})
                             }] \\
    &= \frac{1}{\beta} \exp(F/\beta)\exp[\ln(\kappa\beta/kv)]
                       \exp\p\{{-\exp(F/\beta)\exp[\ln(\kappa\beta/kv)]}\} \\
    &= \frac{1}{\beta} \frac{\kappa\beta}{kv} \exp(F/\beta)
                       \exp\p[{-\kappa\beta/kv\exp(F/\beta)}] \\
    &= \frac{\kappa}{kv} \exp(F/\beta)\exp[-\kappa\beta/kv\exp(F/\beta)] \\
    &= \frac{\kappa}{kv} \exp(F/\beta - \kappa\beta/kv\exp(F/\beta)] \\
    &= \frac{\kappa}{kv} \exp(aF - \kappa/akv\exp(aF)] \\
    &= \frac{\kappa}{kv} \exp(aF - b\exp(aF)]
    \propto h(F) \;.
\end{align}
So our unfolding force histogram for a single Bell domain under
constant loading does indeed follow the Gumbel distribution.

% Consolidate with src/sawsim/discussion.tex

\subsection{Saddle-point Kramers' model}

For the saddle-point approximation for Kramers' model for unfolding
(\citet{evans97} Eqn.~3, \citet{hanggi90} Eqn. 4.56c, \citet{vanKampen07} Eqn. XIII.2.2).
\begin{equation}
  \kappa = \frac{D}{l_b l_{ts}} \cdot \exp\p({\frac{-E_b(F)}{k_B T}}) \;,
\end{equation}
where $E_b(F)$ is the barrier height under an external force $F$,
$D$ is the diffusion constant of the protein conformation along the reaction coordinate,
$l_b$ is the characteristic length of the bound state $l_b \equiv 1/\rho_b$,
$\rho_b$ is the density of states in the bound state, and
$l_{ts}$ is the characteristic length of the transition state
\begin{equation}
  l_{ts} = TODO
\end{equation} 

\citet{evans97} solved this unfolding rate for both inverse power law potentials and cusp potentials.

\subsubsection{Inverse power law potentials}

\begin{equation}
  E(x) = \frac{-A}{x^n}
\end{equation}
(e.g. $n=6$ for a van der Waals interaction, see \citet{evans97} in
the text on page 1544, in the first paragraph of the section
\emph{Dissociation under force from an inverse power law attraction}).
Evans then goes into diffusion constants that depend on the
protein's end to end distance, and I haven't worked out the math
yet.  TODO: clean up.


\subsubsection{Cusp potentials}

\begin{equation}
  E(x) = \frac{1}{2}\kappa_a \p({\frac{x}{x_a}})^2
\end{equation}
(see \citet{evans97} in the text on page 1545, in the first paragraph
of the section \emph{Dissociation under force from a deep harmonic well}).

\section{Double-integral Kramers' theory}

The double-integral form of overdamped Kramers' theory may be too
complex for analytical predictions of unfolding-force histograms.
Rather than testing the entire \sawsim\ simulation (\cref{sec:sawsim}),
we will focus on demonstrating that the Kramers' $k(F)$ evaluations
are working properly.  If the Bell modeled histograms check out, that
gives reasonable support for the $k(F) \rightarrow \text{histogram}$
portion of the simulation.

Looking for analytic solutions to Kramers' $k(F)$, we find that there
are not many available in a closed form.  However, we do have analytic
solutions for unforced $k$ for cusp-like and quartic potentials.

\subsection{Cusp-like potentials}


\subsection{Quartic potentials}