Identification of potential plant bioactive as SARS-CoV-2 Spike protein and human ACE2 fusion inhibitors.

The Spike receptor binding domain (S-RBD) from SARS-CoV-2, a crucial protein for the entrance of the virus into target cells is known to cause infection by binding to a cell surface protein. Hence, reckoning therapeutics for the S-RBD of SARS-CoV-2 may address a significant way to target viral entry into the host cells. Herein, through in-silico approaches (Molecular docking, molecular dynamics (MD) simulations, and end-state thermodynamics), we aimed to screen natural molecules from different plants for their ability to inhibit S-RBD of SARS-CoV-2. We prioritized the best interacting molecules (Diacetylcurcumin and Dicaffeoylquinic acid) by analysis of protein-ligand interactions and subjected them for long-term MD simulations. We found that Dicaffeoylquinic acid interacted prominently with essential residues (Lys417, Gln493, Tyr489, Phe456, Tyr473, and Glu484) of S-RBD. These residues are involved in interactions between S-RBD and ACE2 and could inhibit the viral entry into the host cells. The in-silico analyses indicated that Dicaffeoylquinic acid and Diacetylcurcumin might have the potential to act as inhibitors of SARS-CoV-2 S-RBD. The present study warrants further in-vitro and in-vivo studies of Dicaffeoylquinic acid and Diacetylcurcumin for validation and acceptance of their inhibitory potential against S-RBD of SARS-CoV-2.

Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel coronavirus (CoV) leading to CoV disease (COVID-19). SARS-CoV-2 has also been determined to have an enhanced transmission among the human population in comparison to other CoVs [1,2]⁠ Architecturally, the CoV genome is enclosed in a helical capsid, surrounded by a protective envelope [3]⁠. The Spike (S) protein present on the protective envelope has a critical function in the entry of the virus into the host cells. It binds with the angiotensin converting enzyme 2 (ACE2) receptors of the host, giving rise to a cascade of events which ultimately lead to the fusion of the viral and the host cell membranes, thereby ensuring successful viral invasion [4]⁠. Primarily, ACE2 is a membrane protein found in heart, kidneys, lungs, and intestine which assists in the maturation process of angiotensin (peptide hormone) required for controlling blood pressure by vasoconstriction, thus any hindrance in its functions leads to cardiovascular dystrophies [5].

S protein is a densely glycosylated (comprises N-linked glycans) homo-trimeric protein that protrudes from the outer protective envelope of the virus with each monomer comprising of S1 and S2 subunits. The S1 subunit performs the cellular receptor recognition step and attaches itself to the host cell while S2 contains the fusion machinery required to perform the membrane fusion step [[6], [7], [8], [9]]⁠. The S1 subunit further comprises the receptor binding domain (RBD) which interact directly with the peptidase domain (PD) of the ACE2 receptor via hinge like conformational movements [7,8,10]⁠. Prior to the membrane fusion step, the S1 and S2 subunits are non-covalently bonded to each other and the S2 subunit is stabilized by the S1 subunit in the prefusion conformation [8,10]⁠. Once the S protein attaches to the cellular receptors via S1 RBD, the host proteases act on S and cleaves it at the boundary of S1 and S2 upstream of the location of the fusion peptide [4,10]⁠. This cleavage prompts the S2 to transit form the metastable prefusion state to an activated post fusion state through dramatic and irreversible rearrangements in its conformation, required for membrane fusion [8],⁠ eventually shedding the S1 subunit [4]⁠. The residues Thr333-Gly526 of S-RBD of SARS-CoV-2 and ACE2 residues Ser19-Asp615 were identified to interact with each other [9]⁠. Moreover the S-RBD of S1 subunit bind to ACE2 with a very high affinity implying that it is crucial component of S1 required for interaction with host [9,10].

The natural molecules and their derivatives were stated to hold a broad array of biological actions, including neuro-protective, anti-inflammatory, cardio-protective, anti-cancer, anti-microbial, anti-viral, and anti-oxidant effects [[11], [12], [13], [14]]⁠. Recent articles recommended that the natural compounds and their metabolites in various plants showed striking potential to inhibit different SARS-CoV-2 proteins [[15], [16], [17], [18], [19]]⁠. It is challenging to develop new drugs with all clinical trials in this emergency situation. Hence, common use, cost-effectiveness and low risk of side effects of phytochemicals make them prime candidates for antiviral drug discovery. In this study, bioactive molecules from diverse plants (used in in-house developed immunoboosting drink) were selected to find active ingredients for the prevention of COVID-19. The main objectives were to find highly potential drug candidates as fusion inhibitors of S-RBD and ACE2. Next, to verify the docking outcomes and in-depth study of receptor-ligand complexes by conducting (250 ns) long-term molecular dynamics simulations. Finally, employing Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) evaluations to recognize the most effective molecule by measuring the thermodynamic binding free energy.

