Molecular Dynamics

To mimic the aqueous environment the Stochastic Boundary molecular dynamics method (SBMD) has been used. A particular method of solvation will be used due to this fact. We have used a sphere of pre-equilibrated waters with a radius of 24 Å to solvate the system. This sphere was centered on the stereogenic carbon of the substrate (C$_\alpha$). All the crystallographic waters beyond the sphere were removed along with any water which oxygen is closer than 2.5 Å near any heavy atom of the protein. A soft and deformable boundary potential is applied at the edge of the sphere of waters to mimic the effect of the inexistent bulk solvent.

The classical dynamics is carried out partitioning the system into three zones (see section 1.4.4 for more details in SBMD method). The dynamics region which consists of atoms within a distance of 20 Å from the center; the buffer region which contains the atoms surrounding the dynamics region from 20 Å up to 24 Å; and the reservoir region which includes the remainder of the system and is excluded in the explicit dynamics simulation. The sizes of the dynamics and buffer region are large enough to ensure that the active site and all the possible rearrangements in its surroundings during the reaction are adequately modeled.

The partition of the dynamical system is not done atom-wise, that is, all atoms of an aminoacidic residue are included in the dynamics region if any atom of the residue is within 20 Å from the reference point. If none of the atoms of a residue are in the 24 Å sphere the aminoacid will belong to the reservoir region. While the rest of enzyme atoms that may be found between 20 and 24 Å will be contained in the buffer region. The labels for the protein are assigned at the beginning of the simulation and kept all along the simulation. This will not be the case for the water molecules which will be permitted to diffuse between the dynamics and buffer region. The label list containing the waters of every region is updated every five steps during the simulation.

The parameters needed in the SBMD framework are taken from the original publications [186] which have been tested thoroughly in more recent works [187,189]. The trajectory in the dynamics region is propagated using the Newton's equation of motion, while in the buffer region the Langevin equations with a stochastic term are used. The constants for boundary forces are $K_C=1.30$, $K_O=1.22$, $K_N=1.30$ kcal/(mol $\cdot\AA^2$) for the backbone atoms and 0.73 kcal/(mol $\cdot\AA^2$) for the lateral chain heavy atoms. The friction coefficients are 200 ps$^{-1}$ for all protein atoms and 62 ps$^{-1}$ for oxygen atom in water molecules (this last value corresponds to the self-diffusion constant of bulk water at 300 K[183]). Both friction coefficients and boundary forces applied to protein atoms are scaled by a screening function that depends on the distance from the center (see section 1.4.4).

The final dimension of the different regions is:

ATOMS	protein + substrate	waters	Total
QM atoms	63	0	63
MM atoms	5406	2730	8136
Dynamics zone	3276	1242	4518
Buffer zone	1208	1488	2686
Reservoir zone	985	0	985
All	5469	2730	8199

The leapfrog algorithm to integrate the MD equations is used in all cases with a timestep of 1 fs. The hydrogen atoms are constrained using the SHAKE algorithm (page ). The process of heating and equilibrating the system can be outlined as follows:

1. The 24 Å sphere of waters is added on the protein three times at three different orientations. The corresponding remove criterion is applied every time.
2. A short steepest descent minimization is performed into the dynamics and buffer region to avoid bad contacts
3. A 5 ps NVT Langevin dynamics of the solvent is performed in order to equilibrate the waters.
4. The three previous steps are repeated to obtain a good solvation of the system.
5. 10 ps at 100, 200 and 298 K SBMD are performed to heat the system.
6. The final system is equilibrated during 50 ps at 298 K.

Figure 4.3: Schematic representation of the solvated Mandelate Racemase model