What do you mean by Molecular Dynamics (MD) Simulations?

What are the components of MD simulations?

The following components are involved in an MD simulation engine, and a comparison of the available algorithms and options are also mentioned in Table 1. 

1. Molecular system: The molecular system in MD simulations in ICM-Pro can be a protein-ligand, DNA-protein, DNA-ligand, DNA-protein-ligand, RNA-protein, RNA-ligand, RNA-protein-ligand, protein-covalently bound ligand, protein-protein, & protein-peptide complexes.   
2. Force-field: The force field is the mathematical model that describes how atoms interact with one another. It determines the forces acting on every atom at each simulation step. The total potential energy of the system is calculated as the sum of several interaction terms, Bond stretching, Angle bending, Dihedral (torsional) rotations, Van der Waals interactions, & Electrostatic interactions. These energy terms collectively define how the molecular system behaves during the simulation. The force-field employed in OpenMM MD simulations in ICM-Pro is AMBER ff14SB (2). 
3. Solvent model: The solvent model represents the surrounding environment, usually water. Explicit solvent models simulate individual water molecules and ions, providing realistic hydration effects, while implicit solvent models approximate solvent as a continuous medium to reduce computational cost. The solvent strongly influences protein folding, ligand binding, conformational dynamics, and electrostatic interactions. ICM-Pro employs an explicit solvent model with TIP3P type water model (3). 
4. Counter ions and salt: Atomic partial charges describe how electrical charge is distributed across a molecule. These charges are essential for calculating electrostatic interactions, hydrogen bonding, and salt bridges. ICM-Pro OpenMM employs nnAM1-BCC charge model which carries out addition of counterions and buffer ions like Na+ & Cl– (4). 
5. Simulation Box: A simulation box contains the molecular system, including the protein, ligand, solvent molecules, and ions. It also provides space for molecular motion during the simulation. It defines the boundaries of the simulation & enables the application of Periodic Boundary Conditions (PBCs), allowing the system to mimic an effectively infinite bulk environment. It determines the amount of solvent, computational cost, and accuracy of the simulation. Different types of simulation boxes are:
   1. Cubic: Simple geometry, Easy to construct, & contains more solvent than necessary for roughly spherical proteins.
   2. Rectangular (orthorhombic): Suitable for elongated systems & can reduce unnecessary solvent compared to a cube.
   3. Truncated octahedron: Nearly spherical & requires fewer water molecules. More computationally efficient for globular proteins.
   4. Rhombic dodecahedron: Very efficient for approximately spherical proteins & minimizes solvent volume.

   The default simulation box type in ICM-Pro OpenMM is rectangular.

6. Integrator: The integrator numerically solves Newton’s equations of motion. Using the forces calculated from the force field, it updates the positions and velocities of all atoms at every time step, thereby generating the molecular trajectory over time. The integrator used in ICM-Pro OpenMM is Verlet integrator, with a default timestep of 2 fs (5).
7.  Thermostat & Barostat: The thermostat regulates the system temperature by adjusting atomic velocities. It maintains the desired thermal conditions, compensates for numerical energy drift, and ensures sampling of the appropriate statistical ensemble. Proper temperature control is essential for obtaining physically meaningful simulations. The barostat regulates system pressure by adjusting the simulation box dimensions and volume. It enables simulations under constant-pressure conditions (NPT ensemble) (6), allowing the system density and volume to equilibrate to experimentally relevant conditions. Andersen thermostat and Monte Carlo barostat are employed in ICM-Pro OpenMM as default options (7,8).  
8. Ensembles: The ensemble defines which thermodynamic variables remain constant throughout the simulation. Common ensembles include NVE (constant energy), NVT (constant temperature), and NPT (constant pressure). The selected ensemble determines the physical conditions being modelled and the statistical properties sampled during the simulation ICM-Pro OpenMM employs the NPT ensemble for MD simulations.

