Protein-Ligand Docking and Molecular Dynamics Simulations for Antiviral Drug Design Targeting SARS-CoV-2 Spike Protein

1. Introduction

The development of antiviral agents against emerging coronaviruses requires detailed atomic-level understanding of viral entry mechanisms. The spike glycoprotein of SARS-CoV-2 mediates host cell attachment and fusion by engaging the angiotensin-converting enzyme 2 (ACE2) receptor. This interaction is primarily driven by the receptor-binding domain (RBD) of the spike protein. Computational approaches, including protein-ligand docking and molecular dynamics (MD) simulations, have become essential tools for identifying small molecules and peptides that can disrupt this interaction [1, 2, 3]. These methods enable high-throughput virtual screening of compound libraries and provide mechanistic insights into binding kinetics, conformational changes, and resistance mutations.

In veterinary medicine, understanding spike protein dynamics across host species is critical for assessing zoonotic spillover risk and designing interventions for animal reservoirs. This article reviews the application of protein-ligand docking and MD simulations to the SARS-CoV-2 spike RBD, with emphasis on protocols, algorithms, and interpretation of results for antiviral drug design.

2. The SARS-CoV-2 Spike Receptor-Binding Domain as a Drug Target

The spike protein is a homotrimeric class I fusion glycoprotein. Each monomer contains an S1 subunit (including the RBD) and an S2 subunit responsible for membrane fusion. The RBD adopts two main conformations: a "closed" (down) state inaccessible to ACE2 and an "open" (up) state that exposes the receptor-binding motif. The RBD-ACE2 interface involves a large, relatively flat protein-protein interaction surface with key contact residues on both partners. Targeting this interface with small molecules or peptides is challenging due to the buried surface area and lack of deep binding pockets. However, several cryptic pockets and transient cavities have been identified through MD simulations, offering alternative druggable sites [1, 2, 3].

The RBD also exhibits high sequence variability across SARS-CoV-2 variants, particularly at positions 417, 484, 501, and others. These mutations can alter ACE2 binding affinity and immune evasion; computational prediction of such effects is a major focus of structural bioinformatics. This topic is explored further in the article Computational Prediction of Receptor-Binding Domain Mutations in SARS-CoV-2 Using Molecular Dynamics and Free Energy Perturbation.

3. Protein-Ligand Docking Methodologies

Protein-ligand docking predicts the preferred orientation of a small molecule or peptide when bound to a target protein to form a stable complex. The process involves two components: a sampling algorithm to generate candidate poses and a scoring function to rank them. Commonly used docking programs include AutoDock Vina, which employs a stochastic global optimization algorithm (iterated local search) and a semi-empirical scoring function combining steric, electrostatic, and hydrophobic terms. Detailed protocols for AutoDock Vina are provided in AutoDock Vina Receptor-Ligand Docking: Practical Protocols for Protein-Small Molecule Docking.

For peptide ligands, specialized docking approaches such as flexible peptide docking or ensemble docking are required due to the high conformational flexibility of peptides. Sakib et al. [1] conducted a computational screening of 645 antiviral peptides against the SARS-CoV-2 RBD using a combination of docking and MD simulations. They prepared the RBD structure (PDB ID 6M0J) and performed rigid receptor docking with flexible ligand sampling using AutoDock Vina. The top 10 peptides showed binding energies ranging from -9.5 to -10.8 kcal/mol, with key interactions involving hydrophobic contacts and hydrogen bonds with RBD residues such as Tyr453, Asn487, and Gln493 [1]. Figure 1 illustrates the general docking workflow.

Table 1. Representative docking scores (binding free energy in kcal/mol) from computational studies targeting the SARS-CoV-2 RBD.

Elmaaty et al. [2] performed docking of salvianolic acid B and other compounds against both the spike RBD and the main protease (3CLpro). Their docking protocol employed AutoDock 4.2 with Lamarckian genetic algorithm sampling. The study reported that salvianolic acid B formed multiple hydrogen bonds and hydrophobic interactions with RBD residues Arg403, Glu406, and Lys417, suggesting dual inhibition potential [2]. Muhseen et al. [3] evaluated 83 terpenes as attachment inhibitors to the human ACE2 receptor. They used AutoDock Vina to dock compounds against the RBD (PDB ID 6LZG) and identified several with binding affinities below -8.0 kcal/mol, including oleanolic acid and glycyrrhetinic acid. The terpenes occupied the RBD-ACE2 interface, competing directly with the receptor [3].

4. Molecular Dynamics Simulations for Conformational Sampling and Binding Affinity Refinement

Molecular dynamics simulations provide a physics-based description of protein and ligand motion in explicit solvent. For antiviral drug design targeting the spike RBD, MD simulations are used to assess the stability of docked complexes, analyze conformational changes, and calculate binding free energies with higher accuracy than docking alone.

