3. Results and Discussion
3.2. Molecular dynamic simulations of the CTLD of perlucin and MBP-‐A
3.2.1. Principles of molecular dynamic simulations
The very basic principle (apart from sophisticated/efficient implementations of the algorithms) of a molecular dynamic simulation is to evaluate the force acting on each particle of the simulation from a given potential at a given time and calculate from this
forces new positions of all particles. The force acting on a single particle 𝑝𝑝 can be calculated from
𝐹𝐹!(𝑟𝑟! 𝑡𝑡 ) = 𝑚𝑚!⋅𝑑𝑑!𝑟𝑟! 𝑡𝑡 )
𝑑𝑑𝑡𝑡! = −∇ 𝑉𝑉(𝑟𝑟! 𝑡𝑡 ) (3.2.1.)
where ∇ 𝑉𝑉(𝑟𝑟! 𝑡𝑡 ) denotes the gradient of the potential energy at the location 𝑟𝑟!(𝑡𝑡) of the particle 𝑝𝑝.
During a MD simulation the positions and velocities of an evolving system under a given potential function are calculated at small successive time intervals 𝛥𝛥𝑡𝑡. Cuendet and van Gunsteren (Cuendet & van Gunsteren [2007]) give some introductory remarks about two common integration algorithms (see also Allen & Tildesley [1987]): the Verlet (see also Verlet [1967]) and leapfrog algorithm. To give an idea how the calculation of positions and velocities during an MD simulation can be performed in the following the first-‐order leapfrog algorithm is stated. Let 𝛥𝛥𝛥𝛥 be the time step employed during the simulation and set 𝑡𝑡! = 𝑛𝑛 ⋅ 𝛥𝛥𝛥𝛥 the time after 𝑛𝑛 time steps. The acceleration 𝑎𝑎! 𝑡𝑡! of a particle p at position 𝑟𝑟! 𝑡𝑡! and the force – resulting from the potential energy of the particle – are related through Newton’s law
𝑎𝑎! 𝑡𝑡! =𝐹𝐹!(𝑡𝑡!)
𝑚𝑚! = − 1
𝑚𝑚!∇ 𝑉𝑉(𝑡𝑡!) (3.2.2.)
The first order (with respect to the time step) leapfrog algorithm can then be stated as
𝑟𝑟! 𝑡𝑡!!! = 𝑟𝑟! 𝑡𝑡! + 𝑣𝑣! 𝑡𝑡!!!/! ⋅ 𝛥𝛥𝛥𝛥 (3.2.3.)
𝑣𝑣! 𝑡𝑡!!!/! = 𝑣𝑣! 𝑡𝑡!!!/! −∇ 𝑉𝑉 𝑡𝑡!
𝑚𝑚! ⋅ 𝛥𝛥𝛥𝛥 (3.2.4.)
These two equations allow the calculation of the position and velocity of particles in a MD simulation. Note the characteristic shift of 𝛥𝛥𝛥𝛥 2 between the calculated velocities and positions. Examples of higher order approximations can be found in Cuendet et al.
(Cuendet & van Gunsteren [2007]).
The potential energy used during the simulations with the “sander” software module of AMBER can be described as (e.g. Ponder & Case [2003], Cornell et al. [1995], Wang et al. [2000], Duan et al. [2003])
𝑉𝑉 𝒓𝒓 = 𝐾𝐾! 𝑏𝑏 − 𝑏𝑏! !
!"#$%
+ 𝐾𝐾! 𝜃𝜃 − 𝜃𝜃! !
!"#$%&
+ 𝑉𝑉!
!"!!"#$%& 2
(1 + cos(𝑛𝑛 𝜙𝜙 − 𝛿𝛿))
+ 𝐴𝐴!"
𝑟𝑟!"!"−𝐵𝐵!"
𝑟𝑟!"! +𝑞𝑞! 𝑞𝑞! 𝜀𝜀 𝑟𝑟!"
!"!#"!$%$
!"#$% !"; !!!
(3.2.5.)
Here 𝒓𝒓 denotes the 3𝑁𝑁 Cartesian coordinates of a 𝑁𝑁 particle system. The individual terms have the following meanings. The first sum describes the contribution of the covalent bonds to the potential energy. The potential energy is described as a harmonic potential with equilibrium bond length 𝑏𝑏! between two atoms and force constant 𝐾𝐾!. The next sum considers the angle between three covalent bonded atoms. Consider the water molecule as an example. The hydrogen and oxygen atoms lie in a plane. Then the angle between both of the O-‐H bonds would be the angle 𝜃𝜃 mentioned in 3.2.5. Again a harmonic potential is used to model these interactions with the equilibrium angle 𝜃𝜃! and force constant 𝐾𝐾!. The third sum runs over all dihedral angles. Here four consecutive covalently bound atoms are involved in the formation of one dihedral angle (see Fig. 3.2.19. for more information on dihedral angles). In general the potential energy of one dihedral angle is composed of at least one summand of a Fourier series. This is indicated by the integer 𝑛𝑛 appearing as index of the torsional force parameter 𝑉𝑉! and as the factor of the angle variable 𝜙𝜙. Note that the references stating equation 3.2.5. the phase 𝛿𝛿 is not indexed. A glance at the parameter files (parm99.dat and frcmod.ff03 supplied with AMBER 10) used here reveals that the phase is either 0° or 180° – this is also stated by Duan et al. (Duan et al. [2003]) – and depends on 𝑛𝑛 and the atoms involved in the dihedral angle formation.
The last sum contains three different interactions where 𝑟𝑟!" is the distance between to atoms 𝑖𝑖 and 𝑗𝑗 that are not covalently bound to each other. The first two summands constitute the “Lennard-‐Jones” interaction (see e.g. Lennard-‐Jones [1931]). The term proportional to 𝑟𝑟!!" describes the repulsion of two atoms within short distance due to
Pauli’s exclusion principle. The attractive contribution proportional to 𝑟𝑟!! is the van der Waals interaction. It arises from correlated electric dipole fluctuations (see e.g.
Lennard-‐Jones [1931], Leckband & Israelachvili [2001]). The final contribution is the electrostatic or Coulomb interaction between two point charges. Note that the usual factor 1/4𝜋𝜋𝜀𝜀! is omitted here. In many of the references cited here it does not appear.
At least two reasons are possible. On the one hand constant factors depend on the used system of units (see e.g. Jelitto [1994]) and on the other hand they might be omitted for improved readability and considered only in final results or numbers (see e.g.
Appendix A in Deserno & Holm [1998]). In explicit solvent calculations the dielectric constant is set to 𝜀𝜀 = 1.
In the last sum all atom pairs are considered (only once) that are not covalently bonded and are separated at least by three bonds. The Lennard-‐Jones and Coulomb interactions between atoms that are separated by three covalent bonds are scaled by a factor of 1/2.0 and 1/1.2 (see Cornell et al. [1995], Duan et al. [2003]). Usually the Lennard-‐Jones interactions between non-‐bonded atoms are calculated only if their distance is below a certain threshold (here: 9 Å during the heating phase and 10 Å during the unconstrained MD simulations). Reasons may be the reduction of computational expense and obeying the minimum image convention (see e.g. Allen &
Tildesley [1987]). In terms of the Coulomb interaction this cut-‐off value indicates up to which distance the electrostatic interaction is calculated as a direct space sum (see one of the following paragraphs on the PME method).
The entity of necessary parameters in equation 3.2.5. is termed “force field”. In summary equation 3.2.5. describes a “minimalist” (Cornell et al. [1995], p. 5181) potential function that is able to describe biomolecules in conjunction with an appropriate solvation model (Lee & Duan [2004]). It has fixed charges and polarization effects are not included. Polarizable force fields (see e.g. Halgren & Damm [2001], Ponder & Case [2003] for introductory remarks) are not considered further. The used force field is a revision made by Duan et al. (Duan et al. [2003]) of the “parm99” force field of Wang et al. (Wang et al. [2000]). The latter one is based on the “Cornell” force field (Cornell et al. [1995]).
Non-‐bonded interactions – Periodic boundary conditions
As it can be seen in equation 3.2.5. the calculation of the non-‐bonded interactions of a 𝑁𝑁 particle system would require the evaluation of interactions in the order of 𝑁𝑁!. This complexity can be reduced through cut-‐off distances in the sense that the non-‐bonded interaction is only calculated between atoms with a distance below this cut-‐off.
To give an illustrative (neglecting any screening effects) example: given the parameters for a Cα atom the Lennard-‐Jones interaction energy between two Cα atoms at their equilibrium distance 𝑑𝑑! = 3.816 Å is 𝜖𝜖! = −0.1094 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘/𝑚𝑚𝑚𝑚𝑚𝑚. When the distance between these two atoms is 𝑑𝑑 = 10 Å then the interaction energy increased to 𝜖𝜖 ≈ 0.005 ⋅ 𝜖𝜖!. In contrast the Coulomb energy of two particles with +0.5 unit charge (SI units) each, with a distance of 3.8 Å is ≈ 22 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘/𝑚𝑚𝑚𝑚𝑚𝑚. If these particles have a distance of 10 Å then their Coulomb interaction energy is still ≈ 0.38 ⋅ 22 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘/𝑚𝑚𝑚𝑚𝑚𝑚.
Compared to Lennard-‐Jones interactions the Coulomb interaction is classified as a long-‐range interaction. Allen and Tildesley (Allen & Tildesley [1987]) state that an interaction that decays not faster than 𝑟𝑟!𝒟𝒟(𝒟𝒟 is the system dimension) is frequently defined as a long-‐range interaction. Obviously the Coulomb interaction needs special attention. Cheatham and Brooks (section 4 in Cheatham III. & Brooks [1998]) state that a simple cut-‐off scheme of the long-‐range electrostatic interactions can lead to artefacts such as unrealistic pairing distances of equal charged ions that do not occur with different treatments.
In this work here a biomolecule is simulated in a water box with explicit water molecules to mimic a natural aqueous environment. However a single, isolated water box would be surrounded by a vacuum and therefore introduce an artificial interface not present in real aqueous solutions. Cheatham and Brooks (see section 3 in Cheatham III. & Brooks [1998]) describe that at this interface a surface tension can build up that leads to unwanted effects like water ordering and high pressures in the simulation cells.
To circumvent this the (primary) water box with the solute is repeated periodically in the three dimensional space. This introduces periodic boundary conditions (PBC). Now molecules close to the boundary can interact with molecules in neighbouring boxes (see e.g. Allen & Tildesley [1987]) and therefore are in an environment that is more similar to the environment further away from the edges of the primary box. Since
every molecule in the primary box can interact now with its periodic image an artificial periodicity is introduced into the system.
The Coulomb interaction energy of a system comprising 𝑁𝑁 particles and an infinite number of (periodic) images can now be written as (see e.g. Allen & Tildesley [1987], Deserno & Holm [1998])
𝑉𝑉!"#$ = 1 2
𝑞𝑞! 𝑞𝑞! 𝑟𝑟!"+ 𝑛𝑛
~
!
!
!,!
(3.2.6.)
Here the inner sum runs over all vectors 𝑛𝑛 that describe the position of every image cell relative to the primary cell. The tilde indicates that if 𝑛𝑛 = 0 (primary cell) it must be 𝑖𝑖 ≠ 𝑗𝑗 in the outer sum. The calculation of the electrostatic interactions with PBC requires that the primary cell is neutral otherwise the electrostatic energy of the whole system diverges (see e.g. Allen & Tildesley [1987] or Sagui & Darden [1999]). As described for example by Allen and Tildesley (Allen & Tildesley [1987]) the central idea to calculate the electrostatic energy is in a first step to introduce counter charges around the point charges of the atoms. These charge densities can have the shape of a Gaussian, like
𝜌𝜌! 𝑟𝑟 = 𝑞𝑞! 𝛽𝛽!
𝜋𝜋! exp(−𝛽𝛽! 𝑟𝑟!) (3.2.7a.)
where 𝛽𝛽 determines the spatial extent of the charge distribution. These counter charges screen the interactions of the point charges. It is assumed that the Coulomb interactions between the screened point charges decrease sufficiently fast with increasing charge distance. Therefore these interactions are evaluated in the direct space with a cut-‐off. The artificially introduced Gaussian charge distributions must be cancelled by opposite charge distributions so that the original potential is not changed.
The latter charge distributions are summed in the reciprocal space.
One can also motivate the improved treatment of long-‐range electrostatic interactions by the following identity (see e.g. Deserno & Holm [1998], Sagui & Darden [1999])
1 𝑟𝑟 =
erfc(α 𝑟𝑟)
𝑟𝑟 +1 − erfc(α 𝑟𝑟)
𝑟𝑟 (3.2.7b.)
with the complementary error function (erfc). This allows to split the Coulomb interaction energy into faster converging sums. The first summand in 3.2.7b. is evaluated in the direct space when 𝑟𝑟 is smaller than a specified cut-‐off (this is the same cut-‐off of 9 Å and 10 Å respectively introduced above for the Lennard-‐Jones interaction). The second one is evaluated in the reciprocal space.
Such a decomposition scheme due to Ewald (in (Kittel [2006]) in “Anhang B” the Ewald summation is shown for a crystal) results in
𝑉𝑉!"#$ = 𝑉𝑉!"#+ 𝑉𝑉!"#+ 𝑉𝑉!"## (3.2.8.)
The energy contribution calculated in the reciprocal space 𝑉𝑉!"# involves a Fourier transformation of the charge distribution in the primary cell. This can be evaluated with a three dimensional fast Fourier transformation (FFT) when the charges in the primary cell are mapped onto a three dimensional mesh. This is the Particle Mesh Ewald (PME) method (Darden et al. [1993]) that allows the calculation of the long-‐
range Coulomb interaction with a 𝑁𝑁 ⋅ ln 𝑁𝑁 order algorithm implemented in AMBER (Essmann et al. [1995], a flow chart can be found in Darden et al. [1999]). Toukmaji and Board reviewed different Ewald summation techniques in (Toukmaji & Board Jr.
[1996]).
The artificial periodicity of the system necessary for the Ewald methods can lead to artifacts, e.g. hindered rotation of dipoles, but they seem to be negligible for water environments with high dielectric constants (Smith & Pettitt [1996]). In summary the PME method (or similar techniques) for inclusion of the long-‐range electrostatics seems to be necessary for simulations of macromolecules and introduced errors are outweighed (e.g. Cheatham III. et al. [1995], York et al. [1993], Cheatham III. & Brooks [1998]).
Explicit water – TIP3P water model
The aqueous environment incorporating the protein is represented here through water molecules (“explicit water”) and not as a dielectric continuum (“implicit water”) surrounding the solute (see Fig. 3.2.1.A). Following the review on computer simulations on water by Guillot (Guillot [2002]) water models can be roughly
distinguished in terms of rigidity/flexibility of the internal water bonds, polarizability/non-‐polarizability, (point) charge distribution on a water molecule and Lennard-‐Jones interaction parameters. Usually the parameters of water models are adjusted to reflect the thermodynamical (for example density and heat of vaporization), structural (for example oxygen-‐oxygen radial distributions) and dielectric properties of real water during the molecular dynamic simulation.
Figure 3.2.1. The left image A) shows a perlucin molecule immersed in a water box. The three dimensional shape of the water box is that of a truncated octahedron. The right image B) highlights some physico-‐chemical features of the TIP3P (Jorgensen et al. [1983]) water model.
TIP3P is a rigid and non-‐polarizable water model with three point charges centred on the oxygen and hydrogen atoms respectively. Charges are given as fractional charges of the elementary charge. For the computation of the Lennard-‐Jones interaction between a water molecule and any other atom only parameters assigned to the oxygen atom are used. This means that the Lennard-‐Jones parameters for oxygen include the effect of the hydrogen atoms implicitly. To highlight this property a dashed circle is depicted around the oxygen atom incorporating the hydrogen atoms. According to the “parm99” force field supplied with AMBER the value of this van der Waals-‐radius is 𝑟𝑟! = 1.7683 Å. See also III.H. for further information on the term “van der Waals-‐radius”. The molecules are rendered with VMD (Humphrey et al.
[1996] version 1.9.1) and labels are added with Inkscape (http://inkscape.org).
In this thesis the TIP3P water model (see Fig. 3.2.1.B) Jorgensen et al. (Jorgensen et al.
[1983]) was used. The bonds are rigid and it contains three point charges corresponding on the oxygen and hydrogen atoms respectively. TIP3P is a non-‐
polarizable model. Lennard-‐Jones parameters are assigned to the oxygen atom including the hydrogen atoms implicitly. As pointed out by Price and Brooks (Price &
Brooks [2004]) the TIP3P water model was developed under conditions with short spherical cut-‐offs for Lennard-‐Jones and electrostatic interactions. They compared thermodynamic and dielectric properties of TIP3P water obtained from MD simulations using spherical cut-‐offs for Lennard-‐Jones and Coulomb interactions as well as a PME treatment of electrostatics. Since in both cases the determined properties differed from experimental values the authors refined the TIP3P model.
They employed the original and refined TIP3P water models in MD simulations of native protein structures and concluded “[…] that while three proteins do not behave appreciably different from a previous study, there is some subtle reduction in the hydrophilicity of alcohols and presumably other polar compounds […]” (p. 10103). For example they observed the maximal difference between the Cα-‐RMSd values of protein simulations with different water models to be 0.21 Å. The maximal observed SASA difference was around 200 Å!.
Here the original TIP3P water model was used for all MD simulations. This was solely justified by the use of TIP3P by Wang et al. (Wang et al. [2000]) and Duan et al. (Duan et al. [2003]) to test the “parm99” force field and its revised version respectively.
The periodic boundary conditions required for PME put restrictions on the shape of the water box that contains the macromolecule. It must be possible to arrange an infinite number of boxes to cover the whole space without any “vacuum gaps” between them. AMBER offers cuboid shaped boxes or truncated octahedral (imagine an octahedron with chopped off corners) boxes. Here the latter shape was used due to its smaller volume (see Fig. 3.2.1.A). The box volume is chosen in such a manner that at least a distance of 10 Å was preserved between any solutes atom and the box boundary.
Ion parameters
Since the PME method requires charge neutrality inside the water box usually counter ions have to be added to account for a net charge of the protein. Additionally the perlucin model was tested with a different number of associated Ca2+ ions (see also section 3.1.), which usually required charge neutralization. Through the addition of ion
pairs a certain ionic strength of the solvent can be maintained. Effects of ionic strength on the model stability were not investigated here. Sodium and/or chloride ions were used for charge balancing. The ionic strength was not equal in different MD simulation series. The chloride, sodium and calcium ion Lennard-‐Jones parameters in the
“parm99” force field are taken from Smith et al. (Smith & Dang [1994], chloride) and Åqvist (Åqvist [1990], cations). The sodium Lennard-‐Jones parameters in the force field are already adjusted for the TIP3P water model (see e.g. the discussion in Ross &
Hardin [1994]). The adjustments – as explained in the next paragraph – for the calcium ion parameters needed to be done separately with a force field modification file (see Appendix III.I.).
To explain this adjustment one need to know in the first place that the Lennard-‐Jones parameters in the “parm99” force field are not given as the 𝐴𝐴 and 𝐵𝐵 values from equation 3.2.5. Instead the Lennard-‐Jones interaction is formulated terms of the interaction energy minimum (well depth) 𝜖𝜖 and position of this minimum (or equilibrium distance) 𝑟𝑟! between two atoms (see Appendix III.H.). AMBER calculates the equilibrium distance between two different atom types simply as the sum of the individual atomic radii (here also termed van der Waals-‐radii). Now one has to decide which radii and well depth one assigns to the ions depending on which properties the simulated system has to have. Additionally Åqvist (Åqvist [1990]) determined the ion parameters not within the TIP3P water model. The radii of the cations given in the
“parm99” force field are already calculated in conjunction with the parameters of the TIP3P water model and the AMBER mixing rules. As stated by Ross et al. (Ross &
Hardin [1994], p. 6074) these parameters reproduce first peaks of the radial distribution of ion-‐water oxygen distance and hydration free energies. An exemplary calculation for Ca2+ can be found on the AMBER webpage (http://ambermd.org/Questions/vdw.html, last access 02/07/2013).
Thermostats and barostats
Typically biological processes occur at ambient conditions with respect to temperature and pressure. In a molecular dynamics simulation executed in a box of constant volume 𝑉𝑉 and with a constant number 𝑁𝑁 of particles the total energy 𝐸𝐸 would be conserved (NVE conditions). This means that the microstates of the system during the simulation sample the microcanonical ensemble. But conditions with constant temperature 𝑇𝑇 and
constant pressure 𝑃𝑃 are more appropriate to simulate a biological environment or a typical biochemical laboratory experiment. Since the particle number 𝑁𝑁 is usually constant as well these conditions are denoted as NPT conditions (an isothermal-‐
isobaric ensemble is sampled).
Consequently there is a need to include thermostat and barostat algorithms into the MD simulations. Recently Hünenberger reviewed (Hünenberger [2005]) fundamentals of current thermostat algorithms. The following introductory remarks about thermostats are based on this review and/or the book of Allen and Tildesley (Allen &
Tildesley [1987]) if not stated otherwise.
The equipartition theorem (see e.g. Reif [1987]) states that if a system is in thermodynamic equilibrium – at temperature 𝑇𝑇 – every degree of freedom that is quadratic in one generalised coordinate or momentum contributes the average energy of
𝐸𝐸! = 1
2 𝑘𝑘! 𝑇𝑇 (3.2.9.)
to the energy of the whole system (𝑘𝑘! is the Boltzmann constant). This enables to relate the macroscopic temperature 𝑇𝑇 of a system with its the average kinetic energy
𝐸𝐸!"# = 𝜉𝜉
2 𝑘𝑘! 𝑇𝑇 (3.2.10.)
when 𝜉𝜉 degrees of freedom contribute to the kinetic energy. This leads to the definition of a “instantaneous kinetic temperature” 𝑇𝑇 at each time step during the simulation as
𝐸𝐸!"# = 1
2 𝑚𝑚! 𝑣𝑣!!
!!!
!!!
= 𝜉𝜉
2 𝑘𝑘! 𝑇𝑇 (3.2.11.)
whose average corresponds to the desired macroscopic temperature (Allen & Tildesley [1987], p.46-‐47). Consequently a thermostat can be realized by adjusting the velocity distribution of the atoms simulated.
The “sander” software module of AMBER (version 10) offers three different thermostats (each is discussed by Hünenberger Hünenberger [2005]). The “stochastic-‐
coupling” scheme by Andersen (Andersen [1980]) assigns at intervals new velocities to
atoms of a system and these new velocities are sampled from a Maxwell-‐Boltzmann distribution. Another possibility to maintain a constant temperature is the integration of the Langevin equation (with appropriate properties of the stochastic forces and friction coefficients) instead of the Newton equation of motion. For the simulations discussed here the “weak-‐coupling” scheme of Berendsen et al. (Berendsen et al.
[1984]) was used. In this scheme the modified equation of motion for each particle 𝑖𝑖 reads
𝑎𝑎! 𝑡𝑡 = −∇ 𝑉𝑉 𝑟𝑟! 𝑡𝑡
𝑚𝑚! + 1
2 𝜏𝜏! 𝑇𝑇
𝑇𝑇 𝑡𝑡 − 1 𝑣𝑣!(𝑡𝑡) (3.2.12.)
In this equation the last summand implies a velocity scaling through coupling of the system to an external bath with temperature 𝑇𝑇. The coupling strength can be adjusted with the time constant 𝜏𝜏!. A large time constant simulates a weak coupling that vanishes in the limit of 𝜏𝜏! →∞ which corresponds to a simulation without temperature control. A time constant in the order of the time step of the simulation reflects a tight coupling to the heat bath that does not allow temperature fluctuations (Hünenberger [2005]).
Note that as long as the time constant has a value between these two limiting cases – 𝜏𝜏! ≈ 𝛥𝛥𝛥𝛥 or 𝜏𝜏! →∞ as described above – a simulation with the Berendsen algorithm samples configurations from a particular “weak-‐coupling ensemble” (Morishita [2000]) and not from an expected canonical ensemble (see also Hünenberger [2005]).
Another issue might arise when using velocity rescaling schemes for temperature control: the “flying ice cube” effect (Harvey et al. [1998]). The authors exemplify how kinetic energy can be transferred from vibrational motions to centre-‐of-‐mass translations. As remedies they propose occasional velocity reassignment (similar to the
“stochastic-‐coupling” scheme), frequent removal of the centre-‐of-‐mass translation and rotation or ideally calculate the systems temperature as an appropriate time average and not as an instantaneous temperature. Currently the author of this thesis cannot comment in how far these (or other) precautions are established in AMBER 10. No efforts were undertaken to explore possible issues with the weak-‐coupling scheme.
Despite these drawbacks the weak-‐coupling algorithm is used for temperature control since it can also be used for pressure regulation (AMBER 10 offers only the weak-‐
coupling scheme for pressure regulation). Instead of velocity rescaling the pressure
regulation according to Berendsen et al. (Berendsen et al. [1984]) involves coordinate rescaling. For a simple isotropic and cubic system the coordinates and box lengths are scaled by
1 −𝛽𝛽 𝛥𝛥𝛥𝛥
3 𝜏𝜏! 𝑃𝑃 − 𝑃𝑃 𝑡𝑡 (3.2.13.)
where 𝛽𝛽 is the isothermal compressibility. The magnitude of the scaling depends on the difference between the target pressure 𝑃𝑃 and actual pressure 𝑃𝑃(𝑡𝑡) of the system. As in the case of the thermostat the coupling strength of the system to the barostat can be controlled with 𝜏𝜏!.
MD simulation parameters in this thesis
After setting up the system, which means solvation of the biomolecule in the box of explicit water, no temperature could be assigned to the system since there was no velocity information present. The system had to be adapted to the target conditions which were 𝑇𝑇 = 300𝐾𝐾 ≈ 27°𝐶𝐶 and 𝑃𝑃 = 1 𝑏𝑏𝑏𝑏𝑏𝑏 = 1000 ℎ𝑃𝑃𝑃𝑃 for every MD simulation discussed here. So in a first step a simple energy minimization was performed to relax the solute and the solvent molecules. After this minimization velocities were assigned randomly to all atoms in such a manner that the velocity distribution follows a Maxwellian distribution for a given temperature (25 𝐾𝐾 for the perlucin model and 1 𝐾𝐾 for the reference protein MBP-‐A). Afterwards the system was “coupled” to a thermostat and a barostat. The pressure was set to 𝑃𝑃 = 1 𝑏𝑏𝑏𝑏𝑏𝑏 and the temperature was increased in several short MD simulations to 𝑇𝑇 = 300𝐾𝐾. During this heating phase the time constants were kept small to ensure a tight coupling. Harmonic positional restraints were applied to the atoms of the solute to avoid an unfolding of the protein due to fast heating. This means that the positions of the atoms of the solute (not the water molecules) were restrained with a harmonic potential on their positions after the minimization for the first MD simulation. For the following restrained MD simulations the reference coordinates for the restrained atoms were the coordinates obtained from the first MD simulation. After the target temperature of 𝑇𝑇 = 300𝐾𝐾 was established these restraints were gradually released. Table 3.2.1. summarizes the
important parameters of the MD simulations that were performed subsequently on each initial configuration.
The minimization and heating protocol was kindly provided by Prof. Dr. M. Zacharias.
time
temperature
coupling
constant restraint
weight non-‐bonded cut-‐off initial minimization (steepest descent, 1500 steps) and random velocity
assignment from a Maxwellian distribution (25 𝐾𝐾 or 1 𝐾𝐾) 0 to 20 𝑝𝑝𝑝𝑝 100 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 25 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 20 to 40 𝑝𝑝𝑝𝑝 200 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 25 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 40 to 60 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 25 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 60 to 80 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 12 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 80 to 100 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 6 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 100 to 120 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 2 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 120 to 140 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 1 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å 140 to 200 𝑝𝑝𝑝𝑝 300 𝐾𝐾 0.1 𝑝𝑝𝑝𝑝 0.5 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘
𝑚𝑚𝑚𝑚𝑚𝑚 Å! 9 Å
0.2 to 2.2 𝑛𝑛𝑛𝑛 300 𝐾𝐾 2.5 𝑝𝑝𝑝𝑝 -‐ 10 Å
2.2 to 10.2 𝑛𝑛𝑛𝑛 300 𝐾𝐾 2.5 𝑝𝑝𝑝𝑝 -‐ 10 Å
Table 3.2.1. Summary of the subsequent stages of the MD simulations performed for every initial configuration. The first column indicates the time progression of the whole simulation and implies the duration of each stage. Note that the initial energy minimization is no molecular dynamics simulation and therefore no time can be assigned. The second column contains the temperature of the heat bath. At every stage the target pressure is 1 𝑏𝑏𝑏𝑏𝑏𝑏. The next column highlights the coupling constants used for the temperature and pressure regulation.
Note that the same values were chosen for the thermostat and barostat. In the column labelled
“restraint weight” the parameter that indicates the strength of the harmonic potential acting on the restrained atoms. The last column contains the cut-‐off value of the non-‐bonded interactions. Lennard-‐Jones interactions are not calculated between atoms with a distance greater than this value. The electrostatic interactions between atoms up to this distance are calculated in the direct space sum.
After the system was heated to the desired temperature it was allowed to evolve unconstrained for a longer time. Here the duration of the unconstrained simulation
was 10 𝑛𝑛𝑛𝑛 in total (split up in 2 and 8 𝑛𝑛𝑛𝑛). The time step used in every simulation was 𝛥𝛥𝛥𝛥 = 2 𝑓𝑓𝑓𝑓 = 0.002 𝑝𝑝𝑝𝑝.
The following simple estimation points to a possible problem with this large time step (see also Feenstra et al. [1999]). Consider a covalent bond between the oxygen and hydrogen atom of a hydroxyl group. The force constant (or spring constant) for the harmonic potential that models the covalent bond is 𝑘𝑘 = 553 𝑘𝑘𝑘𝑘𝑘𝑘𝑘𝑘/𝑚𝑚𝑚𝑚𝑚𝑚 Å! (“parm99”
force field; bond between atom types “OH” and “HO”). The position of the “heavy”
oxygen shall be fixed. Then the frequency of the oscillating hydrogen atom around the equilibrium distance can be estimated to be 𝜈𝜈! = 1 𝑇𝑇! = 1 2𝜋𝜋 𝑘𝑘 𝑚𝑚 ≈ 7.7 ⋅ 10!" 𝑠𝑠!!. This implies that the time step of 𝛥𝛥𝛥𝛥 = 2 𝑓𝑓𝑓𝑓 corresponds approximately to 0.15 ⋅ 𝑇𝑇!. Feenstra et al. (Feenstra et al. [1999]) give several further examples of characteristic oscillation periods.
These fast bond stretching movements involving hydrogen atoms were removed from the MD simulations with the SHAKE (Ryckaert et al. [1977]) and SETTLE (Miyamoto &
Kollman [1992]) algorithms by imposing distance constraints. The latter one is for rigid water models like TIP3P. This allows the use of the large time step of 𝛥𝛥𝛥𝛥 = 2 𝑓𝑓𝑓𝑓.
This introduction to the principles of molecular dynamic simulations is concluded with a summary of the characteristics of the different simulations whose results will be presented and discussed in the next sections.
feature perlucin MBP-‐A
identifier run09 run21 run22 run02 run07 run10
no. runs 6 3 3 3 3 3
residues 1-‐135 1-‐131 1-‐131 109-‐221 104-‐221 104-‐221 Ca2+ Ca-‐1
Ca-‐2 Ca-‐3 Ca-‐4
-‐
Ca-‐2 -‐
Ca-‐4
-‐ -‐ Ca-‐1
Ca-‐2 Ca-‐3 -‐
-‐
Ca-‐2 -‐
-‐
ions (c) 7 Cl-‐ 3 Cl-‐ 1 Na+ 4 Na+ 4 Cl-‐ -‐
ions (a) -‐ -‐ 3 Na+/Cl-‐ -‐ -‐ 3 Na+/Cl-‐
no. water 8039 4888 5799 4554 5623 7150
restrained
residues 1-‐139 protein and Ca2+
1-‐131 protein only
1-‐131 protein only
protein
and ions protein
and ions protein only