Table 1: Various components of MD simulations along with options available in OpenMM and in ICM-Pro with OpenMM

MD Simulation in ICM Pro

The steps involved in MD Simulations:

MD simulation in ICM-Pro follows all the steps typically involved in an MD run including (Figure 1):

1. Solvation: Addition of TIP3P water molecules. Smooth 10.0 Angstrom cutoff for non-bonded interactions is used with periodic boundary.
2. Addition of ions: Addition of counterions and buffer ions 
3. Energy minimization: Minimization of complex to reduce steric clashes and optimization of geometry. 
4. Equilibration: Equilibration of the complex with respect to atmospheric pressure and physiological temperature.
5. Production: Integrating force equations for all the atoms under study by assigning them flexibilities (9).

Figure 1: Steps in a Molecular Dynamics Protocol in ICM-Pro

How to setup an MD run and edit MD Simulation parameters?

Setting up an MD simulation through the GUI runs a background MD script (Figure 2, Figure 3).

The same script can also be run from the terminal, and the parameters can be edited as per user needs. The following parameters can be edited:

• -ng: don’t use GPU even if detected
• -w: add water
• time=: simulation length in nanoseconds
• confs=: number of snapshots to store in the stack
• rep=: run multiple (number of repeats) trajectories from the provided starting point
• -m: apply icm minimization/relaxation before the MD run
• -rl: restrain parts of the system far from ligand (default is heteroatom molecule, lig=: can be used to select ligand as desired)
• -tl: like -rl but tether to initial positions instead of internal distance restraints
• -E: embed into a pre-built POPC/POPE lipid membrane and run in membrane mode
• center=: in conjunction with -E, center membrane on selected residue
• -M: membrane mode, assuming pre-embedded input. Box must be consistent with periodic membrane
• -B: print A LOT of debug info
• BOX=: define box by two corners (min/max xyz). A custom box can be used but in MD the box must cover the entire object, with some margin. Box is used to create periodic system, where everything repeats in adjacent boxes. In general, there is rarely a good reason to not use automatic default. If one does use custom box, it is to make it bigger, not smaller.
• output=: name of the output object with an embedded stack of snapshots
• ligfile=: run a series of MD simulations, combining each ligand from filename.icb or filename.sdf with target protein
• gpuid=: GPU device index to run simulation on. Default 0

One of the major advantages of performing molecular dynamics (MD) simulations in ICM-Pro is its integration with the OpenMM simulation engine, which is specifically designed to exploit the massive parallel processing capabilities of modern Graphics Processing Units (GPUs).

Unlike traditional CPU-based MD engines, where force calculations are distributed across a relatively small number of processing cores, GPUs contain thousands of lightweight cores capable of performing millions of mathematical operations simultaneously. Since the most computationally intensive part of an MD simulation is calculating bonded and non-bonded interactions between thousands to millions of atoms at every femtosecond time step, these calculations are highly amenable to parallelization.

OpenMM automatically distributes these force calculations across the GPU, dramatically reducing computation time while maintaining the same level of numerical accuracy as CPU implementations. This enables researchers to simulate significantly longer timescales or larger biomolecular systems within practical timeframes.

Advantages of GPU Acceleration

• Substantial performance gains – Depending on the hardware and system size, GPU acceleration can provide 10–100 times faster simulations than conventional CPU-only calculations. 
• Longer simulation times – Nanosecond, and even microsecond-scale MD simulations become far more accessible. 
• Larger biological systems – Protein complexes, membrane proteins, nucleic acids, and large solvated systems can be simulated efficiently. 
• Reduced computational cost – A single high-end GPU can often outperform a multi-core CPU workstation while consuming less power. 
• Interactive research workflow – Faster simulations to iterate more quickly through model refinement, docking validation, and trajectory analysis.

Exhaustive analysis parameters after running MD Simulations:

