sawsim/methods.tex: Plug sec:salt from sec:sawsim:rate:other summary

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

\section{Methods}
\label{sec:sawsim:methods}

% simulation overview
In simulating the mechanical unfolding process, a force curve is
generated by calculating the amount of cantilever bending as the
substrate surface moves away from the tip.  The cantilever bending is
obtained by balancing the tension in the protein polymer and the
Hookean force of the bent cantilever.  The unfolding probability of
the protein molecules in the polymer is then calculated for that
tension, and whether an unfolding event occurs is determined according
to a Monte Carlo method.  The simulation was implemented in
\citetalias{kernighan88}, and there are a number of
\citetalias{python} modules to facilitate running several simulations
in parallel\footnote{
  Source code available at \url{http://blog.tremily.us/posts/sawsim/}.
}.

In the following sections, we'll discuss models used to determine the
tension of a chain composed of several types of ``domains'' (e.g.~one
cantilever, three folded I27 domains, and seven unfolded I27 domains)
(\cref{sec:sawsim:tension}).  We'll also work through a number of
models for calculating the probability that a domain will transition
from one state (e.g.~folded I27) to another (e.g.~unfolded I27)
(\cref{sec:sawsim:rate}).

\subsection{Modeling polymer tension}
\label{sec:sawsim:tension}

% introduce domains and groups.
The fundamental abstraction of the simulation is the ``domain'', which
represents a discrete chunk of the flexible chain between the
substrate and the cantilever holder (\cref{fig:sawsim:domain-chain}).
Each of these domains is assigned a particular state; for example, the
domain representing the cantilever is assigned to the ``cantilever''
state, and the domains representing protein molecules are assigned to
either the ``folded'' or the ``unfolded'' state.  When balancing the
tension along the chain, we assume that the spatial order of domains
along the chain is irrelevant\citep{li00}, so the domains can be
rearranged and grouped by state (\cref{fig:sawsim:domain-states}).  To
determine the tension in the chain and the amount of cantilever
bending when $N$ states are populated, a system of $N+1$ equations
with $N+1$ unknowns must be solved
\begin{align}
  F_i(x_i) &= F_t        \label{eq:sawsim:tension-balance} \\
  \sum_i x_i &= x_t \;,  \label{eq:sawsim:x-total}
\end{align}
where $F$ are tensions, $x$ are extensions, and the subscripts $i$ and
$t$ represent a particular state group and the total chain
respectively (\cref{fig:unfolding-schematic}).  $F(x_t)$ may be
computed from this system of equations using any multi-dimensional
root-finding algorithm.

\begin{figure}
  \begin{center}
    \subfloat[][]{\label{fig:sawsim:domain-chain}
      \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
        \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
        \node[state] (A)              {domain 1};
        \node[state] (B) [below of=A] {domain 2};
        \node[state] (C) [below of=B] {.~.~.};
        \node[state] (D) [below of=C] {domain $N$};
        \node        (S) [below of=D] {Surface};
        \node        (E) [above of=A] {};

        \path[-] (A) edge (B)
          (B) edge node [right] {Tension} (C)
          (C) edge (D)
          (D) edge (S);
        \path[->,green] (A) edge node [right,black] {Extension} (E);
      \end{tikzpicture}}
    \hspace{.25in}%
    \subfloat[][]{\label{fig:sawsim:domain-states}
      \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
        \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
        \node[state] (A)              {cantilever};
        \node[state] (C) [below of=A] {transition};
        \node[state] (B) [left of=C]  {folded};
        \node[state] (D) [right of=C] {unfolded};

        \path (B) edge [bend left] node [above] {$k_1$} (C)
              (C) edge [bend left] node [below] {$k_1'$} (B)
                  edge [bend left] node [above] {$k_2$} (D)
              (D) edge [bend left] node [below] {$k_2'$} (C);
      \end{tikzpicture}}
    \caption{\protect\subref{fig:sawsim:domain-chain} Extending a
      chain of domains.  One end of the chain is fixed, while the
      other is extended at a constant speed.  The domains are coupled
      with rigid linkers, so the domains themselves must stretch to
      accomodate the extension.  Compare with
      \cref{fig:unfolding-schematic}.
      \protect\subref{fig:sawsim:domain-states} Each domain exists in
      a discrete state.  At each timestep, it may transition into
      another state following a user-defined state matrix such as this
      one, showing a metastable transition state and an explicit
      ``cantilever'' domain.\label{fig:sawsim:domains}}
  \end{center}
\end{figure}

\subsubsection{Hooke's law}
\label{sec:sawsim:tension:hooke}

% introduce particular models, and mention parameter aggregation
Inside this framework, we chose a particular extension model
$F_i(x_i)$ for each domain state.  Cantilever elasticity is described
by Hooke's law, which gives
\index{Hooke's law}
\begin{equation}
  F = \kappa_c x_c \;, \label{eq:sawsim:hooke}
\end{equation}
where $\kappa_c$ is the bending spring constant and $x_c$ is the
deflection of the cantilever (\cref{fig:unfolding-schematic}).

\subsubsection{Wormlike chains}
\label{sec:sawsim:tension:wlc}

\index{Wormlike chains}
Unfolded domains are modeled as wormlike chains
(WLCs)\citep{rief97a,carrion-vazquez00}, which treat the unfolded
polymer as an elastic rod of persistence length $p$ and contour length
$L$ (\cref{fig:wlc}).  The relationship between tension $F$ and
extension (end-to-end distance) $x$ is given by Bustamante's
interpolation formula\citep{marko95,bustamante94}.
\nomenclature{WLC}{Wormlike chain, an entropic spring model}
\nomenclature{$p$}{Persistence length of a wormlike chain
  (\cref{eq:sawsim:wlc})).}
\nomenclature{$L$}{Contour length in a polymer tension model
  (\cref{eq:sawsim:wlc,eq:sawsim:fjc})}
\begin{equation}
  F_\text{WLC}(x,p,L) = \frac{k_B T}{p}
      \p[{  \frac{1}{4}\p({ \frac{1}{(1-x/L)^2} - 1 })
            + \frac{x}{L}  }] \;,
      \label{eq:sawsim:wlc}
\end{equation}
where $p$ is the persistence length.  This interpolation formula is
accurate to within 7\% of the exact $F_\text{WLC}$ for
$F_\text{WLC}\approx k_B T/p_u$\citep{marko95}.  Because most unfolded
proteins studied have persistence lengths on the order of the size of
an amino acid
($p_u\approx3.8\U{\AA}$\citep{rief97a,carrion-vazquez99b,carrion-vazquez00}),
this characteristic force works out to be around $11\U{pN}$.  Most
proteins studied using force spectroscopy have unfolding forces in the
hundreds of piconewtons, by which point the interpolation formula is
in it's more accurate high-extension regime.
\nomenclature{\AA}{{\AA}ngstr{\"o}m, a unit of length.
  $1\U{\AA}=1\E{-10}\U{m}$}

For chain with $N_u$ unfolded domains sharing a persistence length
$p_u$ and per-domain contour lengths $L_{u1}$, the tension of the WLC
is determine by summing the contour lengths
\begin{equation}
  F(x, p_u, L_u, N_u) = F_\text{WLC}(x, p_u, N_u L_{u1}) \;.
  \label{eq:sawsim:multi-wlc}
\end{equation}

\begin{figure}
  \begin{center}
    \subfloat[][]{\asyinclude{figures/schematic/wlc-model}%
      \label{fig:wlc-model}}
    \hspace{.25in}%
    \subfloat[][]{\asyinclude{figures/schematic/wlc-extension}%
      \label{fig:wlc-extension}}
    \caption{\protect\subref{fig:wlc-model} The wormlike chain models
      a polymer as an elastic rod with persistence length $p$ and
      contour length $L$.
      \protect\subref{fig:wlc-extension} Force vs.~extension for a WLC
      using Bustamante's interpolation formula.\label{fig:wlc}}
  \end{center}
\end{figure}

\subsubsection{Folded domains}
\label{sec:sawsim:tension:folded}

Short chains of folded proteins, however, are not easily described by
polymer models.  Several studies have used WLC and FJC models to fit
the elastic properties of the modular protein
titin\citep{granzier97,linke98a},
% TODO: check it really is folded domains \& bulk titin
but native titin contains hundreds of folded and unfolded domains.
For the short protein polymers common in mechanical unfolding
experiments (\cref{sec:polymer-synthesis}), the cantilever dominates
the elasticity of the polymer-cantilever system before any protein
molecules unfold.  After the first unfolding event occurs, the
unfolded portion of the chain is already longer and softer than the
sum of all the remaining folded domains, and dominates the elasticity
of the whole chain.  Therefore, the details of the tension model
chosen for the folded domains has negligible effect on the unfolding
forces (\cref{eq:sawsim:x-total}), which was also suggested by
\citet{staple08}.  Force curves simulated using different models to
describe the folded domains yielded almost identical unfolding force
distributions (data not shown).% TODO: show data

As an alternative to modeling the folded domains explicitly or
ignoring them completely, another approach is to subtract the
end-to-end length of the folded protein from the contour length of the
unfolded protein to create an effective contour length for the
unfolding\citep{carrion-vazquez99b}.  This effectively models the
folded domains as WLCs with the same persistence length as the
unfolded domains.

%The chain of $N_f$ folded domains is modeled as a string, free to
%assume any extension up to some fixed contour length $L_f=N_fL_{f1}$
%\begin{equation}
%  F = \begin{cases}
%         0 & \text{if $x_f<L_f$} \;, \\
%         \infty & \text{if $x_f>L_f$} \;,
%      \end{cases}  \label{eq:sawsim:piston}
%\end{equation}
%where $L_{f1}$ is the separation of the two linking points of a folded
%domain, and $x_f$ is the end-to-end length of the chain of folded
%domains.  In this model, any non-zero tension will fully extend these
%folded domains.  Because the range of possible extensions for folded
%domains is so short, the contribution of the folded domains to the
%elastic behavior of the polymer-cantilever system is relatively
%insignificant.

\subsubsection{Other models}
\label{sec:sawsim:tension:other}

\index{FJC}
\index{Kuhn length}
The unfolded polypeptide chain has been shown to follow the WLC model
quite well\citep{rief97a} (\cref{sec:sawsim:tension:wlc}), though
other polymer models have been tried.  One alternative is the
freely-jointed chain
(FJC)\citep{kellermayer97,linke98a,janshoff00,verdier70}, which models
the polymer as a series of $N$ rigid links, each of length $l$ (the
Kuhn length), which are free to rotate about their joints
(\cref{fig:fjc}).
\index{Langevin function}
\begin{equation}
  F_\text{FJC}(x,l,L) = \frac{k_B T}{l} \Langevin^{-1}\p({\frac{x}{L}}) \;,
  \label{eq:sawsim:fjc}
\end{equation}
where $L=Nl$ is the total length of the chain, and
$\Langevin(\alpha)\equiv\coth{\alpha}-\frac{1}{\alpha}$ is the
Langevin function\citep{hatfield99}.
%
\nomenclature{FJC}{Freely-jointed chain, an entropic spring model
  (\cref{eq:sawsim:fjc}).}
\nomenclature{$\Langevin$}{The Langevin function,
  $\Langevin(\alpha)\equiv\coth{\alpha}-\frac{1}{\alpha}$}
\nomenclature{$\coth$}{Hyperbolic cotangent,
  \begin{equation}
    \coth(x) = \frac{\exp{x} + \exp{-x}}{\exp{x} - \exp{-x}} \;.
  \end{equation}
}
\nomenclature{$l$}{Kuhn length in the freely-jointed chain
  (\cref{fig:fjc-model,eq:sawsim:fjc}).}

\begin{figure}
  \begin{center}
    \subfloat[][]{\asyinclude{figures/schematic/fjc-model}%
      \label{fig:fjc-model}}
    \hspace{.25in}%
    \subfloat[][]{\asyinclude{figures/schematic/fjc-extension}%
      \label{fig:fjc-extension}}
    \caption{\protect\subref{fig:fjc-model} The freely-jointed chain
      models the polymer as a series of $N$ rigid links, each of
      length $l$, which are free to rotate about their joints.  Each
      polymer state is a random walk, and the density of states for a
      given end-to-end distance is determined by the number of random
      walks that have such an end-to-end distance.
      \protect\subref{fig:fjc-extension} Force vs.~extension for a
      hundred-segment FJC.  The WLC extension curve (with $p=l$) is
      shown as a dashed line for comparison.\label{fig:fjc}}
  \end{center}
\end{figure}

More exotic models such as elastic WLCs\citep{janshoff00,puchner08},
elastic FJCs\citep{fisher99a,janshoff00}, and freely rotating
chains\citep{puchner08} (FRCs) have also been used to model DNA and
polysaccharides, but are rarely used to model the relatively short and
inextensible synthetic proteins used in force spectroscopy.
\nomenclature{FRC}{Freely-rotating chain (like the FJC, except that
  the bond angles are fixed.  The torsional angles are not
  restricted)}


\subsubsection{Assumptions}
\label{sec:sawsim:tension:assumptions}

% address assumptions & caveats
In the simulation, the protein polymer is assumed to be stretched in
the direction perpendicular to the substrate surface, which is a good
approximation in most experimental situations, because the unfolded
length of a protein molecule is much larger than that of the folded
form.  Therefore, after one molecule unfolds, the polymer becomes much
longer and the angle between the polymer and the surface approaches
$90$ degrees\citep{carrion-vazquez00}.

The joints between domain groups are also assumed to lie along a line
between the surface tether point and the position of the tip
(\cref{eq:sawsim:x-total} is scalar, not vector, addition).  The
effects of this assumption are also minimized due to greater length of
the unfolded domain compared with the other domains (folded proteins
and cantilever deflection).  For example, a $0.050\U{N/m}$ cantilever
under $200\U{pN}$ of tension will bend $x_c=F/\kappa_c=4\U{nm}$.  The
entire end-to-end length of folded domains such as I27 are also around
$5\U{nm}$ (\cref{fig:I27}).  A single unfolded I27, with its 89 amino
acids\citep{improta96}, should have an unfolded contour length of
$89\U{aa}\cdot0.38\U{nm}=33.8\U{nm}$, equivalent to a cantilever and
five folded domains.

\subsubsection{Velocity-clamp example}
\label{sec:sawsim:velocity-clamp}

% introduce constant velocity and walk through explicit example pull
Consider an experiment pulling a polymer with $N$ identical protein
domains at a constant speed.  At the start of an experiment, the
chain is unstretched ($x_t=0$), which means all the domains are
unstretched, the cantilever is undeflected, and the tip is in contact
with the surface.  There is one domain in the cantilever state, $N$ in
the folded state, and none in the unfolded state.  As the surface
moves away from the tip at a constant speed $v$, the chain becomes
more extended (\cref{fig:unfolding-schematic}), such that
\begin{equation}
  x_t = \sum_i x_i = vt  \label{eq:sawsim:const-v} \;.
\end{equation}
The simulation assumes that the pulling takes discrete steps in space
and treats $x_t$ as constant over the duration of one time step
$\Delta t$.  Because of the adaptive time steps discussed in
\cref{sec:sawsim:timesteps}, the space steps $\Delta x_t = v\Delta t$
may have different sizes, but each step will be ``small''.  At each
step, the total extension is calculated using
\cref{eq:sawsim:const-v}, and the tension $F(x_t=vt)$ is determined by
numerically solving \cref{eq:sawsim:tension-balance,eq:sawsim:x-total}
using the models \cref{eq:sawsim:hooke,eq:sawsim:wlc}
%,eq:sawsim:piston}
for known values of the parameters in the various states $(N_u, N_f,
v, \kappa_c,
% L_{f1}, 
L_{u1}, p_u)$.  When one of the molecules in the
polymer unfolds (\cref{sec:sawsim:rate}), there will be
one domain in the unfolded state and $N-1$ in the folded state.  In
the next step, a newly balanced tension between the cantilever and the
polymer is determined by solving for $F(x_t)$ as discussed above, but
with the total extension $x_t$ incremented by $v\Delta t$ and the new
unfolded contour length $L_{u1}$ and folded contour length
$(N-1)L_{f1}$.  The sudden lengthening of the polymer chain results in
a corresponding abrupt drop in the force, leading to the formation of
one sawtooth in the force curve.  As the pulling continues and more
domains unfold, force curves with a series of sawteeth are generated
(\cref{fig:sawsim:sim-sawtooth}).
%
\nomenclature{$v$}{Cantilever retraction speed in velocity-clamp
  unfolding experiments.}

\subsubsection{Equlibration time scales}
\label{sec:sawsim:timescales}

The tension calculation assumes an equilibrated chain, so
consideration must be given to the chain's relaxation time, which
should be short compared to the loading time scale.  The relaxation
time for a WLC\index{WLC!relaxation time} is given by
\begin{equation}
 \tau \approx \eta \frac{k_BT p}{F^2}
%   < \eta \frac{k_BT p}{(k_BTx/pL)^2}
%     Note:  <  because  F > k_BTx/pL
%   < \frac{\eta p^3 L^2}{k_BT x^2}
%   < \frac{\eta p^2 L}{k_BT} % for x/L > \sqrt{p/L}
   \;, \label{eq:sawsim:tau-wlc}
\end{equation}
where $\eta$ is the dynamic viscosity, $F$ is the tension, and $p$ is
the persistence length\citep{evans99}.  For forces greater than
$1\U{pN}$, with $\eta_\text{water}/k_BT=2.45\E{-10}\U{s/nm$^3$}$,
$\tau<2\U{ns}$ for the protein polymer used in the simulation.
%                  eta/(k_BT)        * (k_B     *T  )**2 * p      / F**2
% python -c 'print(2.45e-10*(1e9)**3 * (1.38e-23*300)**2 * 0.4e-9 / (1e-12)**2)'
%                  s/m**3            * (J/K     *K  )**2 * m      / N**2         = s
% 1.68...e-09
Therefore, the polymer chain is equilibrated almost instantaneously
within a time step, which is on the order of tens of microseconds.
The relaxation time of the cantilever can be determined by measuring
the cantilever deflection induced by liquid motion and fitting the
time dependence of the deflection to an exponential
function\citep{jones05}.  For a $200\U{$\mu$m}$ rectangular cantilever
with a bending spring constant of $20\U{pN/nm}$, the measured
relaxation time in water is $\sim50\U{$\mu$/s}$ (data not shown).
% TODO: show data
This relatively large relaxation time constant makes the cantilever
act as a low-pass filter and also causes a lag in the force
measurement.
%
\nomenclature{$\eta$}{Dynamic viscocity (\cref{eq:sawsim:tau-wlc}).}

\subsection{Unfolding protein molecules by force}
\label{sec:sawsim:rate}

In the previous section, we discussed methods for calculating the
tension of a chain composed of several domains in series
(\cref{fig:sawsim:domain-chain}).  Those methods allow us to calculate
the tension of the chain for any given extension.  We use that tension
to calculate transition rates between states
(\cref{fig:sawsim:domain-states}).  In this section, we'll introduce
the Bell model for unfolding (\cref{sec:sawsim:rate:bell}) and mention
a few more exotic models.  We'll wrap up by pointing out some of the
approximations and assumptions we make when we use these simple models
(\cref{sec:sawsim:rate:assumptions}).

\subsubsection{Bell model}
\label{sec:sawsim:rate:bell}

\index{Bell model}
% introduce Bell, probability calculations, and MC comparison
According to the theory developed by \citet{bell78} and extended by
\citet{evans99}, an external stretching force $F$ increases the
unfolding rate constant of a protein molecule\footnote{
  Also in \xref{hummer03}{equation}{1}, the first paragraphs of
  \citet{dudko06} and \citet{dudko07}, and many other SMFS articles.
}
\index{Bell model}
\begin{equation}
  k_u = k_{u0} \exp{\frac{F\Delta x_u}{k_B T}} \;, \label{eq:sawsim:bell}
\end{equation}
where $k_{u0}$ is the unfolding rate in the absence of an external
force, and $\Delta x_u$ is the distance between the native state and
the transition state along the pulling direction.
%
\nomenclature{$\exp{x}$}{Exponential function,
  \begin{equation}
    \exp{x} = \sum_{n=0}^{\infty} \frac{x^n}{n"!}
      = 1 + x + \frac{x^2}{2"!} + \ldots \;.
  \end{equation}
}
\nomenclature{$e$}{Euler's number, $e=2.718\ldots$.}

\begin{figure}
  \asyinclude{figures/schematic/landscape-bell}
  \caption{Energy landscape schematic for Bell model unfolding
    (\cref{eq:sawsim:bell}), which models folded domains as two-state
    systems parameterized by an unforced unfolding rate $k_{u0}$ and a
    distance $\Delta x$ between the folded and transition
    states.\label{fig:bell-landscape}}
\end{figure}

\subsubsection{Monte Carlo transitions}
\label{sec:sawsim:monte-carlo}

We can use the Bell model (or other models, see
\cref{sec:sawsim:rate:other}) to calculate the unfolding rate $k_u$ at
a given force for a single domain.  The probability for that single
protein domain to unfold under applied force is
\begin{equation}
  P_1 = k_u \Delta t \;,  \label{eq:sawsim:prob-one}
\end{equation}
where $\Delta t$ is the time duration for each pulling step, over
which $F$ is constant.  This expression is accurate for $P_1 \ll 1$.
From the binomial distribution, the probability of at least one of a
group of $N_f$ identical domains to unfold in a given time step is
\begin{equation}
  P = 1 - (1-P_1)^{N_f} \approx N_fP_1 \;,  \label{eq:sawsim:prob-n}
\end{equation}
where the approximation is valid when $N_fP_1 \ll 1$.
%
\nomenclature{$k$}{Rate constant for general state transitions
  (inverse seconds)}
\nomenclature{$k_u$}{Unfolding rate constant}
\nomenclature{$k_{u0}$}{Unforced unfolding rate constant}
\nomenclature{$\Delta x_u$}{Distance between a domain's native state
  and the transition state along the pulling direction.}
\nomenclature{$P$}{Probability for at least one domain unfolding in a
  given time step (\cref{eq:sawsim:prob-n}).}

\begin{figure}
  \asyinclude{figures/schematic/monte-carlo}
  \caption{Once the unfolding probability has been caculated, we need
    to determine whether or not a domain should unfold.  We do this by
    generating a random number, and comparing that number to the
    unfolding probability $P$.  The random number determines which of
    the possible paths we should follow for the current simulation.
    Such ``statistical sampling'' is the hallmark of the Monte Carlo
    approach\citep{metropolis87}.  This cartoon translates the idea
    into the more familiar doors (possible paths) and dice (random
    numbers).%
    \label{fig:monte-carlo}}
\end{figure}

To determine if an unfolding event occurs in a particular time step,
the probability calculated using \cref{eq:sawsim:prob-n} is compared
with a randomly generated number uniformly distributed between $0$ and
$1$ (\cref{fig:monte-carlo}).  If $P$ is bigger than the random
number, a domain unfolds, changing the population of each tension
state, and a new balance between the polymer and the cantilever is
determined.  If no unfolding event occurs the pulling continues and
the unfolding probability is calculated again in the next step at a
higher force.  When all the molecules in the polymer have unfolded,
the pulling continues until a pre-determined force level is reached,
where the polymer is assumed to detach from one of the tethering
surfaces.  The cantilever deflection becomes zero after this point.

\subsubsection{Other models}
\label{sec:sawsim:rate:other}

Although the Bell model (\cref{eq:sawsim:bell}) is the most widely
used unfolding model due to its simplicity and its applicability to
various biopolymers\citep{rief98}, other theoretical models have been
proposed to interpret mechanical unfolding data.  For example,
\citet{walton08} uses a stiffness-corrected Bell model.
\citet{schlierf06} used the mechanical unfolding data of the protein
ddFLN4 to demonstrate that Kramers' diffusion model (in the
spatial-diffusion-limited case, a.k.a. the Smoluchowski
limit)\citep{kramers40,hanggi90,evans97,shillcock98,vanKampen07} fit
the measured unfolding force data better than the Bell model for
proteins with broad free energy barriers.
\index{Kramers' model}
\index{Diffusion coefficient}
\index{Free energy}
\index{Unfolding coordinate}
\begin{equation}
  \frac{1}{k_u}
    = \frac{1}{D}
      \integral{-\infty}{\infty}{x}{%
        \exp{\frac{U_F(x)}{k_B T}}
        \integral{-\infty}{x}{x'}{%
          \exp{\frac{-U_F(x')}{k_B T}}}} \;,
  \label{eq:kramers}
\end{equation}
where $D$ is the diffusion coefficient and $U_F(x)$ is the free energy
along the unfolding cordinate $x$ (\cref{fig:landscape:kramers}).
\nomenclature{$D$}{Diffusion coefficient (square meters per second)}
\nomenclature{$U_F(x)$}{Protein free energy along the unfolding
  coordinate $x$ (joules)}

\begin{figure}
  \begin{center}
    \subfloat[][]{\asyinclude{figures/schematic/landscape}%
      \label{fig:landscape}}
    % \hspace{.25in}%
    \subfloat[][]{\asyinclude{figures/schematic/kramers-integrand}%
      \label{fig:kramers:integrand}}
    \caption{\protect\subref{fig:landscape} Energy landscape schematic
      for Kramers integration (compare with
      \cref{fig:bell-landscape}).
      \protect\subref{fig:kramers:integrand} A map of the magnitude of
      Kramers' integrand, with black lines tracing the integration
      region.  The bulk of the contribution to the integral comes from
      the bump in the upper left, with $x$ near the boundary and $x'$
      near the folded state.  This is why you can calculate a close
      approximation to this integral by restricting the integration to
      $x_\text{min}$ and $x_\text{max}$, located a few $k_B T$ beyond
      the folded and transition states respectively.  The restricted
      integral is much easier to calculate numerically than one bound
      by $\pm\infty$.
      (\cref{eq:kramers}).\label{fig:landscape:kramers}}
  \end{center}
\end{figure}

When you are simulating the double integral form of Kramers' model,
you obviously need to parameterize $U_F(x)$ somehow.  There has not
been much research done in this direction, but \citet{schlierf06} used
cubic splines with 15 variable knots.  \citet{shillcock98} used a
cubic free energy with variable coefficients.  The amount of
information you can extract from fitting such a model to your data is
limitless, but you run the risk of over-specifying if you add too many
parameters when you're fitting noisy data.

There are alternative formulations of Kramers' model besides the full
double integral approach.  You can use a Gaussian steepest-descent
approximation (a.k.a. stationary phase method or saddle-point method)
to reduce the integral to a formula that only depends on the free
energy landscape via the curvature $\npderiv{2}{x}{U_F}$ evaluated at
the folded state and transition state\citep{hanggi90}.  This approach
makes sense for sufficiently sharp folded and transition states, where
these two measurements will capture the shape of the large-integrand
region (\cref{fig:kramers:integrand}).  The steepest-descent
formulation has less to say about the underlying energy landscape, but
it may be more robust in the face of noisy data.

How to choose which unfolding model to use?  For proteins with
relatively narrow folded and transition states, the Bell model
provides a good approximation, and it is the model used by the vast
majority of earlier work in the field.  I will use the Bell model in
my analysis of ion-dependent unfolding (\cref{sec:salt}), but
analyzing my unfolding data with a different transition rate model is
just a matter of changing some command line options and rerunning the
\sawsim\ simulations.

\subsubsection{Assumptions}
\label{sec:sawsim:rate:assumptions}

The interactions between different parts of the polymer and between
the chain and the surface (except at the tethering points) are often
ignored.  This is usually reasonable since these interactions should
not make substantial contributions to the force curve at the force
levels of interest, where the polymer is in a relatively extended
conformation.  However, \citet{li00} showed that while the unfolding
properties of I27 are not effected by I28 flankers, I28 \emph{is}
stabilized by neighboring I27.  The unforced Bell model unfolding rate
for I28 in (I28)\textsubscript{8} was $2.8\E{-5}\U{s$^{-1}$}$, while
in (I27-I28)\textsubscript{4} it dropped to
$2.5\E{-6}\U{s$^{-1}$}$\citep{li00}.

\subsection{Choosing the simulation time steps}
\label{sec:sawsim:timesteps}

The demands on the time step vary throughout a simulated pull due to
the non-linear elasticity of the polymer.  Within a specified time
duration (or pulling distance), the force change is small at low force
levels and large at high force levels.  To be efficient, the
simulation algorithm adapts the time step to keep the time steps large
where large time steps have little effect, while shrinking the time
step where smaller steps are necessary.

Within each time step, the total chain extension $x_t$ is treated as a
constant and a force balance is reached very quickly among the various
domains (\cref{sec:sawsim:timescales}).  This force is used to
determine the unfolding probability
(\cref{eq:sawsim:prob-one,eq:sawsim:prob-n}), which determines the
domain state populations in the next time step.  Therefore, the chain
tension must not change appreciably over the course of the time step
($\Delta F < 1\U{pN}$), and the unfolding probability is only
calculated once for the entire step.  The time step must also be short
enough that the probability of unfolding in a single time step is low
($P<10^{-3}$).  Besides ensuring that the approximations made in
\cref{eq:sawsim:prob-one,eq:sawsim:prob-n} are valid, this restriction
makes time steps which should have multiple unfoldings in a single
time step highly unlikely.  Experimentally measured unfolding are
temporally separated, because the unfolding transition is
characterized by multiple, Markovian attempts over a large energy
barrier, where the probability of crossing the barrier in a single
attempt is very low.  A successful attempt quickly extends the chain
contour length, reducing the tension, dramatically reducing the
likelihood of a second escape in that time step.  The time step used
is recalculated for each step so that both of these criteria are
satisfied.