Antiviral peptides against the main protease of SARS-CoV-2: A molecular docking and dynamics study

2.2. Protein preparation

The three-dimensional structure of the SARS-CoV-2 Mpro (PDB ID: 6LU7) was retrieved from the Protein Data Bank (Protein Data Bank, 2019) with a resolution of 2.16 Å. The protein structure was initially prepared in Discovery Studio (San Diego: Accelrys Software Inc., 2012), where the water molecules and heteroatoms were removed. The cleaned protein structure was further minimized in YASARA tools by employing the AMBER14 force field (Land and Humble, 2018). The minimized protein structure was used for further molecular docking and dynamics studies.

2.3. Molecular docking

Peptide and protein docking was implemented as described in our previous study (Mahmud et al., 2021a). The protein structure of Mpro was input as the receptor molecule, and the peptides were used as the ligands in the PatchDock (Schneidman-Duhovny et al., 2005) webserver. The final clustering was selected on the basis of the root-mean-square deviation (RMSD) value. The top ten solutions were taken from PatchDock and further docked using FireDock (Mashiach et al., 2008) tools. The FireDock program optimizes lateral string conformations and rigid body formations to provide larger refinement. These web tools were employed to obtain rigid protein–peptide docking. The best protein–peptide complex was selected from among the top ten conformers based on energy scoring. All the 104 peptides were then further docked using the ClusPro program (Comeau et al., 2004) to obtain more accurate binding energy profiles from the corresponding peptide and protein complexes. Based on the combination of multiple docking programs, four peptide–protein complexes were selected for non-bonded structure analysis. The interaction analysis was conducted using PyMOL (DeLano, 2002) and Discovery Studio (San Diego: Accelrys Software Inc., 2012) software package.

2.4. Dynamics simulation

The docked peptide–protein molecules and their conformational variability were assessed through molecular dynamics simulations in YASARA dynamics software (Krieger et al., 2004). The AMBER14 force field (Maier et al., 2015) was used in this study, and the docked complexes were initially cleaned, along with hydrogen bond orientation and optimization for simulation. The cubic simulation box with a cell size of 110 × 110 × 110 Å³ was filled with water molecules, with a 0.9899 g/cm³  density, and the system was neutralized with 0.9% NaCl (Krieger et al., 2013). The TIP3P water model was used to solvate the complex. The acid dissociation constant value (pKa) was calculated for the amino acids in the complex. The SCWRL algorithm combined with hydrogen bond network optimization was applied to maintain the correct protonation state of each amino acid residue. The particle mesh Ewald method (Darden et al., 1993) was applied to calculate long-range electrostatic interactions, using a cutoff radius of 8 Å. The simulation cell box was 20 Å larger than the peptide–protein complex to allow for maximum free motion, and the periodic boundary condition was maintained (Krieger and Vriend, 2015). The energy of each simulated system was minimized by using the steepest descent gradient approaches (5000 cycles) via simulated annealing methods. An 80-ps position-restrained molecular dynamics simulation was performed for each system at constant pressure (1 atm) and temperature (310 K). Temperature coupling was performed using the Berendsen thermostat, with a temperature coupling constant of 0.1 ps, and a manometer method was used for pressure coupling, with a reference pressure of 1 atm. For the simulation study, a time step of 1.25 fs was used. The chemical bond lengths involving hydrogen bond atoms were fixed using the SHAKE algorithms (Brooks et al., 1983). After reaching an equilibrium state at 1 ns, the simulation was allowed to run for 250 ns, and the simulation trajectories were saved after every 100 ps. The simulation was conducted three times, and the average values from the trajectories were used to analyze the RMSD, root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), and hydrogen bonds (Bappy et al., 2020, Islam et al., 2020b, Khan et al., 2020, Swargiary et al., 2020).

Furthermore, the simulation snapshots were subjected to binding free energy calculations from MM-PBSA approaches by following equations (Mitra and Dash, 2018).

