Introduction to Molecular Dynamics Simulation

Leave a Comment / Computation Biochemistry, Structural Bioinformatics (BDSL) / By

Introduction

Molecular dynamics (MD) simulation is a computational method for investigating the dynamic behavior of molecules, particularly biomolecules such as proteins, nucleic acids, lipids, and their complexes.

In an MD simulation, a molecular system is represented as a collection of atoms whose interactions are described using a predefined mathematical model called a force field. The biomolecule is usually placed in a simulation box containing solvent molecules, such as water, together with appropriate ions.

The forces acting on each atom are calculated from the force field, and the motion of the atoms is simulated by numerically solving Newton’s equations of motion. By repeatedly updating the positions and velocities of atoms over very short time intervals, the simulation generates a trajectory describing how the molecular system changes over time.

Experimental methods such as X-ray crystallography and cryo-electron microscopy often provide structural snapshots of biomolecules. MD simulation complements these experiments by providing information about molecular motion, flexibility, conformational transitions, interactions, and possible mechanisms of biological function.

However, an MD trajectory should not simply be regarded as a molecular movie. Under appropriate conditions and with sufficient sampling, it can also be used to characterize the distribution of conformational states accessible to the molecular system.

Discussion: What does an experimentally determined structure represent?

Is the structure deposited in the PDB necessarily the equilibrium structure of the molecule in solution? Experimental structures are obtained under particular experimental conditions and may represent one or several conformational states.

Classical Molecular Dynamics

Classical Mechanics

Conventional biomolecular MD simulations are based primarily on classical mechanics.

Atoms are treated as classical particles with properties such as mass, position, velocity, and partial electric charge. Their interactions are described using mathematical functions rather than by explicitly calculating electrons and molecular orbitals from quantum mechanics.

The force field therefore provides an effective approximation of the underlying quantum mechanical interactions.

This approximation makes it possible to simulate systems containing hundreds of thousands or even millions of atoms over relatively long timescales.

However, conventional classical MD generally assumes a fixed chemical topology. Consequently, processes involving substantial changes in electronic structure are usually not represented directly.

Examples include:

  • breaking chemical bonds
  • forming new chemical bonds
  • proton transfer
  • electron transfer
  • chemical reactions
  • substantial changes in electronic polarization

Such problems may require quantum mechanical calculations, reactive force fields, or hybrid quantum mechanics/molecular mechanics approaches such as QM/MM.

Discussion: What cannot conventional classical MD simulate?

Consider enzymatic reactions, proton transfer, metal coordination changes, electron transfer, and changes in protonation state.

Material to study

Molecular Representation

Atoms

In classical MD, each atom is represented using a set of parameters.

These commonly include:

  • atomic mass
  • Cartesian coordinates
  • velocity
  • partial electric charge
  • Lennard-Jones parameters
  • atom type

An atom in a force field should therefore not be interpreted as a literal hard sphere containing an explicit nucleus and electrons.

Instead, the parameters assigned to each atom are designed so that the collection of atoms reproduces important molecular properties observed experimentally or obtained from quantum mechanical calculations.

The apparent size of an atom in an MD simulation arises mainly from its non-bonded interaction parameters rather than from a rigid physical boundary.

Molecular Topology

A biomolecule is represented as a set of atoms together with information describing how those atoms are chemically connected.

This information is called the molecular topology.

The topology defines properties such as:

  • which atoms are bonded
  • bond types
  • bond angles
  • dihedral relationships
  • atom types
  • partial charges
  • parameters describing the interactions

The topology usually remains fixed during a conventional MD simulation.

The Force Field

A force field is a collection of mathematical functions and parameters used to calculate the potential energy of a molecular system.

Examples of commonly used biomolecular force fields include:

  • CHARMM
  • AMBER
  • OPLS
  • GROMOS

Different force fields use somewhat different functional forms and parameterization strategies, but they generally describe similar classes of molecular interactions.

The potential energy of the system can broadly be divided into:

  1. bonded interactions
  2. non-bonded interactions

Bonded Interactions

Atoms that are chemically connected require potential-energy terms that maintain realistic molecular geometry.

Bond stretching

A chemical bond is commonly approximated as a harmonic oscillator.

A simplified bond potential can be written as:

