Structure-based virtual screening and molecular dynamics studies to explore potential natural inhibitors against 3C protease of foot-and-mouth disease virus
Sthitaprajna Sahoo
Hak-Kyo Lee1,2* 
Donghyun Shin- 1Department of Agricultural Convergence Technology, Jeonbuk National University, Jeonju, Republic of Korea
- 2Department of Animal Biotechnology, Jeonbuk National University, Jeonju, Republic of Korea
Foot-and-mouth disease (FMD) is a highly infectious animal disease caused by foot-and-mouth disease virus (FMDV) and primarily infects cloven-hoofed animals such as cattle, sheep, goats, and pigs. It has become a significant health concern in global livestock industries because of diverse serotypes, high mutation rates, and contagious nature. There is no specific antiviral treatment available for FMD. Hence, based on the importance of 3C protease in FMDV viral replication and pathogenesis, we have employed a structure-based virtual screening method by targeting 3C protease with a natural compounds dataset (n = 69,040) from the InterBioScreen database. Virtual screening results identified five potential compounds, STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, with a binding affinity of −9.576 kcal/mol, −8.1 kcal/mol, −7.744 kcal/mol, −7.647 kcal/mol, and − 7.778 kcal/mol, respectively. The compounds were further validated through physiochemical properties and density functional theory (DFT). Subsequently, the comparative 300-ns MD simulation of all five complexes exhibited overall structural stability from various MD analyses such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), H-bonds, principal component analysis (PCA), and free energy landscape (FEL). Furthermore, MM-PBSA calculation suggests that all five compounds, particularly STOCK1N-62634, STOCK1N-96109, and STOCK1N-94672, can be considered as potential inhibitors because of their strong binding affinity toward 3C protease. Thus, we hope that these identified compounds can be studied extensively to develop natural therapeutics for the better management of FMD.
Introduction
Foot-and-mouth disease (FMD) is a highly infectious disease that significantly affects global livestock production by affecting cloven-hoofed animals such as cattle, pigs, sheep, goats, and deer (1). FMD outbreaks have occurred in various parts of the world over the years. These outbreaks can have significant economic and agricultural impacts due to trade restrictions imposed on countries affected by FMD (2). Due to the highly infectious nature of the virus and the heightened international trade in animals and animal products, it is interesting to note that no nation is entirely safe from the threat. The World Organization for Animal Health (WOAH) urges global cooperation in addition to independent regional initiatives to combat the disease successfully. The disease is caused by a single-stranded positive RNA virus named foot and mouth disease virus (FMDV), which belongs to the Picornaviridae family and has seven serotypes (O, A, C, Asia1, and SAT1-3) and numerous strains, making it a challenging disease to control and eradicate (3). In most FMD outbreak regions, including Asia, the Middle East, India, Africa, and South America, the O serotype is prevalent (4). The symptoms of FMDV can vary depending on the species and age of the affected animal as well as the serotype and strain of the virus. Generally, severe symptoms are noticed in cattle and reared swine, and they exhibit clinical symptoms 24 to 48 h after FMDV infection. At this point, high viral counts can be observed in the serum of infected animals. Fever, shivering, salivary drooling, and the development of blisters or vesicles on the epithelium of the tongue, nose, coronary bands, and teats are the hallmarks of FMD (5).
The treatment of FMDV primarily focuses on controlling the spread of the virus and managing the clinical symptoms in affected animals. Vaccination against FMDV was vital in maintaining and preventing most FMD outbreaks (6). Countries adopt different vaccination strategies based on the prevailing FMDV strains and epidemiological situations where vaccine-containing inactivated FMDV strains were commonly used (2, 7). Although FMDV vaccines can reduce some FMD outbreaks, vaccination against FMDV does not prevent the transmission of the virus and takes several days to trigger an immune response against FMDV. Additionally, it might be challenging to match the virus vaccine to the multiple serotypes and subtypes. To significantly cope with the genetic variability and reduce clinical signs, the ideal antiviral treatments should be instantly and widely active on diverse virus serotypes. However, it is essential to note that no specific antiviral therapy is available for FMDV, and most efforts are directed toward prevention, vaccination, and supportive care; thus, antiviral therapeutics are still required for quick prevention and management of FMD spread.
FMDV is a small, non-enveloped, spherical virus belonging to the family of Picornaviridae, having a diameter of 25–30 nm. FMDV consists of a single-stranded +ve sense RNA genome having 8,400 nucleotides (nt), which includes 5’UTR (1,300 nt), open-reading frame (ORF) (7,000 nt), and 3’UTR (90 nt) (1). This long ORF of the viral genome encodes a single large polyprotein, which then undergoes proteolytic processing to form 4 structural and 10 non-structural proteins (Figure 1A). Primary co-translation produces four primary individual products, namely, Lpro, P1/2A, P2, and P3. P1/2A is further processed to form four structural proteins, i.e., VP4, VP2, VP3, and VP1. Lpro is a papain-like protease that releases itself from the polyprotein by cleaving its own C-terminal and N-terminal of VP4. The P2 and P3 portions of the polyprotein were further processed to form mature individual proteins, i.e., 2A, 2B, 2C, 3A, 3B1, 3B2, 3B3, 3Cpro, and 3Dpol (8). Upon entry into receptive host cells, the viral genome is immediately translated into a long polyprotein precursor that must then be broken down by virally encoded protease to produce the 14 individual protein products that are needed for RNA replication and for the construction of new viral particles (9). One of the most highly conserved genes in the viral genome, FMDV 3C protease (3Cpro), is the key component of the aforementioned proteolytic activity and is responsible for 10 of the 13 cleavages. It is a trypsin-like serine protease that also cleaves the host proteins present in the infected cells such as eIF4G (10). Additionally, during FMDV infection, FMDV 3Cpro cleaves histone H3, causing a generalized repression of cellular transcription (11). The crucial role of 3Cpro in the viral replication process and cleaving activity makes it a possible antiviral drug target that can be combined with other strategies to prevent and control FMDV. There are antiviral medications available that target the picornavirus protease. For instance, rupintrivir, an inhibitor of human rhinovirus (HRV) 3C protease, demonstrated strong antiviral activity in cultured cells. However, the effectiveness of rupintrivir against a natural HRV infection was insignificant (12, 13).
Figure 1. (A) Genome organization FMDV. (B) The modeled structure of 3C protease is represented in a pale green cartoon structure where catalytic residues (C163, H46, and D84), and active site residue C142 is represented in the stick model.
Research on natural compounds as antiviral therapeutics has gained considerable attention due to the potential benefits of these compounds, including their broad-spectrum activity, low toxicity, and diverse mechanisms of action (14). Many natural compounds isolated from different natural resources such as plants, microorganisms, and marine species have exhibited inhibitory activity against a wide range of pathogens and viruses (15, 16). Since the extraction of natural compounds and antiviral activity evaluation require funding as well as specialized experimental equipment, exploring potent FMDV inhibitors from natural sources is challenging. This strategy can also be integrated with molecular biological approaches to uncover the novel mechanisms of inhibition using structural relationships. In recent years, the integration of structure-based virtual screening (SBVS), MD simulations, and density functional theory (DFT) analysis enhances the drug discovery process by providing a multi-faceted understanding of the molecular interactions between drugs and their targets. This integrated approach enables more informed decision-making in the identification, optimization, and design of novel drug candidates (17–21). Hence, in the current study, we explored potential natural compounds against FMDV 3Cpro through structure-based virtual screening. Moreover, we have also employed a molecular dynamics (MD) approach to analyze the binding affinities and detailed structural interaction of the selected inhibitors with the targeted protein.
Theory and methods
Preparation of target protein
All the crystal structures of FMDV 3Cpro that were available in the RCSB Protein Data Bank1 contained mutations either at active site residues (C142L and C142S) or catalytic residues (C163A and D84E) (22). As serotype O was prevalent in most of the FMD outbreaks, we retrieved the full-length amino acid sequence (213 amino acids) of FMDV 3Cpro of O serotype from UniProt database (Uniprot ID: P03305)2 and modeled the 3D structure of wild type 3Cpro using the most advanced deep learning model AlphaFold2, where structures having PDB ID: 2BHG, 2 J92, 2WV4, and 5HM2 were used as template (Figure 1B) (23). We have also validated the modeled structure of FMDV 3C pro using the Ramachandran plot (Supplementary Figure S1). AlphaFold2 is an advanced deep learning model developed by DeepMind that specializes in protein folding prediction. While it can be used for a variety of protein structure prediction tasks, including de novo folding, one of its significant applications is in homology modeling. The protein preparation wizard of Schrodinger was used to prepare the target protein structure where missing atoms were incorporated, ionization states were generated, and proper bond orders were assigned for optimal docking analysis (24).
Selection of natural compounds
We retrieved 69,040 compounds from the natural compound dataset of the InterBioScreen database.3 The InterBioScreen database is a comprehensive collection of natural compounds and small molecules that are of interest in drug discovery and medicinal chemistry research. It contains a vast array of compounds that come from naturally occurring substances, such as plants, microorganisms, marine organisms, and other biological materials, that assist researchers in the identification and selection of potential lead compounds for drug development. All the natural compounds from this database were downloaded in a structured data file format and considered for virtual screening.
Structure-based virtual screening
Structure-based virtual screening (SBVS) is a computational technique used in drug discovery to identify and prioritize potential drug candidates from large compound libraries or databases with a certain molecular target (25). It involves the use of computer algorithms and molecular modeling methods to predict the binding affinity and biological activity of compounds toward a target protein or biological target along with their drug-like properties using Lipinski’s rule of five. Glide and virtual screening module of Schrodinger was used to perform molecular docking in three stages, i.e., high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP) using OPLS3 force field (26–28). Prior to the molecular docking, a grid box was generated using the catalytic and binding site residues such as C142, S182, D84, C163, and H46 of FMDV 3C protease (29–31). The interaction of the top five screened compounds was visualized using Maestro, Discovery studio, and PyMOL to determine the crucial amino acid residues that were involved in protein-ligand interaction (32, 33).
Density functional theory calculation
Geometry optimization of the selected five compounds was performed through DFT using the DMOL3+ server in BIOVIA Material Studio version 2017 R2. Generalized gradient approximation (GGA) and Perdew–Burke–Ernzerhof (PBE) functions with double numerical plus polarization (DNP) basis set were chosen to calculate the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), energy gap, and electronic chemical potential (μ) (34). The maximum force, energy convergence accuracy, gradient convergence, maximum displacement, displacement convergence, and optimal iterations were set to be 0.004 Ha/A, 2 × 10−5 Ha, 4 × 10−3 Å, 0.3 Å, 5 × 10−3 Å, and 50, respectively.
Molecular dynamics simulation
Molecular dynamics (MD) simulations of the top five screened protein-ligand complexes were performed using the GROMACS 2022.2 package and CHARMM36 force field to analyze the conformational changes and dynamic behavior at the atomistic level (35–37). The ligand topology of the selected five compounds was generated using the CGenFF webserver, and Na+/Cl- counter ions were added to neutralize the complex structures (38). The system was further solvated using a TIP3P water model and a dodecahedron periodic boundary condition. The steepest-descent algorithm was used to minimize the system’s energy followed by conjugate gradient (39). The particle Mesh Ewald (PME) method was used to compute the long-range electrostatic interactions, and the LINear Constraint Solver LINC0 algorithm was used to calculate the Lennard-Jones and Coulomb interactions with a radius of 10 Å cutoff distance (40–42). The temperature coupling was maintained at 300 K using the Berendsen thermostat (V-rescale) coupling algorithm, while pressure coupling was maintained at 1 Pascal bar using the Parrinello–Rahman pressure coupling algorithm (43, 44). After energy minimization, the system was equilibrated for 100 ps each under NVT and NPT conditions to retain the volume, pressure, and temperature. All five systems were subjected to 300 ns MD production run and coordinates were saved at a consistent time step of 2 fs intervals. From the final trajectory information, we studied the root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), and H-bonds analysis using gromac’s modules. The global motion of all complexes was also analyzed with principal component analysis (PCA) considering the last 20 ns trajectory information, and free energy landscape (FEL) was studied from the first two principal components (PCs) using the gmx sham tool, as reported in our previous study (45, 46).
Binding free energy calculation
Molecular mechanics/Poisson–Boltzmann surface area (MM-PBSA) is a computational method used to estimate the binding free energies of biomolecular complexes such as protein–ligand, protein–DNA, or protein–protein complexes (47). It combines molecular mechanics force fields, continuum solvent models, and empirical solvation energy terms to estimate the thermodynamic properties of protein–ligand interactions. MM-PBSA also allows the decomposition of the free energy of binding into different components, such as van der Waals interactions, electrostatic interactions, solvation-free energy, and entropy contributions. This decomposition helps in understanding the individual energetic contributions of different molecular interactions that allow researchers to focus on specific aspects of ligand optimization, such as improving hydrogen bonding, reducing steric clashes, or enhancing electrostatic interactions. MM-PBSA calculations are commonly used in computer-aided drug design and lead optimization studies to evaluate the binding affinities of ligands, compare different ligand poses, and guide the selection of promising drug candidates (48). MM-PBSA results can be compared with experimental binding affinities to validate the accuracy of the computational predictions. This validation is essential for building confidence in the computational approach and its applicability in drug discovery. Mathematically, binding free energy (ΔGbind) can be represented as follows:
Here, ΔGmm represents the molecular mechanics energy by considering van der Waals and electrostatic interactions. ΔGps refers to polar solvation energy, and ΔGnps refers to non-polar solvation energy. TΔS represents total entropic contribution, where T and S are denoted as temperature and entropy, respectively. The binding-free energy of the selected five compounds was calculated with respect to the 3C protease protein from the last 200 ns MD trajectory using the gmx_MMPBSA tool (49, 50).
Results
Identification of potent inhibitors through SBVS
Structure-based virtual screening (SBVS) in drug discovery offers several advantages, including the ability to screen a large number of compounds efficiently, reduced cost and time compared to experimental screening, and the potential to explore a wide chemical space. In the current study, we have considered the natural compounds subset from the InterBioScreen database (n = 69,040) against FMDV 3C protease. Initially, the large natural compounds dataset was screened with the Qikprop module, following Lipinski’s Ro5 filtering, where 48,952 compounds were screened. The high throughput virtual screening module screened 4,895 compounds, and 2,695 natural compounds were filtered for SP docking. In total, 269 compounds were pulled down from SP, and finally, 26 compounds were considered for XP docking for more accurate screening. Subsequently, the top five protein-ligand complexes with maximum glide score (−9.576 to −7.778 kcal/mol) were considered for further analysis (Table 1).
Table 1. Binding affinity of the top five screened compounds along with the amino acid residues involved in interactions.
All the top five screened compounds such as STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570 exhibited the highest binding affinity of −9.576, −8.1, −7.744, −7.647, and − 7.778 kcal/mol, respectively. STOCK1N-62634 exhibited the highest binding affinity of −9.576 kcal/mol by forming four H-bond interactions with H46, D146, D144, and G184 and two pi-alkyl interactions with M143 and A183. Similarly, compound STOCK1N-96109 exhibited the binding affinity of −8.1 kcal/mol by forming four hydrogen bond interactions with E50, D146, M143, and G184 and one pi-cation interaction with the active site residue H46 (Figure 2). Despite having two H-bond interactions with 3C protease, compounds STOCK1N-94672 and STOCK1N-89819 exhibited a strong binding affinity. The 2D structures of the screened compounds and their 2D molecular interaction profile with the 3C protease protein are shown in Supplementary Figure S2.
Figure 2. Interaction of the top five screened compounds (STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570) with the 3C protease protein. Interacting residues of the protein were labeled in the pale green stick model, and interactions were represented in a violet dashed line.
Assessment of drug-likeness through physiochemical properties
The pharmacokinetics and physiochemical properties of all the five screened compounds were evaluated using the QikProp module of Schrodinger. To acquire a safe side and potentially optimal drug-like behavior, a compound needs to fulfill at least four of the five characteristics of Lipinski’s rule of five, i.e., molecular weight < 500 Daltons, ≤5 hydrogen bond donors, and ≤ 10 hydrogen bond acceptors, lipophilicity <5, and molar refractivity between 40 and 130. The physiochemical and drug-likeness properties were evaluated and mentioned in Table 2. The ADME prediction results revealed that all the selected five natural compounds have an acceptable number of H-bond donors and acceptor, molecular weight, and lipophilicity with no violation of Lipinski’s rule of five. Hence, all these compounds can be considered as ideal drug candidates.
Table 2. Predicted ADME properties for top hit lead compounds.
Quantum mechanics
The quantum mechanics studies were performed using geometry optimization-based DFT calculations. DFT analysis results revealed energy values of HOMO, LUMO, energy gap, and electronic chemical potential (μ) and are shown in Table 3. HOMO values represent the electrophilic attack site by calculating its potential to distribute and donate electrons to the receptor, while LUMO represents the nucleophilic attack site by calculating its ability to accept electrons from the receptor (51).
Table 3. Frontier molecular orbital energy of the selected five compounds.
Figure 3 represents the molecular orbital plot of the selected compounds in which the location of the molecules can be identified from HOMO, LUMO energy values, that is represented in Table 3. Blue and yellow isosurfaces exhibit electron build-up region and electron loss regions, respectively. The HOMO energy values range from −0.205 eV in STOCK1N-89819 to −0.169 eV in STOCK1N-62634, while the LUMO energy values ranged from −0.103 eV in STOCK1N-89819 to −0.050 in STOCK1N-94672. The energy gap (HOMO/LUMO Gap) is the energy difference between the LUMO and HOMO levels, which is inversely proportional to the stabilizing energy of the molecule. The calculated energy gap values of all the selected compounds varied from 0.081 eV to 0.126 eV. Electronic chemical potential (μ) is the total amount of energy required to remove one single electron from the molecule (51). The μ values range from −0.06 eV in STOCK1N-94672 and STOCK1N-80570 to −0.04 eV in STOCK1N-62634. Low values of μ in all the five compounds suggest strong binding affinities. Low energy gap values and negative HOMO and LUMO values of all the selected compounds suggest their ability to donate and accept electrons with the receptor for a strong protein–ligand interaction and, thus, high inhibition of 3C protease protein.
Figure 3. Frontier molecular orbital plot of HOMO and LUMO. Blue and yellow isosurfaces illustrate the electron build-up region and electron loss regions, respectively.
Molecular dynamics simulation analysis
The structural stability and dynamic behaviors of all the top five screened compounds with 3C protease protein was investigated through the 300 ns MD simulation study.
Structural deviation analysis through RMSD
The structural stability of all five screened complexes were analyzed through RMSD of the backbone atoms (C, Cα, and N) from 300 ns MD simulation trajectory. The RMSD graph of the protein backbone with reference to the initial conformation with ligands STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570 exhibited an average RMSD of 0.51 nm, 0.36 nm, 0.37 nm, 0.44 nm, and 0.33 nm, respectively (Figure 4A). With all the five compounds, the 3C protease protein is getting stable toward the last 250 ns of MD simulation. These data also suggest that the compounds STOCK1N-96109, STOCK1N-94672, and STOCK1N-89819 are highly stable complexes relative to others. Moreover, the dynamic behaviors of all five compounds were studied by calculating the RMSD of each ligand, and the result revealed that compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570 exhibited an average RMSD of 0.29 nm, 0.18 nm, 0.24 nm, 0.17 nm, and 0.18 nm, respectively (Supplementary Figure S3A). In particular, compounds STOCK1N-94672 and STOCK1N-96109 exhibited higher ligand stability than other compounds.
Figure 4. Stability analysis through 300 ns MD simulation trajectory. (A) RMSD of backbone atoms. (B) Residual flexibility analysis through RMSF. (C) Compactness analysis through Rg. (D) Changes in the number of H-bonds throughout the simulation. Black, green, red, blue, and purple color represent the complexes with 3C protease and compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively.
Residual flexibility analysis through RMSF
The fluctuation of each amino acid in 3C protease protein during the MD simulation was analyzed through RMSF and shown in Figure 4B. The receptor 3C protease after binding to compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570 exhibited an average RMSF of 0.16 nm, 0.15 nm, 0.15 nm, 0.17 nm, and 0.16 nm, respectively. Interestingly, the residual fluctuations are almost similar for all five screened compounds. After binding to all five compounds, the active site and catalytic triad residues such as C142, D84, C163, and H46 have been revealed to be stable throughout the 300 ns trajectory. G78, V124, D144, and A160 of 3C protease exhibited higher fluctuations after binding to ligands.
Compactness analysis through Rg
The compactness of the receptor after binding to ligands was analyzed through Rg. It is the average distance of each atom in the molecule from the center of mass, which provides insights into the overall spatial distribution of atoms within the molecule. The Rg calculation revealed an average value of 1.69 nm, 1.67 nm, 1.67 nm, 1.68 nm, and 1.66 nm for the complexes with compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively (Figure 4C). Almost all the five compounds are getting stable toward the last 150 ns of simulation. In particular, the complexes with STOCK1N-96109 showed a more stable and compact conformations as compared to others.
Solvent accessible surface area analysis
The SASA of all five complexes was computed throughout the simulation and represented in Supplementary Figure S3B. The SASA results revealed an average value of 121.74nm2, 119.73 nm2, 119.29 nm2, 119.82 nm2, and 118.97 nm2 for receptor 3C protease with ligands STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively. All five complexes exhibited similar kinds of SASA values with an average value of 119 nm2, revealing minimal changes after the binding of each compound.
Interaction analysis through hydrogen bonding
The stability of all five complexes and protein–ligand interactions were analyzed through the number of hydrogen bonds between them throughout the 300 ns MD simulations and represented in Figure 4D. Compounds STOCK1N-62634 and STOCK1N-96109 exhibited an average of 0–2 number of H-bonds throughout the simulation, while compounds STOCK1N-89819 and STOCK1N-80570 exhibited an average of 0–1 number of H-bonds. Compound STOCK1N-94672 showed the highest number of H-bonds interacting with the receptor (average of 2–3) as compared to others.
Principal component analysis-based free energy landscape
Principal component analysis was performed to understand the global concentrated motion of all five complexes after ligand binding. This method can help to identify the most significant conformational changes that take place during the simulation. The first 30 eigenvectors were considered to calculate the eigenvalues and percentage changes from the last 200 ns MD simulation (Figure 5A). The first two PCs exhibited larger motions as compared to others with a cumulative variance of 44.97, 52.16, 41.20, 72.88, and 66.14% for STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively (Figure 5B). Furthermore, a 2D plot was generated by assessing PC1 and PC2 for all five complexes and is represented in Figure 5C. Complexes with compounds STOCK1N-62634, STOCK1N-96109, and STOCK1N-94672 showed less overall distribution of protein conformation as compared to others. We have also analyzed the Gibbs free energy landscape (FEL) using the first two PCs in order to visualize the transformation that takes place during MD simulation where deep blue color represents the global minima (Figure 6). FEL can accurately depict the most stable conformational ensembles by analyzing the structural changes in protein–ligand complexes. The receptor with compounds STOCK1N-62634 and STOCK1N-94672 exhibited one major cluster with two global minima, while others exhibited more than one cluster. Compound STOCK1N-89819 and STOCK1N-80570 exhibited dispersed clusters with one and two global minima, respectively. In addition, we have also generated the low energy conformations from all the global minima and aligned to their respective initial conformation (Supplementary Figure S4). Structural alignment results revealed that all the low-energy conformations except the complexes with compound STOCK1N-94672 are observed to be in the binding pocket of 3C protease.
Figure 5. Principal component analysis. (A) Eigenvalues derived from the last 200 ns simulation. (B) Cumulative percentage of first 30 eigenvectors. (C) Projection of first two principal components (PC1 and PC2) for all five complexes. Black, green, red, blue, and purple color represent complexes with 3C protease and compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively.
Figure 6. Gibbs’ free energy landscape using PC1 and PC2. The lower energy systems were illustrated in deep blue color.
Binding free energy analysis
The binding affinity of each ligand with the receptor was calculated from the last 200 ns of MD simulation using molecular mechanics/Poisson–Boltzmann surface area (MM-PBSA) approach. In this study, we have calculated the van der Waals energy (ΔVDW), electrostatic energy (ΔEFL), electrostatic contribution to the solvation free energy calculated by the Poisson–Boltzmann model (ΔEPB), non-polar contributions (ΔENPOLAR), and entropy (TΔS) of each complex and mentioned it in Table 4 and Figure 7. Compound STOCK1N-62634 showed the highest binding affinity of −192.25 kcal/mol with the receptor followed by compounds STOCK1N-96109 (−106.85 kcal/mol), STOCK1N-94672 (−97.86 kcal/mol), STOCK1N-89819 (−24.23 kcal/mol), and STOCK1N-80570 (−38.32 kcal/mol). According to the binding free energy analysis, electrostatic energy components (ΔEFL) and van der Waals energy components (ΔVDW) contribute significantly to the binding of these natural compounds.
Table 4. MM-PBSA based average binding free energy of 3C protease with each compound in kcal/mol.
Figure 7. Free energy components in kcal/mol. (A) Binding energy, ΔH is the summation of ΔVDW, ΔEFL, ΔEPB, and ΔENPOLAR. (B) Total energy calculated from the binding energy and entropy components, i.e., ΔG = ΔH-TΔS. Black, green, red, blue, and purple color represent complexes with 3C protease and compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively.
In addition, we have also analyzed the residual binding energy of each amino acid in the receptor upon ligand binding. This is a key method for determining the hotspot residues that play a major role in the ligand binding during simulation. Energy decomposition results revealed that most of the arginine residues such as R45, R60, R68, R92, R95, R97, R104, R108, R126, R155, and R196 contribute significantly to the ligand binding, with an average energy of −240 kcal/mol (Figure 8). Among active sites and catalytic residues, D84 contributed significantly with an average energy of −90 kcal/mol for all complexes. Residues such as D7, D24, D53, D66, D80, D98, D103, D174, and E213 also played a major role in energy contribution.
Figure 8. Energy decomposition of each amino acid in 3C protease upon ligand binding. Black, green, red, blue, and purple color represent complexes with 3C protease and compounds STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, respectively.
Discussion
FMDV is a highly infectious ssRNA virus that affects many species such as cattle, sheep, goat, and pig, significantly affecting the global livestock economy. Because of its high mutation rate and contagious nature, it infects more than 70 species of animals and has become a major health concern in livestock industries (1). Several nations encountered severe FMD outbreaks that affected the growth of animal husbandry and foreign trade, negatively impacting the global economy (52). FMDV outbreaks can lead to substantial direct economic losses for the livestock industry. These losses primarily stem from animal deaths, decreased productivity, and trade restrictions imposed to control the spread of the disease. As per a research study published in the journal Transboundary and Emerging Diseases, FMDV outbreaks in the United Kingdom between 2001 and 2010 resulted in direct costs ranging from £8 million (approximately $10 million) to £180 million (approximately $225 million) per outbreak (53). It can also lead to trade disruptions, with many countries imposing restrictions on the imports and exports of livestock and meat products from the affected regions. These disruptions can significantly impact the livestock industry’s ability to export and earn revenue. For example, during the FMDV outbreak in the United Kingdom in 2001, the country banned exports of live animals, meat, and dairy products, causing estimated losses of approximately £2.4 billion (approximately $3 billion) to the industry. Different countries adopted several vaccination strategies based on prevailing FMDV strains, but vaccination does not cure the infection and takes several weeks to elicit an immune response. Since no specific antiviral treatments are available for FMD, there is a need to explore new and effective antiviral drugs against FMDV. It is important to highlight that losses incurred by the livestock industry due to FMDV outbreaks can have severe consequences for livelihoods, animal welfare, food security, and the overall economy of affected regions or countries. By quantifying these losses, it highlights the urgency of developing effective antiviral treatments to combat FMDV, reduce the impact of outbreaks, and protect the livestock industry worldwide. The involvement of 3C protease in FMDV replication and proteolytic activity recommends it as one of the critical targets for antiviral drug development (10, 11). Research on natural compounds as antiviral medications has drawn much attention due to their low toxicity, broad-spectrum activity, and diverse mechanisms of action. However, exploring FMDV inhibitors from natural resources is challenging because the extraction of these compounds and the evaluation of antiviral activity demand high funding and specialized experimental equipment (14). The current study uses computational approaches to explore potential natural compounds from an extensive natural compounds database against 3C protease through structure-based virtual screening followed by MD simulation analysis of the top five screened compounds.
The natural product dataset from the InterBioScreen database, having 69,040 compounds, was retrieved and considered for structure-based virtual screening. Furthermore, the top five screened compounds having high binding affinity and interaction with crucial residues were considered for further analysis. The SBVS results revealed that compound STOCK1N-62634 exhibited the highest binding affinity of −9.575 kcal/mol by forming four hydrogen bonds with H46, D146, D144, and G184. Moreover, residues such as M143 and A183 were involved in pi-alkyl interaction with the compound STOCK1N-62634. STOCK1N-96109 interacts with 3C protease by forming four hydrogen bond interactions at E50, D146, M143, and G184 with a binding affinity of −7.878 kcal/mol. Additionally, it formed one pi-cation interaction with H46 and alkyl bonds with C163, A183, and A160. Compounds STOCK1N-94672 and STOCK1N-89819 exhibited the similar binding affinity of −7.674 kcal/mol and − 7.451 kcal/mol, respectively, by forming two hydrogen bonds and one pi-alkyl interaction, each with A160. Compound STOCK1N-94672 formed two H-bonds with H46 and G161, while compound STOCK1N-89819 formed two H-bonds with M143 and G184. Compound STOCK1N-80570 showed a docking score of −7.34 kcal/mol by forming five H-bonds with M143, H181, G184, N186, and Y190 and one pi-alkyl interaction with A160. Mostly, all five compounds formed at least one pi-alkyl interaction with A160. These non-covalent interactions are crucial in medicinal chemistry for designing drugs with high affinity, selectivity, and favorable pharmacokinetic properties. It enables medicinal chemists to optimize the interactions between drugs and their target proteins, ultimately developing effective and safe pharmaceutical agents. According to one of the recent studies by Theerawatanasirikul et al. (54), a similar attempt to identify small compounds against 3C protease and revealed the binding affinity of top hit compounds in the range from −6.6 to −6.8 kcal/mol. Most computational chemists focus on the structure or reactivity of the ligands, including geometry optimization of the molecule as a main component. According to the frontier molecular orbital (FMO) theory of chemical reactivity, the contact between HOMO and LUMO of the reacting species results in the transition of one electron. HOMO values explain the ability to donate electrons, while LUMO explains the ability to accept electrons from the receptor (55). Hence, geometry optimization of the selected five compounds was performed through DFT calculation, and all of them exhibited spatial arrangements of electron density in the HOMO that disperse across the covalent bonds formed between adjacent p orbital overlaps, whereas the electron density in the LUMO was dispersed across the aromatic rings. The negative HOMO and LUMO energy values derived from DFT calculation indicate their strong electron-accepting and electron-donating potentials. Moreover, the small energy gap values in all the compounds indicate the movement of electrons from the HOMO to LUMO, resulting a solid binding interaction with the receptor protein. Furthermore, we have also predicted the physiochemical properties of all five screened compounds, and all of them revealed good drug-like behavior without violating Lipinski’s rule of five. Therefore, a 300 ns MD simulation analysis was performed for 3C protease with all five screened compounds, i.e., STOCK1N-62634, STOCK1N-96109, STOCK1N-94672, STOCK1N-89819, and STOCK1N-80570, for a detailed structural insight and dynamic behavior of protein and protein–ligand complexes.
The dynamic behavior of all five complexes was analyzed through several structural parameters such as RMSD, RMSF, Rg, SASA, and H-bonds. The conformation stability of backbone atoms was calculated through RMSD analysis, and the results revealed that, except for STOCK1N-80570, all the four complexes exhibited stable conformation throughout the trajectory, but compound STOCK1N-80570 exhibited a slight fluctuation after 200 ns. However, it again gets stabilized toward the end. The conformational compactness of all complexes was also calculated through the Rg. The lower the Rg values, the higher the compactness. The calculated Rg values over the 300 ns MD simulation time scale showed overall compact and stable conformations of all five complexes. Moreover, all five complexes exhibited lower Rg values toward the last 200 ns of simulation, resulting in more compact and structural stability. In addition, we have calculated the RMSF, SASA, and H-bonds throughout the 300 ns MD trajectory, and interestingly, no residual fluctuations were observed in the active site and catalytic residues of 3C protease after ligand binding. The SASA analysis also supports the RMSD and Rg analysis by suggesting stable and minimal changes after ligand binding. All five screened compounds’ interaction with 3C protease was analyzed through H-bond analysis throughout the trajectory, and the results indicate that compound STOCK1N-94672 exhibited the highest average number of H-bonds with the receptor. In comparison, others formed at least one H-bond with the receptor throughout the trajectory, which indicates a strong interaction between all five natural compounds and the receptor. The global dynamic behavior of all complexes was analyzed through PCA and Gibb’s free energy landscape (FEL), and the results derived from the first two PCs indicate that 3C protease with all five compounds exhibited fewer dynamics. Mainly, compounds STOCK1N-96109 and STOCK1N-94672 exhibited less overall relative motion than others. Furthermore, we have calculated the binding free energy and each residual binding energy using the MM-PBSA method from the last 200 ns of MD simulations. MM-PBSA is a widely used approach to calculate the free energy in biomolecular complexes such as the protein–ligand complex, protein–DNA complex, or protein–protein complex (47, 49). MM-PBSA calculations are commonly used in computer-aided drug design and lead optimization studies to evaluate the binding affinities of ligands, compare different ligand poses, and guide the selection of promising drug candidates. MM-PBSA analysis suggests that compound STOCK1N-62634 showed the highest binding affinity of −192.25 kcal/mol with the receptor, followed by STOCK1N-96109 (−106.85 kcal/mol), STOCK1N-94672 (−97.86 kcal/mol), STOCK1N-89819 (−24.23 kcal/mol), and STOCK1N-80570 (−38.32 kcal/mol), and these natural compounds, particularly STOCK1N-62634, STOCK1N-96109, and STOCK1N-94672, can be considered as a potential lead for the inhibition of 3C protease (48, 56). From the above results, we believe that these natural compounds could be potential antiviral inhibitors for the 3C protease of FMDV.
Conclusion
Foot-and-mouth disease virus (FMDV) is a highly contagious virus that infects many species such as cattle, sheep, goats, and pigs and has become a significant health concern in livestock industries. Currently, no particular antiviral treatment is available for FMD, making exploring new inhibitors against FMDV critically important. Many countries adopted different vaccine strategies based on the prevailing FMDV strains, but it is hard to cure the infection, and it takes several days to elicit an immune response with the vaccine. Given the importance of 3C protease in viral replication and pathogenesis, we employed a structure-based virtual screening and MD simulation approach to explore potential inhibitors from a large natural compound dataset. The integration of computer-based drug discovery (CBDD) and MD simulations in drug discovery for viral diseases holds the potential to revolutionize the field by enabling more efficient, targeted, and personalized approaches to drug development. These technologies contribute to a better understanding of the underlying biology, accelerate the identification of potential drug candidates, and improve the overall success rate of drug discovery pipelines. The results identified five potential lead antiviral inhibitors, particularly STOCK1N-62634, STOCK1N-96109, and STOCK1N-94672, with effective drug-likeness properties, which can be further evaluated using clinical and experimental approaches for future management of FMDV infection.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.
Author contributions
SS: Conceptualization, Formal analysis, Investigation, Writing – original draft. H-KL: Funding acquisition, Supervision, Writing – review & editing. DS: Conceptualization, Formal analysis, Funding acquisition, Supervision, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was supported by the Ministry of Education of the Republic of Korea and the National Research Foundation of Korea (NRF-2022R1A2C4002510).
Acknowledgments
The authors thank Jeonbuk National University for providing necessary facilities to carry out this study. The authors greatly thank Center of Excellence in Green and Efficient Energy Technology (CoE-GEET), Central University of Jharkhand, Ranchi-835205, Jharkand, India for the help in DFT analysis. The authors also thank Mahesh Samantaray for the technical assistance and help in virtual screening analysis.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2023.1340126/full#supplementary-material
Abbreviations
FMD, foot-and-mouth disease; FMDV, foot-and-mouth disease virus; 3Cpro, 3C protease; HRV, human rhinovirus; ORF, open reading frame; SBVS, structure-based virtual screening; MD, molecular dynamics; SP, standard precision; XP, extra precision; RMSD, root mean square deviation; RMSF, root mean square fluctuation; Rg, radius of gyration; SASA, solvent accessible surface area; H-bond, hydrogen bond; PCA, principal component analysis; FEL, free energy landscape; MM-PBSA, molecular mechanics Poisson–Boltzmann surface area; PC, principal component; DFT, density functional theory; HOMO, highest occupied molecular orbital; LUMO, lowest unoccupied molecular orbital.
Footnotes
References
1. Jamal, SM, and Belsham, GJ. Foot-and-mouth disease: past, present and future. Vet Res. (2013) 44:116. doi: 10.1186/1297-9716-44-116
2. Knight-Jones, TJD, and Rushton, J. The economic impacts of foot and mouth disease – what are they, how big are they and where do they occur? Prev Vet Med. (2013) 112:161–73. doi: 10.1016/j.prevetmed.2013.07.013
3. King, AMQ, Adams, MJ, Carstens, EB, and Lefkowitz, EJ. Order – Picornavirales In: K AMQ, editor. Virus Taxonomy. San Diego: Elsevier (2012). 835–9.
4. Klein, J. Understanding the molecular epidemiology of foot-and-mouth-disease virus. Infect Genet Evol. (2009) 9:153–61. doi: 10.1016/j.meegid.2008.11.005
5. Wong, CL, Yong, CY, Ong, HK, Ho, KL, and Tan, WS. Advances in the diagnosis of foot-and-mouth disease. Front Vet Sci. (2020) 7:477. doi: 10.3389/fvets.2020.00477
6. Diaz-San Segundo, F, Medina, GN, Stenfeldt, C, Arzt, J, and De Los Santos, T. Foot-and-mouth disease vaccines. Vet Microbiol. (2017) 206:102–12. doi: 10.1016/j.vetmic.2016.12.018
7. Theerawatanasirikul, S, Thangthamniyom, N, Kuo, C-J, Semkum, P, Phecharat, N, Chankeeree, P, et al. Natural phytochemicals, Luteolin and Isoginkgetin, inhibit 3C protease and infection of FMDV. In Silico In Vitro Viruses. (2021) 13:2118. doi: 10.3390/v13112118
8. Rodríguez Pulido, M, and Sáiz, M. Molecular mechanisms of foot-and-mouth disease virus targeting the host antiviral response. Front Cell Infect Microbiol. (2017) 7:252. doi: 10.3389/fcimb.2017.00252
9. Curry, S, Roqué-Rosell, N, Sweeney, TR, Zunszain, PA, and Leatherbarrow, RJ. Structural analysis of foot-and-mouth disease virus 3C protease: a viable target for antiviral drugs? Biochem Soc Trans. (2007) 35:594–8. doi: 10.1042/BST0350594
10. Li, W, Ross-Smith, N, Proud, CG, and Belsham, GJ. Cleavage of translation initiation factor 4AI (eIF4AI) but not eIF4AII by foot-and-mouth disease virus 3C protease: identification of the eIF4AI cleavage site. FEBS Lett. (2001) 507:1–5. doi: 10.1016/S0014-5793(01)02885-X
11. Falk, MM, Grigera, PR, Bergmann, IE, Zibert, A, Multhaup, G, and Beck, E. Foot-and-mouth disease virus protease 3C induces specific proteolytic cleavage of host cell histone H3. J Virol. (1990) 64:748–56. doi: 10.1128/JVI.64.2.748-756.1990
12. Patick, AK, Binford, SL, Brothers, MA, Jackson, RL, Ford, CE, Diem, MD, et al. In vitro antiviral activity of AG7088, a potent inhibitor of human rhinovirus 3C protease. Antimicrob Agents Chemother. (1999) 43:2444–50. doi: 10.1128/aac.43.10.2444
13. Binford, SL, Maldonado, F, Brothers, MA, Weady, PT, Zalman, LS, Meador, JW, et al. Conservation of amino acids in human rhinovirus 3C protease correlates with broad-Spectrum antiviral activity of Rupintrivir, a novel human rhinovirus 3C protease inhibitor. Antimicrob Agents Chemother. (2005) 49:619–26. doi: 10.1128/aac.49.2.619-626.2005
14. Zitterl-Eglseer, K, and Marschik, T. Antiviral medicinal plants of veterinary importance: a literature review. Planta Med. (2020) 86:1058–72. doi: 10.1055/a-1224-6115
15. El Sayed, KA. Natural products as antiviral agents. Stud Nat Prod Chem. (2000) 24:473–572. doi: 10.1016/S1572-5995(00)80051-4
16. Musarra-Pizzo, M, Pennisi, R, Ben-Amor, I, Mandalari, G, and Sciortino, MT. Antiviral activity exerted by natural products against human viruses. Viruses. (2021) 13:828. doi: 10.3390/v13050828
17. Cavasotto, CN, Adler, NS, and Aucar, MG. Quantum chemical approaches in structure-based virtual screening and Lead optimization. Front Chem. (2018) 6:188. doi: 10.3389/fchem.2018.00188
18. Duay, SS, Yap, RCY, Gaitano, AL, Santos, JAA, and Macalino, SJY. Roles of virtual screening and molecular dynamics simulations in discovering and understanding antimalarial drugs. Int J Mol Sci. (2023) 24:9289. doi: 10.3390/ijms24119289
19. Samantaray, M, Pattabiraman, R, Murthy, TPK, Ramaswamy, A, Murahari, M, Krishna, S, et al. Structure-based virtual screening of natural compounds against wild and mutant (R1155K, A1156T and D1168A) NS3-4A protease of hepatitis C virus. J Biomol Struct Dyn. (2023) 1–18. doi: 10.1080/07391102.2023.2246583
20. Sethi, G, Hwang, JH, and Krishna, R. Structure based exploration of potential lead molecules against the extracellular cysteine protease (EcpA) of Staphylococcus epidermidis: a therapeutic halt. J Biomol Struct Dyn. (2023) 1–17. doi: 10.1080/07391102.2023.2250455
21. Pathak, RK, Seo, Y-J, and Kim, J-M. Structural insights into inhibition of PRRSV Nsp4 revealed by structure-based virtual screening, molecular dynamics, and MM-PBSA studies. J Biol Eng. (2022) 16:4. doi: 10.1186/s13036-022-00284-x
22. Berman, HM, Westbrook, J, Feng, Z, Gilliland, G, Bhat, TN, Weissig, H, et al. The Protein Data Bank. Nucleic Acids Res. (2000) 28:235–42. doi: 10.1093/nar/28.1.235
23. Jumper, J, Evans, R, Pritzel, A, Green, T, Figurnov, M, Ronneberger, O, et al. Highly accurate protein structure prediction with AlphaFold. Nature. (2021) 596:583–9. doi: 10.1038/s41586-021-03819-2
25. Lavecchia, A, and Di Giovanni, C. Virtual screening strategies in drug discovery: a critical review. Curr Med Chem. (2013) 20:2839–60. doi: 10.2174/09298673113209990001
26. Halgren, TA, Murphy, RB, Friesner, RA, Beard, HS, Frye, LL, Pollard, WT, et al. Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening. J Med Chem. (2004) 47:1750–9. doi: 10.1021/jm030644s
27. Repasky, MP, Shelley, M, and Friesner, RA. Flexible ligand docking with glide. Curr Protoc Bioinformatics. (2007) 18:8.12.1–8.12.36. doi: 10.1002/0471250953.bi0812s18
28. Friesner, RA, Murphy, RB, Repasky, MP, Frye, LL, Greenwood, JR, Halgren, TA, et al. Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein−ligand complexes. J Med Chem. (2006) 49:6177–96. doi: 10.1021/jm051256o
29. Curry, S, Roqué-Rosell, N, Zunszain, PA, and Leatherbarrow, RJ. Foot-and-mouth disease virus 3C protease: recent structural and functional insights into an antiviral target. Int J Biochem Cell Biol. (2007) 39:1–6. doi: 10.1016/j.biocel.2006.07.006
30. Zunszain, PA, Knox, SR, Sweeney, TR, Yang, J, Roqué-Rosell, N, Belsham, GJ, et al. Insights into cleavage Specificity from the crystal structure of foot-and-mouth disease virus 3C protease complexed with a peptide substrate. J Mol Biol. (2010) 395:375–89. doi: 10.1016/j.jmb.2009.10.048
31. Birtley, JR, Knox, SR, Jaulent, AM, Brick, P, Leatherbarrow, RJ, and Curry, S. Crystal structure of foot-and-mouth disease virus 3C protease: new insights into catalytic mechanism and cleavage specificity *. J Biol Chem. (2005) 280:11520–7. doi: 10.1074/jbc.M413254200
32. BIOVIA, Dassault Systèmes. Discovery studio visualizer, [v21.1.0.20298]. San Diego: Dassault Systèmes (2020).
34. Delley, B. DMol3 DFT studies: from molecules and molecular environments to surfaces and solids. Comput Mater Sci. (2000) 17:122–6. doi: 10.1016/S0927-0256(00)00008-2
35. Hirano, Y, Okimoto, N, Fujita, S, and Taiji, M. Molecular dynamics study of conformational changes of Tankyrase 2 binding subsites upon ligand binding. ACS Omega. (2021) 6:17609–20. doi: 10.1021/acsomega.1c02159
36. Bekker, H, Berendsen, H, Dijkstra, E, Achterop, S, Vondrumen, R, Vanderspoel, D, et al. GROMACS - a PARALLEL computer for molecular-dynamics simulations: 4th International Conference on Computational Physics (PC 92). Physics Computing’ 92, pp. 252–256. (1993).
37. Huang, J, and MacKerell, AD. CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data. J Comput Chem. (2013) 34:2135–45. doi: 10.1002/jcc.23354
38. Vanommeslaeghe, K, and MacKerell, AD. Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. J Chem Inf Model. (2012) 52:3144–54. doi: 10.1021/ci300363c
40. Darden, T, York, D, and Pedersen, L. Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. J Chem Phys. (1993) 98:10089–92. doi: 10.1063/1.464397
41. Hess, B. P-LINCS: a Parallel linear constraint solver for molecular simulation. J Chem Theory Comput. (2008) 4:116–22. doi: 10.1021/ct700200b
42. Hess, B, Bekker, H, Berendsen, HJC, and Fraaije, JGEM. LINCS: a linear constraint solver for molecular simulations. J Comput Chem. (1997) 18:1463–72. doi: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H
43. Bussi, G, Donadio, D, and Parrinello, M. Canonical sampling through velocity rescaling. J Chem Phys. (2007) 126:14101. doi: 10.1063/1.2408420
44. Parrinello, M, and Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. J Appl Phys. (1981) 52:7182–90. doi: 10.1063/1.328693
45. Sahoo, S, Samantaray, M, Jena, M, Gosu, V, Bhuyan, PP, Shin, D, et al. In vitro and in silico studies to explore potent antidiabetic inhibitor against human pancreatic alpha-amylase from the methanolic extract of the green microalga Chlorella vulgaris. J Biomol Struct Dyn. (2023) 1–11. doi: 10.1080/07391102.2023.2244592
46. Sahoo, S, Son, S, Lee, H-K, Lee, J-Y, Gosu, V, and Shin, D. Impact of nsSNPs in human AIM2 and IFI16 gene: a comprehensive in silico analysis. J Biomol Struct Dyn. (2023) 1–13. doi: 10.1080/07391102.2023.2206907
47. Genheden, S, and Ryde, U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discovery. (2015) 10:449–61. doi: 10.1517/17460441.2015.1032936
48. Poli, G, Granchi, C, Rizzolio, F, and Tuccinardi, T. Application of MM-PBSA methods in virtual screening. Molecules. (2020) 25:1971. doi: 10.3390/molecules25081971
49. Valdés-Tresanco, MS, Valdés-Tresanco, ME, Valiente, PA, and Moreno, E. gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS. J Chem Theory Comput. (2021) 17:6281–91. doi: 10.1021/acs.jctc.1c00645
50. Miller, BRI, TD, MG, Swails, JM, Homeyer, N, Gohlke, H, and Roitberg, AE. MMPBSA.Py: an efficient program for end-state free energy calculations. J Chem Theory Comput. (2012) 8:3314–21. doi: 10.1021/ct300418h
51. Uwakwe, KJ, Okafor, PC, Obike, AI, and Ikeuba, AI. Molecular dynamics simulations and quantum chemical calculations for the adsorption of some imidazoline derivatives on iron surface. Global J Pure Appl Sci. (2017) 23:69–80. doi: 10.4314/gjpas.v23i1.8
52. Molecular mechanisms of suppression of host innate immunity by foot-and-mouth disease virus. Veterinary. Vaccine. (2023) 2:100015. doi: 10.1016/j.vetvac.2023.100015
53. Arzt, J, Baxt, B, Grubman, MJ, Jackson, T, Juleff, N, Rhyan, J, et al. The pathogenesis of foot-and-mouth disease II: viral pathways in swine, small ruminants, and wildlife; Myotropism, chronic syndromes, and molecular virus–host interactions. Transbound Emerg Dis. (2011) 58:305–26. doi: 10.1111/j.1865-1682.2011.01236.x
54. Theerawatanasirikul, S, Lueangaramkul, V, Pantanam, A, Mana, N, Semkum, P, and Lekcharoensuk, P. Small molecules targeting 3C protease inhibit FMDV replication and exhibit Virucidal effect in cell-based assays. Viruses. (2023) 15:1887. doi: 10.3390/v15091887
55. Rozhenko, AB. Density functional theory calculations of enzyme–inhibitor interactions in medicinal chemistry and drug design In: L Gorb, V Kuz’min, and E Muratov, editors. Application of computational techniques in pharmacy and medicine. Challenges and advances in computational chemistry and Physics. Dordrecht: Springer Netherlands (2014). 207–40.
Keywords: virtual screening, MD simulation, FMDV, DFT, 3C protease
Citation: Sahoo S, Lee H-K and Shin D (2024) Structure-based virtual screening and molecular dynamics studies to explore potential natural inhibitors against 3C protease of foot-and-mouth disease virus. Front. Vet. Sci. 10:1340126. doi: 10.3389/fvets.2023.1340126
Edited by:
Rajesh Kumar Pathak, Chung-Ang University, Republic of KoreaReviewed by:
Jyoti Kant Choudhari, Maulana Azad National Institute of Technology, IndiaJoycee Jogi, Nanaji Deshmukh Veterinary Science University, India
Shubham Pant, Central Electrochemical Research Institute (CSIR), India
Copyright © 2024 Sahoo, Lee and Shin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hak-Kyo Lee, breedlee@jbnu.ac.kr; Donghyun Shin, sdh1214@gmail.com