BindingEnergy=EpotRecept+EsolvRecept+EpotLigand+EsolvLigand-EpotComplex-EsolvComplexThe YASARA macro was used to calculate the MM-PBSA binding energy where positive energy indicates the better bindings (Srinivasan and Rajasekaran, 2016).

2.5. Physiochemical properties

The physiochemical properties of the four best peptides were evaluated using ProtParam (Gasteiger et al., 2005) tools. This tool provides relevant properties, including molecular weight, net charge at pH 7, volume, peptide properties, stability, and charge.

3.2. Molecular dynamics

The binding interaction and the conformational variations that occur during binding can be understood using a docking study to provide insight into the binding sites to a limited extent. Further changes in the ligand or the bound complex can be studied by varying the temperature and pressure in a molecular dynamics simulation, which can minimize the cost of performing such experiments in the laboratory (Arfin et al., 2018, Durrant and McCammon, 2011, Nair and Miners, 2014). A combined docking and molecular dynamics approach can be used to validate the docking-derived results (Gioia et al., 2017).

Molecular dynamics simulations were implemented for the peptide and SARS-CoV-2 Mpro complexes to confirm the structural rigidity and validate the docking outcomes for the complexes. The RMSD values of the C-alpha atoms were explored to understand the structural rigidity. Fig. 3 (a) demonstrated that the DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804 complexes showed initial RMSD increases due to instability. The increase in RMSD was higher for the DRAMP00877 peptide complex than for the other complexes. The upward trend in RMSD was maintained until 40 ns, after which all four complexes achieved stability, which was maintained throughout the remainder of the simulation period. Although the DRAMP00877 complex had a higher RMSD trend during the initial phase, the RMSD profile decreased from 40 to 100 ns. The DRAMP02333 complex had a lower RMSD than the other complexes, indicating increased structural integrity. The average RMSD profile of the complexes formed with DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804 were 1.88 Å, 1.58 Å, 1.82 Å, and 1.81 Å, respectively. The overall RMSD trend did not exceed 2.5 Å, which indicates overall structural stability as higher RMSD profiles for C-alpha atoms are indicators of low stability (Mahmud et al., 2021b).

Open in a separate window

Fig. 3

The molecular dynamics simulation of the peptide and main protease complex, here (a) root mean square deviation of the c-alpha atoms, (b) solvent accessible surface area, (c) radius of gyration, (d) hydrogen bonding of the complexes, (e) root mean square fluctuation of the complexes to understand the flexibility of the amino acid residue.

The SASA values of the complexes were analyzed to evaluate changes in the surface of the SARS-CoV-2 Mpro in response to binding with the peptide molecules. The DRAMP00877, DRAMP02333, and DRAMP03804 complex had similar SASA trends throughout the entire simulation trajectory, indicating no change in protein volume. However, the DRAMP02669 complex showed an initial upward trend in SASA until 50 ns, indicating an expansion of the surface area following initial complex formation [Fig. 3 (b)]. The complex then stabilized and showed a trend similar to that of the DRAMP02333 peptide complex. The average SASA profiles of the four complexes were 15663.42 Å², 16302.79 Å², 16066.2Ų and 15590.45 Å² for DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804, respectively. The SASA profiles of the peptide–protein complexes had similarly stable trends, with little fluctuation, indicating a lack of expansion or contraction for the protein complexes during the simulation time (Rakib et al., 2021).

The compactness of the protein complexes was evaluated using Rg, for which higher values indicate a more labile nature and a lower value indicates the firmness of the protein. The DRAMP02669 complex displayed fluctuations in the Rg value until 20 ns, indicating an initially loose packaging system for the protein [Fig. 3 (c)]. However, this complex maintained stability for the remainder of the simulation time. The DRAMP03804 complex had a larger Rg trend than all other complexes but was associated with a lower degree of deviation, demonstrating the achievement of a steady state. The comparative higher Rg of this peptide may be responsible for relatively more flexible state and mobile nature of the complexes. The average Rg values were 22.504, 22.527, 22.501, and 22.867 Å for DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804, respectively. The higher Rg profile of the DRAMP03804 complex than other complexes may defines the more mobile nature and folding and unfolding of the protein complexes. In addition, hydrogen bond formation between the protein and peptide in complex plays a vital role in maintaining complex stability. Fig. 3 (d) demonstrated that the DRAMP02333 complex contained the highest number of hydrogen bonds throughout the entire simulation time. All of the evaluated peptide and Mpro complexes were stable, and high levels of fluctuation were generally absent for these complexes.

