You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
The :class:`AbsoluteBindingProtocol <.AbsoluteBindingProtocol>` calculates the absolute binding free energy,
8
+
which is the free energy difference between a ligand in solution and the ligand bound to a protein.
9
+
10
+
The absolute binding free energy is calculated through a thermodynamic cycle.
11
+
In this cycle, the interactions of the molecule are decoupled, meaning turned off,
12
+
using a partial annihilation scheme (see below) both in the solvent and in the complex phases.
13
+
14
+
Restraints are required to keep the weakly
15
+
coupled and fully decoupled ligand in the binding site region and thereby reduce the phase
16
+
space that needs to be sampled. In the :class:`AbsoluteBindingProtocol <.AbsoluteBindingProtocol>`
17
+
we apply orientational, or Boresch-style, restraints, as described below.
18
+
19
+
The absolute binding free energy is then obtained via summation of free energy differences along the thermodynamic cycle.
20
+
21
+
.. figure:: img/abfe-cycle.png
22
+
:scale:50%
23
+
24
+
Thermodynamic cycle for the absolute binding free energy protocol.
25
+
26
+
Scientific Details
27
+
------------------
28
+
29
+
Orientational restraints
30
+
~~~~~~~~~~~~~~~~~~~~~~~~
31
+
32
+
Orientational, or Boresch-style, restraints are automaticallly (unless manually specified) applied between three
33
+
protein and three ligand atoms using one bond, two angle, and three dihedral restraints.
34
+
Reference atoms are picked based on different criteria, such as the root mean squared
35
+
fluctuation of the atoms in a short MD simulation, the secondary structure of the protein,
36
+
and the distance between atoms, based on heuristics from Baumann et al. [1]_ and Alibay et al. [2]_.
37
+
Two strategies for selecting protein atoms are available, either picking atoms that are bonded to each other or that can span multiple residues.
38
+
This can be specified using the ``restraint_settings.anchor_finding_strategy`` settings.
39
+
40
+
Partial annihilation scheme
41
+
~~~~~~~~~~~~~~~~~~~~~~~~~~~
42
+
43
+
In the :class:`.AbsoluteBindingProtocol` the coulombic interactions of the molecule are fully turned off (annihilated).
44
+
The Lennard-Jones interactions are instead decoupled, meaning the intermolecular interactions are turned off, keeping the intramolecular Lennard-Jones interactions.
45
+
46
+
The lambda schedule
47
+
~~~~~~~~~~~~~~~~~~~
48
+
49
+
Molecular interactions are turned off during an alchemical path using a discrete set of lambda windows.
50
+
For the transformation in the binding site, the following steps are carried out, starting with the ligand fully interacting in the binding site.
51
+
52
+
1. Restrain the ligand using orientational restraints.
53
+
2. Turn off the electrostatic interactions of the ligand.
54
+
3. Decouple Lennard-Jones interactions of the ligand.
55
+
4. Release the restraints of the now dummy ligand analytically.
56
+
57
+
The lambda schedule in the solvent phase is similar to the one in the complex, except that no restraints are applied.
58
+
A soft-core potential is applied to the Lennard-Jones potential to avoid instablilites in intermediate lambda windows.
59
+
The soft-core potential function from Beutler et al. [3]_ is used by default.
60
+
The lambda schedule is defined in the ``lambda_settings`` objects ``lambda_elec``, ``lambda_vdw``, and ``lambda_restraints``.
61
+
62
+
Simulation overview
63
+
~~~~~~~~~~~~~~~~~~~
64
+
65
+
The :class:`.ProtocolDAG` of the :class:`.AbsoluteBindingProtocol` contains :class:`.ProtocolUnit`\ s from both the complex and solvent transformations.
66
+
This means that both legs of the thermodynamic cycle are constructed and run concurrently in the same :class:`.ProtocolDAG`.
67
+
This is different from the :class:`.RelativeHybridTopologyProtocol` where the :class:`.ProtocolDAG` only runs a single leg of a thermodynamic cycle.
68
+
If multiple ``protocol_repeats`` are run (default: ``protocol_repeats=3``), the :class:`.ProtocolDAG` contains multiple :class:`.ProtocolUnit`\ s of both complex and solvent transformations.
69
+
70
+
Simulation steps
71
+
""""""""""""""""
72
+
73
+
Each :class:`.ProtocolUnit` (whether complex or solvent) carries out the following steps:
74
+
75
+
1. Parameterize the system using `OpenMMForceFields <https://github.com/openmm/openmmforcefields>`_ and `Open Force Field <https://github.com/openforcefield/openff-forcefields>`_.
76
+
2. Equilibrate the fully interacting system using a short MD simulation using the same approach as the :class:`.PlainMDProtocol` (including rounds of NVT and NPT equilibration).
77
+
3. Create an alchemical system.
78
+
4. Add orientational restraints to the complex system.
79
+
5. Minimize the alchemical system.
80
+
6. Equilibrate and production simulate the alchemical system using the chosen multistate sampling method (under NPT conditions).
81
+
7. Analyze results for the transformation.
82
+
83
+
.. note:: Three different types of multistate sampling (i.e. replica swapping between lambda states) methods can be chosen; HREX, SAMS, and independent (no lambda swaps attempted).
84
+
By default the HREX approach is selected, this can be altered using ``solvent_simulation_settings.sampler_method`` or ``complex_simulation_settings.sampler_method`` (default: ``repex``).
85
+
86
+
Simulation details
87
+
""""""""""""""""""
88
+
89
+
Here are some details of how the simulation is carried out which are not detailed in the :class:`.AbsoluteBindingSettings`:
90
+
91
+
* The protocol applies a `LangevinMiddleIntegrator <https://openmmtools.readthedocs.io/en/latest/api/generated/openmmtools.mcmc.LangevinDynamicsMove.html>`_ which uses Langevin dynamics, with the LFMiddle discretization [4]_.
92
+
* A MonteCarloBarostat is used in the NPT ensemble to maintain constant pressure.
93
+
94
+
Getting the free energy estimate
95
+
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
96
+
97
+
The free energy differences are obtained from simulation data using the `MBAR estimator <https://www.alchemistry.org/wiki/Multistate_Bennett_Acceptance_Ratio>`_ (multistate Bennett acceptance ratio estimator) as implemented in the `PyMBAR package <https://pymbar.readthedocs.io/en/master/mbar.html>`_.
98
+
Both the MBAR estimates of the two legs of the thermodynamic cycle, and the overall absolute binding free energy (of the entire cycle) are obtained,
99
+
which is different compared to the results in the :class:`.RelativeHybridTopologyProtocol` where results from two legs of the thermodynamic cycle are obtained separately.
100
+
101
+
In addition to the estimates of the free energy changes and their uncertainty, the protocol also returns some metrics to help assess convergence of the results, these are detailed in the :ref:`multistate analysis section <multistate_analysis>`.
102
+
103
+
See Also
104
+
--------
105
+
106
+
**Setting up AFE calculations**
107
+
108
+
* :ref:`Defining the Protocol <defining-protocols>`
.. [1] Broadening the Scope of Binding Free Energy Calculations Using a Separated Topologies Approach, H. Baumann, E. Dybeck, C. McClendon, F. Pickard IV, V. Gapsys, L. Pérez-Benito, D. Hahn, G. Tresadern, A. Mathiowetz, D. Mobley, J. Chem. Theory Comput., 2023, 19, 15, 5058–5076
133
+
.. [2] Evaluating the use of absolute binding free energy in the fragment optimisation process, I. Alibay, A. Magarkar, D. Seeliger, P. Biggin, Commun Chem 5, 105 (2022)
134
+
.. [3] Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations, T.C. Beutler, A.E. Mark, R.C. van Schaik, P.R. Greber, and W.F. van Gunsteren, Chem. Phys. Lett., 222 529–539 (1994)
135
+
.. [4] Unified Efficient Thermostat Scheme for the Canonical Ensemble with Holonomic or Isokinetic Constraints via Molecular Dynamics, Zhijun Zhang, Xinzijian Liu, Kangyu Yan, Mark E. Tuckerman, and Jian Liu, J. Phys. Chem. A 2019, 123, 28, 6056-6079
Below are the settings which can be tweaked in the protocol. The default settings (accessed using :meth:`AbsoluteBindingProtocol.default_settings`) will automatically populate settings which we have found to be useful for running binding free energy calculations. There will however be some cases (such as when calculating difficult to converge systems) where you will need to tweak some of the following settings.
0 commit comments