Vbond=12kb(r−r0)2V_{\mathrm{bond}} = \frac{1}{2}k_b(r-r_0)^2

where:

  • (r) is the current bond length,
  • (r_0) is the equilibrium bond length,
  • (k_b) is the bond force constant.

The energy therefore increases when a bond is stretched or compressed away from its preferred length.

Angle Bending

Three connected atoms define a bond angle.

Deviation from the preferred angle is also commonly represented using a harmonic potential:

Vangle=12kθ(θ−θ0)2V_{\mathrm{angle}} = \frac{1}{2}k_\theta(\theta-\theta_0)^2

This term helps maintain realistic local molecular geometry.

Dihedral Angles

Four sequentially connected atoms define a dihedral angle.

Rotation around chemical bonds is extremely important because it allows biomolecules to adopt different conformations.

Dihedral potentials describe the energetic preference for particular rotational states.

For example, rotation around many single bonds is possible but not equally favorable at all angles.

Double bonds, aromatic systems, and other chemically restricted structures are represented using parameters that strongly favor particular geometries.

Classical force fields do not normally describe explicit sigma and pi electrons. Instead, the underlying electronic properties of these chemical bonds are incorporated indirectly into the force-field parameters.

Improper Dihedrals

Improper dihedral terms are commonly used to maintain properties such as:

  • planarity
  • chirality
  • appropriate geometry around particular functional groups

CHARMM36 Bonded Terms

For example, the CHARMM36 force field includes several bonded interaction terms such as:

  • bond stretching
  • angle bending
  • Urey-Bradley interactions
  • proper dihedrals
  • improper dihedrals

Discussion: Why are bonds, angles, and dihedrals represented separately?

What structural property is controlled by each term, and why would bond lengths alone be insufficient to maintain a realistic molecular structure?

Non-Bonded Interactions

Atoms also interact with atoms to which they are not directly bonded.

The two major non-bonded interactions in conventional biomolecular MD are:

  • van der Waals interactions
  • electrostatic interactions

Lennard-Jones Potential

Van der Waals interactions are commonly modeled using the Lennard-Jones potential.

A typical form is:

VLJ(r)=4ϵ[(σr)12+(σr)6]V_{\mathrm{LJ}}(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} +\left(\frac{\sigma}{r}\right)^6 \right]

At intermediate distances, the attractive term approximates dispersion interactions.

At very short distances, the potential becomes strongly repulsive.

This repulsive term prevents atoms from unrealistically occupying the same region of space and approximates short-range quantum mechanical effects associated with overlapping electron clouds.

It should not be interpreted simply as an elastic collision between rigid spheres.

Electrostatic Interactions

Atoms in classical force fields are assigned partial electric charges.

The electrostatic interaction between two charged atoms is approximately described by Coulomb’s law:

VCoulomb(r)=14πϵ0qiqjrV_{\mathrm{Coulomb}}(r) = \frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r}

Electrostatic interactions are particularly important for:

  • salt bridges
  • hydrogen-bonding environments
  • protein-ligand recognition
  • protein-DNA interactions
  • ion interactions
  • solvent organization

Because electrostatic interactions are long-range, their calculation requires special treatment in periodic MD simulations.

Particle Mesh Ewald (PME) is commonly used for this purpose.

Discussion: How are long-range interactions calculated efficiently?

Topics to explore:

  • interaction cutoffs
  • neighbor lists
  • periodic boundaries
  • Ewald summation
  • Particle Mesh Ewald

Potential Energy, Force, and Motion

The force field defines the potential energy (V) of the molecular system.

The force acting on atom (i) is calculated from the gradient of this potential energy:

𝐅i=−∇iV\mathbf{F}_i = -\nabla_i V

The force is then related to acceleration through Newton’s second law:

𝐅i=mi𝐚i\mathbf{F}_i = m_i\mathbf{a}_i

Thus, the conceptual sequence of an MD simulation is:

atomic coordinates→potential energy→forces→acceleration→new positions and velocities\text{atomic coordinates} \rightarrow \text{potential energy} \rightarrow \text{forces} \rightarrow \text{acceleration} \rightarrow \text{new positions and velocities}

This procedure is repeated many times throughout the simulation.

Kinetic Energy and Temperature

