Introduction

COVID-19 caused by the severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2) has been declared a pandemic in early 2020 by the WHO. This new RNA virus has several folds higher transmission among humans and was first reported in December 2019 in the Wuhan province of China [1]. A wide range of symptoms were observed in people infected with SARS-CoV-2, which has a higher transmission rate among humans [2]. As of September 30, 2021, the WHO reported more than 23 crore confirmed cases of COVID-19 worldwide, including more than 48 lac deaths [3].

SARS-CoV-2 is a betacoronavirus that is the seventh human coronavirus. SARS-CoV-2 uses a two-step invasion strategy that involves attachment of the receptor-binding domain (RBD) to the angiotensin-converting enzyme-2 (ACE2) receptor followed by fusion to the cell membrane [4,5,6]. Viral attachment and entry are of particular interest among possible therapeutic targets because they represent the first steps in the replication cycle and take place at a relatively accessible extracellular site [7, 8]. Structurally, SARS-CoV-2 has four main structural proteins, including the spike (S), small envelope (E), membrane (M) glycoprotein, and nucleocapsid (N), as well as several accessory proteins. Human spike protein facilitates binding of envelope viruses to host cells by human angiotensin-converting enzyme 2 (hACE2), followed by proteolytic activation by human proteases. Hence, blocking of this protein–protein interaction can disrupt infection efficiency and pathogenicity. The S protein, which was seen at the atomic level using cryo-electron microscopy, is responsible for viral attachment and membrane fusion. The total length of the SARS-CoV-2 S protein is 1273 amino acids long and consists of a signal peptide (amino acids 1–13) located at the N-terminus, the S1 subunit (14–685 residues), and the S2 subunit (686–1273 residues). The spike protein has an extracellular N-terminus, a transmembrane (TM) domain anchored in the viral membrane, and a short intracellular C-terminal segment with a molecular weight of 180–200 kDa [9,10,11,12,13,14].

Even though COVID-19 can be prevented using vaccines and already accessible medications, it is still important to create new therapies that are more patient-friendly or less restrictive [15]. Peptides are more selective and effective than small molecule drugs, and since they are also more tolerated, they are an appealing alternative to small molecule inhibitors. Peptide synthesis can also be readily established and fine-tuned. On the other hand, small molecules frequently necessitate the creation of new time-consuming and costly synthesis techniques. Large-scale peptide synthesis, on the other hand, might be costly, depending on the peptide’s nature. In general, peptide inhibitors must be bioavailable, powerful, stable, and safe in order to be used in clinical trials. Peptides are highly selective and tolerated due to their chemical composition, but they have low (oral) bioavailability and short half-lives due to their comparatively large size and proteolytic breakdown. Shortening peptide sequences, on the other hand, can increase activity and stability. Peptide size is an important factor as a therapeutic target. The larger peptide size prevents it from passing through the cell membrane. The extracellular membrane fails to adopt the correct conformation for proteins. In peptides, hotspot residues, which are typically no more than three or four residues, are known to contribute significantly to target recognition and receptor activation [16,17,18].

Computational approaches have been used to search for potential therapeutics against the S protein of SARS-CoV-2 [19]. Molecular docking is a well-known in silico structure-based approach used in drug discovery. Docking makes it possible to find new drugs with therapeutic potential by predicting molecular interactions between ligands and targets [20]. A computational method called molecular dynamics makes it possible to forecast how molecular systems will change over time. The identification of cryptic or allosteric binding sites is just one of the many functions that molecular dynamics simulations can perform in the drug discovery process. It improves on traditional virtual-screening techniques and allows for more accurate prediction of ligand binding energies [14].

Small molecules with a high affinity for the S protein of SARS-CoV-2 were discovered through a similar screening of prospective medications. Unfortunately, the majority of these substances fail to bind to the RBD-ACE2 complex’s binding surface. Although hesperidin was supposed to be present on the surface of RBD, it only covered part of the interface. Short peptide inhibitors were researched in the early attempts to block SARS-CoV, and amino acid changes were made to the S protein of SARS-CoV. The proposed peptide, however, was unable to preserve secondary structure, making it unable to completely block the SARS-CoV binding surface. Cyclodextrins and broad-spectrum antiviral nanoparticles were developed, simulated, and used to prevent additional viruses [21,22,23,24]. They are category 2 or 3 inhibitors, although it is unclear if they will work against SARS-CoV-2. Promising COVID-19 treatments may include proteins or stiff peptides with particular (multivalent) binding domains and conformations that match RBD.