The RMSF of the peptide and Mpro complexes were explored to understand the flexibility across the amino acid region of Mpro. Fig. 3 (e) demonstrated that almost all residues had low RMSF, except Ser1 (helix-strand), Glu47 (beta-turn), Asp48 (beta-turn), Arg60 (helix-strand), Lys97 (beta-turn), Lys137 (beta-turn), Arg222 (beta-turn), Phe294 (helix-strand), Ser301 (beta-turn), Gly302 (beta-turn), Val303 (beta-turn), Thr304 (beta-turn), Phe305 (beta-turn), and Gln306 (beta-turn).

The pre- and post-molecular dynamics docked structures were superimposed to explore variations among these structures. The RMSD of the DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804 peptides were found to be 2.076, 1.658. 1.793, and 1.778 Å, respectively, indicating limited flexibility among the structures (Fig. 4 ). Moreover, the binding free energy of the peptide-protein complexes were calculated from the MM-PBSA method (Fig. 5 ) where more positive score indicates the better binding of the complexes. The average free energy of the DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804 peptide complexes were −245.612, −243.441, 19.51237, −51.5827KJ/mol respectively. The DRAMP02669 peptide complex exhibited more favorable binding than the other complexes.

Open in a separate window

Fig. 4

The superimposition of the Pre and Post MD structure; (A) Pre and Post MD structure of DRAMP00877 peptide, (B) Pre and Post structure of DRAMP02333 peptide, (C) superimposed structure of DRAMP02669 peptide, (D) superimposed structure of DRAMP03804 peptide. The figure was analyzed in Discovery Studio software package.

Open in a separate window

Fig. 5

The binding free energy of the top 4 peptide and protein complexes which were calculated by MM-PBSA methods.

The flexibility of the amino acid residues within the peptides was also assessed through RMSF analysis, and all four peptide molecules were found to display RMSF values below 2.5 Å, except Lys23 in the DRAMP00877 peptide; Tyr29, Asp33, and Val34 in the DRAMP02333 peptide; Glu13, Gly17, Thr18, Tyr21, and Gly22 in the DRAMP02669 peptide; and Val1, Val2, and Arg16 in the DRAMP03804 peptide. Lower RMSF values among the amino acid sequence indicate a less flexible and more stable nature of the peptide molecule (Fig. S6). The secondary structures of the peptide molecules were assessed to explore changes in the secondary structures, included the formation of α helices, turns, random coils, 310 helices, and pi helices (Fig. S7). The DRAMP00877 peptide displayed stable proportions of α-helices, β-sheets, turns, random coils, and pi helices, although some deviations were observed for 310 helices; however, 310 helices represented a lower percentage of secondary structure contents for the DRAMP00877 peptide. Random coils and helices were primarily noted for the DRAMP02333 peptide, which displayed a stable trend throughout the entire simulation times, with minor deviations in the proportions of 310 helices and turns. Random coils were primarily noted for the DRAMP02669 peptide, and stable profiles were observed for all secondary structures observed for this peptide throughout the simulation. Turns and random coils were primarily observed for the DRAMP03804 peptide, which featured stable trends, with some aberrations for α helix and 310 helix formation.

