Step 8
The first step in solvating a system is to determine the size of the solvent box and the number of counterions, if any, that are to be added. This is achieved with the following simple program:
# . Get the system.
system = Unpickle ( "step7.pkl" )
system.Summary ( )
# . Determine solvation parameters.
DetermineSolvationParameters ( system )
The function DetermineSolvationParameters prints out information about possible (orthorhombic) solvation boxes for the system, including the numbers of counterions that it is necessary to add to have a solution of a particular molarity. It is usual to have a solvent layer of sufficient thickness such that proteins within adjacent boxes interact as little as possible. For PKA, we choose to have a solvent buffer of 10 Ă… minimum between the protein and the edge of the box (i.e. at least 20 Ă… between protein molecules) which means having a box with approximate dimensions 75 Ă— 60 Ă— 60 Ă…3.
The total charge of the PKA system is +3 so 3 monovalent anions would need to be added to have a simulation system that is neutral overall. However, we also would like to mimic the molarity of the intracellular environment within which the protein normally works. If this is taken to be 0.154 M, then approximately 24 monovalent cations and 27 monovalent anions are required for a box of the specified volume.
A solvent box of the appropriate size could be constructed using a program similar to that of Example 27 in the pDynamo distribution which uses a Monte Carlo refinement procedure. Here, however, a geometry optimization/molecular dynamics procedure is employed. The program is straightforward and is:
# . Parameters.
_Density = 1000.0 # . Density of water (kg m^-3).
_Refine = True
_Steps = 10000
# . Box sizes.
_XBox = 75.0
_YBox = 60.0
_ZBox = 60.0
# . Define the solvent MM and NB models.
mmModel = MMModelOPLS.WithParameterSet ( "bookSmallExamples" )
nbModel = NBModelCutOff.WithDefaults ( )
# . Define the solvent molecule.
molecule = ImportSystem ( "water.mol" )
molecule.Summary ( )
# . Create a symmetry parameters instance with the correct dimensions.
symmetryParameters = SymmetryParameters ( )
symmetryParameters.SetCrystalParameters ( _XBox, _YBox, _ZBox, 90.0, 90.0, 90.0 )
# . Create the basic solvent box.
solvent = BuildSolventBox ( CrystalSystemOrthorhombic ( ), symmetryParameters, molecule, _Density )
solvent.label = "Water Box"
solvent.DefineMMModel ( mmModel )
solvent.DefineNBModel ( nbModel )
solvent.Summary ( )
solvent.Energy ( )
# . Refine the system using minimization and dynamics.
if _Refine:
# . Minimization.
ConjugateGradientMinimize_SystemGeometry ( solvent ,
maximumIterations = 200 ,
logFrequency = 1 ,
rmsGradientTolerance = 10.0 )
# . Define a normal deviate generator in a given state.
normalDeviateGenerator = NormalDeviateGenerator.WithRandomNumberGenerator ( RandomNumberGenerator.WithSeed ( 714717 ) )
# . Dynamics.
LangevinDynamics_SystemGeometry ( solvent ,
collisionFrequency = 25.0 ,
logFrequency = 100 ,
normalDeviateGenerator = normalDeviateGenerator ,
steps = _Steps ,
temperature = 300.0 ,
timeStep = 0.001 )
# . Final energy.
solvent.Energy ( )
# . Calculate and print the final density.
logFile.Paragraph ( "Solvent density = {:.2f} kg m^-3.".format ( SystemDensity ( solvent ) ) )
# . Save the water box.
Pickle ( "waterBox.pkl", solvent )
The function BuildSolventBox works by placing solvent molecules evenly, but in random orientations, throughout a box of appropriate shape and dimension. This is followed up by a geometry optimization to obtain a configuration of reasonable energy (with few or no intermolecular overlaps), and then a molecular dynamics simulation to more fully equilibrate the system.
The third stage in the solvation process is to add counterions to the protein system. The full program is:
# . Parameters.
# . Box sizes.
_XBox = 75.0
_YBox = 60.0
_ZBox = 60.0
# . Number and type of ions to add.
_NNegative = 27
_NPositive = 24
_NegativeIon = "chloride"
_PositiveIon = "potassium"
# . Reorient option.
_Reorient = False
# . Get the system.
system = Unpickle ( "step7.pkl" )
system.Summary ( )
# . Reorient the system if necessary (see the results of GetSolvationInformation.py).
if _Reorient:
masses = Array.FromIterable ( [ atom.mass for atom in system.atoms ] )
system.coordinates3.ToPrincipalAxes ( weights = masses )
# . Get the positive and negative ions.
if _NNegative > 0:
anion = ImportSystem ( _NegativeIon + ".mol" )
anion.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )
anion.Summary ( )
if _NPositive > 0:
cation = ImportSystem ( _PositiveIon + ".mol" )
cation.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )
cation.Summary ( )
# . Add the counterions.
newSystem = AddCounterIons ( system, _NNegative, anion, _NPositive, cation, ( _XBox, _YBox, _ZBox ) )
# . Save the combined system.
Pickle ( "step8_a.pkl", newSystem )
After some initial parameter definitions, the program starts by reading in the protein system and then, optionally, reorienting the protein coordinates. This is not done here as the results of the program GetSolvationInformation suggest that this is not necessary to obtain a simulation box of reasonable size. After this, systems corresponding to a single anion (chloride) and a single cation (potassium) are defined and the function AddCounterIons is invoked which places the required number of anions and cations within a certain box around the protein. Finally, the new, merged protein/ion system is saved for subsequent use.
The function AddCounterIons places ions at random positions around the protein. It is, in principle, possible to use algorithms which use an energetic criterion for placement (i.e. choose positions with the lowest energy) but similar results can be obtained, if necessary, by carrying out Langevin molecular dynamics or Monte Carlo simulations on the combined protein/ion system. This would require using an implicit solvent model (e.g. an NB model with a dielectric appropriate to that of water – about 80), fixing the protein atoms, performing a simulation and then selecting energetically favorable ion configurations.
The final step in the solvation procedure is to merge the protein/ion system and the water box. As before it would be possible to employ a Monte Carlo refinement procedure similar to that of Example 28 in the pDynamo distribution. Here, though, geometry optimization and molecular dynamics is again used. The program is:
# . Parameters.
# . Refine options.
_Optimize = True
_Refine = True
_Steps = 2000
# . Define the NB model.
nbModel = NBModelCutOff.WithDefaults ( )
# . Get the solute system.
solute = Unpickle ( "step8_a.pkl" )
solute.Summary ( )
# . Get the solvent system.
solvent = Unpickle ( "waterBox.pkl" )
solvent.DefineNBModel ( nbModel )
solvent.Summary ( )
# . Create the solvated system.
solution = SolvateSystemBySuperposition ( solute, solvent, reorientSolute = False )
solution.label = "Solvated PKA"
solution.DefineNBModel ( nbModel )
solution.Summary ( )
# . Refinement.
if _Refine:
# . Check energies.
solute.Energy ( )
for system in ( solvent, solution ):
system.DefineNBModel ( nbModel )
system.Energy ( )
# . Get selections for the protein atoms in the A and I chains.
protein = ( AtomSelection.FromAtomPattern ( solution, "A:*:*" ) | AtomSelection.FromAtomPattern ( solution, "I:*:*" ) )
hydrogens = Selection.FromIterable ( [ index for ( index, atom ) in enumerate ( solution.atoms ) if atom.atomicNumber == 1 ] )
heavies = protein - hydrogens
# . Fix all protein atoms.
solution.freeAtoms = protein.Complement ( upperBound = len ( solution.atoms ) )
solution.Summary ( )
# . Optimization.
if _Optimize:
ConjugateGradientMinimize_SystemGeometry ( solution ,
maximumIterations = 200 ,
logFrequency = 10 ,
rmsGradientTolerance = 10.0 )
# . Define a normal deviate generator in a given state.
normalDeviateGenerator = NormalDeviateGenerator.WithRandomNumberGenerator ( RandomNumberGenerator.WithSeed ( 517152 ) )
# . Dynamics loop - fixed atoms then tether restraints.
for forceConstant in ( None, 500.0, 100.0, 20.0, 4.0 ):
# . Harmonically restrain protein heavy atoms.
tethers = None
if ( forceConstant is not None ):
reference = Clone ( solution.coordinates3 )
tetherEnergyModel = RestraintEnergyModel.Harmonic ( 0.0, forceConstant )
tethers = RestraintModel ( )
tethers["Tethers"] = RestraintMultipleTether.WithOptions ( energyModel = tetherEnergyModel ,
reference = reference ,
selection = heavies )
solution.DefineRestraintModel ( tethers )
# . Dynamics.
LangevinDynamics_SystemGeometry ( solution ,
collisionFrequency = 25.0 ,
logFrequency = 100 ,
normalDeviateGenerator = normalDeviateGenerator ,
steps = _Steps ,
temperature = 300.0 ,
timeStep = 0.001 )
# . Unfix atoms and remove tethers.
solution.freeAtoms = None
solution.DefineRestraintModel ( None )
# . Save the system.
Pickle ( "step8_b.pkl", solution )
As in Example 28, the program employs the function SolvateSystemBySuperposition to perform the solvation. This is followed by a refinement of the water and ion coordinates, keeping the protein atoms fixed, using geometry optimization, followed by a series of molecular dynamics simulations in which the protein heavy atoms are harmonically tethered by force constants of decreasing value.
The strategy used by SolvateSystemBySuperposition is to overlap the solute and solvent systems and then remove those solvent molecules which overlap with solute atoms. This is a straightforward method to implement but the solvated systems it produces depend crucially upon the initial configurations of the solute and solvent systems. This means, for example, that certain regions of the solute, such as internal cavities or nooks and crannies on its surface, may be insufficiently solvated. This problem can be alleviated by performing several cycles of solvation (by superposition) and, optionally, refinement using different solvent configurations. There also exist more specialized solvation algorithms, such as those that search for solvatable regions within or on the surface of a solute. These may be implement in future versions of pDynamo.