Remove silly tex/ directory.

[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.
\Cref{fig:sawsim:sim-hist} shows the distribution of the unfolding forces,
\ie, the highest force in each peak (except the last peak in a force
curve), from a total of $400$ force curves ($3200$ force values).  The
unfolding forces have an average of $281\U{pN}$ with a standard
deviation of $25\U{pN}$.
% DONE: discuss noise? no.

\begin{figure}
\vspace{-1in}
\begin{center}
\subfloat[][]{\asyinclude{figures/sim-sawtooth/sim-sawtooth}%
  \label{fig:sawsim:sim-sawtooth}%
}\\
\subfloat[][]{\asyinclude{figures/sim-hist/sim-hist}%
  \label{fig:sawsim:sim-hist}%
}
\caption{(a) 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.  (b) The
  distribution of the unfolding forces from $400$ simulated force
  curves ($3200$ data points) such as that shown in (a). 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 ($3200$, for this
  histogram).\label{fig:sawsim:sim-all}}
\end{center}
\end{figure}

\subsection{Dependence of the unfolding force on the unfolding order
and polymer length}
\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, leading to a decrease in the
unfolding force.  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 depending 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
\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 $\sigma^2 = (\pi k_BT/\Delta x_u)^2/6$,
and and average
\begin{equation}
  \avg{F(i)} = \frac{k_BT}{\Delta x_u}
               \p[{-\ln\p({\frac{N_fk_{u0}k_BT}{\kappa v\Delta x_u}})
                   -\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$.

From \cref{fig:sawsim:order-dep}, it can be seen 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.  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.  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.  \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}

The zero-force unfolding rate $k_{u0}$ and the distance $\Delta x_u$
from the native state to the transition state are the two kinetic
parameters obtainable for mechanical unfolding experiments by matching
the simulated data with measured results.  \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{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 discussed in \cref{sec:sawsim:results-scaffold}.  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.

In most published literature, determination of the values of $k_{u0}$
and $\Delta x_u$ was mostly done by carrying out simulations using a
handful of possible unfolding parameters and selected the best fit by
eye%
%\citep{us,CV1999}
.  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.
% DONE: mention $D_\text{JS}$ features to explains selection over $\chi^2$? no.
\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{(a) 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.  (b) The dependence of standard deviation of the unfolding
    force distribution on the pulling speed for the simulation data
    shown in (a), 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 (a), 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}