Remove silly tex/ directory.

[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
C\footnote{Source code available at
  \url{http://www.physics.drexel.edu/~wking/sawsim/}}.

\subsection{Generating force curves}
\label{sec:sawsim:methods-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.  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}, and therefore, the domains can be rearranged
and grouped by state.  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}).  From this $F(x_t)$ may be
computed using any multi-dimensional root-finding algorithm.

% 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
\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}).
Unfolded domains are modeled as WLCs (\cref{sec:tension:wlc}).

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.  As discussed in \cref{sec:tension:folded}, the
contribution of the folded domains to the elastic behavior of the
polymer-cantilever system is relatively insignificant.

% 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 assumed to lie along a line between the surface tether
point and the position of the tip (\cref{eq:sawsim:x-total}).  The
effects of this assumption are also minimized due to greater length of
the unfolded domain.  Finally, the interactions between different
parts of the polymer and between the chain and the surface (except at
the tethering points) are not considered.  This is 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.

% introduce constant velocity and walk through explicit example pull
Consider an experiment of pulling a polymer with $N$ identical protein
molecules 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:methods-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:methods-unfolding}), 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}).

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 timescale.  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.
% python -c 'print 2.45e-10*(1e9)**3 * (1.38e-23*300)**2 * 0.4e-9 / (1e-12)**2'
% 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.

\subsection{Unfolding protein molecules by force}
\label{sec:sawsim:methods-unfolding}

% 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
\begin{equation}
  k_u = k_{u0} \exp\p({\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.  The probability for
a protein molecule to unfold under an 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$.

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

% address unfolding models
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{schlierf06} used the mechanical unfolding data of the protein
ddFLN4 to demonstrate that Kramers' diffusion model fit the measured
unfolding force data better than the Bell model for proteins with
broad free energy barriers.  For proteins with relatively narrow
folded and transition states, the Bell model provides a good
approximation.

\subsection{Choosing the simulation time steps}
\label{sec:sawsim:methods-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 (see \cref{sec:sawsim:methods-tension} for equilibration 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.