The SARS-CoV-2 Mpro consists of three domains; Domain 1, Domain 2, and Domain 3. Domain 1 contains residues 8–101; Domain 2 is formed by residues 102–184, which features an antiparallel β-barrel structure; and Domain 3 consists of residues 201–303 (Jin et al., 2020). Domain III (residues 201–303) comprises of five α- helices and is liable for the dimerization of enzymes exhibiting involvement in dimer formation (Chou et al., 2004, Kneller et al., 2020). In the Mpro dimerization, the troublesome the foreword of domain III has been adjusted by fragment elimination demonstrating that a decollated enzyme inexistent domain III stays as a monomer which is catalytically passive (Suárez and Díaz, 2020). Thus, Mpro constructs a functioning dimer via intermolecular interactivity, predominantly mesne the helical domains eliciting indispensable accosts of Domain III in the protease functionality (Kneller et al., 2020). The substrate-binding position is residing in the cleft interim domains I and II and the protomers, that cohere with every other via N-terminus residues 1–7, which are lying middle domains II and III with playing roles in the formation of the substrate-binding site. Four subsites namely; S1′, S1, S2, and S4 comprised the substrate-binding cleft (Mengist et al., 2021). The S1 subsite is formed with Phe140, Gly143, His163, His172, Cys145 and Glu166 while S2 comprised with Thr25, His41, and Cys145 amino acid residues; chiefly inlaid in electrostatic and hydrophobic interactions. The perfunctory subsites S3-S5 consist of Met49, His41, Met165, Glu166, and Gln189 amino acid residues. These shallow subsites can endure varied performance (Jin et al., 2020, Yang et al., 2005).

The binding interactions of the docked complexes were analyzed after taking snapshots at the 100, 200, and 250 ns time points during the simulation to understand interaction stability and deviations that occur across the simulation time. The DRAMP00877 peptide had several stable interactions at Domain 3 (Val 297, Phe294, His246, Glu240, Ile249, and Asp248), and equal stable interactions were found with Domain 2 (Gln107, Ile106, Pro108, Val104, Asp153, and Gln110) after 100 ns of simulation (Table S2, Fig. S2). Meanwhile, the DRAMP00877 peptide had the most stable interactions at Domain 3 (Phe294, Glu240, Ile249, His246, Asp248, and Asp245), and two stable interactions were observed with Domain 2 (Gln110 and Asp153) after 200 ns of simulation. Howbeit, the DRAMP00877 peptide had two stable interactions at Domain 2 (Gln110 and Asp153) and diverse stable interactions at Domain 3(Phe294, Ile249, Glu240, His246, and Asp245).

