[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
      \fref{figure}{2.a} and \iref{figure}{2.b}, see
      \cref{fig:sawsim:kappa-sawteeth}.\label{tab:sawsim:model}}
  \end{center}
\end{table}

Because the unfolding behaviors 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:plugins}), 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}

\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}

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}.
%
\nomenclature[sr ]{$N_f$}{The number of folded domains in a protein
  chain (\cref{sec:sawsim:results:scaffold}).}
\nomenclature[sr ]{$N_u$}{The number of unfolded domains in a protein
  chain (\cref{sec:sawsim:results:scaffold}).}

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{\frac{F\Delta x_u}{k_B T}}
    = \frac{1}{\rho}\exp{\frac{F-\alpha}{\rho}} \;,
\end{align}
where $N_f$ is the number of folded domain, $\kappa$ is the spring
constant of the cantilever-polymer system, $\kappa v$ is the force
loading rate\footnote{
  \begin{equation}
    \deriv{t}{F} = \deriv{x}{F}\deriv{t}{x} = \kappa v \;.
  \end{equation}
  Alternatively,
  \begin{align}
    F &= \kappa x = \kappa vt \\
    \deriv{t}{F} &= \kappa v \;.
  \end{align}
  See the text before \xref{evans97}{equation}{11} or
  \xref{dudko06}{equation}{4} for similar explanations.
}, 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)$.
We can approximate $\kappa$ as a series of Hookean springs,
\begin{equation}
  \kappa=\p({\frac{1}{\kappa_c}+\frac{N_u}{\kappa_\text{WLC}}})^{-1} \;,
  \label{eq:kappa-system}
\end{equation}
where $\kappa_\text{WLC}$ is the effective spring constant of one
unfolded domain, assumed constant for a particular polymer/cantilever
combination.

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{\frac{F-\alpha}{\rho}
                                           -\exp{\frac{F-\alpha}{\rho}}
                                           } \;.  \label{eq:sawsim:gumbel}
\end{equation}
The distribution has a mode $\alpha$, 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\citep{NIST:gumbel}.  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[sr ]{$r_{uF}$}{Unfolding loading rate (newtons per
  second).}
\nomenclature[sg a ]{$\alpha$}{The mode unfolding force,
  $\alpha\equiv-\rho\ln(N_f k_{u0}\rho/\kappa v)$
  (\cref{eq:sawsim:gumbel}).}
\nomenclature[sg r ]{$\rho$}{The characteristic unfolding force,
  $\rho\equiv k_BT/\Delta x_u$ (\cref{eq:sawsim:gumbel}).}
\nomenclature[sg ce ]{$\gamma_e$}{Euler--Macheroni constant,
  $\gamma_e=0.577\ldots$.}
\nomenclature[sg s ]{$\sigma$}{Standard deviation.  For example,
  $\sigma$ is used as the standard deviation of an unfolding force
  distribution in \cref{eq:sawsim:gumbel}.  Not to be confused with
  the photodiode sensitivity $\sigma_p$.}

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{li00}.  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.