In the proposed study, small peptide inhibitors were designed based on the SARS-CoV-2 spike protein and the virus binding domain (amino acids) of hACE2, which were found in the crystal structure of the novel coronavirus spike receptor-binding domain complex with its receptor ACE2. The main objectives of the current study were to design small peptide inhibitors based on size and hACE-2 structure and to investigate their binding affinity against SARS-CoV-2 by using molecular docking approach. The binding affinity of designed peptides was compared with known antiviral drugs. Molecular dynamics study was further performed to understand the stability of the most potent compound. Also, the ADME properties of the designed peptides were evaluated.

Materials and methods

Design of inhibitors

As therapeutic medications, the size of the peptide is always a major issue; hence, we used a modest residue inhibitor in the presented manuscript. Furthermore, while target identification can occur with as few as a few residues, other data suggest that even extensive binding grooves can be bound by peptides with short residue lengths [25].

Three-residue peptides were designed on the basis of hACE2-spike protein and SARS-CoV-2 interactions [26]. There are four critical binding components viz. α1, α2, β1, and β2 comprises of critical amino acids from 19 to 393 which binds with SARS-CoV-2. Small peptide inhibitors were designed based on the critical amino acids of hACE-2, which bind to SARS-CoV-2. Understanding the X-ray crystal structures of the novel coronavirus spike receptor-binding domain complex with its receptor ACE2 (PDB ID: 6LZG), the active amino acids of hACE2 that interact with spike protein were studied and used to design small peptide inhibitors [27].

Preparation of protein

The X-ray crystal structures of the novel coronavirus spike receptor-binding domain complex with its receptor ACE2 (PDB ID: 6LZG) were downloaded from the RCSB PDB (Protein Data Bank) under the criteria of resolution 2.5 Å. It consists of two chains: chain A, angiotensin-converting enzyme 2, and chain B, spike protein. Moreover, the obtained structure was also combined with water molecules and co-factors. Water molecules were removed manually once the PDB was imported into Molegro Virtual Docker (MVD version 6.0 software) [27, 28].

The charge was determined by MVD and allocated to the models throughout the preparation process. Additionally, potential missing bonds were assigned, and plausible explicit hydrogens were established. During this preparation process, molecules were assigned bond order, bonds, charge, hybridization, flexible torsion, and explicit hydrogens in ligands [29].

Ligand preparation

The structures of peptides and small-molecule antiviral agents were drawn using Chemdraw ultra-8.0, copied to Chem3D ultra 8.0 to create a 3D model, and finally subjected to energy minimization using molecular mechanics (MM2) and the molecular orbital package (MOPAC). All 3D structures were saved in the *.mol format for further docking process. The minimization was executed until the value of the root mean square gradient reached a value smaller than 0.001 kcal/mol [29].

Docking validation

The MVD scoring algorithm needed to be evaluated first for the crystal structures to ensure that docked ligands indicated a valid score and correct binding with receptor. In view of this, the re-docking of the incorporated cocrystallized ligand (2-acetamido-2-deoxy-beta-D-glucopyranose) was performed. The given program is considered to perform effectively if it can return a pose that is less than a predetermined root mean square deviation (RMSD) value from the known conformation [30].

Molecular docking study (Molegro Virtual Docker)

Molecular docking of designed small peptides and small-molecule antiviral agents in ongoing clinical trials for COVID-19 was performed using the software Molegro Virtual Docker (MVD) v 6.0. The software MVD was used to generate the grid, calculate the dock score, and evaluate conformers in the current molecular docking study [22, 31]. The Piecewise Linear Potential (PLP) algorithm is used to calculate the various scoring functions. The active binding site region was defined as a spherical region that encompasses all proteins within 15.0 Å bound crystallographic ligand atoms with selected coordinates of X, Y, and Z axes, respectively. During the docking process, the active site center was set to X =  −12.34, Y =  −16.21, and Z = 64.32, and the grid resolution was set to 0.30 Å. Default parameters were used for all calculations, which considered of maximum iteration of 1,500 and the highest population size of 50. To improve docking accuracy with the differential evolution algorithm, 10 independent runs were performed, followed by 10 reranked results. The MVD relies on a differential evolution algorithm in which it measures the MolDock score, rerank score, and Hbond interactions. Here, Eq. (1) demonstrates the MolDock score energy or total Escore measurements in which Einter stands for the energy of interaction between ligand and receptor. Whereas, Eintra represents the ligand’s internal energy. Aside from that, the Einter can be calculated by Eq. (2) [29, 30, 32, 33].

