Alchemical Transfer Method
The Alchemical Transfer Method (ATM) is a really interesting technique originally developed by the Gallichio lab at CUNY.
In this post, we will show you how use BioSimSpace to set up and run a Relative Binding Free Energy (RBFE) calculation using ATM on a pair of ligands bound to Tyrosine kinase 2 (TYK2).
In ATM, a relative binding free energy calculation is implemented by swapping a pair of ligands that are simulated in a single box. For more details about the theory we refer you to the documentation from the Gallichio lab.
System setup and visualisation
Let’s first load our inputs
The protein used in this case is tyrosine kynase 2, or TYK2, a common benchmark system for binding free energy calculations. You can verify this by visualising the loaded protein input.
Unlike in other RBFE methods, a system used in ATM calculation includes both ligands simultaneously, with one ligand bound to the protein and the other free in the solvent.
We will use BioSimSpace to create an ATM compatible system from the inputs loaded above, choosing which ligand is to be bound and which is to be free.
One decision that needs to be made prior to creating an ATM compatible system is the choice of “rigid core atoms” within each ligand. These atoms are used to define a series of aligning forces that will be used to maintain the relative orientation of the two ligands during minimisation, equilibration and production simulations.
Maintaining the orientation of the ligands is important to the stability of ATM simulations as, in the second half of production (i.e. windows with a lambda value greater than 0.5) the free and bound ligands will be swapped. If the free ligand is out of alignment during this swap it will cause direct overlaps with protein atoms, causing the simulations to immediately explode.
BioSimSpace provides a helper function that visualises the choice of these atoms.
Making a good choice of rigid core atoms can be complex. A good starting point is to find a common core between your two ligands and choose three well-separated atoms within this common core.
Now that a sensible choice of rigid core atoms has been made we can now create the system using the BioSimSpace.FreeEnergy.ATMSetup class.
System preparation returns a pair of objects, the first is a system containing the protein and two ligands, the second is a dictionary containing information on the setup of the system that will be used to ensure that consistent settings are used throughout all simulations.
By default, BioSimSpace attempts to find a best-fit vector along which to translate the free ligand relative to the bound, and then translates the free ligand along this vector by a distance of 20 Angstroms. The translation vector can also be set manually by passing a displacement vector to atm_setup.prepare.
Now lets visualise the system to make sure that the free ligand is far enough away from the protein
After that we can solvate our system using BioSimSpace.Solvent
Minimisation and Equilibration
We are now ready to minimise and equilibrate the system.
ATM simulations require a significant number of restraints, as well as the gradual introduction of the ATMForce itself. As such standard minimisation and equilibration protocols are not appropriate, and a series of custom protocols is needed.
In this section we will cover a full minimisation and equilibration protocol for the ATM system we created above.
The simulations we will run in this section are significantly cut down, real production simulations should be minimised and equilibrated much more thoroughly than we will here, with at least 10000 minimisation steps and 100ps of runtime for each equilibration.
ATM simulations are best run with a series of restraints:
- core alignment; these are the the rigid core restraints that are applied to the atoms we found earlier.
- positional restraints; these are flat-bottom restraints used to keep specific atoms in place, these are generally applied to the alpha carbons of the protein.
- centre of mass distance restraints; these are used to maintain the relative positions of the protein and ligands and are applied to the centre off mass of the atoms listed in the
datadictionary (these were found automatically when we prepared the system, but they can be set manually if the ones found are not appropriate).
Now let’s find the alpha carbons that we will be restraining using BioSimSpace search functionality and setup our minimisation.
We will now run our first equilibration, including all of the forces we used in the minimisation.
The system has now been minimised and equilibrated as a vanilla system, that is without the ATMForce present.
We now need to introduce the ATMForce to the system; this introduction needs to be gradual, otherwise our simulations will crash.
This gradual introduction can be done by annealing the system to a lambda value of 0.5. Starting from lambda=0, the annealing protocol simulates a series of windows in which the value of lambda is gradually increased. In this case we will use 10 annealing cycles with a total runtime of 1ps, meaning that the value of lambda is increased by a value of 0.05 every 0.1ps
Finally, we need to perform a post-annealing equilibration. This equilibration step is performed at a lambda value of 0.5, and is designed to introduce the free ligand to the protein-ligand complex without actually performing a swap. This step dramatically increases the stability of production simulations in which `ligand_free’ is in the binding site (i.e. those with a lambda value greater than 0.5) by decreasing the probability of atom overlaps occurring when the ligand is placed in the binding site.
When running full production simulations this is the step that ytou should look to first if you experience instability. It is usually advisable to run this post-annealing equilibration for at least 500ps to give the system enough time to settle.
Production and Analysis
Now that the system has been properly minimised and equilibrated it is ready for production.
Before running production simulations we need to decide how many lambda windows to run. Keep in mind that when performing an ATM calculation we are running the equivalent of all 4 legs (bound and free, forward and reverse) of a standard RBFE calculation simultaneously, so we should expect to have to run a larger number of windows.
A good middle ground for most systems is 22 windows, and this is the default chosen by BioSimSpace, but some systems might need more sampling.
For the TYK2 simulations we are running in this tutorial we will need to change some of the soft-core settings, as the BioSimSpace defaults are not appropriate.
Once production simulations are complete output_directory will be populated by a series of folders, one for each lambda window. Within each of these folders should be a file called openmm.csv which contains the information we need to calculate the ΔΔG value for our system. In this case we will be analysing the pre-prepared outputs that we downloaded earlier.
In-built BioSimSpace analysis will give us a relative binding free energy value in kcal/mol, along with the error, also in kcal/mol.
If you want to try this out for yourself, check out our detailed alchemical transfer tutorial by logging in to our JupyterHub server at try.openbiosim.org (GitHub account required) and following the instructions in the file TUTORIALS.txt.
This tutorial was written by Matthew Burman