Atoms possess velocity and therefore have kinetic energy.

For a single atom:

Ek=12mv2E_k = \frac{1}{2}mv^2

The total kinetic energy of the molecular system is related statistically to its temperature.

At finite temperature, atoms continuously move and exchange energy through their interactions with surrounding atoms.

The total energy of a system can broadly be expressed as:

Etotal=Ekinetic+EpotentialE_{\mathrm{total}} =E_{\mathrm{kinetic}} + E_{\mathrm{potential}}

In an isolated idealized system, total energy should be conserved.

However, many biomolecular simulations intentionally exchange energy or volume with external algorithms in order to maintain a desired temperature or pressure.

Discussion: Potential energy, kinetic energy, and total energy

How does temperature relate to kinetic energy? What should happen to total energy in an NVE simulation compared with an NVT simulation?

Numerical Integration and Time Steps

The equations of motion cannot generally be solved analytically for a large molecular system.

Instead, MD programs calculate the trajectory using numerical integration.

Common integration algorithms include:

  • Verlet
  • velocity Verlet
  • leap-frog

The system is advanced in very small discrete time intervals called time steps.

For atomistic biomolecular simulations, a typical time step is approximately:

1−2 fs1-2\ \mathrm{fs}

where:

1 fs=10−15 s1\ \mathrm{fs}=10^{-15}\ \mathrm{s}

The time step must be short enough to accurately describe the fastest important molecular motions.

Bond vibrations involving hydrogen atoms are particularly fast. Constraints such as LINCS can be used to restrict selected bond vibrations, allowing somewhat longer time steps.

Discussion: Why are MD time steps so short?

Explore the relationship between molecular vibration frequencies, hydrogen-containing bonds, numerical stability, and constraint algorithms such as LINCS.

Simulation Environment

Solvent

Biomolecules normally function in an aqueous environment.

A protein is therefore typically placed inside a simulation box containing explicit water molecules.

Common water models include:

  • TIP3P
  • TIP4P
  • SPC
  • SPC/E

The choice of water model is usually made together with the selected force field.

Ions may also be added to:

  • neutralize the total charge of the system
  • approximate physiological ionic strength

For example, a simulation may contain approximately 0.15 M NaCl.

Periodic Boundary Conditions

A simulation can contain only a finite number of water molecules.

If the edge of the simulation box represented a physical wall, molecules near that boundary would experience an artificial environment.

This problem is commonly avoided using periodic boundary conditions.

The simulation box is conceptually repeated in all directions.

When a molecule leaves one side of the box, an equivalent image appears from the opposite side.

This creates an approximately bulk-like environment without requiring an infinitely large number of solvent molecules.

Discussion: Periodic boundary conditions

What does it mean for a simulation box to be periodically repeated? What artifacts can occur if the box is too small?

Temperature Control

Most biomolecular simulations are performed at a specified temperature, such as:

  • approximately 298 K for room temperature
  • approximately 310 K for physiological temperature

A thermostat is used to generate or maintain the desired temperature distribution.

A thermostat should not simply be regarded as an external environment that manually adds energy whenever the system becomes cold.

Rather, it is a mathematical algorithm designed to control the statistical distribution of kinetic energy.

Common thermostats include:

  • velocity-rescaling thermostat
  • Nosé-Hoover thermostat

Discussion: What does a thermostat really do?

A thermostat does not necessarily force the instantaneous temperature to equal the reference temperature. Temperature naturally fluctuates. The important question is whether the thermostat generates the appropriate statistical ensemble.

Pressure Control

Biomolecules are commonly studied at approximately atmospheric pressure.

A barostat is used to control pressure.

A barostat typically changes the volume and dimensions of the simulation box so that the system samples states consistent with the specified pressure.

Examples include:

  • Berendsen barostat
  • C-rescale barostat
  • Parrinello-Rahman barostat

Different barostats have different properties and are not necessarily appropriate for every phase of a simulation.

Statistical Ensembles

MD simulations can be performed under different thermodynamic conditions.

NVE Ensemble

The NVE ensemble maintains:

  • N: number of particles
  • V: volume
  • E: total energy

No thermostat or barostat is used in an ideal NVE simulation.

NVT Ensemble