$${E}_{\mathrm{score}}= {E}_{\mathrm{inter}}+ {E}_{\mathrm{intra}}$$
(1)
$${E}_{\mathrm{inter}}= \sum\nolimits_{i=\mathrm{ligand}}\sum\nolimits_{j=\mathrm{protein}}[{E}_{\mathrm{PLP} }({r}_{ij}) + 332.0 \frac{{q}_{i}{q}_{j}}{4{r}_{ij}^{2}}]$$
(2)

Molecular docking study (AutoDock)

To compare the binding energies, the docking was also performed on AutoDock. The empirical free energy function was built using the regular approach and default parameters of the molecular docking AutoDock 4.2 software. In AutoDock Tools, only polar hydrogens were added to the protein file, and all water molecules were deleted. For these ligands, the docking technique was used to find many different conformations. The one with the best pose and the lowest binding energy was chosen [26, 34].

Molecular dynamics

Molecular dynamics simulation (MDS) was carried out using GROMACS. Optimization and parameter simulations were used as previously reported. AMBERFF14SB and the AMBER force field were used to handle each ligand–protein complex, which was solvated within a cubic box (100 100 100), allowing for a minimum of 10 marginal distance between proteins. To eliminate steric overlap, each system underwent many steps of steepest descents energy minimization. Following that, all of the systems went through a two-step equilibration phase at 300 K called NVT (constant number of particles, volume, and temperature) and NPT (constant number of particles, pressure, and temperature). NVT ensemble of 50 ps to stabilize the temperature of the system, followed by an NPT ensemble for 1 ns to stabilize the pressure of the system by relaxing the system and keeping the protein restrained. The simulation of the top two compounds with the highest moldock score was made at 150 ns. The Parrinello-Rahman algorithm and the V-rescale thermostat algorithm were used to control the pressure and temperature, respectively. After every 10 ps, the output data was collected. Using gmx rmsd and gmx rmsf, the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) analyses were computed. RMSF is a calculation of individual residue flexibility, or how much a given residue moves (fluctuates) during a simulation, and it reflects positional changes between complete structures over time [35, 36].

Calculation of binding energy through MM-PBSA

The g-mmpbsa tool was used to compute the ligand’s binding free energy. The g-mmpbsa used the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA) method to compute the binding free energy of compounds. For the last 5 ns of the MDS trajectory, the potential energy of solvation (electrostatic + van der Waals contact) and the free energy of solvation (polar + non-polar solvation energies) were calculated [37].

Swiss ADME web tool filtration

The free SwissADME web tool, available from the Swiss Institute of Bioinformatics (SIB), was used to calculate the physicochemical descriptors as well as to predict the ADME parameters, pharmacokinetic properties, drug-like nature, and medicinal chemistry friendliness of the selected compounds after clustering [38]. The compounds’ structures were converted to SMILES notations and then submitted to the online server for calculation. Accessing http://www.swissadme.ch in a web browser displays directly the submission page of SwissADME, where molecules’ SMILES are input [39].

Ligand energy inspector

Information about energy interactions for a ligand was evaluated by using the ligand energy inspector tool of the molegro virtual docker. It is to be taken into account that the ligand energy inspector evaluates the energy of the ligand (or pose) when invoked.

Results and discussion

Design of peptide inhibitors

The size of the peptide is always a strong consideration for therapeutic drugs, so in the given manuscript, we have considered a small residue inhibitor. Also, some reports suggest that even large binding grooves can be bound by small residue length peptides, despite the fact that target identification can occur with as few as a few residues [40, 41]. In the structure of ACE2 and SARS-CoV-2 (PDB: 6LZG), we analyzed the interacting amino acids of ACE2 with SARS-CoV-2. In total, 23 amino acids of ACE2 interact with SARS-CoV-2, viz. 19 (S), 24 (Q), 27 (T), 28 (F), 30 (D), 31 (K), 34 (H), 35 (E), 37 (E), 38 (D), 41 (Y), 42 (Q), and 45 (L) are in α1; 79 (L), 82 (M), 83 (Y) are in α2; 325 (Q), 329 (E), 330 (N), 353 (K), 354 (G), 355 (D), 357 (R), and 393 (R) are in β1 + β2. Nine inhibitors with three residues have been designed. Combinations of peptides were designed, which are depicted in Table 1.

