3.2 Molecular dynamics
Molecular dynamics (MD) simulation techniques directly solving Newton’s
equations of motion are commonly used to model systems of biomolecules
and biomaterials because they can track individual atoms and, therefore,
answer questions to specific material properties [21, 22]. In MD
simulations, the starting point is defining the initial coordinates and
initial velocities of the atoms characterizing the model system, for
example, the desired biomolecule plus the biologically relevant
environment, i.e., water molecules or other solvent and/or membranes.
The coordinates of the desired biomolecule can usually be found as
structural data (X-ray or NMR) deposited into the protein data bank
(PDB) [23] (www.pdb.org); otherwise, it
is possible to derive initial geometry and coordinate data from model
building techniques, including homology methods. This step also
typically includes the placement and positioning of the environment of
the molecules (solvation, ionic strength, etc.). The initial velocities
are typically derived from the Maxwell-Boltzmann distributions at the
desired temperature of the simulation. The potential of interactions of
each of the atoms is calculated using a force field function, which
parameterizes the non-bonded and bonded interaction terms of each atom
depending on its constituent atom connectivity: bond terms, angle terms,
dihedral terms, improper dihedral terms, non-bonded Lennard-Jones terms,
and electrostatic terms. The potential interactions are summed across
all the atoms contained in the system, to compute an overall potential
energy:
\(U\left(\overset{}{R}\right)=\sum_{\text{bonds}}{K_{b}\left(b-b_{0}\right)^{2}}+\sum_{\text{angles}}{K_{\theta}\left(\theta-\theta_{0}\right)^{2}+}\sum_{\text{dihedrals}}{K_{\chi}(1+\cos{(\eta\chi-\delta))}+\sum_{\text{impropers}}{K_{\phi}\left(\phi-\phi_{0}\right)^{2}+\sum_{\text{nonbonded}}{\left(\varepsilon_{\text{ij}}\left[\left(\frac{R_{\min_{\text{ij}}}}{r_{\text{ij}}}\right)^{12}-\left(\frac{R_{\min_{\text{ij}}}}{r_{\text{ij}}}\right)^{6}\right]\right)+\ \frac{q_{i}q_{j}}{\text{ϵr}_{\text{ij}}}}}}\)(Eq. 11)
Taking the derivative of the potential energy function yields the force,
and from Newton’s second law, this is equal to mass times acceleration.
Although the process seems simple, the derivative function results in a
set of 3N-coupled 2nd order ordinary differential
equations that must be solved numerically. The solution consists of a
numerical recipe to advance the positions and the velocities by one
timestep. This process is repeated over and over again to generate MD
trajectories of constant energy. Constant temperature dynamics are
derived by coupling the system to a thermostat using well-established
formulations such as the Langevin dynamics or the Nose-Hoover
methodologies [24]. Application of MD simulations to biomolecules is
facilitated by several popular choices of force fields such as CHARMM27
[25] (www.charmm.org), AMBER [26]
(www.ambermd.org), and GROMOS [27]
(www.gromacs.org), as well as dynamic
simulations packages and visualization/analysis tools such as NAMD
[28]
(www.ks.uiuc.edu/Research/namd/)
and VMD [29]
(www.ks.uiuc.edu/Research/vmd/).
MD simulations for commonly modeled molecules such as proteins, nucleic
acids, and carbohydrates that have well-established force fields can be
performed directly using a favorite software package such as LAMMPS,
GROMACS or HOOMD-blue
(http://glotzerlab.engin.umich.edu/hoomd-blue).