The NVT ensemble maintains:

  • N: number of particles
  • V: volume
  • T: temperature

A thermostat is used to control temperature while the box volume remains fixed.

NPT Ensemble

The NPT ensemble maintains:

  • N: number of particles
  • P: pressure
  • T: temperature

A thermostat and barostat are used.

For biomolecules in solution, NPT conditions are commonly used because they allow the solvent density and simulation-box volume to adjust at the desired temperature and pressure.

Typical MD Workflow

A conventional protein MD simulation commonly includes the following stages.

1. Structure Preparation

The starting molecular structure is prepared.

This may include:

  • removing or retaining crystallographic molecules
  • adding missing atoms
  • resolving alternate conformations
  • assigning protonation states
  • adding hydrogen atoms
  • parameterizing ligands
  • checking unusual residues

The quality of the starting structure strongly affects the simulation.

2. Solvation

The biomolecule is placed in a simulation box and surrounded by water molecules.

3. Ion Addition

Counterions are added to neutralize the system.

Additional ions may be introduced to reproduce an approximate experimental ionic strength.

4. Energy Minimization

The initial molecular system may contain unfavorable contacts or steric clashes.

Energy minimization adjusts atomic coordinates to reduce these severe interactions before dynamic simulation begins.

Energy minimization is not itself an MD simulation because it does not represent physical motion over time.

5. NVT Equilibration

The system is gradually brought to the desired temperature while maintaining a fixed simulation-box volume.

Position restraints may initially be applied to the biomolecule while the solvent reorganizes around it.

6. NPT Equilibration

Pressure control is introduced.

The simulation box and solvent density can then adjust toward conditions appropriate for the target temperature and pressure.

7. Production MD

After the system has been sufficiently prepared and equilibrated, a production simulation is performed.

The resulting trajectory is used for subsequent scientific analysis.

A common simplified workflow is therefore:

\text{Preparation} \rightarrow \text{Energy minimization} \rightarrow \text{NVT} \rightarrow \text{NPT} \rightarrow \text{Production MD}

Trajectories and Molecular Ensembles

During a simulation, atomic coordinates are periodically saved.

The resulting collection of structures is called a trajectory.

A trajectory can be visualized as a molecular movie, but this interpretation is incomplete.

At finite temperature, a biomolecule does not possess only one structure. Instead, it occupies an ensemble of conformational states.

Some states may be highly populated, whereas others may occur only rarely.

MD simulation attempts to sample these states over time.

Consequently, many biologically meaningful quantities are better described as distributions or ensemble averages rather than as properties of one particular structure.

Examples include:

  • average distances
  • probability of a hydrogen bond
  • distribution of dihedral angles
  • frequency of loop opening
  • conformational populations
  • residue flexibility

Discussion: Trajectory versus ensemble

A trajectory is one realization of molecular motion over time. Under what assumptions can observations from a finite trajectory be interpreted as properties of the molecular ensemble?

Sampling and Timescale

One of the major limitations of MD simulation is sampling.

Molecular processes occur over a wide range of timescales.

Examples may include:

  • bond vibrations: femtoseconds
  • side-chain rotations: picoseconds to nanoseconds
  • loop rearrangements: nanoseconds to microseconds or longer
  • domain motion: microseconds to milliseconds
  • protein folding: microseconds to seconds or longer

A simulation can only observe processes that occur within the accessible simulation time and are sampled sufficiently frequently.

Therefore, the absence of a conformational transition in a simulation does not demonstrate that the transition is biologically impossible.

Similarly, observing a transition only once may not be sufficient to estimate its probability.

Multiple independent simulations, longer trajectories, or enhanced-sampling methods may therefore be required.

Discussion: Simulation time versus biological time

Does a 100-ns simulation adequately sample a protein’s conformational behavior? The answer depends on the molecular process being investigated.

Assessing a Simulation

Before interpreting biological mechanisms from an MD trajectory, the quality and behavior of the simulation should be examined.

Several different concepts should be distinguished:

  • numerical stability
  • thermodynamic stability
  • structural equilibration
  • sampling
  • convergence

These are related but are not equivalent.

Temperature, Pressure, and Energy

Basic simulation properties should first be inspected.

These may include:

  • temperature
  • pressure
  • density
  • potential energy
  • kinetic energy
  • total energy
  • simulation-box volume