Table 1 Design of small peptides with three residues amino acids

Validation of docking

The reliability of molecular docking was assessed by the validation of docking methods. In each complex, the RMSD between the docked structure and the incorporated co-crystalized ligand was computed. The prediction is considered successful if the RMSD of the docked posture is less than or equal to 1.0 from the experimentally observed conformation. With an RMSD of 0.97, Redock compound binds to the native co-crystallized ligand at the same position and orientation (Fig. 1).

Fig. 1
figure 1

Validation of docking with RMSD of 0.97 Å of co-crystallized (2-acetamido-2-deoxy-beta-D-glucopyranose) and redocked ligand

Molecular docking study

The designed inhibitors’ binding affinity with the spike protein was investigated. The docking results were given in terms of moldock score, rerank score, H-bond energy, and interacting amino acid residues present at the predicted active site of the protein. Ligands that possess the lowest or minimal energy create better interactions with proteins [15]. All the designed small peptide inhibitors have a binding energy of less than −120 kcal/mol and bind at the active site of protein. Among all, HisGluAsp showed the highest binding affinity to the predicted active site of the protein, with a moldock score of −152.38 kcal/mol and a rerank score of −120.39 kcal/mol. From all the directed non-covalent interactions between biomolecules, hydrogen bonding is considered the most important one, because it facilitates protein–ligand binding. With each additional hydrogen bond, the binding affinity of the ligand increases by one order of magnitude [42]. Moreover, HisGluAsp formed two H-bonds with active amino acids Lys 458 and one with Gln 474 residues present at the predicted active site of the protein (Fig. 2a). Compounds with the highest moldock score are inserted into the cavity of the protein’s activity site (the best 03 compounds are shown in Fig. 2b).

Fig. 2
figure 2

a Hydrogen bond interactions of HisGluAsp. b Best three docked compounds in the activity site of protein. Atom color, HisGluAsp; yellow, PheAsnLys; blue, GlyAspArg

In contrast, SerGluThr showed the lowest binding affinity to the predicted active site of protein with a moldock score of −124.58 kcal/mol, a rerank score of −82.39 kcal/mol, and an H bond interaction of −7.20. It was observed that SerGlyThr forms no H-bond with active amino acids. PheAsnLys showed considerable binding affinity in comparison with HisGluAsp with a moldock score of −147.35 kcal/mol, a rerank score of −110.13 kcal/mol, and an H bond interaction value of −8.13. It forms two hydrogen bonds with active amino acids Gln 474 and Lys 458. The results of docking performed on molegro was compared with AutoDock. Interestingly, it was observed that binding energy of AutoDock concord with MolDock score (Table 2). Table 3 elicits the distance, energy, and interaction of atoms from ligand and protein in hydrogen bond interactions of HisGluAsp. It has been observed that the distance of all amino acids is at significant range that is less than 3.5 Å.

Table 2 MolDock score, rerank score, hydrogen bond interactions, and active amino acids of designed peptides
Table 3 Hydrogen bond interaction of most potent designed peptide with its distance, energy, and participation of atom from ligand and protein

Molecular docking was performed for five small-molecule antiviral drugs remdesivir, tenofovir, ribavirin, favipiravir, and galidesivir. These drugs are in ongoing clinical trials for COVID-19. It has been observed that remdesivir showed highest moldock score among all five antiviral agents with moldock score value of −133.09 kcal/mol, rerank score −95.3 kcal/mol, and Hbond interaction of −12.84 kcal/mol. Other than this, all antiviral drugs displayed moldock scores of less than −70 kcal/mol. All antiviral agents showed less value of moldock score in comparison to seven small designed peptides.

Steric interactions (van der Waal) play an important role and can influence the binding of protein and ligand. The primary force behind the development of docking techniques has been steric complementarity, frequently with the addition of physicochemical complementarity. Figure 3 depicted HisGluAsp forms prominent steric interactions with various amino acids viz. Lys 458, Ile 422, Arg 454, Arg 457, Ser 459, Ser 469, Glu 471, Tyr 473, and Gln 474. Aside from hydrogen bond interactions, the strength and distance of steric interactions play an important role in ligand–protein binding (Table 4).

Fig. 3
figure 3

Steric and electrostatic interactions of HisGluAsp with amino acid

Table 4 Steric interactions of most potent designed peptide with its distance and strength

Molecular dynamics