Materials and methods

Datasets

X-ray crystallographic co-ordinates for the tertiary structure of the S-RBD (PDB id: 6M0J) [9]⁠ of SARS-CoV-2 were downloaded from the Protein Data Bank. In this crystal structure, only S1 subunit of the S glycoprotein bound to the ACE2 was diffracted at a resolution of 2.45 Å [9]⁠. The molecules employed in this analysis comprised of a wide range of bioactive (33) molecules (from plants used in in-house developed immunoboosting drink) that were collected from Pubchem database [20]⁠ which were also reported in earlier studies [13,[21], [22], [23]]⁠. These molecules were listed in Table S1 with their plant of origin. Prior to docking, the structure of the protein is prepared via the “prepare protein” methodology [24]⁠. Subsequently, the ligand molecules were geometrically optimized via the minimization protocols of Gaussian16's DFT [25].

Protein-protein molecular dynamics simulation

The crystal structure obtained from PDB comprised of the S-RBD of SARS-CoV-2 bound with the human ACE2 protein. This structure was subjected to protein-protein molecular dynamics (MD) simulation of 125 ns to identify the binding site residues that are required for interaction between the two proteins. The simulation was performed by GROMACS under the GROMOS96 43a1 force field [26,27]⁠. Here, the two proteins in the structure were prepared via the pdb2gmx module and then set as the initiation points of the simulation. The entire simulation was performed inside a simple point charge/extended (SPC/E) water model solvated simulation box [28]⁠. This box was stabilized by the add-on of counter ions via the “gmx genion” script of GROMACS [29,30]⁠. Furthermore, steepest descent method was used for energy minimization of the system followed by NPT and NVT equilibration simulations. During the simulation the two proteins interacted with each other forming new long range and short range interactions. These interactions are restrained via the LINCS algorithm [31]⁠. For the complete simulation experiment, the temperature was maintained at 300K with the help of Berendsen Thermostat [32],⁠ whereas the pressure was maintained at 1 bar with the help of Parrinello-Rahman barostat [33]⁠. After completion of MD run, end-state thermodynamics free energy was evaluated by MM-PBSA approach [34]⁠. The binding free energy was further decomposed into per residue contribution energy to identify residues with most favorable energies. The identified residues of S-RBD involved in interactions with ACE2 determined through protein-protein MD simulations were accepted as the ligand-binding site for docking experiments.

Molecular docking

Next we characterized the interaction pattern followed by the S-RBD with the 33 chosen bioactive molecules through molecular docking experiments. This experiment was performed on the Accelrys Discovery Studio 2018 software [35]⁠ that utilizes the CDOCKER algorithm to meet this purpose. In this protocol of docking, conformation of the protein receptor is kept stagnant while the ligand molecule is movable and flexible so that it is able to adapt the conformational changes occurring in the receptor. This is stated as flexible-ligand and protein-fixed docking. All the CDOCKER parameters for this experiments were kept default.

As there was no ligand binding site defined in the crystal structure, so the residues identified from the formerly performed protein-protein MD simulation were selected as the ligand binding site. Since, the residues of S-RBD (Thr333, Arg346, Arg355, Lys356, Arg357, Lys378, Lys386, Arg403, Arg408, Lys417, Lys424, Lys444, Arg454, Arg457, Lys458, Lys462, Arg466, and Arg509) were involved in interaction with the ACE2 protein of humans. The shape of this region is crucial for the accurate simulations of interactions between the protein-ligand. Given that CDOCKER is a grid-based approach [36,37],⁠ a grid box with a spacing of 0.50 and a grid angle of 90° was also defined. The X, Y, and Z coordinates of the binding site within the grid box are (−33.550696, 38.344658, 5.080347) bearing a radius of 12.98 Å. The molecules were then docked onto the chosen ligand binding site and favorable poses for the protein-ligand complexes were generated. High temperature kinetics was administered during the docking process to ensure that the ligand molecules occupy a flexible conformational space in the protein, and simulated annealing verified the docking accuracy [36]⁠. Furthermore, energy minimization of the docked complexes was handled by the CHARMm force field minimizing their total energy to 0.01 kcal/mol/Å thereby stabilizing the conformations [37]⁠. CHARMm also refines the complexes by removing any hetero-atoms and bound inhibitors present in them, that aren't requisite for analysis as well as inspecting them for any structural issues [36,37]⁠. Pose generation was proceeded by computation of CDOCKER interaction energy that was used as a parameter to rank the generated poses, where a pose having the lowest interaction energy was considered to have the most favorable association between the protein and ligand [36]⁠. Top 2 molecules with best CDOCKER interaction energy scores were selected for molecular dynamics simulation analysis. Interactions formed in a protein-ligand complex were visualized via PLIP server [38]⁠.