4.1. Simulation Setup and Force Fields

Standard MD protocols involve solvating the protein-ligand complex in a water box (e.g., TIP3P water model), adding counterions to neutralize the system, and performing energy minimization followed by equilibration in the NVT and NPT ensembles. The AMBER or CHARMM force fields are commonly used for the protein; ligand parameters are generated using tools such as GAFF or CGenFF. Simulations are typically run for 50-500 ns using programs such as GROMACS or AMBER. Detailed guidance is available in GROMACS Molecular Dynamics: Setting Up, Simulating, and Analyzing Protein-Water Systems.

Sakib et al. [1] followed their docking with 50 ns MD simulations in GROMACS using the CHARMM36 force field. They analyzed root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), and hydrogen bond occupancy. The peptide-RBD complexes maintained stable RMSDs below 0.3 nm, and the peptides induced changes in RBD flexibility near the binding interface. Elmaaty et al. [2] also performed 50 ns MD simulations on the top docked complexes using GROMACS with the AMBER99SB-ILDN force field. Their analysis confirmed stable binding of salvianolic acid B to the RBD, with persistent hydrogen bonds to Arg403 and Asp405. Muhseen et al. [3] conducted 50 ns MD simulations of the top terpene-RBD complexes and calculated binding free energies using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method. The MM-PBSA approach combines molecular mechanics energies, solvation free energies, and entropy estimates to predict binding affinities. This methodology is discussed in Free Energy Perturbation Calculations in Drug Discovery.

4.2. Conformational Dynamics and Cryptic Binding Sites

MD simulations can reveal conformational changes in the spike RBD that are not captured by static docking. The RBD transitions between the open and closed states on microsecond timescales. Simulations of the RBD in complex with small molecules or peptides may stabilize or destabilize specific conformations. For instance, Sakib et al. [1] observed that certain antiviral peptides induced an opening of the RBD, potentially facilitating interactions with ACE2, a counterproductive effect that highlights the importance of evaluating functional consequences.

The article Molecular Dynamics Simulations of Viral Envelope Proteins for Drug Docking and Design provides a broader overview of these techniques.

5. Free Energy Perturbation and Binding Affinity Prediction

While docking scores correlate roughly with experimental binding affinities, more accurate predictions require alchemical free energy methods such as free energy perturbation (FEP) or thermodynamic integration. FEP calculations involve gradually mutating the ligand to a dummy molecule in both the bound and free states to compute the relative binding free energy. These methods require extensive sampling and careful setup but can achieve accuracy within 1 kcal/mol of experimental values.

For spike-RBD targeting, FEP has been applied to assess the impact of mutations on ACE2 binding and to optimize lead compounds. The article Computational Prediction of Receptor-Binding Domain Mutations in SARS-CoV-2 Using Molecular Dynamics and Free Energy Perturbation elaborates on this topic.

None of the three provided references explicitly used FEP; however, the MM-PBSA approach employed by Muhseen et al. [3] and the MM-GBSA approach (used by Elmaaty et al. [2] for the 3CLpro but not reported for spike) represent end-point free energy methods that are computationally less expensive but also less accurate than FEP. The calculated binding free energies from MD-based methods are generally more reliable than docking scores for ranking compounds.

6. Case Studies: Computational Screening Against Spike RBD

6.1. Peptide Library Screening (Sakib et al. 2021)

Sakib et al. [1] curated a library of 645 antiviral peptides from the Antimicrobial Peptide Database (APD3) and the PeptideRadar tool. They docked each peptide against the RBD and selected the top 10 based on docking score, favorable interactions, and absence of unfavorable contacts. Subsequent MD simulations confirmed stable binding for three peptides (APD-001, APD-027, and APD-058). These peptides formed networks of hydrogen bonds and salt bridges with RBD residues critical for ACE2 recognition, such as Lys417, Tyr453, and Asn487. The computational pipeline used by Sakib et al. is representative of a typical structure-based virtual screening workflow.

6.2. Natural Products: Salvianolic Acid B (Elmaaty et al. 2022)

Elmaaty et al. [2] evaluated salvianolic acid B, a polyphenolic compound from Salvia miltiorrhiza, as a dual inhibitor of the spike RBD and 3CLpro. Docking against the RBD (PDB ID 6M0J) revealed that salvianolic acid B occupies a groove near the ACE2 binding interface, interacting with Arg403 and Lys417 via hydrogen bonds and with Phe490 and Leu492 via hydrophobic contacts. MD simulations showed that the RBD-salvianolic acid B complex remained stable throughout 50 ns, with minimal conformational changes in the RBD loops. These findings support further experimental validation for this compound as a potential attachment inhibitor [2].