MolSoft’s ICM-Pro offers an exhaustive analysis after completion of an atomistic MD simulation run, allowing a holistic idea of the stability of the complex under study.

• Root-mean square fluctuations (RMSF):

RMSF is a residue- or atom-based measure that quantifies the average positional fluctuation of each atom or residue around its mean position during a molecular dynamics (MD) simulation. The RMSF of atom or residue i is calculated as:

• Radius of Gyration (Rg):

Rg is a global structural parameter used in molecular dynamics (MD) simulations to measure the compactness of a biomolecule. It describes how the atoms of a molecule are distributed around its center of mass (COM). The Radius of Gyration is defined as:

• Ligand contacts: 

Ligand contacts monitor the interactions between a ligand and its binding site residues throughout the simulation. It provides detailed information about which amino acid residues remain in contact with the ligand, how frequently these contacts occur, and the nature of the interactions. The types of ligand contacts include Hydrogen Bonds, Hydrophobic contacts, Salt-bridge interactions, ϖ-ϖ stacking, ϖ-cation and ϖ-anion interactions. Additionally, in presence of a covalently bound ligand, the formation of covalent bond throughout the trajectory frames can also be analyzed.    

Molecular Mechanics Energy (EMM) includes Van der Waals interactions, Electrostatic interactions, Bond stretching, Angle bending, & Dihedral energies. Generalized Born Energy (GB) represents the polar contribution to solvation. Surface Area Energy (SA) represents the nonpolar contribution to solvation. It accounts for Hydrophobic effects, Cavity formation, & Solvent exclusion.

• RTCNN score:

RTCNN score is a is a neural network-based rescoring function developed by MolSoft for protein–ligand binding affinity prediction. It is used to rescore snapshots generated from MD simulations. It evaluates each snapshot independently using a trained neural network. 

• Cluster Tree:

Cluster Tree Analysis is used to group similar molecular conformations into clusters based on their 3D coordinates. Cluster analysis identifies representative conformational states that the biomolecule adopts during the simulation. 

How to perform Stack calculation for MD parameters using the GUI in ICM Pro? 

 2. RMSD to reference:

To calculate the RMSD for an ICM object, firstly the reference object (docking object or the PDB structure or Modelled object) and MD object must be superimposed. To calculate this parameter, select the reference ICM object (For Ligand RMSD, select ligand in reference ICM object, for protein RMSD, select protein object in reference ICM object, for complex RMSD, select the entire complex in the reference ICM object) for calculation. Additionally, all these individual graphs can also be plotted together on a single plot (Figure 5). There are three options to calculate RMSD for an ICM object:

1. Static: Every atom in the current frame is compared directly to the corresponding atom in the reference structure without any structural alignment. It measures both internal conformational changes as well as translational and rotational changes.
2. Superimpose: An optimal least-squares alignment between the selected atoms of the trajectory frame and the selected atoms of the reference structure is performed before RMSD calculation. It removes the translation & rotation motions. It measures only conformational changes while ignoring the global movements in the complex.
3. Super Other: Instead of aligning the atoms whose RMSD is being measured, ICM aligns using a different selection (For eg: The protein backbone or a specific protein chain), then computes RMSD on the selection of interest.

Figure 5: A combined RMSD plot for Ligand, Protein, and Complex (MC) showing RMSD (in Ǻ) vs Time (in ns)

3. RMSF:

To calculate the RMSD for an ICM object, firstly the reference object (docking object or the PDB structure or Modelled object) and MD object must be superimposed. To calculate this parameter, the amino acid residues must be selected, or the entire protein residues will be used for calculation (Figure 6). Additionally, the RMSF can also be specifically calculated for residues within 5 Å vicinity from the bound ligand. There are two options to calculate RMSD for an ICM object:

     a) Superimpose: This option superimposes the structure onto the reference. After removing global translation and rotation, it calculates the fluctuation of the selected residues around its average position.

     b) Static: No alignment is performed. The coordinates are used exactly as stored in every trajectory frame. Therefore, RMSF contains contributions from internal flexibility, protein translation, & protein rotation. The output is a plot of amino acid residues vs RMSF (in Angstroms)