Protein-ligand molecular dynamics simulation and binding free energy calculation

Lastly, to analyze the strength of the interactions between the chosen protein-ligand complexes as well as to reaffirm their overall conformational stability. We have taken the complexes to undergo protein-ligand MD simulation analysis of 250 ns. This simulation was executed by the GROMOS96 43a1 force field [26,27],⁠ where the ligand topology development was carried out by the PRODRG server [39]⁠ and the protein topology development was carried out by pdb2gmx module of GROMACS. Complexes were kept inside the statistical process control water model (solvated simulation box), stabilized by the add-on of counter ions via the “gmx genion” script of GROMACS. Furthermore, steepest descent method was used for energy minimization of the system followed by NPT and NVT equilibration simulations. LINCS and SETTLE algorithms were applied to restrain the bond length and water molecules in the complexes respectively [40,41]⁠. Particle Mesh Ewald (PME) approach was used for the handling of the long range electrostatic interaction, where the Fourier grid spacing was kept 0.15 [42]⁠. A time step of 2fs was allowed via the Shake algorithm, which constrains the H-bond lengths, whereas the coulombic and van der Waal interactions were clipped at 1.0 nm. For the complete simulation experiment, the temperature was maintained at 300K with the help of Berendsen Thermostat [32], whereas the pressure was maintained at 1 bar with the help of Parrinello-Rahman barostat [33]. In-built GROMACS scripts were used to analyze the results and produce trajectories of root mean square deviation (RMSD) through “gmx_rms” script and Hbond through “gmx_hbond” script [29]⁠. Snapshots of the simulated complexes were saved using Coot package [43]⁠. Dynamically stable complexes topologies obtained from simulation [44]⁠ were further put through Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations by applying the g_mmpbsa script [30,34]⁠, which provided the binding free energy trajectories for the selected complexes.

Results and discussion

Identification of binding site

The SARS-CoV-2 utilizes the S-RBD to interact with the human cell receptor ACE2. Thus, we conducted protein-protein MD simulation study of 125 ns to find out crucial residues involved in the binding of S-RBD and ACE2. The binding free energy for S-RBD and ACE2 was −918.744 kJ/mol. This energy was further decomposed to per residue contribution energy to identify crucial residues responsible for interaction between S-RBD and ACE2. Based on lower per residue contribution energy, S-RBD interface residues (Table S2) were selected to define the binding site for molecular docking study. The binding site was defined with a 12.98 Å radius containing all the residues selected from per residue contribution energy.

Binding interaction analysis

We adopted molecular docking, the most commonly accepted method in structure-based drug design to examine the nature of binding affinity and interaction between selected compounds and target protein. This method can predict the ligand interaction within the target binding site with a significant degree of precision. A library of bioactive molecules from different medicinal plants used in in-house developed immunoboosting drink were used to predict their binding capability with S-RBD of SARS-CoV-2. Based on best interaction energy, Diacetylcurcumin and Dicaffeoylquinic acid were selected for further analysis. The interaction energies of these molecules were presented in Table S3. The docking poses of top two molecules were shown in Fig. 1 , while the 2-D interaction poses of the following three molecules (Pipercide, Guineensine, and Demethoxycurcumin) were provided in Fig. S1. In the docked conformations, Diacetylcurcumin formed one hydrogen bond and salt bridge with residue Lys417, π stacking with residue Tyr489, and a hydrophobic interaction with residue Phe456. On the other hand, Dicaffeoylquinic acid formed seven hydrogen bonds with the residues Lys417, Glu484, Arg457, Tyr473, Gln493, and Lys458. Additionally, residue Lys417 also formed a salt bridge. The residues Phe456 and Tyr489 were also involved in double hydrophobic interactions. The visual analysis of docked complexes revealed that the Dicaffeoylquinic acid displayed maximum number of hydrogen bond with the residues responsible for significant role in binding. Previously published studies reported that the residues Lys417, Phe456, Gln493, Glu484, Tyr473, and Tyr489 of S-RBD interface majorly contributed in binding with ACE2 [7,9,[45], [46], [47], [48], [49], [50], [51], [52]]⁠. Moreover, in various in-silico studies, many phytochemicals (Bis-demethoxycurcumin, compound-4, compound-2, hesperidin, emodin, chrysin, nabiximols, pectolinarin, epigallocatechin gallate, rhoifolin, andrographolide, artemisinin, apigenin, berberine, emodin, capsaicin, glabridin, colchicine, harmaline, kaempferol, harmine, niranthin, phyllanthin, oleic acid, rosavin, withanolide A, withaferin A, curcumin, apigenin, and chrysophanol) were also recommended as S-RBD inhibitors on the basis of molecular docking studies [[53], [54], [55], [56]]⁠.

