Skip to content

Commit

Permalink
preparing V1 release
Browse files Browse the repository at this point in the history
  • Loading branch information
jmichel80 committed Dec 14, 2023
1 parent 455bfc6 commit 0a151af
Show file tree
Hide file tree
Showing 7 changed files with 60 additions and 40 deletions.
6 changes: 4 additions & 2 deletions LIVECOMS/02_funnel_metad/livecoms.tex
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,9 @@ \subsubsection{Introduction}\label{Introduction}}
and unbound phases, as well as increase the rate of convergence of the binding free energy, the exploration of the ligand in 3D space is limited by funnel-shaped
restraints.

\emph{fun-metaD} was originally implemented by Limogelli \emph{et al}~\cite{Limongelli2013}. This tutorial is based on the implementation described by Rhys \emph{et al}~\cite{Evans2020}, and Saleh \emph{et al}~\cite{Saleh2017}. The main difference between the original and later implementations is the functional form of the funnel restraints: the original \emph{fun-metaD} relied on a cone and a cylinder joined to make a funnel using a step function, while the new implementation uses a single
Readers that lack familiarity with metadynamics may wish to study first a simple \href{https://biosimspace.openbiosim.org/tutorials/metadynamics.html}{BioSimSpace metadynamics tutorials} for a alanine dipeptide molecule in the gas phase.

\emph{fun-metaD} was originally implemented by Limogelli et al.~\cite{Limongelli2013}. This tutorial is based on the implementation described by Rhys et al~\cite{Evans2020}, and Saleh et al~\cite{Saleh2017}. The main difference between the original and later implementations is the functional form of the funnel restraints: the original \emph{fun-metaD} relied on a cone and a cylinder joined to make a funnel using a step function, while the new implementation uses a single
sigmoid function. The Limogelli implementation also requires the protein to be realigned with a reference structure to keep the funnel strictly in place over the binding site, which negatively affects performance. The implementation described here allows the funnel to move with the protein.

A challenge with using \emph{fun-metaD} for ABFE calculations is that it can be difficult to select optimised funnel parameters, and tedious to write multiple input files for PLUMED. To facilitate use of \emph{fun-metaD} by non-experts BioSimSpace implements automated setup algorithms for selecting funnel parameters, and producing input files.
Expand Down Expand Up @@ -55,7 +57,7 @@ \subsubsection{Theory}\label{the-theory}}
In order to have a good separation between the bound and unbound phases for the ABFE estimations, the funnel needs to point \emph{out}, with the narrow end in the solvent,
and away from any protein residues. The BSS automated funnel assignment code does this generally well. Typically the vectors picked for defining the $P_{0}$ and $P_{1}$ points offer an unobstructed exit path for the ligand. It is nonetheless still a good idea to visually check the proposed funnel, especially for new protein systems. This can be done through use of the funnel visualisation functionality within BSS (see below).

The size of the funnel radius has to be determined on a case-by-case basis. The metadynamics code tracks the exploration of the center of mass of the ligand, so the volume that the small molecule will explore is much larger than implied by the visualization of the funnel. There is usually only one binding site and the funnel should enclose only it, excluding other protein features, by setting a small value of $h$. This helps accelerate convergence by preventing the ligand from exploring irrelevant regions in the free energy surface (FES).
The size of the funnel radius has to be determined on a case-by-case basis. The metadynamics code tracks the center of mass the ligand, therefore ligand atoms are able to sample positions outside the visualised funnel, provided the center of mass remains inside the funnel. the volume that the small molecule will explore is much larger than implied by the visualization of the funnel. There is usually only one binding site and the funnel should enclose only it, excluding other protein features, by setting a small value of $h$. This helps accelerate convergence by preventing the ligand from exploring irrelevant regions in the free energy surface (FES).

\hypertarget{setupfunmetad}{%
\subsubsection{Setting up funnel metadynamics simulations}\label{setupfunmetad}}
Expand Down
11 changes: 6 additions & 5 deletions LIVECOMS/03_steered_md/livecoms.tex
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
\subsubsection{Introduction}
Many relevant biological processes, such as transmembrane permeation or transitions between active and inactive protein conformations, occur on a timescale of $\mu$ s-s~\cite{Zwier2010,Choy2017,Wells2007}. However, even with GPU acceleration, the timescales accessible via MD simulations are only a few hundred ns/day~\cite{HecBioSim_benchmark}. One of the methods to get around this limitation is steered molecular dynamics (sMD). sMD involves applying a harmonic restraint to bias the system towards a conformation defined through one or more collective variables (CVs):
Many relevant biological processes, such as transmembrane permeation or transitions between active and inactive protein conformations, occur on a timescale of microseconds to seconds ~\cite{Zwier2010,Choy2017,Wells2007}. However, even with GPU acceleration, the timescales accessible via MD simulations are only a few hundred ns/day ~\cite{HecBioSim_benchmark}. One of the methods to get around this limitation is steered molecular dynamics (sMD). sMD involves applying a harmonic restraint to bias the system towards a conformation defined through one or more collective variables (CVs):

\begin{equation}
V(\vec{s},t) = \frac{1}{2} \kappa(t) ( \vec{s} - \vec{s}_0(t) )^2,
\label{eq:sMD}
\end{equation}

