Molecular Dynamics Simulations of Viral Envelope Proteins for Drug Docking and Design

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-05-10

Key Takeaways

• Molecular Dynamics (MD) simulations are critical for understanding the dynamic conformational ensembles of viral envelope proteins, which are essential for host cell entry via receptor recognition and membrane fusion. These simulations reveal transiently accessible binding sites and allosteric pockets that are often not apparent in static structural data.
• Accurate MD simulations rely on well-parameterized force fields (e.g., AMBER, CHARMM) and appropriate system setup, including explicit solvent models and lipid bilayer embedding for membrane-bound proteins, to capture realistic protein dynamics and interactions.
• Ensemble docking, where multiple protein conformations from MD trajectories are used for ligand screening, significantly improves the identification of potential drug candidates by accounting for protein flexibility, a limitation of traditional rigid-receptor docking.
• Advanced computational methods like Markov State Models (MSMs) and free energy calculations (MM/GBSA, FEP) are integrated with MD and docking to more accurately predict binding affinities and rank potential inhibitors, facilitating rational drug design against viral targets.
• MD simulations have been successfully applied to diverse viral families, including coronaviruses (e.g., SARS-CoV-2 spike protein), influenza viruses (hemagglutinin), flaviviruses, and filoviruses, to identify inhibitors targeting fusion mechanisms, receptor binding, or ion channels.
• Future directions involve enhancing sampling methods to capture slow protein motions, improving force field accuracy for glycosylated regions, and integrating machine learning for accelerated discovery of antibodies and small-molecule inhibitors, with increasing relevance for veterinary virology.

Introduction

Viral envelope proteins mediate the initial steps of host cell infection, including receptor recognition, attachment, and membrane fusion [<a href="#ref-1">1</a>, <a href="#ref-2">2</a>, <a href="#ref-3">3</a>]. These glycoproteins exist in multiple conformational states, often transitioning from a metastable pre-fusion form to a post-fusion form upon binding to host receptors or exposure to low pH [<a href="#ref-4">4</a>, <a href="#ref-5">5</a>, <a href="#ref-6">6</a>]. The dynamic nature of these proteins presents both opportunities and challenges for structure-based drug design. Molecular dynamics (MD) simulations have become an indispensable tool for capturing the conformational ensembles of envelope proteins, identifying cryptic binding pockets, and estimating binding free energies for small molecules, peptides, and antibodies [<a href="#ref-7">7</a>, <a href="#ref-8">8</a>, <a href="#ref-9">9</a>]. In veterinary virology, MD simulations are employed to study envelope proteins of pathogens such as canine coronavirus, avian influenza virus, and white spot syndrome virus (WSSV) [<a href="#ref-2">2</a>, <a href="#ref-10">10</a>, <a href="#ref-11">11</a>]. This article provides an exhaustive technical review of the principles and applications of MD simulations in the context of drug docking and design targeting viral envelope proteins.

Force Fields and Simulation Setup

All-atom MD simulations resolve the motions of every atom in a protein, solvent system over time, typically on nanosecond to microsecond scales [<a href="#ref-5">5</a>, <a href="#ref-12">12</a>, <a href="#ref-13">13</a>]. The accuracy of these simulations depends critically on the choice of force field, which defines the potential energy function for bonded and non-bonded interactions. Commonly employed force fields for envelope protein simulations include AMBER and CHARMM families, both extensively parameterized for proteins, lipids, and explicit solvent models such as TIP3P [<a href="#ref-5">5</a>, <a href="#ref-13">13</a>, <a href="#ref-14">14</a>]. Simulations of viral glycoproteins are frequently performed with the AMBER ff14SB or CHARMM36 force fields, as these parameters adequately represent the conformational flexibility of loops and glycosylated regions [<a href="#ref-4">4</a>, <a href="#ref-12">12</a>]. The preparation of a simulation system involves solvation in a periodic water box, neutralization with counterions, and energy minimization followed by equilibration under constant temperature (typically 310 K) and pressure (1 bar) using algorithms such as the Langevin thermostat and Berendsen barostat [<a href="#ref-6">6</a>, <a href="#ref-13">13</a>]. Long-range electrostatic interactions are treated with the Particle Mesh Ewald (PME) method, and hydrogen bond lengths are constrained with the SHAKE or LINCS algorithm to permit integration timesteps of 2 fs [<a href="#ref-12">12</a>, <a href="#ref-13">13</a>].