Figure 6: RMSF plot showing RMSF (in Ǻ) vs Time (in ns)

• PCA analysis:

PCA is a statistical technique used to identify the dominant collective motions of a biomolecular system during an MD simulation. Rather than analyzing the motion of individual atoms independently, PCA extracts the most significant, large-scale conformational changes from the trajectory. The PCA analysis can be run by pasting the following ICM script as it is all at once in the terminal:

l_commands = no
N = Nof( a_ conf )
s_at = “c,ca,n”   # backbone (can be just ‘ca’)
M = Nof( a_//$s_at )*3
MM = Matrix( N, M )
load conf a_ 1
xyz1 = Xyz( a_//$s_at )
MM[1,?] = Rarray(xyz1)

for i=2,N
  load conf a_ i
  xyz = Xyz( a_//$s_at )
  superimpose  xyz1 xyz  # align to first frame
  MM[i,?] = Rarray(xyz)
endfor

TT = Table(MM)
make pca TT
make plot TT “x=PC1;y=PC2;xHash=yes;color=PC3” 

Furthermore, the PC1 vs PC2 can be plotted to understand the dominant conformations in the generated MD stack. 

• Water-mediated hydrogen bonds:

The key Hydrogen bonding interactions mediated through the water molecules can also be analyzed throughout the entirety of the trajectory, by pasting the following ICM script all at once in the terminal:

myHbond = Iarray(0)
confunum = Iarray(0)
for i=1, 2002
load conf a_dock_pose1_MD. i
show hbond a_dock_pose1_MD.water a_dock_pose1_MD.a
myHbond  //= i_out
confnum //= i
endfor

add column t confnum myHbond

NOTE: The parts of the above script highlighted in orange should be edited before pasting. “2002” should be changed to the number of trajectory frames in the generated stack output. “dock_pose1_MD” should be changed to the name given by the user to their respective MD object.

Advantages of OpenMM integrated with ICM-Pro:

• Easy MD setup for both solution and membrane complexes.
• Offers exhaustive analysis for atomistic MD simulations.
• Additional parameters like VLS Binding Score and RTCNN score for each trajectory frame.
• Advanced parameters like DCCM, PCA, Water-mediated hydrogen-bonding analysis can also be performed.
• User-friendly visualization of trajectory, frame by frame analysis of interactions also enabled.
• Dynamic and attractive graphical options for MD parameter plots, with ability for multi-plot integration, plot labelling options and colouring styles.
• Another key feature is the option to utilize GPU for an MD run, thus saving tremendous computational time and cost.

References:

1. Comput Sci Eng. 2010 Feb 05;12(4):34–39. doi: 10.1109/MCSE.2010.27
2. J Chem Theory Comput. 2015 Jul 23;11(8):3696–3713. doi: 10.1021/acs.jctc.5b00255
3. J. Phys. Chem. A. 2001 Oct 06; 105 (43): 9954–60. doi: https://doi.org/10.1021/jp003020w
4. J. Comput. Chem. 2002 Oct 18;23(16):1623-41. doi: 10.1002/jcc.10128  
5. Entropy 2014 Sep 19; 16(1), 138-62; https://doi.org/10.3390/e16010138
6. Comput. Phys. Comm. 2006 Feb 15; 174(4):263-9 doi: https://doi.org/10.1016/j.cpc.2005.10.011
7. Comm. Pure Appl. Math., 61: 96-136. https://doi.org/10.1002/cpa.20198
8. J Comput Chem. 2021 Dec 18;43(6):431–4. doi: 10.1002/jcc.26798
9. Neuron. 2018 Sep 19;99(6):1129–43. doi: 10.1016/j.neuron.2018.08.011