Molecular dynamics simulation is a method for determining the conformational behavior, stability of proteins, structural dynamics, and protein–ligand complexes. With the introduction of MDS, the field of structure-based drug discovery has evolved [43]. The MDS investigation was carried out to better understand the dynamic perturbations that occur in the novel coronavirus spike receptor-binding domain complex with its receptor during ligand binding. It also allowed us to forecast the protein–ligand complexes’ stability (HisGluAsp and PheAsnLys). The RMSD and RMSF of protein–ligand complexes were determined in this study. After 30 ns, the entire trajectory was stabilized, according to RMSD analysis. As a result, all of the above parameters were estimated for the last 20 ns of the trajectory.

The root mean square deviation (RMSD) measures the divergence of the protein backbone from its initial structural conformation to its final conformation. The variations generated during the simulation of the protein can be used to evaluate its stability with respect to its conformation. Protein backbone variations are less in stable protein structures, and vice versa. Thus, high RMSD values would therefore be associated with considerable instability as a result of changes in the studied molecule’s structure.

Figure 4a shows the graph of RMSD (nm) and time (ns) for two most potent inhibitors HisGluAsp and PheAsnLys. Fluctuation in the RMSD was observed from 20 to 40 ns for both the ligands, whereas after 50 ns, the trajectories were in stable state. It is worth noting that the average RMSD values, throughout the plateau MD simulation interval (50–150 ns), were higher for PheAsnLys compared to HisGluAsp (0.25 nm vs. 0.17 nm). Throughout the plateau phase, the latter differential dynamic activity provides a more stabilized and confined accommodation for HisGluAsp within the binding site. The HisGluAsp protein trajectory was seen to remain constant from 50 to 150 ns till the MD simulation ended. HisGluAsp showed lower RMSD value indicates the marked stability of the compound with the protein [44, 45].

Fig. 4
figure 4

a Time dependent root mean square deviation plot of HisGluAsp and PheAsnLys after molecular dynamics simulation of 150 ns. b The RMSF of HisGluAsp and PheAsnLys

For each ligand-bound protein compared to the protein apo-state, the per residue root-mean-square fluctuation (RMSF) profile was calculated in order to acquire more insights into the stability of the complex binding site. The GROMACS “gmx rmsf” command line was used to calculate the specific backbone RMSF of each protein. The contribution of specific protein residues to the structural variations of the ligand/protein complex is shown by this flexibility validation criterion. Within the reduced starting structures, RMSF calculates the average deviation over time for each residue from its reference point. For assessing the significant change within structural motions, an RMSF cut-off value of 0.30 nm was appropriate, and residues with values greater than 0.30 were considered to have decreased mobility [46, 47].

The RMSF of loosely organized structures like loops, twists, and coils is higher, but the RMS fluctuation is lower in well-structured regions like the -helix and -sheets. The RMSF value was determined to predict the structural changes in protein structure caused by ligand binding. Figure 4b shows the plot of RMSF for HisGluAsp and PheAsnLys. The average RMSF value for HisGluAsp and PheAsnLys was found to be 0.15 and 0.24, respectively. HisGluAsp showed less RMSF value implies that it might act as a potential inhibitor. Positive RMSF values were observed for both ligands, with HisGluAsp showing lower RMSF values than PheAsnLys from N-terminal free residues to vicinal residues (0 to 150). It was observed that RMSF value declined for PheAsnLys near to C-terminal residues. Since the hACE2 pocket residues are close to the protein's C-terminal side, our results show that HisGluAsp and protein complexes are more stable than those of PheAsnLys because of lower value of RMSF.

During MD modeling, the interaction with residues Arg 454, Lys 458, Glu 471, and Gln 474 exhibited the most interactions with the ligand. The hydrogen bond with the ligand was formed by amino acid residues Arg 454, Lys 458, Glu 471, and Gln 474, which accounted for 90%, 100%, 90%, and 100% of the simulation time in the 50 ns trajectory, respectively. During the MD simulation, the majority of the interactions between the protein and ligand that were seen during docking were consistently maintained (Fig. 5).

Fig. 5
figure 5

The residues that interact with the ligand in each trajectory frame

MM-PBSA calculation