For membrane-bound envelope proteins such as influenza hemagglutinin or coronavirus spike trimers, the protein is embedded in a lipid bilayer (e.g., POPC or a complex raft mixture) using tools such as CHARMM-GUI [<a href="#ref-5">5</a>, <a href="#ref-6">6</a>]. The membrane environment imposes anisotropic forces that stabilize certain helical bundles and influence the exposure of hydrophobic fusion peptides [<a href="#ref-6">6</a>, <a href="#ref-15">15</a>]. Simulations of the SARS-CoV-2 envelope (E) ion channel have demonstrated how membrane embedding affects conformational dynamics and drug binding [<a href="#ref-6">6</a>]. The setup and production phases are typically followed by extensive validation: monitoring root-mean-square deviation (RMSD) of the backbone, radius of gyration, and solvent-accessible surface area (SASA) to ensure equilibration [<a href="#ref-5">5</a>, <a href="#ref-13">13</a>].

Capturing Conformational Dynamics and Cryptic Pockets

Envelope proteins often sample multiple conformations that are not visible in a single X-ray or cryo-EM structure [<a href="#ref-5">5</a>, <a href="#ref-14">14</a>]. MD simulations can reveal transient openings of loops or rearrangements of secondary structure elements that create cryptic or allosteric binding sites [<a href="#ref-5">5</a>, <a href="#ref-7">7</a>, <a href="#ref-16">16</a>]. For example, simulations of the SARS-CoV-2 spike trimer identified a potential allosteric site that stabilizes the receptor-binding domain (RBD) in the “down” conformation, thereby reducing ACE2 binding affinity [<a href="#ref-14">14</a>]. Similarly, ensemble-based docking to MD-generated conformations of the spike protein uncovered hidden pockets suitable for small-molecule inhibitors [<a href="#ref-5">5</a>, <a href="#ref-16">16</a>].

Cryptic pocket detection often employs methods such as solvent mapping (FTMap), principal component analysis (PCA) of trajectory data, and Markov state models (MSMs) that capture slow conformational transitions [<a href="#ref-5">5</a>, <a href="#ref-14">14</a>]. MSMs, in particular, can partition the conformational space into metastable states and compute transition probabilities, enabling the identification of low-population, druggable states [<a href="#ref-5">5</a>, <a href="#ref-13">13</a>]. In the context of flavivirus envelope proteins, MD simulations of dengue virus (DENV) pre-fusion envelope protein were used to screen inhibitors that bind to a pocket near the hinge region, preventing the pH-induced conformational change required for fusion [<a href="#ref-7">7</a>]. The combination of long timescale simulations (microseconds) with free energy calculations has proven effective for such targets.

Integration with Docking and Virtual Screening

MD simulations provide an ensemble of receptor conformations that more accurately reflect the dynamic binding landscapes encountered by ligands in vivo [<a href="#ref-5">5</a>, <a href="#ref-7">7</a>, <a href="#ref-17">17</a>]. Traditional rigid-receptor docking can miss interactions that require protein rearrangements. To overcome this, an ensemble docking protocol is employed: multiple snapshots from an MD trajectory are extracted, and each is used as a receptor structure for docking libraries of small molecules or peptides [<a href="#ref-7">7</a>, <a href="#ref-17">17</a>, <a href="#ref-18">18</a>]. The docking scores are then averaged or ranked based on consensus across frames [<a href="#ref-17">17</a>, <a href="#ref-19">19</a>]. This approach has been applied to identify inhibitors of the respiratory syncytial virus (RSV) fusion protein, where MD-derived conformations of the pre-fusion (Pre-F) protein were targeted with benzimidazole-based libraries [<a href="#ref-17">17</a>, <a href="#ref-19">19</a>].