\citet{benedetti11} have since proposed an alternative
parameterization for \cref{eq:kappa-system}, using
\begin{equation}
  \kappa = \p({\frac{1}{\kappa_c}
               + \frac{N_f}{\kappa_f} + \frac{N_u}{\kappa_u}})^{-1}
    \equiv \frac{\kappa'}{1 - A N_f} \;,
  \label{eq:kappa-system-benedetti}
\end{equation}
where $\kappa'$ is the spring constant of the completely unfolded
chain and $A$ is a correction term for the supramolecular scaffold.
This is effectively a first order Taylor expansion for $\kappa^{-1}$
about $N_f=0$, but the remaining analysis is identical.
\begin{align}
  f(N_f) \equiv \kappa^{-1}
    &= \frac{1}{\kappa_c} + \frac{N_f}{\kappa_f} + \frac{N - N_f}{\kappa_u} \\
    &= f(0) + \left.\deriv{N_f}{f}\right|_{N_f=0} N_f + \order{N_f^2} \\
    &\approx \p({\frac{1}{\kappa_c} + \frac{N}{\kappa_u}}) +
             \p({\frac{1}{\kappa_f} - \frac{1}{\kappa_u}}) N_f
  \label{eq:kappa-system-taylor}
\end{align}
In the case where the wormlike chain stiffnesses $\kappa_f$ and
$\kappa_u$ are fairly constant over the unfolding region, there are no
higher order terms and the first order expansion in
\cref{eq:kappa-system-taylor} is exact.  Comparing
\cref{eq:kappa-system-benedetti,eq:kappa-system-taylor}, we see
\begin{align}
  \kappa' &= \frac{1}{\kappa_c} + \frac{N}{\kappa_u} \\
  -\kappa' A &= \frac{1}{\kappa_f} - \frac{1}{\kappa_u} \\
  A &= \frac{\frac{1}{\kappa_u} - \frac{1}{\kappa_f}}
            {\frac{1}{\kappa_c} + \frac{N}{\kappa_u}}
\end{align}
By focusing on the $A=0$ case (\ie~$\kappa_f=\kappa_u$),
\citet{benedetti11} avoid running Monte Carlo simulations when
modeling unfolding histograms.  This simplification does not hold for
our simulated data (\cref{fig:sawsim:order-dep}), but for some
experimental analysis the loss of accuracy may be acceptable in return
for the reduced computational cost.

\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 \;.  \label{eq:sawsim:p_m}
\end{equation}
%
\nomenclature[sr ]{$D_\text{JS}$}{The Jensen--Shannon divergence
  (\cref{eq:sawsim:D_JS}).}
\nomenclature[sr ]{$D_\text{LK}$}{The Kullback--Leibler divergence
  (\cref{eq:sawsim:D_KL}).}
\nomenclature[sr ]{$p_m(i)$}{The symmetrized probability distribution
  used in calculating the Jensen--Shannon divergence
  (\cref{eq:sawsim:D_JS,eq:sawsim:p_m}).}
% 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\citep{NIST:chi-square},
\begin{equation}
  D_{\chi^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$.
%
\nomenclature[sg x2 ]{$\chi^2$}{The chi-squared distribution.}
\nomenclature[sr ]{$D_{\chi^2}$}{Pearson's $\chi^2$ test
  (\cref{eq:sawsim:X2}).}

\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}[commandchars=\\\{\}]
#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 \citetalias{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{minted}[samepage]{console}
$ 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{minted}
That's a bit of a mouthful, so let's break it down.  Without the
\sawsim\ template (\imint{sh}|-f ...|), we can focus on the comparison
options:
\begin{minted}[samepage]{console}
$ sawsim_hist_scan.py ... -r '[1e-5,1e-3,50],[0.1e-9,1e-9,50]' --logx histograms.txt
\end{minted}
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 \imint{sh}|%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[text ]{MPI}{Message passing interface, a parallel
  computing infrastructure.}
\nomenclature[text ]{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.}

\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 (distributed with \sawsim) that compares simulated unfolding
force histograms with analytical histograms for a number of situations
where solving for the analytical histogram is possible.  In the
following subsection, I'll work out the theoretical unfolding force
distribution for a number of tractable cases.  The sawsim test suite
generates simulated unfolding curves for these tractable cases
(e.g. single domain Bell model unfolding with a constant loading
rate), and compares the simulated unfolding force histograms with the
expected theoretical distribution.  The simulated histograms match the
theoretical distributions for each combination of models regardless of
the parameters you feed into the models, so we can be confident that
\sawsim\ correctly implements at those models.

The instantaneous likelyhood of a protein unfolding is given by
$\deriv{F}{N_u}$, and the unfolding histogram is merely this function
discretized over a bin of width $W$\footnote{
  This is similar to \xref{dudko06}{equation}{2}, remembering that
  $\dot{F}=\kappa v$, that their probability density is not a
  histogram ($W=1$), and that their probability density function is
  normalized to $N=1$
}.
\begin{equation}
  h(F) \equiv \deriv{\text{bin}}{F}
    = \deriv{F}{N_u} \cdot \deriv{\text{bin}}{F}
    = W \deriv{F}{N_u}
    = -W \deriv{F}{N_f}
    = -W \deriv{t}{N_f} \deriv{F}{t}
    = \frac{W}{\kappa v} N_f k_u \label{eq:unfold:hist}
\end{equation}
Solving for theoretical histograms is merely a question of taking your
chosen $k_u$, solving for $N_f(F)$, and plugging into
\cref{eq:unfold:hist}.  We can also make a bit of progress solving for
$N_f$ in terms of $k_u$ as follows:
\begin{align}
  k_u &\equiv -\frac{1}{N_f} \deriv{t}{N_f} \\
  -k_u \dd t \cdot \deriv{t}{F} &= \frac{\dd N_f}{N_f} \\
  \frac{-1}{\kappa v} \integral{0}{F}{F'}{k_0(F')}
    &= \left. \ln(N_f(F')) \right|_0^F
    = \ln\p({\frac{N_f(F)}{N_f(0)}})
    = \ln\p({\frac{N_f(F)}{N}}) \\
  N_f(F) &= N\exp{\frac{-1}{\kappa v}\integral{0}{F}{F'}{k_u(F')}} \;,
  \label{eq:N_f}
\end{align}
where $N_f(0) = N$ because all the domains are initially folded.
%
\nomenclature[sr ]{$W$}{Bin width of an unfolding force histogram
  (\cref{eq:unfold:hist}).}

\subsubsection{Constant unfolding rate}

In the extremely weak tension regime, the protein's unfolding rate is
independent of tension, so we can simplify \cref{eq:N_f} and plug into
\cref{eq:unfold:hist}.
\begin{align}
  N_f &= N\exp{\frac{-1}{\kappa v}\integral{0}{F}{F'}{\colA{k_u(F')}}}
     = N\exp{\frac{-\colA{k_{u0}}}{\kappa v}\colB{\integral{0}{F}{F'}{}}}
     = N\exp{\frac{-k_{u0} \colB{F}}{\kappa v}} \\
  h(F) &= \frac{W}{\kappa v} N_f k_u
     = \frac{W k_{u0} N}{\kappa v} \exp{\frac{-k_{u0} F}{\kappa v}} \;.
\end{align}
A constant unfolding-rate/hazard-function gives exponential decay.
This is not an 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.

\subsubsection{Bell model}

Stepping up the intensity a bit, we come to Bell's model for unfolding
(\cref{sec:sawsim:rate:bell}).  We can simplify the following
calculation by parametrizing with the characteristic force $\rho$
defined in \cref{sec:sawsim:results:scaffold} and the similar
single-domain mode $\alpha'\equiv-\rho\ln(k_{u0}\rho/\kappa v)$.  With
these substitutions, \cref{eq:sawsim:bell} becomes
\begin{equation}
  k_u = k_{u0} \exp{\frac{F}{\rho}} \;.
\end{equation}
The unfolding histogram is then given via \cref{eq:N_f,eq:unfold:hist}.
\begin{align}
  N_f &= N\exp{\frac{-1}{\kappa v}\integral{0}{F}{F'}{\colA{k_u}}}
     = N\exp{\frac{-1}{\kappa v}
       \integral{0}{F}{F'}{\colAB{k_{u0}}{\colA{\exp{\frac{F'}{\rho}}}}}}
     = N\exp{\frac{-\colB{k_{u0}}}{\kappa v}
       \colA{\integral{0}{F}{F'}{\exp{\frac{F'}{\rho}}}}}
     = N\exp{\frac{\colB{-}k_{u0}\colA{\rho}}{\kappa v}
       \colAB{\p({\exp{\frac{F}{\rho}}-1})}} \\
     &= N\exp{\colA{\frac{k_{u0}\rho}{\kappa v}}
       \colB{\p({1 - {\exp{\frac{F}{\rho}}}})}}
     = N\exp{\colAB{\exp{\frac{-\alpha'}{\rho}}}
       \colB{\p({1 - {\exp{\frac{F}{\rho}}}})}}
     = N\exp{\colB{\exp{\frac{-\alpha'}{\rho}} -
       \exp{\frac{F-\alpha'}{\rho}}}} \\
  h(F) &= \frac{W}{\kappa v} \colA{N_f} \colB{k_u}
     = \frac{W}{\kappa v}
       \colA{N\exp{\exp{\frac{-\alpha'}{\rho}} - \exp{\frac{F-\alpha'}{\rho}}}}
       \colB{k_{u0}\exp{\frac{F}{\rho}}}
     = \frac{W N \colAB{k_{u0}}}{\colA{\kappa v}}
       \exp{\colB{\frac{F}{\rho}} - \exp{\frac{F-\alpha'}{\rho}} +
         \exp{\frac{-\alpha'}{\rho}}} \\
     &= \frac{W N}{\colA{\rho}}
       \exp{\frac{F \colA{-\alpha'}}{\rho} - \exp{\frac{F-\alpha'}{\rho}} +
         \colB{\exp{\frac{-\alpha'}{\rho}}}}
     = \frac{W N}{\rho}
       \exp{\frac{F-\alpha'}{\rho} - \exp{\frac{F-\alpha'}{\rho}}}
       \colB{\exp{\exp{\frac{-\alpha'}{\rho}}}} \\
     &= \frac{W N \exp{\exp{\frac{-\alpha'}{\rho}}}}{\rho}
       \exp{\frac{F-\alpha'}{\rho} - \exp{\frac{F-\alpha'}{\rho}}}
  \label{eq:unfold:bell_pdf}
\end{align}
which matches \cref{eq:sawsim:gumbel} except for a constant
prefactor due to the range\footnote{
  The Gumbel distribution in \cref{eq:sawsim:gumbel} is normalized for
  the range $-\infty < F < \infty$, but \cref{eq:unfold:bell_pdf} is
  normalized for the range $0 \le F < \infty$.  This distinction will
  alter the analytical mean and variance listed after
  \cref{eq:sawsim:gumbel}, but with the experimental unfolding
  histograms showing few zero-force unfolding events, the effective
  difference will be negligible.
}.
%
\nomenclature[sg a' ]{$\alpha'$}{The mode unfolding force for a single
  folded domain, $\alpha'\equiv-\rho\ln(k_{u0}\rho/\kappa v)$
  (\cref{eq:unfold:bell_pdf}).}

\subsubsection{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}
  k_u = \frac{D}{l_b l_{ts}} \cdot \exp{\frac{-U_b(F)}{k_B T}} \;,
    \label{eq:kramers-saddle}
\end{equation}
where $U_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.
%
\nomenclature[sr ]{$U_b(F)$}{The barrier energy as a function of force
  (\cref{eq:kramers-saddle}).}
\nomenclature[sr ]{$l_b$}{The characteristic length of the bound state
  $l_b \equiv 1/\rho_b$ (\cref{eq:kramers-saddle}).}
\nomenclature[sg r_b ]{$\rho_b$}{The density of states in the bound
  state (\cref{eq:kramers-saddle}).}
\nomenclature[sr ]{$l_{ts}$}{The characteristic length of the
  transition state (\cref{eq:kramers-saddle}).}

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

\section{Review of current research}

There is a long history of protein unfolding and unbinding
simulations.  Early work by \citet{grubmuller96} and
\citet{izrailev97} focused on molecular dynamics (MD) simulations of
receptor-ligand breakage.  Around the same time, \citet{evans97}
introduced a Monte Carlo Kramers' simulation in the context of
receptor-ligand breakage.  The approach pioneered by \citet{evans97}
was used as the basis for analysis of the initial protein unfolding
experiments\citep{rief97a}.  However, none of these earlier
implementations were open source, which made it difficult to reuse or
validate their results.
%
\nomenclature[text ]{MD}{Molecular dynamics simulation.  Simulate the
  physical motion of atoms and molecules by numerically solving
  Newton's equations.}

Within the Monte Carlo simulation approach, there are two main models
for 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 its
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.

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

There have been some tangential attempts towards even fancier models:
\citet{dudko03} attempted to reduce the restrictions of the
single-unfolding-path model and \citet{hyeon03} attempted to measure
the local roughness using temperature dependent unfolding.  However,
further work on these lines has been slow, because the Bell model fits
the data well despite its simplicity.  For more complicated models to
gain ground, we need larger, more detailed datasets that expose
features which the Bell model doesn't capture.