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 can be constructed using a
program very similar to that of Example 27 in the pDynamo distribution.
The program is straightforward and is:

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 two
constant-volume Monte Carlo calculations, the first to obtain a
configuration of reasonable energy (with few or no
intermolecular overlaps) and the second 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:

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. This involves a procedure that is
very similar to that of Example 28 in the pDynamo distribution. The
program is:

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 both
geometry optimization and
molecular dynamics simulation.

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.