Virtual screening pipelines integrate MD with ligand docking and scoring functions such as AutoDock Vina, Glide, or GOLD [<a href="#ref-7">7</a>, <a href="#ref-20">20</a>, <a href="#ref-21">21</a>]. More advanced workflows include free energy perturbation (FEP) calculations or MM/GBSA (Molecular Mechanics Generalized Born Surface Area) methods to compute binding affinities more accurately than empirical scoring functions [<a href="#ref-7">7</a>, <a href="#ref-22">22</a>, <a href="#ref-23">23</a>]. For instance, MM/GBSA was used to rank the binding of pinoresinol (an olive-derived lignan) to the spike RBD of Omicron variants, confirming stable interaction [<a href="#ref-23">23</a>]. The synergy between MD and docking also facilitates the design of peptide inhibitors targeting heptad repeat domains, as demonstrated for MERS-CoV and canine coronavirus spike proteins [<a href="#ref-1">1</a>, <a href="#ref-2">2</a>]. The overall workflow is summarized in Figure 1.

Analysis of Protein, Ligand Interactions

Post-simulation analysis of protein, ligand complexes involves monitoring RMSD, hydrogen bond occupancy, salt bridge stability, and per-residue decomposition of binding free energies [<a href="#ref-4">4</a>, <a href="#ref-22">22</a>, <a href="#ref-24">24</a>]. Root-mean-square fluctuation (RMSF) identifies flexible regions of the envelope protein that become rigidified upon ligand binding [<a href="#ref-4">4</a>, <a href="#ref-12">12</a>]. For antibody-based inhibitors, simulations can assess the stability of the complementarity-determining regions (CDRs) and the effect of mutations in the viral epitope [<a href="#ref-9">9</a>, <a href="#ref-25">25</a>]. For example, computational affinity dynamics between SARS-CoV-2 spike variants and ACE2 were characterized using MD to compute binding free energies and identify escape mutations [<a href="#ref-12">12</a>, <a href="#ref-26">26</a>]. In veterinary applications, the design of a peptide inhibitor targeting canine coronavirus spike-mediated fusion was guided by MD simulations of the HRC domain, validating stable helical interactions with the HR1 domain [<a href="#ref-2">2</a>].

Applications to Specific Viral Envelope Proteins

Coronaviruses

Coronavirus spike proteins are class I fusion glycoproteins that undergo large conformational changes from pre-fusion to post-fusion states [<a href="#ref-1">1</a>, <a href="#ref-5">5</a>, <a href="#ref-14">14</a>]. MD simulations have been instrumental in mapping the flexibility of the RBD and identifying druggable pockets in the S2 subunit. Allosteric inhibitors that stabilize the RBD down state have been discovered through ensemble docking of spike trimer simulations [<a href="#ref-14">14</a>]. Simulations of the spike from Omicron variants revealed altered hydrogen bonding networks that enhance ACE2 binding despite escape from neutralizing antibodies [<a href="#ref-23">23</a>, <a href="#ref-26">26</a>]. The SARS-CoV-2 envelope (E) ion channel has also been simulated in a membrane environment to study inhibitor binding that blocks ion conductance, which is critical for viral pathogenesis [<a href="#ref-6">6</a>].

Influenza A Virus

Influenza hemagglutinin (HA) is another class I fusion protein. MD simulations of HA have been used to study the stability of the fusion peptide and to screen inhibitors that block the low-pH-induced conformational change [<a href="#ref-10">10</a>, <a href="#ref-15">15</a>]. The M2 proton channel, though not an envelope protein per se, is a viroporin that cooperates with HA and has been targeted by virtual docking of approved drugs such as pamiparib [<a href="#ref-10">10</a>]. For drug-resistant M2 mutants (V27A/S31N), MD-guided design of novel inhibitors has been reported [<a href="#ref-15">15</a>]. Avian influenza H3N2 inhibitors were also identified by docking and MD validation of natural products targeting HA [<a href="#ref-20">20</a>].

Flaviviruses and Other Enveloped Viruses

MD simulations of dengue virus envelope protein have enabled the discovery of potent inhibitors that bind to the pre-fusion pocket, preventing the conformational rearrangements required for endosomal fusion [<a href="#ref-7">7</a>]. For West Nile virus (WNV), phytocompounds were screened against the envelope glycoprotein to block host cell attachment and membrane fusion [<a href="#ref-3">3</a>]. The E8 surface protein of monkeypox virus was also investigated using MD to identify potential therapeutic agents [<a href="#ref-27">27</a>]. In the field of invertebrate virology, the envelope protein of white spot syndrome virus (WSSV) was targeted with peptide inhibitors using computational docking and MD, offering a strategy for managing outbreaks in shrimp aquaculture [<a href="#ref-11">11</a>].

Filoviruses and Paramyxoviruses