These measurements help identify major problems such as numerical instability or inappropriate equilibration.

However, stable temperature and pressure do not necessarily mean that the biomolecule itself has fully equilibrated.

Root Mean Square Deviation

One of the most widely used structural measurements is the root mean square deviation (RMSD).

RMSD measures the structural difference between two sets of atomic coordinates, generally after appropriate structural alignment.

For example, the protein structure at each simulation time point can be compared with the initial structure.

RMSD can indicate whether the protein:

  • remains close to its initial conformation
  • relaxes into a different structural state
  • undergoes large conformational transitions
  • shows possible structural instability

However, RMSD must be interpreted carefully.

A plateau in RMSD does not prove that the simulation has reached thermodynamic equilibrium or fully converged.

Likewise, a large RMSD does not necessarily indicate a failed simulation.

A large change may represent:

  • normal relaxation from the experimental starting structure
  • loop rearrangement
  • domain motion
  • ligand-induced conformational change
  • transition between functional states
  • unfolding
  • simulation artifact

The biological context and structural details should therefore be examined.

Discussion: Does a stable RMSD mean convergence?

Usually not. A stable RMSD shows only that the overall deviation from the selected reference remains within a certain range.

Stability, Equilibration, and Convergence

These concepts should be distinguished.

Stability

A simulation is stable when it does not show obvious numerical or structural failure.

Examples of instability include:

  • rapidly increasing energy
  • unrealistic temperature
  • severe structural distortion
  • numerical integration failure

Equilibration

A system can be considered equilibrated with respect to a particular property when that property no longer shows a strong systematic dependence on the initial preparation.

Different components may equilibrate on different timescales.

For example:

  • temperature may equilibrate rapidly
  • solvent density may equilibrate relatively rapidly
  • side-chain conformations may take longer
  • large protein domains may require much longer sampling

Thus, there is no single universal equilibration time.

Convergence

Convergence is a stronger concept.

A measured property can be considered reasonably converged when additional simulation does not substantially change its estimated value or distribution.

Examples include convergence of:

  • conformational populations
  • average distances
  • free-energy profiles
  • principal components
  • contact probabilities

Convergence is often difficult to establish from a single trajectory.

Independent replicate simulations can provide important evidence about reproducibility and sampling.

Structural and Dynamic Analyses

Once the simulation has been evaluated sufficiently, the trajectory can be analyzed in greater detail.

The appropriate analysis depends on the biological question.

RMSD

Root mean square deviation

Measures overall structural deviation relative to a reference structure.

Useful for:

  • monitoring global conformational changes
  • comparing different trajectories
  • identifying large structural transitions

RMSF

Root mean square fluctuation

Measures how strongly individual atoms or residues fluctuate around their average positions.

Useful for identifying:

  • flexible loops
  • rigid structural regions
  • ligand-induced changes in flexibility

Radius of Gyration

Measures the overall spatial compactness of a molecule.

It can help identify:

  • expansion
  • compaction
  • unfolding-related changes

Distance and Angle Analysis

Distances or angles between selected structural features can be monitored directly.

Examples include:

  • distance between two domains
  • catalytic residue distance
  • ligand-protein distance
  • opening angle of a protein domain
  • loop opening and closing

These measurements are often more biologically interpretable than global RMSD.

Hydrogen-Bond Analysis

Hydrogen bonds can be monitored over time.

Questions may include:

  • How frequently is a hydrogen bond formed?
  • Is an interaction stable throughout the trajectory?
  • Does ligand binding introduce new hydrogen bonds?

Contact Analysis

Residue-residue or protein-ligand contacts can be monitored.

A contact map can reveal how interactions change across different conformational states.

Secondary Structure

The formation and loss of:

  • alpha helices
  • beta sheets
  • turns

can be monitored throughout the trajectory.

Solvent Accessible Surface Area

Solvent-accessible surface area (SASA) can be used to assess:

  • protein exposure to solvent
  • burial of hydrophobic regions
  • accessibility of ligand-binding regions

Clustering

Trajectory structures can be grouped according to structural similarity.

Clustering can help identify representative conformational states rather than examining thousands of individual frames.

Essential Dynamics