where $\kappa$ is the force constant, $\vec{s}_0$ is the expected CV value at a specific timestep, and $\vec{s}$ is the actual CV value at that timestep\cite{Isralewitz2001,Tribello2014}.
where $\kappa$ is the force constant, $\vec{s}_0$ is the expected CV value at a specific timestep, and $\vec{s}$ is the actual CV value at that timestep \cite{Isralewitz2001,Tribello2014}.

This section of the tutorial summarizes how to use BioSimSpace to set up and run sMD simulations. BSS prepares input files for PLUMED, which is the software that works together with MD engines such as AMBER and GROMACS to add the restraint in eq~\ref{eq:sMD}.

We use protein tyrosine phosphatase 1B (PTP1B) as the system of choice for this tutorial. It is a negative regulator of insulin signalling~\cite{sMD_ptp1b-diabetes} and is an attractive target for type II diabetes~\cite{sMD_Wiesman}. The function of PTP1B depends on the conformation of its WPD loop, which can be closed (active) or open (inactive) (Figure~\ref{fig:ptp1b}). The WPD loop of PTP1B opens and closes on a $\mu$s timescale~\cite{Choy2017}, and therefore this transition is not observed on conventional computational timescales.
We use protein tyrosine phosphatase 1B (PTP1B) as the system of choice for this tutorial. It is a negative regulator of insulin signalling ~\cite{sMD_ptp1b-diabetes}, and is an attractive target for type II diabetes ~\cite{sMD_Wiesman}. The function of PTP1B depends on the conformation of its WPD loop, which can be closed (active) or open (inactive) (Figure~\ref{fig:ptp1b}). The WPD loop of PTP1B opens and closes on a $\mu$s timescale ~\cite{Choy2017}, and therefore this transition is not observed on conventional computational timescales.

\begin{figure}[htp]
\includegraphics[width=\linewidth]{LIVECOMS/03_steered_md/open-close.png}
Expand Down Expand Up @@ -40,8 +40,9 @@ \subsubsection{sMD trajectory analysis}

The notebook also illustrates a "failed" steered MD trajectory, where the steering duration and force were insufficient to reach the target CV value.

\subsubsection{Markov State Models}
While the information provided here focuses on running sMD simulations with BSS, there are multiple potential applications, such as studying membrane permeability~\cite{Wells2007} or ligand residence time\cite{Potterton2019}. Another use is for the additional exploration of conformational space for predicting allosteric modulation using Markov State Models (MSMs). MSMs give the probability of protein conformations and therefore can be used to model how a ligand affects the conformation ensemble of a target (e.g. whether it decreases the active state probability and therefore is an allosteric inhibitor). There is a lot to consider when building MSMs, and the method is not covered in this tutorial. Here the Python library \href{http://emma-project.org/latest/}{PyEMMA} was used, which has extensive examples and documentation\cite{Wehmeyer_2019}. The integration of sMD in this allosteric modulation prediction workflow is illustrated in Figure \ref{fig:ensemble-protocol}. \href{https://github.com/OpenBioSim/BioSimSpaceTutorials/blob/main/03_steered_md/02_trajectory_analysis.ipynb}{02-trajectory-analysis} shows how the sMD trajectory can be sampled to extract a range of protein conformations. Hardie \emph{et al} report a detailed study of allosteric modulators of PTP1B using this sMD/MSM methodology~\cite{Hardie2023}.
\subsubsection{Example application - combining steered MD with Markov State Modeling}
While the information provided here focuses on running sMD simulations with BSS, there are multiple potential applications, such as studying membrane permeability~\cite{Wells2007} or ligand residence time\cite{Potterton2019}. Here we briefly highlight one application of sMD simulations enabled by BioSimSpace in the AMMo software project. AMMo ("Allostery in Markov Models") was developed to evaluate the allosteric effects of protein mutations or ligand binding by combining sMD with Markov State Models (MSMs). MSMs are used to give the probability of protein conformations and therefore can be used to model how a ligand affects the conformation ensemble of a target (e.g. whether the presence of a ligand decreases the active state probability and therefore is an allosteric inhibitor).
There is a lot to consider when building MSMs, and the method is not covered in this tutorial. AMMo uses the Python library \href{http://emma-project.org/latest/}{PyEMMA} was to implement MSMs. Extensive examples and documentation for PyEMMA are available elsewhere ~\cite{Wehmeyer_2019}. The integration of sMD with MSM in this allosteric modulation prediction workflow is illustrated in Figure \ref{fig:ensemble-protocol}. The notebook \href{https://github.com/OpenBioSim/BioSimSpaceTutorials/blob/main/03_steered_md/02_trajectory_analysis.ipynb}{02-trajectory-analysis} shows how a sMD trajectory can be sampled to extract a range of protein conformations suitable for inputs to an MSM workflow. Hardie et al. report a detailed study of allosteric modulators of PTP1B using this sMD/MSM methodology ~\cite{Hardie2023} and notebooks for the PTP1B case study are available on the \href{https://github.com/michellab/AMMo/tree/main/examples/ptp1b}{GitHub} page for AMMo.

\begin{figure}[htp]
\includegraphics[width=\linewidth]{LIVECOMS/03_steered_md/ensemble-md-protocol.png}
Expand Down
Loading

0 comments on commit 0a151af

Please sign in to comment.