In case of the DRAMP02333 peptide molecule, there were diverse interactions in the simulation environment at Domain 3 (Lys236, Asn238, Leu287, Tyr237, Asp289 and Leu286) and four stable interactions at Domain 2 (Lys137, Cys128, Tyr126 and Ser123 in conjunction with an interaction at Domain 1 (Ala7). Domain2 and Domain 3 of Mpro are connected by a long loop region, and DRAMP02333 formed stable contacts with Asp197 within the loop region. Furthermore, the DRAMP02333 peptide molecule exhibited several interactions at Domain 2 (Val125, Gly124, Ser139, Lys137, Ser123, and Tyr126) and three stable interactions at Domain 3 (Asn238, Lys236, and Leu286) in tandem with an interaction at Domain 1 (Glu14 and Ala7) in the simulation environment after 200 ns of simulation (Table S3, Fig. S3), The DRAMP02333 constituted stable contacts with Asp197 within the loop region that connects Domain 2 and Domain 3 of Mpro. Howbeit, after 250 ns of simulation, the DRAMP02333 peptide molecule demonstrated various interactions in the simulation environment, including nine static interactions were noticed at Domain 2 (Val125, Lys137, Ser139, Gly124, Gly138, Tyr126, Pro122, Ser123 and Glu166) and six stable interactions at Domain 3 (Asn238, Lys236, Tyr237, Leu286, Leu287, and Met235) in conjunction with two rigid interactions at Domain 1 (Glu14, and Ala7). Here, the DRAMP02333 formed stable contacts with Thr199 within the loop region connecting Domain2 and Domain 3 of Mpro simulation.

For the DRAMP02669 peptide, sole one interaction audited at Domain 1 (Lys100), in addition with six interactions, which observed at Domain 2 (Lys102, Ile106, Val104, Gln110, Pro108, Asp153) and five interactions observed at Domain 3 (Phe294, Val297, Pro252, Asp248, Asp245) after the 100 ns simulation (Table S4, Fig. S4). Additionally, four interactions noticed at Domain 2 (Lys102, Val104, Gln107, Asp153), with five interactions observed at Domain 3 (Asp245, His246, Asp248, Phe294, Val297) after the 200 ns simulation. Besides, after the 250 ns simulation only three interactions were seen at Domain 2 (Lys102, Val104, Asp153), with notably seven interactions observed at Domain 3 (Pro252, His246, Phe294, Val297, Asp248, Ile249, Asp245).

There are six interactions observed for the DRAMP03804 peptide were formed with Domain 1 (His41, Ser46, Leu50, Thr25, Thr24, Met49), with six interactions observed at Domain 2 (Asn142, Met165, Pro168, Glu166, Phe140, Leu141) after the 100 ns simulation (Table S5, Fig. S5). Moreover, afterwards the 200 ns simulation, five interactions were observed at Domain 1 (Ser46, Met49, Leu50, Thr25, Thr24), including seven interactions observed at Domain 2 (Asn142, Met165, Pro168, Leu141, Phe140, Glu166, His172). Furthermore, only three interactions were observed at Domain 1 (Cys44, Ser46, Leu50), with six interactions observed at Domain 2 (Asn142, Glu166, Met165, Pro168, Leu141, His163) after the 250 ns simulation and this peptide constructed static contacts with GLN 189, ALA 191 among the loop region that appends Domain 2 and Domain 3 of Mpro. The interactions formed by the DRAMP03804 peptide primarily involve the active residues of the SARS-CoV-2 Mpro, which may interfere with the Mpro function. The co-crystalized ligand obtained with the Mpro structure binds in the substrate-binding pocket in extended conformations, and the inhibitor backbone atoms from an antiparallel sheet that interacts with residues 164–168 (Jin et al., 2020). Interestingly, the DRAMP03804 peptide molecules also featured three conserved interactions within these regions, at Met165, Pro168 and Glu166. It is cognizant that, DRAMP03804 peptide has Asn142, Gln189, Met165, Pro168, Leu141 and Glu166 residues that elicit resemblance with sundry of the abuzz residues of SARS-CoV-2 Mpro (Jin et al., 2020).

The peptide molecule DRAMP00877 interacted with the SARS-CoV-2 Mpro through conserved residues Phe294, His246, Glu240, Ile249, Asp153, and Gln110, whereas the DRAMP02333 peptide interaction involved Lys137, Tyr126, Ser123, Lys236, Ala7 and Leu286 throughout the simulation at 0 ns, 100 ns, 200 ns and 250 ns. Moreover, in all respects of the entire simulation trajectory, the amino acid residues Lys102, Phe294, Asp153, Val297, Val104, Asp248 and Asp245 of DRAMP02669 were involved in the interaction with SARS-CoV-2 Mpro and the Asn142, Gln189, Met165, Pro168, Leu141, Ser46, Glu166 and Leu50 residues contributed to the interaction between Mpro and DRAMP03804.

3.3. Physiochemical properties

We assessed the physicochemical properties of the best four peptide molecules (Table 4 ). The largest peptide contained 44 amino acids, whereas the smallest peptide only contained 30 amino acids. The molecular weights and theoretical isoelectric point (pI) were higher for the DRAMP02333 and DRAMP03804 peptides than for the other peptides. The peptides DRAMP00877, DRAMP02333, DRAMP02669, and DRAMP03804 contained 1,5,1, and 1 negatively charged residues, respectively. The total number of atoms included in the entire peptide molecule ranged from 400 to 600, except for DRAMP02333.