Proteins often contain many atomic fluctuations, but only a smaller number of collective motions may dominate their behavior.

Principal component analysis (PCA), also called essential dynamics in the context of molecular simulation, can be used to identify these dominant motions.

PCA can reveal patterns such as:

  • domain opening and closing
  • collective loop movements
  • large-scale twisting
  • hinge motion

These motions may be related to biological function.

Dynamic Correlation

Different regions of a protein may move together or in opposite directions.

Dynamic cross-correlation analysis can be used to examine relationships between these motions.

It may help identify:

  • mechanically coupled residues
  • communication between distant domains
  • residues associated with allosteric regulation

However, correlated motion alone does not prove direct causality or molecular communication.

Discussion: Correlation versus causation

If two residues move together, does one residue control the other? Additional structural and experimental evidence is generally required.

From Motion to Biological Mechanism

The objective of MD analysis should not simply be to calculate as many structural measurements as possible.

Instead, analyses should be selected according to a biological hypothesis.

For example, suppose a protein is believed to become active through opening of an active-site loop.

A possible analysis strategy could be:

  1. Measure the distance or angle describing loop opening.
  2. Determine the distribution of open and closed conformations.
  3. Identify residues whose interactions change during opening.
  4. Use PCA to determine whether loop opening represents a dominant collective motion.
  5. Use correlation analysis to identify residues whose motions are associated with the loop.
  6. Compare wild-type and mutant simulations.
  7. Test whether the computational prediction is consistent with experimental observations.

In this way, MD becomes a tool for testing mechanistic hypotheses rather than simply generating molecular animations.

Limitations of Classical MD

Although molecular dynamics is powerful, several important limitations should always be considered.

Force-Field Approximation

The accuracy of the simulation depends on the quality of the force field.

Force fields are approximations of molecular interactions and may not reproduce every physical property equally well.

Fixed Chemical Topology

Conventional MD generally does not allow bonds to form or break.

Chemical reactions therefore require other approaches.

Electronic Polarization

Many commonly used biomolecular force fields use fixed atomic partial charges.

Real electron distributions respond dynamically to their environment.

This polarization is only approximated indirectly in standard fixed-charge force fields.

Sampling Limitations

The simulation may be too short to sample important conformational transitions.

A trajectory that appears stable may therefore represent only one region of a much larger conformational landscape.

Initial Structure Dependence

Simulation results may depend on the starting conformation.

Independent replicas can help determine whether observations are reproducible.

Model Dependence

Results can be affected by choices such as:

  • force field
  • water model
  • ion parameters
  • protonation states
  • simulation-box dimensions
  • thermostat
  • barostat
  • simulation length

These choices should therefore be documented and justified.

Key Conceptual Summary

A molecular dynamics simulation can be summarized conceptually as:

\boxed{ \text{Molecular structure} \rightarrow \text{Force field} \rightarrow \text{Potential energy} \rightarrow \text{Forces} \rightarrow \text{Atomic motion} \rightarrow \text{Trajectory} \rightarrow \text{Conformational ensemble} \rightarrow \text{Biological interpretation} }

The force field describes how the molecular system stores potential energy.

The gradient of this energy determines the forces acting on atoms.

Newtonian mechanics determines how those forces change atomic positions and velocities.

Repeated calculations generate a molecular trajectory.

The trajectory can then be analyzed to understand the conformational states and motions of the molecule.

The major scientific challenge is not simply generating a trajectory, but determining whether the simulation has sampled the molecular behavior sufficiently to support a reliable biological interpretation.

Topics for Further Study

Resources

MD Tutorials

  • http://www.mdtutorials.com/
  • GROMACS tutorials developed by Justin A. Lemkul
  • A useful practical introduction to preparing, running, and analyzing molecular dynamics simulations with GROMACS.

GROMACS Documentation

The official GROMACS documentation is particularly useful when connecting these theoretical concepts to practical simulation parameters, including:

  • integrators
  • constraints
  • thermostats
  • barostats
  • electrostatics
  • cutoffs
  • trajectory output
  • analysis methods

ChatGPT

This document was originally written by Seungpyo Hong and subsequently substantially edited with the assistance of ChatGPT.

Views: 27

Leave a Comment

Your email address will not be published. Required fields are marked *

Scroll to Top