From ce839c063b0fd57a0f3dc50a32e960333d42d50d Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 21 May 2013 16:41:39 -0400 Subject: [PATCH] sawsim/discussion.tex: Rework Bell testing to match scaffold phrasing This fills in the derivation of the Gumbel pdf from the exponential hazard function. In the scaffold section, we rephrased the hazard function and then cited NIST:gumbel to make the transition to the pdf. Maybe this derivation would be better of in an appendix? Or dropped entirely? I'll leave it for now. --- src/apparatus/cantilever-calib.tex | 1 + src/sawsim/discussion.tex | 333 +++++++++++++---------------- src/sawsim/methods.tex | 5 +- 3 files changed, 159 insertions(+), 180 deletions(-) diff --git a/src/apparatus/cantilever-calib.tex b/src/apparatus/cantilever-calib.tex index 961630b..4fb942e 100644 --- a/src/apparatus/cantilever-calib.tex +++ b/src/apparatus/cantilever-calib.tex @@ -33,6 +33,7 @@ over a long time interval. \begin{equation} \avg{A} \equiv \iLimT{A} \;. \end{equation}} +\nomenclature{$\equiv$}{Defined as (\ie\ equivalent to)} To calculate the spring constant $\kappa$ using \cref{eq:equipart}, we need to measure the buffer temperature $T$ and the thermal vibration diff --git a/src/sawsim/discussion.tex b/src/sawsim/discussion.tex index 5a5c7a9..ae438df 100644 --- a/src/sawsim/discussion.tex +++ b/src/sawsim/discussion.tex @@ -162,7 +162,18 @@ The rate of unfolding events with respect to force is \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, and $k_u$ is the unfolding rate constant +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, @@ -183,9 +194,10 @@ the scale and location parameters, respectively\citep{hummer03} -\exp{\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 +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} \;, @@ -209,7 +221,12 @@ unfolding force histograms produced under such a variable loading rate. % \nomenclature{$r_{uF}$}{Unfolding loading rate (newtons per second)} -\nomenclature{$\gamma_e$}{Euler-Macheroni constant, $\gamma_e=0.577\ldots$} +\nomenclature{$\alpha$}{The mode unfolding force, + $\alpha\equiv-\rho\ln(N_f k_{u0}\rho/\kappa v)$ + (\cref{eq:sawsim:gumbel}).} +\nomenclature{$\rho$}{The characteristic unfolding force, + $\rho\equiv k_BT/\Delta x_u$ (\cref{eq:sawsim:gumbel}).} +\nomenclature{$\gamma_e$}{Euler--Macheroni constant, $\gamma_e=0.577\ldots$} \nomenclature{$\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 @@ -594,208 +611,113 @@ 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. -\section{Review of current research} - -There are two main approaches to modeling protein domain unfolding -under tension: Bell's and Kramers'\citep{schlierf06,hummer03,dudko06}. -Bell introduced his model in the context of cell -adhesion\citep{bell78}, but it has been widely used to model -mechanical unfolding in -proteins\citep{rief97a,carrion-vazquez99b,schlierf06} due to 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. - -\subsection{Evolution of unfolding modeling} - -Evans introduced the saddle-point Kramers' approximation in a protein unfolding context in 1997 (\citet{evans97} Eqn.~3). -However, early work on mechanical unfolding focused on the simpler Bell model\citep{rief97a}.%TODO -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}.%TODO -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}.%TODO - -There have been some tangential attempts towards even fancier models. -\citet{dudko03} attempted to reduce the restrictions of the single-unfolding-path model. -\citet{hyeon03} attempted to measure the local roughness using temperature dependent unfolding. - -\subsection{History of simulations} - -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{MD}{Molecular dynamics simulation. Simulate the - physical motion of atoms and molecules by numerically solving - Newton's equations.} - -\subsection{History of experimental AFM unfolding experiments} - -\begin{itemize} - \item \citet{rief97a}: -\end{itemize} - -\subsection{History of experimental laser tweezer unfolding experiments} - -\begin{itemize} - \item \citet{izrailev97}: -\end{itemize} - -\section{Single-domain proteins under constant loading} - -TODO: consolidate with \cref{sec:sawsim:results:scaffold}. - -Let $x$ be the end to end distance of the protein, $t$ be the time -since loading began, $F$ be tension applied to the protein, $N_f$ be -the surviving population of folded proteins. Make the definitions -\begin{align} - v &\equiv \deriv{t}{x} && \text{the pulling velocity} \\ - \kappa &\equiv \deriv{x}{F} && \text{the loading spring constant} \\ - N_{f0} &\equiv N_f(t=0) && \text{the initial number of folded proteins} \\ - N_u &\equiv N_{f0} - N_f && \text{the number of unfolded proteins} \\ - k_u &\equiv -\frac{1}{N_f} \deriv{t}{N_f} && \text{the unfolding rate} -\end{align} -\nomenclature{$\equiv$}{Defined as (\ie\ equivalent to)} -The proteins are under constant loading because -\begin{equation} - \deriv{t}{F} = \deriv{x}{F}\deriv{t}{x} = \kappa v\;, -\end{equation} -a constant, since both $\kappa$ and $v$ are constant (\citet{evans97} -in the text on the first page, \citet{dudko06} in the text just before -\fref{equation}{4}). - 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$ (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$). +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}{vk} N_f\kappa \label{eq:unfold:hist} + = \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 +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} \int k_0 \dd F &= \ln(N_f) + c \\ - N_f &= C\exp{\frac{-1}{\kappa v}\integral{}{}{F}{k_u}} \;, + \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 $c \equiv \ln(C)$ is a constant of integration scaling $N_f$. +where $N_f(0) = N$ because all the domains are initially folded. +% +\nomenclature{$W$}{Bin width of an unfolding force histogram + (\cref{eq:unfold:hist}).} -\subsection{Constant unfolding rate} +\subsubsection{Constant unfolding rate} -In the extremely weak tension regime, the proteins' unfolding rate is independent of tension, we have +In the extremely weak tension regime, the proteins' unfolding rate is +independent of tension, so we can simplify \cref{eq:N_f} and plug into +\cref{eq:unfold:hist}. \begin{align} - P &= C\exp{\frac{-1}{kv}\integral{}{}{F}{\kappa}} - = C\exp{\frac{-1}{kv}\kappa F} - = C\exp{\frac{-\kappa F}{kv}} \\ - P(0) &\equiv P_0 = C\exp{0} = C \\ - h(F) &= \frac{W}{vk} P \kappa - = \frac{W\kappa P_0}{vk} \exp{\frac{-\kappa F}{kv}} + 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} -So, a constant unfolding-rate/hazard-function gives exponential decay. -Not the most earth shattering result, but it's a comforting first step, and it does show explicitly the dependence in terms of the various unfolding-specific parameters. +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. -\subsection{Bell model} +\subsubsection{Bell model} Stepping up the intensity a bit, we come to Bell's model for unfolding -(\citet{hummer03} Eqn.~1 and the first paragraph of \citet{dudko06} and \citet{dudko07}). +(\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} - \kappa = \kappa_0 \cdot \exp{\frac{F \dd x}{k_B T}} - = \kappa_0 \cdot \exp{a F} \;, + k_u = k_{u0} \exp{\frac{F}{\rho}} \;. \end{equation} -where we've defined $a \equiv \dd x/k_B T$ to bundle some constants together. -The unfolding histogram is then given by +The unfolding histogram is then given via \cref{eq:N_f,eq:unfold:hist}. \begin{align} - P &= C\exp{\frac{-1}{kv}\integral{}{}{F}{\kappa}} - = C\exp{\frac{-1}{kv} \frac{\kappa_0}{a} \exp{a F}} - = C\exp{\frac{-\kappa_0}{akv}\exp{a F}} \\ - P(0) &\equiv P_0 = C\exp{\frac{-\kappa_0}{akv}} \\ - C &= P_0 \exp{\frac{\kappa_0}{akv}} \\ - P &= P_0 \exp{\frac{\kappa_0}{akv}\p({1-\exp{a F}})} \\ - h(F) &= \frac{W}{vk} P \kappa - = \frac{W}{vk} P_0 - \exp{\frac{\kappa_0}{akv}\p({1-\exp{a F}})} \kappa_0 \exp{a F} - = \frac{W\kappa_0 P_0}{vk} - \exp{a F + \frac{\kappa_0}{akv}\p({1-\exp{a F}})} \;. + 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} -The $F$ dependent behavior reduces to -\begin{equation} - h(F) \propto \exp{a F - b\exp{a F}} \;, -\end{equation} -where $b \equiv \kappa_0/akv \equiv \kappa_0 k_B T / k v \dd x$ is -another constant rephrasing. - -This looks similar to the Gompertz / Gumbel / Fisher-Tippett -distribution, where -\begin{align} - p(x) &\propto z\exp{-z} \\ - z &\equiv \exp{-\frac{x-\mu}{\beta}} \;, -\end{align} -but we have -\begin{equation} - p(x) \propto z\exp{-bz} \;. -\end{equation} -Strangely, the Gumbel distribution is supposed to derive from an -exponentially increasing hazard function, which is where we started -for our derivation. I haven't been able to find a good explaination -of this discrepancy yet, but I have found a source that echos my -result (\citet{wu04} Eqn.~1). TODO: compare \citet{wu04} with -my successful derivation in \cref{sec:sawsim:results-scaffold}. - -Oh wait, we can do this: -\begin{equation} - p(x) \propto z\exp{-bz} = \frac{1}{b} z'\exp{-z'}\propto z'\exp{-z'} \;, -\end{equation} -with $z'\equiv bz$. I feel silly... From -\href{Wolfram}{http://mathworld.wolfram.com/GumbelDistribution.html}, -the mean of the Gumbel probability density -\begin{equation} - P(x) = \frac{1}{\beta} \exp{\frac{x-\alpha}{\beta} - -\exp{\frac{x-\alpha}{\beta}}} - \label{eq:sawsim:gumbel-x} -\end{equation} -is given by $\mu=\alpha-\gamma\beta$, and the variance is -$\sigma^2=\frac{1}{6}\pi^2\beta^2$, where $\gamma=0.57721566\ldots$ is -the Euler-Mascheroni constant. Selecting $\beta=1/a=k_BT/\dd x$, -$\alpha=-\beta\ln(\kappa\beta/kv)$, and $F=x$ we have -\nomenclature{$\mu$}{The mean of a distribution (e.g. the Gumbel - distribution, \cref{eq:sawsim:gumbel-x}).} -\begin{align} - P(F) - &= \frac{1}{\beta} \exp{\frac{F+\beta\ln(\kappa\beta/kv)}{\beta} - -\exp{\frac{F+\beta\ln(\kappa\beta/kv)}{\beta}}} \\ - &= \frac{1}{\beta} \exp{F/\beta}\exp{\ln(\kappa\beta/kv)} - \exp{-\exp{F/\beta}\exp{\ln(\kappa\beta/kv)}} \\ - &= \frac{1}{\beta} \frac{\kappa\beta}{kv} \exp{F/\beta} - \exp{-\kappa\beta/kv\exp{F/\beta}} \\ - &= \frac{\kappa}{kv} \exp{F/\beta}\exp{-\kappa\beta/kv\exp{F/\beta}} \\ - &= \frac{\kappa}{kv} \exp{F/\beta - \kappa\beta/kv\exp{F/\beta}} \\ - &= \frac{\kappa}{kv} \exp{aF - \kappa/akv\exp{aF}} \\ - &= \frac{\kappa}{kv} \exp{aF - b\exp{aF}} - \propto h(F) \;. -\end{align} -So our unfolding force histogram for a single Bell domain under -constant loading does indeed follow the Gumbel distribution. - -% Consolidate with src/sawsim/discussion.tex +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$. +}. +% +\nomenclature{$\alpha'$}{The mode unfolding force for a single folded + domain, $\alpha'\equiv-\rho\ln(k_{u0}\rho/\kappa v)$ + (\cref{eq:unfold:bell_pdf}).} -\subsection{Saddle-point Kramers' model} +\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). @@ -821,9 +743,10 @@ $l_{ts}$ is the characteristic length of the transition state \nomenclature{$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. +\citet{evans97} solved this unfolding rate for both inverse power law +potentials and cusp potentials. -\section{Double-integral Kramers' theory} +\subsubsection{Double-integral Kramers' theory} The double-integral form of overdamped Kramers' theory may be too complex for analytical predictions of unfolding-force histograms. @@ -836,3 +759,55 @@ portion of the simulation. Looking for analytic solutions to Kramers' $k(F)$, we find that there are not many available in a closed form. However, we do have analytic solutions for unforced $k$ for cusp-like and quartic potentials. + +\section{Review of current research} + +There are two main approaches to modeling protein domain unfolding +under tension: Bell's and Kramers'\citep{schlierf06,hummer03,dudko06}. +Bell introduced his model in the context of cell +adhesion\citep{bell78}, but it has been widely used to model +mechanical unfolding in +proteins\citep{rief97a,carrion-vazquez99b,schlierf06} due to 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. + +\subsection{Evolution of unfolding modeling} + +Evans introduced the saddle-point Kramers' approximation in a protein unfolding context in 1997 (\citet{evans97} Eqn.~3). +However, early work on mechanical unfolding focused on the simpler Bell model\citep{rief97a}.%TODO +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}.%TODO +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}.%TODO + +There have been some tangential attempts towards even fancier models. +\citet{dudko03} attempted to reduce the restrictions of the single-unfolding-path model. +\citet{hyeon03} attempted to measure the local roughness using temperature dependent unfolding. + +\subsection{History of simulations} + +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{MD}{Molecular dynamics simulation. Simulate the + physical motion of atoms and molecules by numerically solving + Newton's equations.} + +\subsection{History of experimental AFM unfolding experiments} + +\begin{itemize} + \item \citet{rief97a}: +\end{itemize} + +\subsection{History of experimental laser tweezer unfolding experiments} + +\begin{itemize} + \item \citet{izrailev97}: +\end{itemize} diff --git a/src/sawsim/methods.tex b/src/sawsim/methods.tex index 3f54d56..951c327 100644 --- a/src/sawsim/methods.tex +++ b/src/sawsim/methods.tex @@ -402,7 +402,10 @@ approximations and assumptions we make when we use these simple models % 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 +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} -- 2.26.2