Explicit Solvent

Building an Explicit Solvent Model

In MELD, an explicit solvent simulation can be started either from a protein sequence or from pre-generated amber topology and coordinate files. As of yet, explicit solvent is not set to work from any other input files including PDB.

From Sequence

A simple explicit solvent model can be built in the setup.py script by first generating a protein sequence and setting up an explicit solvent SystemBuilder that defines the system:

p = subsystem.AmberSubSystemFromSequence("NALA ALA CALA")
b = builder.AmberSystemBuilder(solvation="explicit")

The system is then built from the molecular sequence according to the specifications of the builder:

self.system = b.build_system([p])

This set of commands will instruct MELD to use tleap to generate a solvent box with 6 Å between any atom in the protein and the edge of the box, and produce the necessary amber topology (system.top) and coordinate files (system.mdcrd). By default MELD uses the TIP3P water model, however 'spce', 'spceb', 'obc', 'tip4pew' are all allowed values for solvent_forcefield:

b = builder.SystemBuilder(explicit_solvent=True, solvent_forcefield="tip3p")

The distance between the protein and the edge of the box (in Å) can also be specified:

b = builder.SystemBuilder(explicit_solvent=True, solvent_distance=10)

The box is not automatically neutralized, which therefore requires the user to determine the charge of the molecule in order to accurately neutralize. Positive and negative ions can be added to the water box:

b = builder.SystemBuilder(explicit_solvent=True, explicit_ions=True, p_ioncount=5, n_ioncount=5)

By default the positive ions are Na+, and the negative ions are Cl-. Allowed values for p_ions include: Na+, K+, Li+, Rb+, Cs+, Mg+. The allowed values for n_ions include: Cl-, I-, Br-, F-. They can be specified by:

b = builder.SystemBuilder(explicit_solvent=True, explicit_ions=True, p_ion="K+", p_ioncount=5, n_ion="F-", n_ioncount=5)

As not all forcefields include parameters for all ion types, make sure to select ions that are supported by your forcefield.

From Amber Files

An explicit solvent system can also be generated from pre-made amber topology (.top) and coordinate (.crd) files. An example tleap script to generate a 4-mer peptide in a 10 Å TIP3P water box is included below:

xleap
        source leaprc.protein.ff14SB
        foo = {ALA ALA ALA ALA}
        source leaprc.water.tip3p
        solvatebox foo TIP3PBOX 10.0
        saveamberparm foo system.top system.crd
        quit

The system is then built in the setup.py script as:

s = system.builder.load_amber_system('system.top', 'system.crd')

When generating the amber files, keep in mind that MELD can only use the following specifications:

Parameter	Allowed Values
solvent_forcefield	‘tip3p’, ‘spce’, ‘spceb’, ‘obc’, ‘tip4ew’
p_ions	‘Na+’, ‘K+’, ‘Li+’, ‘Rb+’, ‘Cs+’, ‘Mg+’
n_ions	‘Cl-‘, ‘I-‘, ‘Br-‘, ‘F-‘

and is limited to rectangular box shapes.

Setting Options

Once an explicit solvent system is built either with the builder or from pre-generated amber files, it must be declared to the runtime options:

options = system.RunOptions(solvation='explicit')

Explicit Solvent and Replica Exchange

For a simple system, MELD’s regular Hamiltonian replica exchange algorithm is likely to be sufficient. Keep in mind that, due to the greater size of the system, this will require a larger number of replicas and a longer run time, with careful monitoring of replica exchanges. The exact parameters are system specific. For larger systems (such as proteins) it is advisable to use replica exchange with solute tempering (REST2). A full explanation of the method can be found at:

Wang, Lingle, Richard A. Friesner, and B. J. Berne. “Replica exchange with solute scaling: a more efficient version of replica exchange with solute tempering (REST2).” The Journal of Physical Chemistry B 115.30 (2011): 9431-9438. DOI : 10.1021/jp204407d.

In short, REST2 works by dividing the system into protein and water, and scaling the intramolecular potential energy function at a constant temperature. This is physically similar to changing the temperature of the protein while keeping the surroundings at the target temperature as we climb the temperature ladder. Therefore the acceptance probability scales only with the size of the biomolecule, and not with the number of water molecules.

In MELD, we employ REST2 by first declaring a temperature scaler for the target temperature of the system:

s.temperature_scaler = system.ConstantTemperatureScaler(300.)

We then define the REST2 scaler that will be applied to the potential energy function:

rest2_scaler = system.GeometricTemperatureScaler(alpha_min=0.5, alpha_max=1, temperature_min=300., temperature_max=450.)

This scaler can be geometric or linear. In the above example, the minimum and maximum temperatures will be turned into a scaling factor that when applied to the potential energy function achieves a physically similar result to scaling the temperature of the replicas from alpha_min to alpha_max across the ladder from temperature_min to temperature_max.

In the options, we declare that we are using REST2 by:

options.use_rest2 = True

And set the REST2 scaler:

options.rest2_scaler = system.REST2Scaler(300., rest2_scaler)

Periodic Boundary Conditions

Explicit solvent simulations in MELD make use of periodic boundary conditions (PBC’s) through OpenMM. In OpenMM, if a periodic box is enforced then the center of every molecule is translated so that it lies in the same periodic box. This means that unconnected molecules, say a peptide and a protein in a complex that are not bound together, could be translated / reimaged differently. In OpenMM all atoms involved in each bond are treated as a single molecule. Therefore, in MELD we group together all of the atom pairs that are restrained by: distance, hyperbolic distance, torsions, Gaussian mixture models (GMM’s), distance profiles, and torsion profiles. This creates a single “molecule” that will be reimaged when PBC’s are enforced. For example, by placing a distance restraint on a protein and the peptide it is complexed with, they become a single molecule in terms of PBCs. Confinement and Cartesian restraints work across periodic boundaries through OpenMM’s periodicdistance() function in the CustomExternalForce class.