The total of all non-bonded interactions is the binding free energy [48]. It was calculated using the MM-PBSA method for HisGluAsp and PheAsnLys. For the last 5 ns of the MD trajectory, interaction energies such as van der Waal’s energy, electrostatic energy, SASA energy, polar solvation energy, and binding energy were determined and presented in Table 5. As can be seen, van der Waals (−238.54 ± 15.32) dominates and favors complex interactions, while electrostatic interactions (−48.32 ± 8.31) contribute far less. In the formation of receptor-peptides, the polar solvation energy (115.21 ± 10.35) was shown to be more unfavorable. HisGluAsp and PheAsnLys have binding free energies of −198.31 ± 15.24 kJ/mol and −182.24 ± 12.84 kJ/mol, respectively. According to the results, these compounds have a considerable binding affinity; however, HisGluAsp has a larger binding affinity than PheAsnLys (Table 5).

Table 5 Van der Waal’s, electrostatic, polar salvation, SASA, and binding energy in kJ/mol for top two complexes

Physicochemical properties of peptides

Peptide-based and small-molecule drugs fulfill different chemical spaces due to variations of native amino acids in their side chains, putting peptides outside the established descriptor cutoffs. The molecular weight of the most potent compound in relation to binding affinity follows Lipinski’s molecular weight rule with a value of 396.33. The number of hydrogen bond donors and acceptors is a fundamental molecular descriptor to predict the oral bioavailability of small drug candidates. It is generally assumed that hydrogen bond donors and acceptors impact passive diffusion across cell membranes, a fundamental event during drug absorption and distribution. In the case of Lipinski hydrogen bond acceptors (LHBAs), the HisGluAsp does not follow the Lipinski rule, having a value of 13. The Lipinski hydrogen bond donors (LHBDs) are determined by counting the numbers of OH and NH bonds in each molecule. HisGluAsp follows the trend of LHBD with a value of 1. Previous reports have revealed that while physical models of permeability for conventional drug-like molecules are well established, the descriptors of membrane diffusion for molecules are not completely understood, hindering the rational design of cell-permeable, “beyond-Rule- of-Five” (bRo5) molecules as therapeutic agents. Other properties validated for HisGluAsp include molar refractivity (40–130), rotatable bonds (1–10), total polar surface area (140), and total atoms (20–70) (Table 6) [49,50,51,52].

Table 6 Physicochemical properties of designed peptides

Energy inspector of peptide

Ligand energy inspector allows to get detailed information about the energy interactions of a given ligand with a given protein. In Table 7, the energy contribution of each atom of HisGluAsp with maximum binding affinity is depicted. The short-range interactions promote short-range conformational order, i.e., short segments of the chain have higher probabilities of assuming the same local conformations as in the native conformation. In Table 7, ligand energy inspector of the most potent molecule, HisGluAsp, is depicted and both short- and long-range interactions are in the negative range, which suggests promising interactions between ligand and protein and also the interaction of an atom with protein [53].

Table 7 Ligand energy inspector of most potent peptide (HisGluAsp)

BOILED-Egg model

The BOILED-Egg model is proposed as an accurate predictive model that works by computing the lipophilicity and polarity of small molecules. The BOILED-Egg model delivers a rapid, intuitive, easily reproducible yet statistically unprecedented method to predict the passive gastrointestinal absorption and brain access of small molecules useful for drug discovery and development. The white region is the physicochemical space of molecules with the highest probability of being absorbed by the gastrointestinal tract, and the yellow region (yolk) is the physicochemical space of molecules with the highest probability of permeating the brain. Yolk and white areas are not mutually exclusive. In Fig. 6a, the radar graph represents the value of various physicochemical properties, and in Fig. 6b, the most potent molecule has been shown to be highly absorbable in the gastrointestinal tract [30].

Fig. 6
figure 6

a Radar plot of ADME properties of most active molecule. b BOILED-Egg model of most active molecule

Conclusion

In summary, using molecular docking study, it was shown that peptide inhibitors designed from ACE2 provide promising targets against SARS-CoV-2. A total of nine small peptides were designed from the α1 section with 41–42-45 number of amino acids showed the most potent binding affinity, followed by β1 + β2 section. Maximum of designed peptide showed potential binding affinity in comparison to antiviral drugs. HisGluAsp was found to interact with COVID-19 spike protein by exhibiting highest binding affinity through forming a hydrogen bond with active amino acid Lys 458 and Gln 74. Throughout the all-atom 150 ns MD simulation, HisGluAsp depicted superior stability. The molecule has drug-like properties and has the best chance of being absorbed by the gastrointestinal tract. Van der Waal and electrostatic energy showed significant accessibility of the pocket for small drug-like peptides. The results of this study could aid in the development of potent drugs to combat against SARS-CoV-2.