Connecting the role of Protein-Protein Interactions to relevant disease pathways can assist in drug-design and discovery studies. In this blog, we have explored the binding dynamics of Protein-Protein Interactions for a reported complex by employing Protein-Protein Global FFT dock module of ICM-MolSoft and further analyzing the structural dynamics with OpenMM-based MD Simulation module of ICM-MolSoft.
What is Protein-Protein Interaction?
Protein-protein Interactions (PPIs) are interactions between two or more proteins bind non-covalently and specific to carry out niche biological functions. PPIs are so common that there exist about 130,000–650,000 PPIs in the human body. An entire set of PPIs within a human is called as ‘Interactome’ for that human. These interactions are quintessential for adequate functioning of multiple biological processes involving transcription, replication, signal transduction, eliciting immunological responses, as well as maintaining the cellular assembly. Disruptions of these PPIs often results in abnormal cell behaviour, leading to disease phenotypes. The PPIs are functionally involved in a variety of biological functions like modification of kinetic properties of the enzymes, substrate channelling, inactivation or suppression of protein expression, formation of binding sites for small biological molecules, change substrate specificity for a particular protein structure, and regulation of upstream or downstream processes. Given, their intricate functional involvement in many critical pathways, they represent a key target for therapeutic intervention in any disease condition. Inhibition of such critical pathways can completely alter the disease mechanisms. Electrostatic surface representations are important for understanding the molecular basis of protein–protein interactions because they reveal the distribution of charged regions across the protein surface. Positively and negatively charged patches help guide proteins toward each other through long-range electrostatic attraction, a process often referred to as electrostatic steering. The red regions typically indicate negatively charged surfaces (electron-rich; acidic residues such as Asp and Glu), while blue regions indicate positively charged surfaces (basic residues such as Lys and Arg). White or faintly coloured areas generally correspond to neutral or hydrophobic regions (Figure 1).
Figure 1: Electrostatic surfaces for ACE2 (red) and Spike Protein (blue)
What are the challenges while studying Protein-Protein interactions In Silico?
The Computational drug-design and discovery encircling PPIs are often hindered by following factors:
- Very high dimensionality of the conformational space
- Inherent protein flexibilities
- Identifying probable ‘hot spot’ amino acid residues
How does ICM-MolSoft overcome these challenges?
ICM-MolSoft employs two key technologies in its framework:
- ICM: Internal Coordinate Mechanics – ICM-MolSoft computes internal coordinates (ICs) instead of Cartesian coordinates for any molecular modelling operation. ICs are composed of bond lengths (b), bond angles (ω), torsion angles (ϕ), and phase dihedral angles (φ). This accounts for inherent protein flexibilities.
- BPMC: Biased Probability Monte Carlo – ICM-MolSoft employs BPMC algorithm which allows efficient sampling of the conformational space, making it explore complete set of probable conformations with high accuracy, reduced time-length and lesser computational expense.
- ODAs: Optimal Docking Areas – Additionally, ICM-MolSoft also integrates an algorithm called ‘ODA’, Optimal Docking Areas, which predicts the ‘hot spot’ residues, allowing the user to make data-guided decisions for protein-protein docking rather than scanning through the entire 3D structures for both the proteins.
- Enhanced Sampling Precision – 4 different options to choose from Quick Test Mode, Coarse, Regular, and Fine in the increasing power of sampling precision.
- Scoring and Energy parameters – Score calculated as initial score resulting from grid-based energy scoring function, and Total Energy parameter composed of Electrostatic energy, non-bonded interatomic pairwise interactions van der Waals energy, and surface energy term allowing more accurate binding pose selection.
- Integrated powerful processing capacity – Ease of submitting a job with access to 1, 2, 4 and 8 processors enabling faster process completion.
ICM-MolSoft Methodology for Protein-Protein Docking:
The Protein-Protein docking is a rigid-docking protocol which is carried out by Fast-Fourier Transform (FFT-based) algorithm. It comprises of 2 stages:
Stage 1: Translation and rotational-based conformational search
A) FFT-based translation search and initial scoring function
The first stage uses a simplified scoring function representing steric fit and hydrophobic/hydrophilic contact matching. FFT is then used for translational search using systematic search of rotations from 60 x 27 (coarse) to 256 x 125 (fine) orientations. The scoring function is a combination of two functions,
E(rt)=∫∫∫Φ(r)Ψ(r-rt) dr
where Φ(r) and Ψ(r) represent atomic coordinates of protein 1 and protein 2 respectively. These functions correspond to the van der Waals interaction reflecting shape complementarity and lipophilicity.
B) Systematic rotational search
Rotational degrees of freedom in a system of two interacting macromolecules are explored through a systematic search strategy. To reduce the computational cost associated with grid potential evaluations and forward transformations, rotations are applied to both interacting partners. One molecule (or sub-assembly) undergoes global rotation, sampling the full three-dimensional Euler angle space (φ, ψ, θ), whereas the second molecule is subjected to finer rotational adjustments within a limited region surrounding its initial orientation. The global rotational sampling employs a polyhedral framework to define the φ and ψ angles, combined with uniformly spaced values across the 0–360° range for θ. In contrast, the local rotational refinement follows a structured three-dimensional grid defined by increments iδ, jδ, kδ, where the indices i, j, and k vary from –n to n, and δ represents the angular step size. Although this method introduces some non-uniformity in sampling density, it significantly reduces the number of potential grids that require computation and transformation. For example, at a representative sampling resolution, 256 grid sets are generated for one molecule and 125 for the other, enabling the evaluation of a total of 32,000 rotational configurations.
The algorithm has 3 options to run FFT docking protocol- coarse, medium and fine search. At the coarse sampling level, global rotational exploration employs an icosahedral scheme, yielding 12 distinct (φ, ψ) orientation pairs, combined with a θ increment of 72°. Concurrently, local rotational refinement of the second molecule is performed using a 3 × 3 × 3 angular grid with a step size δ = 30°. At the intermediate (medium) resolution, global sampling is based on a dodecahedral distribution providing 20 (φ, ψ) pairs, with the θ angle sampled at 60° intervals. The corresponding local refinement uses a 5 × 5 × 5 grid with δ = 15°. While, the highest (fine) resolution, global rotations follow a hybrid dodecahedral–icosahedral scheme, generating 32 (φ, ψ) orientation pairs, while θ is incremented by 45°. Local rotational sampling remains on a 5 × 5 × 5 grid but with a finer angular step of δ = 9°, enabling more detailed exploration of nearby conformational space.
Stage 2: Scoring refinement for the docked poses
The top few thousand docking poses, typically between 3,000 and 20,000 are re-evaluated with an optimized Solvent Accessible Surface Area (SASA)-based solvation term which is known to improve ranking of the near-native docking solutions generated by Monte-Carlo sampling. It gives a score comprising of van der Waals, electrostatics, and solvation terms. Score is calculated as follows:
Score = E vw, hb + w el . E el + w sf . Score sf
where E vw, hb is the combined van der Waals and hydrogen bonding terms calculated according to the ECEPP/328 force-field, Eel is the electrostatic term, and Scoresf (surface) is the solvation term.
The technology behind Optimal Docking Areas (ODAs):
Optimal Docking Areas (ODAs) represents a computationally efficient and accurate approach for mapping protein surfaces to identify regions that exhibit favorable energetic contributions upon burial during protein–protein association. The method detects contiguous surface patches characterized by minimal desolvation energy, calculated using atom-based solvation parameters specifically calibrated for protein–protein docking applications. Validation of the approach is performed on unbound structures of 66 structurally diverse, non-homologous proteins participating in non-obligate heteromeric complexes with experimentally resolved architectures. In approximately half of the cases, distinct surface regions exhibiting significantly reduced docking-associated desolvation energy are identified. Notably, the ODA-derived “hot spots” correspond to experimentally observed binding interfaces in nearly 80% of the proteins when mapped onto their X-ray structures in the unbound state. These energetically favorable surface regions likely play a critical role during the early stages of complex formation. Their preferential burial upon association can stabilize transient encounter complexes by biasing interactions toward near-native interfacial arrangements, even in the presence of otherwise non-specific relative orientations.
Internal Case Study:
Here, we look at an example of Angiotensin-converting enzyme 2 (ACE2) in complex with SARS-Cov-2 spike protein with the PDB ID (6m0j). We use herein, the powerful Protein-Protein Docking module of ICM-MolSoft to accurately replicate the binding pose of both these proteins in complex.
The adopted Methodological Workflow in ICM-MolSoft:
Step 1: Load the PDB IDs 1r42 and 7l4z
Figure 2: ICM-MolSoft GUI with loaded PDB IDs 1r42
Figure 3: ICM-MolSoft GUI with loaded PDB IDs 7l4z
Step 2: Convert the XR crystal structures to an ICM objects
Right-click on the object, choose the ‘Convert PDB’ option and select the options mentioned below (Figure 4):
Figure 4: Pop-up window and relevant options using Convert PDB option
Step 3: Identify the ODAs for each protein to be docked
To identify the ODAs for each protein structure, Go to Tools à 3D Predict otein interface by ODA. In the resulting pop-up window, choose Atomic envelope and Create residue table options, and click on OK. This step was carried out individually for both the proteins under this study. (Figures 5 & 6).
Figure 5: Identification of ODAs for ACE2 protein for the PDB ID 1r42
Figure 6: Identification of ODAs for SARS-CoV-2 Spike RBD for the PDB ID 7l4z
The light and dark-red coloured spheres indicate high probability of binding with these residues. Whereas the blue coloured spheres bear extremely low probability.
Step 4: Setting up the docking protocol
Select the ligand residues in orange selection by clicking on the Spike Protein object à Select à By Res. Numbers à Arrange the AverageODA values in descending order to get most probable ‘hot spot’ residues à Enter the top 5 residue numbers in the pop-up window, click on OK. Click on the Orange-Green selection button to colour this selection as orange. Similarly, follow the same protocol to make green selection in the receptor object which is the ACE2 protein structure. (Figures 7 & 8).
Figure 7: Orange selections for the amino acid residues in the ligand object
Figure 8: Green selection for the amino acid residues in the receptor object
Set up the docking protocol by clicking on Docking à Protein-protein à Run Global FFT-dock. Select ACE2 protein as the receptor object and Spike Protein as the Ligand object. Choose ‘Fine’ method to allow the algorithm to explore a complete chemical space (Figure 9).
Figure 9: Pop-up window to run FFT Protein-Protein docking by selection of ODAs
Results for Protein-Protein Docking of ACE2 and SARS2-Spike Protein RBD:
The docking results identified the best binding pose of the Protein-protein complex by identifying 50 best docking hits. The first entry in the results table had a total energy of -68.82. Further, Refine FFT solutions option was selected and the best pose from the obtained results was loaded onto the Graphical window (Figure 10).
Figure 10: Protein-Protein docking table results employing the ODAs method
Furthermore, this obtained binding pose was subjected to MD simulations using the OpenMM module integrated with ICM-MolSoft to analyze the stability of the complex over a run time of 100 ns across 2000 trajectory frames (Figure 11).
Figure 11: Pop-up window to run Protein-Protein docked complex for MD Simulation (100 ns)
Analysis of Protein-protein refined complex by ODAs method after MD run for 100 ns:
A) RMSD graphs and their interpretations
After an initial sharp equilibration phase, all three curves stabilize, indicating that the system reaches a relatively stable conformational state. The ACE2 graph (red) gradually drifts downward, suggesting tightening or slight compaction over time. Whereas RMSD graph of the complex (green) shows good stability with minimal fluctuations, while the Spike RBD graph (blue) exhibits higher deviations, indicating greater flexibility or conformational rearrangements. Overall, the complex remains stable, but with asymmetric flexibility between the interacting partners (Figure 12).
Figure 12: Combined RMSD plots for ACE2, Spike Protein and Complex
B) RMSF graphs and their interpretation
The RMSF plot of the ACE2 protein highlights residue-wise flexibility throughout the 100 ns simulation. Most residues exhibit low fluctuations, indicating a largely stable conformation, while several sharp peaks correspond to highly flexible loop or terminal regions. These spikes suggest localized mobility rather than global instability, which is typical for surface-exposed segments. The relatively stable baseline implies that the core of ACE2 remains structurally rigid during binding. Overall, the pattern supports a stable protein with dynamic peripheral regions contributing to interaction adaptability. Similarly, the RMSF profile of the Spike protein RBD shows overall low-to-moderate fluctuations, indicating a relatively stable core during the simulation. The initial high peak corresponds to flexible terminal residues, which quickly stabilize along the sequence. Most residues maintain low RMSF values, suggesting that the binding interface remains structurally rigid upon interaction with ACE2. However, increased fluctuations toward the C-terminal region indicate flexible loop segments that may contribute to conformational adaptability. Overall, the RBD appears stable with localized flexibility at terminal and loop regions (Figure 13 & 14).
Figure 13: RMSF graph for ACE2 protein
Figure 14: RMSF graph for Spike Protein
C) Radius of Gyration graph and its interpretation
After a brief initial adjustment, the Rg stabilizes, suggesting that the complex quickly reaches an equilibrated state. A gradual increase in Rg over time points to slight expansion or loosening of the complex, possibly due to conformational rearrangements at the interface. However, the fluctuations remain within a narrow range, indicating no major unfolding events. Overall, the complex maintains structural integrity with minor flexibility-driven expansion (Figure 15).
Figure 15: Radius of Gyration graph for the complex
D) Interaction contacts between SARS2 Spike Protein and ACE2 protein
Persistent vertical bands indicate key residues—typically polar and charged ones such as Lys, Arg, Asp, and Glu—forming stable hydrogen bonds and salt bridges, which anchor the complex. Hydrophobic residues like Leu, Phe, and Tyr contribute through van der Waals and π–π interactions, visible as continuous but less intense contact patterns. Intermittent signals suggest transient hydrogen bonds and water-mediated interactions, reflecting dynamic rearrangements at the interface. Overall, the binding is stabilized by a combination of strong electrostatic interactions and supportive hydrophobic contacts involving hotspot residues in the RBD (Figure 16).
Figure 16: Amino acid residue interactions between SARS2 Spike Protein and ACE2 protein
E) Calculations for contact areas and its interpretation:
The contact area graph for the ACE2-Spike protein RBD complex displays a relatively stable binding interface over the 100 ns simulation, with values primarily fluctuating between 900 and 1100 Ų. A noticeable increase in contact area occurs between 40 and 55 ns, suggesting a period of enhanced structural packing or a conformational shift that momentarily strengthened the protein-protein interaction. Following this peak, the system equilibrates back to a steady state, indicating that the complex maintains high structural integrity and a consistent binding footprint. This sustained contact area is a strong indicator of a stable, favorable interaction between the viral RBD and the human receptor (Figure 17).
Figure 17: Contact area (sq. Ǻ) vs time (ns) between SARS2 Spike Protein and ACE2 protein
Conclusion:
The development of MolSoft ICM-Pro is centered on the work of Ruben Abagyan and Maxim Totrov. They changed molecular modeling by using Internal Coordinates (IC) instead of standard Cartesian coordinates. This approach significantly reduces complexity by focusing on torsion angles rather than the X, Y, Z positions of every atom.
- The 1994 Nature Breakthrough
In 1994, Abagyan and Totrov published a landmark paper in Nature Structural Biology. They successfully predicted a lysozyme-antibody complex with 1.6 Angstrom accuracy. This was achieved using Biased Probability Monte Carlo (BPMC) optimization, which remains the core engine of ICM software today.
- The 1996 Blind Docking Challenge
The method was validated in a 1996 “blind” challenge involving beta-lactamase. The researchers used a two-step process: first, a rigid-body docking search, followed by full side-chain refinement. This proved that ICM could accurately handle the flexibility of protein surfaces.
- Success in CAPRI Competitions
Since the early 2000s, ICM-Pro has been a top performer in CAPRI (Critical Assessment of Predicted Interactions), often referred to as the “Olympics” of protein docking.
- Key Technical Innovations
Internal Coordinates: Reducing the degrees of freedom to make sampling more efficient.
REBEL Method: A fast boundary element algorithm for calculating accurate electrostatic potentials.
Soft Docking: Using adjusted energy potentials to account for minor movements in the protein backbone during the initial docking stages.
ICM-DISCO: Introduced in 2003, this method allowed for global energy optimization while keeping sidechains fully flexible.
ODA (Optimal Docking Area): Added in 2005, this tool identifies surface patches with optimal desolvation energy to predict binding “hotspots.”
In this internal case study, ICM-MolSoft is efficiently able to predict the binding pose of the ACE2-SARS CoV-2 Spike protein complex within the RMSD cut-off range of 2 Ǻ with the Global FFT, Protein-Protein docking module by selecting the ODAs. Furthermore, MD simulations with OpenMM integrated with ICM-MolSoft enabled accurate determination of complex stability with various parameters like RMSD, RMSF and Radius of Gyration.
Recommended workflow:
- Load the PDB structure
- Predict the ODAs for both the protein structures
- Perform Protein-Protein docking by specifying the hotspot residues as selected ODAs
- Run Global FFT dock and Refine results obtained
- Select the best Protein-Protein docked complex with least total energy
- Run MD simulations with OpenMM
- Analyze the Protein-Protein complex with various parameters like RMSD, RMSF and Radius of gyration.
References:
- Akbarzadeh, S., Coşkun, Ö., & Günçer, B. (2024). Studying protein–protein interactions: Latest and most popular approaches. Journal of Structural Biology, 216(4), 108118. https://doi.org/10.1016/j.jsb.2024.108118
- Alam, A., Fernandez Cabezudo, M. J., Hamoudi, R., & Al-Ramadi, B. K. (2026). Targeting protein protein interactions and their modulators to enable new therapeutic strategies for human diseases. Npj Systems Biology and Applications. https://doi.org/10.1038/s41540-026-00698-z
- Nooren, I. M. A. (2003). NEW EMBO MEMBER’S REVIEW: Diversity of protein-protein interactions. The EMBO Journal, 22(14), 3486–92. https://doi.org/10.1093/emboj/cdg359
- Rao, V. S., Srinivas, K., Sujini, G. N., & Kumar, G. N. S. (2014). Protein-Protein Interaction Detection: Methods and Analysis. International Journal of Proteomics, 2014, 1–12. https://doi.org/10.1155/2014/147648
- Okereke, M., David, K. B., & Adedeji, O. J. (2023). Computational Approaches to Targeting Protein–Protein Interactions in Cancer: a Pathway to Drug Discovery. Future Drug Discovery, 5(1). https://doi.org/10.4155/fdd-2023-0003
- Arnautova, Y. A., Abagyan, R., & Totrov, M. (2018). Protein-RNA Docking Using ICM. Journal of Chemical Theory and Computation, 14(9), 4971–84. https://doi.org/10.1021/acs.jctc.8b00293
- Fernandez‐Recio, J., Totrov, M., Skorodumov, C., & Abagyan, R. (2005). Optimal docking area: A new method for predicting protein–protein interaction sites. Proteins: Structure, Function, and Bioinformatics, 58(1), 134–43. https://doi.org/10.1002/prot.20285