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).