Molecular docking and molecular dynamics simulation of glycyrrhizic acid in multitarget agents as potential inhibitors of respiratory influenza viruses

Retrieving herbal compounds and preparation of drug library

The 3D structures of GA from liquorice were retrieved from the PubChem database. We used the Structure File Generator (https://cactus.nci.nih.gov/translate/) and Online SMILES Translator to draw and convert compounds with missing.mol/.sdf files [8,9].

The way molecules interact in living things is not fully captured by molecular docking because it is based on the static structures of proteins and ligands. This simplification can affect the accuracy of predicted binding affinities and poses a challenge in interpreting the results as definitive.

Target protein retrieval and reparation

Based on the literature survey, we found that various targets, such as the ATP-bound state of BiP (5E84), main protease (Mpro) (6LU7), spike receptor binding domain (6LZG), RNA-dependent RNA polymerase (6M71), spike glycoprotein (6VSB), NSP15 endonuclease (6VWW), Nsp9 RNA-binding protein (6W4B), papain-like protease (6W9C), and neurominidase from the H1N1 (5NZ4) are major targets to study antiviral activity against respiratory influenza viruses, including SARS-CoV-1 and SARS-CoV-2 [10–14]. Therefore, we retrieved the fasta sequences of the aforementioned targets for antiviral activity against respiratory influenza viruses, such as SARS-CoV-1 and SARS-CoV-2, from the National Centre for Biotechnology information server. We also used the Basic Local Alignment Search tool to search for similar biological sequences on the Protein Data Bank. We then sorted the top 5 to 10 selected sequences based on their improved query coverage, percentage identity, and E-value. We used the PDB databank to obtain the three-dimensional X-ray crystallographic structure of various targets, such as the ATP-bound state of BiP (5E84), main protease (Mpro) (6LU7), spike receptor binding domain (6LZG), RNA-dependent RNA polymerase (6M71), spike glycoprotein (6VSB), NSP15 endonuclease (6VWW), Nsp9 RNA-binding protein (6W4B), papain-like protease (6W9C), and neurominidase from the H1N1 (5NZ4). We validated each target’s accession number using parameters including resolution, mutation, wwPDB validation, co-crystal ligand, and Ramachandran plot [15,16]. Table 1 presents a comparison between the standard values and the recovered protein, which serves as validation for the docking study protein.

Table 1. Comparison of standard values with recovered protein for validating docking study protein. [Click here to view]

The angles from a Ramachandran plot can be used to verify the solution to a crystal structure as well as determine the role of an amino acid in a secondary structure. It also assists with constraining structure prediction simulations and defining energy functions. It displays the protein structure regions that energetically allow backbone dihedral angles against amino acid residues. Table 2 below displays Ramachandran plots for several study targets. The PROCHECK web tool determined the three-dimensional geometry of the protein model, calculated the Ramachandran plot, and generated results for residues in various colored regions, namely red (favorite), yellow (additionally allowed), pale yellow (generously allowed), and whitish yellow (disallowed) [17,18].

Table 2. Ramachandran Plots Analysis of Protein Structures Using PROCHECK: comparative assessment for proteins (a) ATP-bound state of BiP (5E84), (b) Main protease (6LU7), (c) Spike receptor-binding domain(6LZG), (d) RNA-dependent RNA polymerase (6M71), (e) Spike glycoprotein (6VSB), (f) NSP15 Endoribo-nuclease (6VWW), (g) Papain-like protease (6W4B), (h) Papain-like protease (6W9C), and (i) Neuraminidase from H1N1 (5NZ4). [Click here to view]

Every molecule has the potential to be either a macromolecule or a micromolecule; however, it is necessary to optimize and minimize the latter before conducting a docking study. To confirm the binding pocket, we used the protein data bank’s PDB [sum server], a pictorial library of 3D structures for conventional inhibitor interactions with proteins. We completed the proteins with missing residues and synthesized side chains using CHIMAERA v1.16. Subsequently, they underwent optimization and minimization before being included in the docking study. We achieve the optimization of proteins by configuring 1,000 iterations of the Steepest Descent algorithm with a step size of 0.1 A0, followed by 100 iterations of the Conjugate Gradient algorithm with a size of 0.1 A0. All hydrogen atoms, including the ones with slower rates of addition, were included. Protonation statuses were assigned to the histidine residues. Additional fees were applied for both conventional (using the AMBER ff14SB force field) and unconventional (using the AM1-BCC force field) residues. The net charges for all non-standard entities were stabilized to enable the computation of their atomic partial charges using ANTECHAMBER charges. The protein was purified by removing any nonstandard residues, such as water molecules, cocrystal ligands, and superfluous chains, using Biovia Discovery Studio visualizer V21.1.0.20298 after optimization and minimization [19,20].

Grid generation

Auto-Dock Tools, Chimera, and Maestro were used for receptor grid identification. The workspace showed many targets, such as BiP (5E84) bound to ATP, Mpro (6LU7), the spike receptor binding domain (6LZG), RNA-dependent RNA polymerase (6M71), spike glycoprotein (6VSB), NSP15 endonuclease (6VWW), Nsp9 RNA-binding protein (6W4B), papain-like protease (6W9C), and neurominidase from the H1N1 protein (5NZ4). The grid’s volume was determined by utilizing the pocket’s dimensions [21]. The size of the enclosing box was minimized to ensure its compatibility with the protein’s active site and the predicted ligands for docking. Computer-based target identification primarily encompasses the identification of disease-related targets, the identification of binding sites, and the evaluation of drug ability. The ability of the ligand and target protein to interact can be ascertained using the binding site [20]. The grid parameters specified in Table 4 were used to determine the active sites in the docking studies. The problem of low utility or druggability commonly encountered in clinical trials can be effectively mitigated through early druggability assessment of proteins. Moreover, it is necessary to ascertain whether the targeted protein is suitable for therapeutic purposes. Finding binding sites and testing druggability are important for figuring out how proteins work, understanding how cells work, doing molecular docking, and coming up with smart drugs [22–24]. Table 3 below illustrates the active sites for binding and grid generation.

Table 3. Analysis of active site binding and grid generation parameters for docking studies. [Click here to view]

Table 4. Grid parameters utilized for determining active sites in docking studies. [Click here to view]

The size of the enclosing box was minimized to ensure its compatibility with the protein’s active site and the anticipated ligand for docking.

Molecular docking of target protein with ligands

Once the ligands and proteins were obtained, their structures were converted to the pdbqt format using an internal bash script created with AutoDock Tools 1.5.6 for ligands and ADFRsuit for proteins. This script allowed all the rotatable bonds of the ligands to rotate freely while considering the receptor as stiff [26]. For docking studies, we utilized the AutoDock Vina 1.2.3, with 0.375 A spacing between grid points. The grid box was precisely positioned at the active site of the enzyme with great accuracy, enabling the program to explore potential interaction sites between the ligands and the receptor. Alternative arrangements were deemed the standard [24,27]. The XYZ center has coordinates (X*Y*Z), and the grid box has dimensions of 20 * 20 * 20 A0. The CPU parameter was set to 23, the exhaustiveness parameter was set to 32, the number of modes parameter was set to 9, and the energy range parameter was set to 3. The redockings were executed using identical settings as the previous dockings.

Simulation set-up

We examined 6LZG’s dock complexes with standard Ribavirin and test drug GA using MD simulations on a Desmond 2020.3 computer [25]. We implemented this system in a 10 × 10 × 10-period boundary salvation box using the OPLS-2005 force field [29–32] and an explicit solvent model with SPC water molecules [33]. Neutralizing the charge with 0.15 M NaCl solutions created an environment akin to living organisms. The protein-ligand complexes underwent retraining after equilibration with an NVT ensemble for 10 ns. After the previous phase, an NPT ensemble performed a brief (12 ns) equilibration and minimization. The NPT ensemble was established using the Nose-Hoover chain coupling approach [34] with the same 27°C temperature, 1.0 ps relaxation period, and 1 bar of pressure employed throughout all of the simulations. We set the timestep to 2 fs. To regulate pressure, we used a barostat approach based on the Martyna–Tuckerman–Klein chain coupling system with a relaxation duration of 2 ps [35]. The particle mesh Ewald technique [36] was used to calculate long-range electrostatic interactions with a 9-coulomb interaction radius. The RESPA integrator was used to compute the bonded forces of each trajectory at a time step of 2 fs. The final manufacturing run time for 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA was 100 ns. The final production run was finished in 100 ns per unit. To check how consistent the MD was, we found the root mean square deviation (RMSD), the radius of gyration (Rg), the root mean square fluctuation (RMSF), and the hydrogen amount (H-bonds). The stability was monitored using simulations.

Binding free energy analysis

Molecular mechanics and generalized The Born surface area technique (MM-GBSA) was used to calculate the binding free energies of 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA complexes. In the last fifty frames of the simulation trajectory, the MM-GBSA binding free energy was estimated using the Python script thermal mmgsba.py and the VSGB solvation model, as well as the OPLS_2005 force field and a one-step sample size. To calculate the binding free energy of MM-GBSA (kcal/mol), the distinct energy modules of columbic, covalent, hydrogen bond, van der Waals, self-contact, lipophilic, and ligand and protein solvation were aggregated using the additivity concept. You may figure out how much Gbind costs by inserting the following numbers into the equation [34].

? G_bind = Δ G_MM + Δ G_Solv – Δ G_SA

where

ΔG_bind designates the binding free energy,

ΔG_MM designates the difference between the free energies of ligand-protein complexes and the total energies of protein and ligand in isolated form,

ΔG_Solv designates the difference in the GSA solvation energies of the ligand-receptor complex and the sum of the solvation energies of the receptor and the ligand in the unbound state,

ΔG_SA designates the difference in the surface area energies for the protein and the ligand.

The accuracy of MD simulations is dependent on the force field parameters used to model the molecular interactions. Although widely validated force fields are employed, they still represent approximations and may not capture all the nuances of real molecular behavior.

Target prediction

The GA was investigated for molecular interactions with antiviral targets (Table 6). It has the highest binding affinity shows hydrogen bonds with ten residues and shows Hydrophobic Interactions with one and two residues, respectively.

Table 6. The interaction and distance between inhibitors and critical amino acids of nine different targets in docking complexes were assessed using the LIGPLOT software. [Click here to view]

Stability of GA-protein complexes

Simulations based on MDs and experiments were employed to investigate the stability and convergence of the Cα-backbone of 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA protein-ligand complexes. Figure 2A represents the RMSD of three different proteins: 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA, over a time span of 100 nanoseconds (ns). RMSD is a common metric in structural biology and MDs, used to measure the average distance between the atoms of superimposed molecules. High RMSD values typically indicate significant conformational changes from a reference structure, while lower RMSD values suggest that the structure is relatively similar to the reference. 6LZG_GA exhibits an increase in RMSD till 75 ns, which then stabilizes till the end of the simulation with an average RMSD of 2.51 Å. This could indicate a significant structural change or rearrangement that then stabilizes. The RMSD values for the Apo state were found stable throughout the simulation with an average RMSD of 2.62 Å. While 6LZG_Ribavirin shows a gradual increase over time of 52 ns and then stabilizes till the end of 100 ns with an average RMSD of 2.47 Å, suggesting continuous conformational changes or deviations from the reference structure and then get stabilizes. All three Protein states show RMSD below 3 Å indicating stable states while comparing, the Apo state displays slightly higher fluctuations in RMSD values than standard Ribavirin and GA throughout the simulation. The stability and dynamics of each combination can be inferred from their RMSD trajectories. The 6LZG_Ribavirin, and 6LZG_GA, after an initial drift from the reference structure, seem to maintain a consistent structure. The Apo state might be the less stable among the three due to its dynamic nature, with significant conformational changes. The complexes 6LZG_Ribavirin, and 6LZG_GA demonstrate similar behaviors in terms of their structural deviations from a reference. 6LZG_Apo state undergoes consistent changes, while 6LZG_Ribavirin and 6LZG_GA shows late stabilization. Figure 2B presents the RMSF values of residues for three protein-ligand complexes: 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA. The RMSF is a measure of the average deviation of a set of atoms (in this case, protein residues) over time from a reference position, typically used in MDs to determine the flexibility or variability of specific residues in a protein structure. 6LZG_Apo displays variable fluctuations throughout the residue indices. Notably, there are sharp peaks around residue 116–120, indicating specific regions of high flexibility. Overall, 6LZG_Ribavirin has moderately high fluctuations with pronounced peaks at residues 116–121, 137–138, 272, and 408–410, suggesting these regions might be particularly dynamic or flexible when interacting with Ribavirin. The complex 6LZG_GA shows significant fluctuations, especially around residues 83 and 117–120. The higher the RMSF value for a specific residue, the greater its deviation from the reference, indicating increased flexibility or movement. From graph 1B, all three complexes show increased flexibility at 116–120, which is the region of the protein that is naturally more mobile because of the presence of loop regions. After critical observation, it has been found that fluctuations happen at residues in the loop regions which are not found in the binding site. Figure 2C illustrates the Rg for three complexes: 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA, plotted against simulation time in nanoseconds (ns). The Rg measures the distribution of a set of points in space, and in the context of MDs, it offers insights into the compactness or spatial distribution of a protein molecule over time. Range of RG values 6LZG_Apo, 6LZG_Ribavirin, and 6LZG_GA are 24.33–28.79, 24.41–33.34, and 24.46–28.32, respectively, which indicates protein with GA has a more compact structure than the other two states while protein with Ribavirin has less compact structure. A larger Rg generally implies a more expanded or less compact protein conformation, whereas a smaller Rg suggests a more compact or folded state. All three complexes demonstrate dynamic behavior throughout the simulation, as evidenced by changes in their Rg values. These fluctuations provide insights into the protein’s conformational changes, potentially influenced by the compounds’ binding. The results can be a foundation for understanding the structural effects of these compounds on 6LZG. Figure 2D represents the number of hydrogen bonds formed between the protein 6LZG and two different ligands (Ribavirin, and GA) over a time course of 100 ns. Hydrogen bonds play a crucial role in molecular recognition, protein-ligand interactions, and stabilizing the molecular complex. The strength and number of hydrogen bonds can be indicative of the ligand’s binding affinity to the protein. The number of hydrogen bonds in 6LZG_Ribavirin shows the range of 0–6 throughout the 100 ns simulation time. The relatively low count suggests either transient or weak bonding or perhaps the molecule interacts through other non-hydrogen bond mechanisms. 6LZG_GA complex reveals a higher and more consistent hydrogen bond count than the 6LZG_Ribavirin interaction. The bond count fluctuates mainly between 0 and 8, indicating a stronger and more stable interaction between 6LZG and Glycyrrhizic acid. The stability and duration of a protein-ligand complex can be gauged by the consistency and number of hydrogen bonds. A higher count suggests a potentially stronger interaction, which might result in higher binding affinity. The data provides valuable insights into the ligand binding behavior, which could be instrumental in understanding their potential therapeutic effects or drug design strategies.

Figure 2. MD simulation analysis of 100 ns trajectories of (A) Cα backbone of 6LZG_Apo (black), 6LZG_Ribavirin (red), and 6LZG_GA (blue), (B) RMSF of Cα backbone of 6LZG_Apo (black), 6LZG_Ribavirin (red), and 6LZG_GA (blue). (C) The radius of Rg of Cα backbone of 6LZG_Apo (black), 6LZG_Ribavirin (red), and 6LZG_GA (blue). (D) Formation of hydrogen bonds in 6LZG_Ribavirin (red), and 6LZG_GA (blue). [Click here to view]

Molecular mechanics generalized born surface area (MM-GBSA) calculations

The Table 7 showcases the binding free energy in the form of MM-GBSA components for protein-ligand complexes 6LZG_Ribavirin, and 6LZG_GA, as calculated by MD-simulation. The table represents the binding energies (in kcal/mol) of two ligands (Ribavirin and GA) to a protein (6LZG). The energy values are categorized based on various contributions such as lipid, van der Waals (vdW), Coulombic, hydrogen bond, solvent effect (SolvGB), and covalent energies. Overall binding energies (ΔG_bind) of 6LZG_GA have the most negative binding energy, suggesting the strongest binding affinity among the two ligands. Ribavirin has somewhat comparable energy, but weaker binder than GA. Lipophilic binding energies (ΔG_bindLipo) for all ligands show negative values, indicating favorable lipophilic interactions. Ribavirin exhibits the most significant lipophilic interaction. van der Waals binding energies (ΔG_bindvdW) shows positive values for both the ligands suggesting weaker or repulsive vdW interactions. Coulombic Binding Energies (ΔG_bindCoulomb) for all ligands demonstrate negative values, indicating favorable electrostatic interactions. GA shows the strongest Coulombic interaction. Hydrogen bonding energies (ΔG_bindHbond) for all ligands exhibit strong hydrogen bond contributions. GA shows the most significant hydrogen bonding energy. The positive values of the solvent effect (ΔG_bindSolvGB) indicate an unfavorable contribution from the solvent effect upon binding. Both ligands have the same unfavorable solvent contribution. Covalent binding energies (ΔG_bindCovalent) were found to be negative for all ligands indicating favorable covalent interactions However, GA has stronger covalent interactions than Ribavirin. GA seems to have the strongest overall binding affinity for 6LZG, based on its most negative overall binding energy. This strong binding can be attributed to its significant lipophilic, Coulombic, hydrogen bonding, and covalent interactions

Table 7. Binding free energy components for the 6LZG_Ribavirin, and 6LZG_GA, calculated by MM-GBSA. [Click here to view]