Ebola virus glycoprotein (GP) and Marburg virus VP40 matrix protein have been subjected to MD simulations for drug design. Fragment-based design of monoterpenoid inhibitors targeting Ebola GP employed QSAR and MD to optimize binding [<a href="#ref-8">8</a>]. For Marburg VP40, machine learning and MD identified natural compounds with antiviral activity [<a href="#ref-28">28</a>]. Similarly, Sudan ebolavirus VP40 was studied with in silico peptide lead identification [<a href="#ref-29">29</a>]. In paramyxoviruses, human metapneumovirus (HMPV) fusion protein was targeted with dimeric catechins that stabilize the pre-fusion state; MD simulations elucidated the molecular logic behind this stabilization through galloylation-driven anchoring [<a href="#ref-4">4</a>]. Machine learning and MD also revealed potential inhibitors of HMPV fusion protein [<a href="#ref-30">30</a>]. RSV fusion protein inhibitors derived from benzimidazoles were validated through MD simulations [<a href="#ref-19">19</a>].

Challenges and Future Directions

Despite its power, MD-based drug design faces several challenges. The timescales accessible to conventional all-atom MD (microseconds) may be insufficient to capture slow domain motions of envelope proteins, such as the large-scale hinge movements of spike trimer subunits [<a href="#ref-5">5</a>, <a href="#ref-13">13</a>]. Enhanced sampling methods (e.g., replica exchange, metadynamics, and accelerated MD) are often required to explore the full conformational landscape [<a href="#ref-5">5</a>, <a href="#ref-14">14</a>]. Force field inaccuracies, particularly for glycans and lipid interactions, can bias simulations of heavily glycosylated envelope proteins [<a href="#ref-5">5</a>, <a href="#ref-7">7</a>]. The glycan shield can be explicitly modeled, but carbohydrate parameterization remains an active area of development.

Integration with machine learning is a growing trend. Deep learning models can predict high-affinity antibodies against envelope proteins, as shown for Zika virus, where MD validated the predicted interactions [<a href="#ref-9">9</a>]. Reinforcement learning has been applied to design ab initio inhibitors for RSV Pre-F from natural fragment libraries [<a href="#ref-17">17</a>]. AlphaFold predictions of envelope protein structures are increasingly used as starting points for MD simulations, though the dynamic ensembles generated by MD remain essential for drug docking [<a href="#ref-9">9</a>, <a href="#ref-22">22</a>]. For veterinary applications, the development of species-specific force fields or scaling models for animal body temperatures (e.g., 38-42 °C for poultry) could improve relevance [<a href="#ref-2">2</a>, <a href="#ref-11">11</a>]. Continued advances in hardware (GPU acceleration) and software (e.g., GROMACS, NAMD, Amber) will enable longer and more complex simulations, including multi-protein systems and full viral particles.

Conclusion

Molecular dynamics simulations have become a cornerstone of computational virology, providing atomic-level insights into the conformational dynamics of viral envelope proteins and guiding the rational design of inhibitors. From identifying cryptic pockets and allosteric sites to validating docking poses and estimating binding affinities, MD simulations bridge the gap between static structural biology and dynamic drug binding. As veterinary medicine increasingly adopts computational approaches, MD-based workflows will play a vital role in developing antivirals for emerging zoonotic and animal-specific pathogens.

The following Mermaid diagram summarizes the integrated workflow of MD simulations for drug docking and design.

flowchart TD
 A["Obtain envelope protein structure (X-ray, cryo-EM, or AlphaFold)"] --> B["Prepare system: force field, solvation, membrane embedding"]
 B --> C["Equilibration and production MD simulation"]
 C --> D["Conformation ensemble analysis: RMSD, RMSF, PCA, MSM"]
 D --> E["Detect cryptic/allosteric pockets via solvent mapping or MD snapshots"]
 E --> F["Ensemble docking of small-molecule/peptide libraries"]
 F --> G["Scoring and ranking: docking scores, MM/GBSA, FEP"]
 G --> H{"Experimental validation?"}
 H -->|"Yes"| I["Lead optimization with additional MD rounds"]
 H -->|"No"| J["Re-run simulation with modified candidate"]
 I --> K["Final candidate for in vitro/vivo testing"]

References

Related Clinical & Scientific Guides

• A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
• Computational Immunology: Modeling the Immune System
• How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices