This page provides a brief explanation of the workflow to perform APR calculations with pAPRika. For more detail users are recommended to go through the tutorials, which details further on how to setup and run APR simulations from start to finish.
Aligning the host-guest complex
The starting structure for the APR simulation can be configured with paprika. The APR calculation is most efficient in a rectangular box with the long axis parallel to the pulling axis (reduces the number of water molecules in the system). To make this easy to set up, the Align
module provides functions to shift and orient a structure. For example, we can translate a structure to the origin and then orient the system to the z-axis by running
Adding dummy atoms
To provide something to pull "against," we add dummy atoms that are fixed in place with strong positional restraints. These dummy atoms can be added to a the host-guest structure using the Dummy Atoms
module in paprika.
We will need the mol2
and frcmod
files for the dummy atoms, which we will need to generate the AMBER topology
Building the topology
Finally, we can use the tleap
wrapper to combine all of these components to generate the topology and coordinate files.
Solvating the structure
The TLeap
wrapper also provides an API for choosing water models when the user wants to solvate their structure. tleap
provides a number of models for both 3-point and 4-point water models. There are also sub-types for each water model, e.g. for TIP3P we can choose to use the ForceBalance optimized variant called TIP3P-FB or the host-guest binding optimized model called Bind3P. To choose a water model for solvation we use the set_water_model
method of the TLeap
wrapper. The method requires the user to specify the water model and optionally the sub-type as the model_type
attibute. The supported water models are:
- spc: None (SPCBOX), "flexible" (SPCFWBOX), "quantum" (QSPCFWBOX)
- opc: None (OPCBOX), "three-point" (OPC3BOX)
- tip3p: None (TIP3PBOX), "flexible" (TIP3PFBOX), "force-balance" (FB3BOX)
- tip4p: None (TIP4PBOX), "ewald" (TIP4PEWBOX), "force-balance" (FB4BOX)
Below is an example for solvating a system with 2000 TIP3P water molecules with ForceBalance
optimized parameters.
In APR calculations we apply restraints on the host (or protein) and the guest molecules. The restraints can be grouped into four categories: (1) static restraints, (2) varying restraints, (3) wall restraints and (4) positional restraints. The equilibrium target values and force constants can be specified as either a float or Pint quantity through the openff-units wrapper.
(1) Static Restraints
Static restraints do not change during the whole APR process and do not affect the free energy. We apply static restraints on the host (or protein) molecule to orient the host/protein degrees of freedom. The static restraints are composed of distance, angle, and torsional (DAT) restraints based on the choice of anchor atoms. For host-guest systems, we need to define three anchor atoms [H1,H2,H3]
and combined with three dummy atoms [D1,D2,D3]
, we apply a total of six static restraints on the host molecule (three for the translation and three for orientation).
To generate static restraints we use the function static_DAT_restraints
. As an example, to apply a distance restraint on D1
and H1
with a force constant of 5 kcal/mol/Å2 we call
The equilibrium target for the harmonic restraint is estimated from the ref_structure
.
(2) Varying Restraints
As the name suggests, these restraints change during the APR process. During the attach and release phases, the force constants of these restraints changes. In the pull phase, varying restraints can have their equilibrium position change, and this can be used as the restraint to pull the guest molecule out of the host molecule.
To generate varying restraints, we use the DAT_restraint
class. The code below shows a restraints r that starts from 6.0 Å to 24 Å in the pull phase and stays restrained at 24 Å during the release phase.
(3) Wall Restraints (optional)
Wall restraints are half-harmonic potentials that is useful for preventing guest molecules from leaving the binding site (for weak binding) or preventing the guest molecule from flipping during the attach phase. We still use the DAT_restraint
class to generate the restraints but will use the custom_restraint_values
method to generate the half-harmonic potential.
Below is an example for generating a "lower wall" restraint that prevents the angle of [D2,G1,G2]
from decreasing below 91 degrees.
(4) Positional Restraints
Positional restraints in APR simulations are applied to the dummy atoms. Together with static restraints, this provides a laboratory frame of reference for the host-guest complex. Different MD programs handles positional restraints differently. For example, in AMBER
you can define positional restraints in the input configuration file using the ntr
keyword (Chapter 19 in the AMBER20 manual). For other programs like GROMACS
and NAMD
that uses Plumed
, positional restraints can be applied using the method add_dummy_atom_restraints()
.
Creating the APR windows and saving restraints to file
To create the windows for the APR calculation we need to parse a varying restraint to the utility function create_window_list
. This function will return a list of strings for the APR protocol
It may also be useful to save both the windows list and the restraints to a JSON file so you do not need to redefine again. The restraints can be saved to a JSON file using the utility function save_restraints
.
Extending/adding more windows
Sometimes it may be necessary to add more windows in the APR calculation due to insufficient overlap between neighboring windows. For convenience we can add the windows at the end of the current list instead of inserting them in order. For example, let's say that we have a defined a restraint that spans from 8.4 to 9.8 Å and we want to add three windows between 8.6 and 9.0 Å.
We will just need to append the target_list of this dictionary and reinitialize the restraints
We can save the updated restraints to a new file and pass it to the analysis script. The fe_calc
class will take care of the window ordering thus there is no need to manually order the windows.
paprika provides wrappers with the Simulate
module for a number of MD engines enabling us to run the simulations in python.
Once the simulation is complete, the free energy can be obtained using the Analysis
module, which will also estimate the uncertainties using the bootstrapping method. There are three types of methods that you can do with the Analysis
module: (1) thermodynamic integration with block-data analysis ("ti-block"), (2) multistate Benett-Acceptance-Ratio with block-data analysis ("mbar-block"), and (3) multistate Benett-Acceptance-Ratio with autocorrelation analysis ("mbar-autoc").
We can also estimate the free energy cost of releasing the restraints on the guest molecule semianalytically. To do this we need to extract the restraints that is specific to the guest molecule. The extract_guest_restraints
function from the restraints
module and pass this to the analysis object.
The results are stored in the variable results
as a python dictionary and you can save this to a JSON file.
The processed simulation data can also be saved to a JSON file so that you do not need to re-read the MD trajectories if you need to do further analysis.