6.3. Terpenes as Attachment Inhibitors (Muhseen et al. 2020)

Muhseen et al. [3] investigated a set of 83 terpenes from dietary and medicinal plants. Docking identified six compounds with binding affinities below -7.5 kcal/mol, including oleanolic acid, ursolic acid, and glycyrrhetinic acid. MD simulations and MM-PBSA calculations confirmed that these terpenes maintained stable interactions with RBD residues Ser494, Tyr505, and Gln493. The authors proposed that these natural compounds could serve as lead scaffolds for developing spike-ACE2 interface inhibitors [3].

7. Integrating Computational and Experimental Data

Validation of computational predictions against experimental data is essential. Binding affinities obtained from surface plasmon resonance (SPR), biolayer interferometry, or fluorescence polarization can be compared with docking scores or MD-derived binding free energies. Additionally, structural data from X-ray crystallography or cryo-electron microscopy can validate predicted binding modes. The article Structural and Evolutionary Dynamics of Coronavirus Spike Protein: Integrating Cryo-EM, Molecular Dynamics, and Phylogenetic Surveillance discusses such integration.

For the studies reviewed, experimental validation has been reported separately: Sakib et al. [1] did not perform in vitro assays but provided rationale for future testing; Elmaaty et al. [2] referenced their own experimental work for 3CLpro but not for spike; Muhseen et al. [3] relied solely on computational predictions. This highlights a common gap in the literature where docking and MD predictions outpace experimental confirmation.

8. Workflow for Computational Antiviral Design Targeting Spike RBD

flowchart TD
    A[Target Selection: SARS-CoV-2 Spike RBD], > B[Structure Preparation<br>PDB download, protonation, energy minimization]
    B, > C[Ligand Library Preparation<br>Peptides / Small molecules / Natural products]
    C, > D[Protein-Ligand Docking<br>AutoDock Vina / Glide]
    D, > E[Pose Selection and Scoring<br>Top candidates by binding energy]
    E, > F[Molecular Dynamics Simulations<br>GROMACS / AMBER, 50-500 ns]
    F, > G{Stable Complex?}
    G, >|Yes| H[Binding Free Energy Calculation<br>MM-PBSA / MM-GBSA / FEP]
    G, >|No| D
    H, > I[Hit Prioritization<br>Rank by affinity, stability, synthetic accessibility]
    I, > J[Experimental Validation<br>SPR, ITC, pseudovirus neutralization]
    J, > K[Lead Optimization<br>Iterative docking and MD cycles]

This workflow integrates docking for initial screening, MD for stability assessment, and free energy calculations for quantitative ranking. It is applicable to both small-molecule and peptide inhibitors targeting the spike RBD.

9. Conclusion

Protein-ligand docking and molecular dynamics simulations provide a powerful computational framework for antiviral drug design targeting the SARS-CoV-2 spike RBD. Studies on peptide inhibitors [1], natural polyphenols [2], and terpenes [3] demonstrate the utility of these methods in identifying candidates that disrupt RBD-ACE2 interactions. Docking algorithms such as AutoDock Vina enable high-throughput virtual screening, while MD simulations refine binding poses and predict conformational stability. Free energy methods further improve affinity predictions.

In veterinary contexts, these computational approaches can be extended to spike proteins from animal coronaviruses, aiding in the design of species-specific inhibitors and in assessing cross-species transmission risks. Continued integration of computational predictions with experimental validation will accelerate the development of antiviral interventions for both human and animal health.

References

[1] Sakib MMH, Nishat AA, Islam MT, et al. Computational screening of 645 antiviral peptides against the receptor-binding domain of the spike protein in SARS-CoV-2. Comput Biol Med. 2021. URL: https://pubmed.ncbi.nlm.nih.gov/34403938/

[2] Elmaaty AA, Darwish KM, Khattab M, et al. In a search for potential drug candidates for combating COVID-19: computational study revealed salvianolic acid B as a potential therapeutic targeting 3CLpro and spike proteins. J Biomol Struct Dyn. 2022. URL: https://pubmed.ncbi.nlm.nih.gov/33928870/

[3] Muhseen ZT, Hameed AR, Al-Hasani HMH, et al. Promising terpenes as SARS-CoV-2 spike receptor-binding domain (RBD) attachment inhibitors to the human ACE2 receptor: Integrated computational approach. J Mol Liq. 2020. URL: https://pubmed.ncbi.nlm.nih.gov/33041407/ *** Disclaimer: This article is for educational and informational purposes only. It is not intended to substitute for professional veterinary advice, diagnosis, treatment, or regulatory guidance. Always consult a licensed veterinarian or qualified specialist regarding animal health, disease diagnosis, and therapeutic decisions.