Fig. 1

3D interactions of S-RBD of SARS-CoV-2 protein with bioactive molecules (a) Diacetylcurcumin, (b) Dicaffeoylquinic acid.

Analysis of post-dynamics trajectories

A crucial but often overlooked perspective to be acknowledged in molecular docking is the flexibility of the target binding site [57]⁠. The protein may experience conformational fluctuations during protein-ligand interactions. Structural changes caused by the dynamic activity of a protein's structural components can also affect the biological functioning [57]⁠. Several obscure biological roles of proteins and their intense dynamic mechanisms can be discovered by observing their internal changes [17]. Molecular docking, although effective, only presents a static pose of protein-ligand interaction. The structural perturbations during protein-ligand interactions can be dealt with robust MD simulations. Accordingly, 250 ns of MD simulations were conducted for two best-docked complexes. MD simulations were used to analyze the dynamic properties of the selected complexes required for the inhibition mechanism and structural changes.

The dynamic function of whole simulated systems was analyzed in-depth, by employing various MD-driven parameters, such as RMSD, ensemble clustering, and H-bond analysis. To validate the interaction poses generated by molecular docking, RMSD of backbone C-α atoms of both the selected complexes was calculated. The RMSD is an established measure of the conformational stability of the protein-ligand complexes [58]⁠. RMSDs of both the complexes were interpreted in a graphical mode to analyze the stability of protein-ligand complexes (Fig. 2 ). The average RMSD value for complexes with Diacetylcurcumin and Dicaffeoylquinic acid were ~0.37 nm and ~0.5 nm respectively. The complex with Diacetylcurcumin showed a gradual increase in RMSD value for the first ~30 ns, while the RMSD trajectory of Dicaffeoylquinic acid increased rapidly till ~60 ns. Afterwards, both the RMSD trajectories equilibrated and followed stable trajectories till the end of simulation. The protein-ligand complexes containing potential phytochemicals suggested by recent in-silico studies [[59], [60], [61], [62], [63], [64]]⁠ also showed average RMSD values within ~0.5 nm. The low RMSD values of complexes validated the stability of molecular docking poses. The low degree of fluctuations and small variations in the average RMSD values directed equilibration of trajectories, suitable for further computational interpretations.

Fig. 2

RMSD of backbone C-α atoms of S-RBD of SARS-CoV-2 protein complexes with bioactive molecules.

The RMSD clustering simulation strategy is an adequate and dynamic process of examining the structural flexibility of proteins [65]⁠. The well equilibrated (last 5ns) simulation trajectories of the complexes were collected and subjected for the interpretation of central conformations obtained from the average structure of each complex during the sampling time. We clustered the examined configurations trajectories based on the Daura's approach [66]⁠. The atomic motion was analyzed and adequately described the conformational changes during the simulation (Fig. 3 ). In complex with S-RBD, Diacetylcurcumin formed 51 clusters with an average RMSD of 0.19 nm, while Dicaffeoylquinic acid formed 51 clusters with an average RMSD of 0.18 nm. The binding of selected molecules had minimal effect on the secondary structures of S-RBD. All the helices and sheets showed minimal fluctuations and retained native conformations throughout the sampling period. Moreover, the selected molecules displayed intermediate structural fluctuations within the binding pocket by completely occupying it; hence possibly ceasing its interaction with ACE2. Further, to access the stability of protein-ligand interactions, we carried out analysis of hydrogen bonds between the selected molecules and S-RBD of SARS-CoV-2.

Fig. 3

Pictorial representation of conformational flexibility using cluster analysis for S-RBD in complex with (a) Dicaffeoylquinic acid, (b) Diacetylcurcumin.

* Number of clusters.

Hydrogen bonding is one of the main element in biological systems responsible for molecular interactions. Hydrogen bonds render the basis for molecular recognition by allowing explicitness and directionality to molecular interactions [67,68]⁠. We analyzed the formation of hydrogen bonds for selected complexes during the entire run of MD simulations (Fig. 4 ). The molecule Dicaffeoylquinic acid formed eight average hydrogen bonds, while Diacetylcurcumin formed an average of four hydrogen bonds with S-RBD of SARS-CoV-2. The formation of hydrogen bonds between protein and ligands continuously persisted throughout the simulation with minor fluctuations. These results confirmed that the molecules formed a significant number of hydrogen bonds with S-RBD during the simulations and showed the stability of protein-ligand interactions. Interactions of inter-molecular hydrogen bonds limits the internal movement of the binding site and enhances its stability [69]⁠. The bioactive molecules were capable of maintaining strong interactions with the binding pocket of S-RBD during the simulation period. These results showed the potential of both the selected molecules to be developed as fusion inhibitors to tackle SARS-CoV-2.

Fig. 4

Hydrogen bond profiles of S-RBD of SARS-CoV-2 protein complexes with bioactive molecules.

Thermodynamic binding free energy analysis

The MM-PBSA is an advanced and broadly adopted method that couples the continuum solvent forms and molecular mechanics to measure binding free energy between two molecules [34]⁠. From the equilibrated MD simulation trajectories, all the conformations of the complexes were inspected using the MM-PBSA method. The binding free energies for the Dicaffeoylquinic acid and Diacetylcurcumin were −193.74 kJ/mol and −187.60 kJ/mol, as depicted in Fig. 5 . The binding free energy suggested that Dicaffeoylquinic acid binds more tightly to the S-RBD of SARS-CoV-2 than Diacetylcurcumin. The electrostatic and van der Waals energies were the major driving energies behind binding free energy of Dicaffeoylquinic acid and S-RBD of SARS-CoV-2. In contrast, the polar solvation energy negatively contributed to the binding free energy, as displayed in Table 1 . In the already published literature, bioactive and their derivatives (O-Demethyldemethoxycurcumin, Tellimagrandin-II, 2,3‐Dihydrowithaferin A, and Asparoside D) in complex with S-RBD bestowed −12.48, −32.08, −87.60, and −66.49 kJ/mol binding free energies [60,61,70]⁠. However, our suggested molecules (Diacetylcurcumin and Dicaffeoylquinic acid) with S-RBD showed more favorable (−193.74 and −187.60 kJ/mol) binding free energies. The MM-PBSA analysis stated that the Dicaffeoylquinic acid and Diacetylcurcumin could be adopted as a lead molecules to design inhibitors against S-RBD of SARS-CoV-2.

Fig. 5

Graphical representation of the Delta_E_Binding free energy in kJ/mol showing S-RBD of SARS-CoV-2 protein complexes with bioactive molecules, Dicaffeoylquinic acid (green), Diacetylcurcumin (red).

Table 1

MM-PBSA calculations of all components of binding free energy for both the selected complexes.

S. No.	Bioactive molecules	Binding energy (kJ/mol)	Polar solvation (kJ/mol)	Electrostatic (kJ/mol)	SASA (kJ/mol)	Van der Waal (kJ/mol)
1.	Diacetylcurcumin	−187.6	89.79	−27.94	−17.81	−231.65
2.	Dicaffeoylquinic acid	−193.75	179.55	−170.6	−17.07	−185.63

Conclusion

The present study aimed to discover small bioactive molecules that could inhibit SARS-CoV-2 by binding to the S-RBD and prohibiting the virus from penetrating into the host cells. A potent binding site was identified by per residue contribution energy obtained by protein-protein MD simulations. Bioactive molecules from different plants were docked on the predicted binding site. The top two molecules (Diacetylcurcumin and Dicaffeoylquinic acid) with best interaction energies were selected for further computational analysis. To validate our selection, the selected complexes were subjected to long-term MD simulations. The MD-driven data suggested that these potential molecules obstructed the S-RBD-ACE2 interactions by tightly adhering to the binding site. Our results assured the testing of Dicaffeoylquinic acid and Diacetylcurcumin by in-vitro and in-vivo examinations against S-RBD of SARS-CoV-2 to inhibit the viral entry into host cells.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Declaration of competing interest

All authors hereby declare that they have no competing interests related to this work.