1 Introduction

The prolific spread of COVID-19 caused by a novel coronavirus (SARS-CoV-2) from its epicenter in Wuhan, China, to every nook and cranny of the world after December 2019, jeopardize the prevailing health system in the world and has raised serious concerns about human safety [1]. The World Health Organization (WHO), on March 20, 2020, declared COVID-19 as a pandemic due to its worldwide penetration in a few months. The updates of August 07, 2020 have been reported 19.2 million infections and 717,754 death cases. The symptoms manifested by most patients with COVID-19 include fever, coughing, dyspnea, myalgia, shortness of breath, and radiological evidence of ground-glass lung opacities compatible with atypical pneumonia. However, some patients with asymptomatic or mildly symptomatic cases have also been reported [2,3,4].

The members of the subfamily Orthocoronavirinae (family Coronaviridae) are distributed within four genera: alpha (α), beta (β), gamma (γ), and delta (δ) [5, 6]. The host of alpha and beta-coronaviruses are mammals, whereas the host of gamma and delta coronaviruses are birds [7, 8]. The major epidemics of 2003, 2012, and 2019 (caused by SARS-CoV-1, MERS, and SARS-CoV-2, respectively) revealed the vulnerability of humanity before the growing threat and offensive of the members of beta-coronaviruses [9, 10]. The SARS-CoV-1, MERS-CoV, and SARS-CoV-2 transmitted from bats to palm civets, dromedary camels, and pangolin. These coronaviruses became lethal due to their ability of transmission from animals to humans as well as from humans to humans [11]. CFR (case fatality ratio: calculated by dividing the number of deaths from the disease by the number of diagnosed cases of disease and multiplied by 100) of SARS-CoV-1, MERS-CoV, and SARS-CoV-2 were reported 10%, 35%, and 5%, respectively. The phenomenal spread of SARS-CoV-2 across the globe exposed everyone to this lethal pathogen, and investigations are underway to use various approaches such as novel vaccine development, drug repurposing, and novel drug development to mitigate risks and threats associated with SARS-CoV-2 [12].

The SARS-CoV-2 has a genome of ~ 30 kb and is responsible for encoding the whole proteome of the pathogen. Both the structural and non-structural protein-coding region and accessory protein-coding region are the major parts of the intact coding RNA [13]. The four structural proteins constituting the main envelope of the virus, and eight accessory proteins are encoded by the genes located on the 3′-terminus. The structural protein encodes for the spike surface glycoprotein (S), a small envelope protein (E), membrane protein (M), and nucleocapsid protein (N), whereas the eight accessory proteins codes for 3a, 3b, p6, 7a, 7b, 8b, 9b, and ORF14 [14]. The non-structural proteins constitute the replication/transcription complex (RTC) that are encoded by two overlapping genes, ORFIa and ORF1ab, located on the 5′-terminus of the viral genome. The OFR1a and OFR1ab encode two long polypeptides, pp1a and pp1ab, respectively. These polypeptides, pp1a and pp1ab, are directly translated and then proteolytically processed by two main viral proteases, papain-like protease (PLpro) and 3-chymotrypsin-like protease (3CLpro) also known as main proteases (Mpro) [15]. The former one is responsible for the cleavage of non-structural proteins (nsp) 1, 2, and 3, while the latter (3CLpro) cleaves the polyprotein at 11 discrete sites downstream of nsp4 to produce different non-structural proteins that play a pivotal role in the viral life cycle. On account of these multi-faceted roles, the 3CLpro is considered an attractive target for anti-coronaviruses drug development [16].

Using in silico approaches to design peptides is widely practiced by researchers. For instance, Neeraj et al. designed antimicrobial peptides to target the PBP5 protein [17]. Similarly, Junaid et al. used extensive residues scan approach to design potent peptide inhibitors from the cagA and ASPP2 interfaces [12, 18,19,20,21]. Using a similar approach, Carolina et al. designed small peptides of 15-mers in length to target the Hepatitis E virus (HEV) [22]. Taking the advantage of these successful implementations, herein using peptide inhibitor of SARS to perform residue scan-based methodology to potentially design more active peptides against the SARS-CoV-2.

Earlier studies suggested an important role of 3CLpro (Mpro) in cell replication and maturation; therefore, inhibition of this target will greatly contribute in controlling the COVID-19. 3CLpro is composed of three domains I (8–101), II (102–184), and III (201–303). Domain II and III three are connected through a loop (185–200) [23]. In this study, in silico mutagenesis-based remodeling of SARs-CoV-1 peptide (ATLQAIAS) was performed to inhibit the activities of 3CLpro of SARS-CoV-2 for the curtailment of COVID-19. The activity was confirmed by molecular dynamics simulation and free energy calculations. This study provides a basis for designing peptide inhibitors against the SARS-CoV-2 to confront this pandemic.

2 Materials and Methods

2.1 Structure Retrieval and Modeling of Peptide (ATLQAIAS)

The crystal structure of the 3CLpro complex was retrieved from the protein databank (http://www.rcsb.org/) and assessed for Structural deformities [24]. The structure of the SARS-CoV-2 main protease (3CLpro) contains N3 ligand which was removed prior to peptide docking. The peptide sequence (ATLQAIAS) was obtained from the previous literature and modeled using PEP-FOLD3 (https://bioserv.rpbs.univ-paris-diderot.fr/services/PEP-FOLD3/) [25], which is an online web server for peptides modeling. Previous computational and experimental studies reported that this peptide, along with others, remained active against the SARS virus and potentially inhibited the virus replication [26, 27].

2.2 Protein–Peptide Docking

To test the activity of the (ATLQAIAS) peptide against the SARS-CoV-2, molecular docking of the peptide and 3CLpro from SARS-CoV-2 was carried out to understand the interaction pattern. For this purpose, ZDOCK (http://zdock.umassmed.edu/) [28] was used with a blind docking option to increase and test the reliability of the server. Among multiple docked conformations, only those conformations (10 conformations) were selected and analyzed, which make interactions with active sites key residues such as His41 and Cys145 which form a catalytic dyad and do the proteolysis. However, the other conformations which are docked far from the active sites were discarded.

2.3 Interface Analysis and Peptide Library Construction

Understanding the binding interface of two interacting proteins or peptide is vital for designing small peptide inhibitors. The use of small peptides derived from the native peptide has been widely used in the designing of short peptide derivatives. There is a similar strategy to manipulate the interface residues of the reference-native-peptide-3CLpro and to design decoy peptides against SARS-CoV-2 using the machine learning protocol implemented in MOE [29]. The study of Zhang et al. 2006 suggested that ATLQAIAS peptide forms interaction with the key active site residues when docked with the 3CLpro. After this confirmation, information regarding the key residues in the interface is important for mutants modeling. For this purpose, the Alanine scanning strategy was exercised to calculate the impact of each residue in the interaction with the 3CLpro. To grasp this strategy, the dAffinity and dStability parameters were monitored in the ASM (alanine scanning module) of MOE. The dAffinity and dStability values show the relative change in binding energy when a particular amino acid is changed into alanine. These scores reveal important information regarding the importance of a particular residue, after the confirmation of key residues selection of the hotspot residues for residue scan to be replaced by the 19 amino acids. For this purpose, the residue scan module of MOE was used using the Unary Quadratic Optimization (UQO) parameter under the LowModeMD ensemble. This residue scan generated a database of mutant peptides with their respective scores. The detailed mechanism of this alanine scanning mutagenesis and residue scan approaches has been discussed previously [30].

2.4 Molecular Dynamics Simulation

A comprehensive protein dynamics technique was used to evaluate the dynamic features of the 3CLpro and the designed peptides via the residue scan method. Rationally developed peptides with fierce dAffinity and dStability scores were obtained using the Amber18 package and subjected to molecular dynamics simulation (MDS) with the ff14SB force field [31]. Using the TIP3P water model with box dimension 10.0 Å, each structure was correctly solved and neutralized by incorporating Na + counter ions. Two steps of energy minimization were performed. For the first energy minimization, 6000 steps, while for the second 3000 steps of conjugate gradient minimization, were completed. Following the minimization, each system was heated using default parameters (300 K for 200 ps). Weak restraint was used for density equilibration for 2 ns, while the whole system at a constant pressure for 2 ns. A 100 ns MD under constant pressure was performed. For the temperature control, Langevin thermostat (1 atm, 300 K) was used [32]. Particle Mesh Ewald (PME) algorithm was used to compute long-range interactions [33, 34]. The cutoff distances were set to 10 Å. For the covalent bonds involving hydrogen, the SHAKE algorithm was used [35]. GPU accelerated simulation using (PMEMD.CUDA) was used for all the processes

$${\text{RMSD}} = \sqrt {\frac{{\mathop \sum \nolimits_{{i = 0}}^{N} \left[ {m_{i} *\left( {X_{i} *~Y_{i} } \right)^{2} } \right]}}{M}} ,$$

where N is the number of atoms, mi is the mass of atom i, Xi is the coordinate vector for target atom i, Yi is the coordinate vector for reference atom i, and M is the total mass. If the RMSD is not mass-weighted, all mi = 1 and M = N.

2.5 Post-simulation Analyses

Post-simulation analyses such as RMSD (root-mean-square deviation), RMSF (root-mean-square fluctuation), and radius of gyration (Rg) were calculated using CPPTRAJ and PTRAJ modules of Amber [36].

2.6 Binding Affinity Calculations

To estimate the binding free energy of the designed peptides towards the 3CLpro MMPBSA.PY script was used [37]. This method has been previously used by various studies for free energy calculations [20, 38,39,40]. Using 5000 frames from the 100 ns trajectory, the free energy of binding was estimated using the following equation:

$$\Delta G_{{{\text{bind}}}} = \Delta G_{{{\text{complex}}}} - \left[ {\Delta G_{{{\text{receptor}}}} + \Delta G_{{{\text{ligand}}}} } \right].$$

In the above equation, the ∆Gbind represents the total BFE while the other components in the equation represent BFE of the complex, the protein, and ligand molecules. Each term in the binding free energy was estimated using the following equation:

$$G = G_{{{\text{bond}}}} + G_{{{\text{ele}}}} + G_{{{\text{vdW}}}} + G_{{{\text{pol}}}} + G_{{{\text{npol}}}} - {\text{TS}}{\text{.}}$$

The above equation represents the polar, nonpolar, electrostatic, solvent-accessible surface area SASA, and van der Waals interactions, respectively.

2.7 Per-Residue Energy Decomposition Analysis

Per-residue energy decomposition is the best approach to understand the energetic contribution at the residues level significantly. Herein, to understand the impact of each substitution on the binding of the designed peptides, we used per-residue energy decomposition analysis.

2.8 Clustering of MD Trajectories Using PCA

To comprehend the motion of MD trajectories, an unsupervised learning method known as Principal Component Analysis (PCA) [41, 42] was performed to acquire knowledge regarding the internal motion of the system. For this purpose, an Amber module known as CPPTRAJ was used. The spatial covariance matrix was determined for eigenvector and their atomic co-ordinates. Using the orthogonal coordinate transformation, a diagonal matrix of eigenvalues was generated. Based on the eigenvectors and eigenvalues, the principal components were extracted. Using these PCs, the dominant motions during the simulation were plotted [43, 44].

3 Results

3.1 Structure Retrieval and Preparation

The crystal structure of the SARS-CoV-2 3CLpro was retrieved from RCSB using accession number 6LU7. The structure was prepared and minimized using Galaxyweb server [45]. The peptide sequence (ATLQAIAS) was retrieved from previous literature, which reported that ATLQAIAS inhibits the SARs-CoV-1 by targeting the main protease 3CLpro. The peptide was modeled and prepared for docking. 3CLpro has three domains with a connecting loop, which connects domains II and III. The structure of 3CLpro and the peptide (ATLQAIAS) is given in Fig. 1.

Fig. 1
figure 1

The 3D structure of 3CLpro in brown color while the active site is shown in balls in blue color. The Magenta color shows the ribbon representation of the octapeptide (ATLQAIAS). The interaction of the octapeptide with the 3CLpro key residues such as Thr26, Thr45, Ser46, Asn142, Gly143, Cys145, and Gln189 is given in the panel (A) right side

Full size image

3.2 The Interaction Analysis of the Wild (ATLQAIAS) and 3CLpro Structure

The ZDOCK server predicted different conformations for the peptide docked at different positions. The conformation surrounds the important catalytic dyad (His and Cys45) was selected. Using ZDOCK the initial interaction pattern was obtained. It can be seen that the native peptide ATLQAIAS formed eight hydrogen bonds with the active site residues. It can be seen that Thr26, Thr45, Ser46, Asn142, Gly143, Cys145, and Gln189 are involved in hydrogen bonding with the octapeptide. Among these residues, Cys145 acts as a catalytic dyad and is vital for the function of the 3CLpro. Also, Gln189 from the active site formed two hydrogen bonds. Hence, this confirms that the octapeptide acts as an inhibitor for the 3CLpro and ultimately blocks the SARS-CoV-2 protease function. The interaction pattern of the native peptide (ATLQAIAS) is given in Fig. 1.

3.3 Residue Scan to Design a Peptide Library

The interface residues of the peptide and the 3CLpro were identified using the residue scan approach, an integrated module in the Molecular Operating Environment (MOE). In the native peptide, which was taken as reference peptide, the key hotspot residues include Ala1, Thr2, Leu3, and Ala5, which forms significant interactions with the 3CLpro; however, the other residues were involved as supplementary and were found out to be less significant at the interface to make interactions. In the first step, we used the alanine scanning approach to select the most appropriate residue for fixed amino acid substitution to increase the binding affinity of the native peptide. To scan for the most important residues, alanine scanning reported the dAffinity of each residue substitution, which calculates the relative binding affinity changes upon substitution. The alanine scanning analysis suggested that five residues are preferable for substitution, which could impact the binding of the peptide. These residues include Thr2, Leu3, Ala5, and Ala7, and could be replaced by other residues to increase the binding affinity of the native reference peptide. Nineteen different amino acids replaced each selected residue, and the binding affinity was calculated for each. This residue scan analysis was observed via dAffinity, which was measured as a difference between the binding energies of the wild-type residue and the corresponding mutated amino acid residue. A mutational space of 95 mutations was designed and searched. By residue scanning, four different peptides were generated based on Affinity and dAffinity scores given in Table 1. The interface and the designed peptide library structure are given in Fig. 2.

Table 1 Final peptides’ library generated from residue scan approach based on the total affinity and dAffinity given in kcal/mol
Full size table
Fig. 2
figure 2

A The interface analysis of the wild peptide and the 3CLpro. B The designed peptide library, and each amino acid replacement is colored differently

Full size image

3.4 Interaction Analysis of the Designed Decoy Peptides with 3CLpro

As given in Table 1, it can be seen that the substitution at position (AWLQAIAS) T2W increased the binding affinity from − 9.4 kcal/mol to − 11.70 kcal/mol. The substitution of small polar threonine by large aromatic Tryptophan increased the binding affinity by forming an extra hydrogen bond with Gln189 of the 3CLpro. As given in Fig. 3A an extra h-bond is also formed with His41, which cannot be seen in the wild type. In the case of (AYLQAIAS) T2Y substitution, the binding affinity was observed to be − 11.56 kcal/mol, while the dAffinity was reported to be -1.00 kcal/mol. This substitution of Tyrosine, an aromatic amino acid, favors an extra interaction with Thr190 along with those reported in T2W. The replacement of (ATRQAIAS) L3R shifts the interaction paradigm and favors significant interactions through the binding energy is relatively lower than the first two peptides. Ten H-bonds were formed with the key active site residues, including the catalytic dyad His41 and Cys145. Furthermore, a single hydrogen bond was formed with Thr26, Asn142, and Glu166 was reported, while two H-bonds with Meth165 and three H-Bonds with Gln189 were formed. The affinity for L3R was reported to be − 11.64 kcal/mol. In the case of A5W (ATLQWIAS) substitution, extra interactions with the key residue His41 were observed. Here, the main interacting residues include Thr25, Thr26, Asn142, Gly143, Glu166, and Gln189. Hence, the substitution makes a more favorable environment for interaction with the 3CLpro and suggests that the designed mutant peptides efficiently interact with the 3CLpro and block the key residues which are required for the function of the protein. The interaction pattern of the designed decoy peptides is given in Fig. 3.

Fig. 3
figure 3

Interaction pattern of the designed decoy inhibitory peptides of SARS-CoV-2 3CLpro. A AWLQAIAS, B AYLQAIAS, C ATRQAIAS, and D ATLQWIAS

Full size image

3.5 Dynamics Stability of 3CLpro-Peptides Complexes

To investigate the dynamic properties of these designed peptides in complex with 3CLpro molecular dynamics simulation of all the complexes was performed. The root-mean-square deviation as a function of time was calculated to quantify the stability of each system. Herein, in the case of the wild type, the system showed invariant behavior at a different interval. Initially, the RMSD was found to be ~ 1.5Å until 40 ns, but then later on, the RMSD converged at different intervals until 100 ns. On the other hand, the T2W system relatively showed higher RMSD, but no major convergence was observed like the wild type. The average RMSD for the T2W was ~ 2.0Å. In addition, the T2Y complex remained more stable when compared to the wild and T2W. The RMSD remained ~ 1.8Å with a small convergence between 58 and 60 ns, but then until 100 ns, the system entered the equilibration state again. In the case of L3R, the average RMSD remained ~ 2.0Å, with an acceptable convergence between 55 and 60 ns. On the other hand, the A5W system also remained stable, with a convergence between 87–88 ns. The average RMSD for A5W was also reported to be ~ 2.0Å. These results show that the mutant systems relatively remained stable during the course of simulation than the wild type and imply that the mutant decoy peptides bind more tightly than the wild type and favor the stability of the whole complex. All the RMSD graphs are given in Fig. 4.

Fig. 4
figure 4

The dynamic stability of each system as RMSD. A different color represents each complex. The x-axis is showing time in nanoseconds, while the y-axis is showing RMSD in Angstrom

Full size image

3.6 Residual Flexibility of the Complexes

The RMSF graphs revealed that, in each formulated peptide, substitution has a specific impact on the residues' flexibility. As given in Fig. 5, the wild type showed relatively high residual flexibility, which is in uniformity with the RMSD results. The wild type was reported to be relatively unstable; hence, the residual flexibility also remained higher than the others. It can be seen that the other systems exhibited a more similar pattern of residual flexibility. In the case of the T2W and T2Y, the flexibility increased at the substitution point, while the L3R and A5W flexibility decreased. These results reveal that the designed decoy peptides bind to the protein target more efficiently than the wild type and affect the residual dynamics by stabilizing the flexibility as compared to the wild type. Figure 5 shows the flexibility of each system in a different color.

Fig. 5
figure 5

The residual flexibility of each system. The x-axis is showing the number of residues, while the y-axis is showing RMSF in Angstrom

Full size image

3.7 The Compactness of the Peptide Bound Systems

The Gyration radius (RoG) was calculated to understand the compactness of each peptide bound system. Herein, the Rg calculation was performed to understand how the peptides remained intact with the protein target during the course of the simulation. This compactness also determines the stability of the system. As shown in Fig. 6, the wild-type system showed highly converged Rg. Initially, the Rg value was 22.2Å, but after reaching 20 ns, the Rg of the wild type converged all of a sudden, and the Rg value increased to 23Å. The Rg value fell between 20 and 40 ns but remained higher until 100 ns. This shows that the wild-type system is not that stable, and the peptide is probably released and the water entered the active pocket and lessened the peptide binding.

Fig. 6
figure 6

Rg of all the systems. The x-axis is showing the number of frames, while the y-axis is showing Rg in Angstrom. Wild (black), T2W (blue), T2Y (Magenta), and L3R (navy), while the A5W is shown in red color

Full size image

On the other hand, the Rg of the T2W remained stable, and the system showed more compactness. A little convergence between 75 and 80 ns was observed, but overall the system remained highly compact, and the average Rg value was reported to be 22.2Å. In the case of T2Y, the system was comparatively more compact than the wild type and T2W. The average Rg value for T2Y was 22.2Å. The Rg value for T2Y decreased after 50 ns, which shows the tight binding of the peptide to the receptor. In the case of L3R, the systems remained compact, but convergence between 80 and 100 ns was observed, which is in comparison with the RMSD, which shows a little unstable behavior. The average Rg for L3R was observed to 22.2Å. In addition, the Rg value for the A5W was reported to be 22.1Å. No major convergence was observed except for a fluctuation at 40 ns and 80 ns. These results confirmed that the designed decoy peptides remained intact with the active site residues during the simulation time period. These results implied that the decoy peptides efficiently bind and inhibit 3CLpro. Rg(s) of all the systems are given in Fig. 6.

3.8 Binding Free Energy Calculation

A common technique, MM-GBSA, has been used to accurately predict the binding free energy of all the complexes. The Gibbs free energy determines the binding affinity between interacting molecules. Each energy term, including vdW, electrostatic energy, polar solvation, solvent-accessible surface area, and total binding free energy of all the systems, are given in Table 2. As given in the total binding energy for the wild reference peptide is − 59.13 kcal/mol. The vdW for wild reference peptide was reported to be − 72.73 kcal/mol, electrostatic − 110.34 kcal/mol, while SASA − 8.56 kcal/mol. For T2W, the total binding energy was observed to be -70.04, while the vdW was reported to be − 81.01 kcal/mol. The electrostatic energy for T2W was reported to be − 111.12 kcal/mol, while SASA was observed to be − 9.6 kcal/mol. For the decoy peptide T2Y, the total binding energy was reported to be − 66.47 kcal/mol. The vdW and electrostatic energies for T2Y were observed to be − 79.08 kcal/mol and − 106.72 kcal/mol, respectively. The SASA for T2Y was reported to be − 9.55 kcal/mol. Furthermore, the total binding energy for L3R was reported to be − 67.47 kcal/mol with − 80.73 kcal/mol vdW, − 225.82 kcal/mol electrostatic, and -9.52 kcal/mol SASA.

On the other hand, the total binding energy for A5W was reported to be − 65.55 kcal/mol. In the case of A5W, both the vdW and electrostatic energies increased when compared to the wild type. The vdW for A5W was reported to be − 76.45 kcal/mol, while electrostatic for A5W was observed to be − 124.13 kcal/mol. These results significantly suggest that the designed decoy peptides possess strong inhibitory potential than the wild reference peptide. The contribution of the substituted residues increases the vdW and electrostatic energies significantly and eventually increases the total binding energy. The total binding energy and all the related terms are given in Table 2.

Table 2 The total binding energy of all the systems, including the wild type and the designed decoy peptides
Full size table

3.9 Per-Residue Energy Decomposition Analysis

Furthermore, to understand the energetic contribution at residues’ level, we used per-residue energy decomposition approach. This analysis further clarified the impact of these substation and its contribution to the total energy. In case of the wild-type Threonine at position 2, − 1.19 kcal/mol. Upon mutation to Tryptophan, the energy increased up to − 7.52 kcal/mol. At the same position when Threonine was replaced by Tyrosine, the energy observed was − 5.2 kcal/mol. With this substitution, the energy of the Glutamine at position 4 also increased to − 4.35 kcal/mol. In case of the L3R peptide the energy paradigm changed by at different residues. On the other hand, the substitution at position 5, replacing alanine by Tryptophan, the energy changed from − 0.1 kcal/mol to − 2.92 kcal/mol. These results significantly justify the impact of the substitution on the total binding affinity of the designed decoy peptides. The bar graphs representing the residual contribution of each peptide are given in Fig. 7.

Fig. 7
figure 7

Per-residue energy decomposition analysis of all the designed peptides. The graphs are given in wild type, T2W, T2Y, L3R, and A5W sequence

Full size image

3.10 Principal Motions of the Wild Type and Decoy Designed Peptides

PCA was performed and plotted to pinpoint the principal structural variations exhibited by each system enacted by the binding of the designed peptides. The objective of the PCA analysis was to extract the knowledge on the conformational states of the Wild Type and the designed peptides complexes, using the trajectory of 100 ns. Significant motion in all the systems for the first three eigenvectors was observed. The first three eigenvectors indicated significant fluctuations, while the remaining eigenvectors showed localized fluctuations in each complex (Fig. 8). In the case of wild-type peptide complex, the first three eigenvectors contributed 72% variance to the total observed motion, while in T2W 72%, in T2Y 42%, in L2R and A5W accounted for 33% variance in motion. This behavior may explain the structural rearrangement due to the peptides binding.

Fig. 8
figure 8

Fraction of the first ten eigenvectors. The (%) contribution of each eigenvector obtained from the covariance matrix plotted against the corresponding eigenvector indices constructed from the MD trajectory

Full size image

The first two eigenvectors were plotted against each other to gain conceivable attributed motions. The continuous representation of the red-to-blue color shows the switching along simulation time from one conformation to another. The dots represent each frame, beginning with red and ending in blue. In the wild-type and T2Y peptide complexes, more periodic jumps and continuous overlapping can be observed, while in the case of T2W, L3R and A5W showed localized fluctuations (Fig. 9). Together, all these observations imply that mutations significantly affect the structure and change the internal dynamics of the complexes.

Fig. 9
figure 9

Principal component analysis (PCA) of wild and designed decoy peptides bound to the 3CLpro. The first PC1 and second PC2 from the PCA of the backbone carbon were used. The sequence is wild type, T2W, T2Y, L3R, and A5W

Full size image

4 Discussion

The prolific spread of COVID-19 caused by a novel coronavirus (SARS-CoV-2) from its epicenter in Wuhan, China, to every nook and cranny of the world after December 2019, jeopardize the prevailing health system in the world and has raised serious concerns about human safety [1]. The World Health Organization (WHO), on March 20, 2020, declared COVID-19 as a pandemic due to its worldwide penetration in a few months. The updates of August 07, 2020 have been reported 19.2 million infections and 717,754 death cases. The symptoms manifested by most patients with COVID-19 include fever, coughing, dyspnea, myalgia, shortness of breath, and radiological evidence of ground-glass lung opacities compatible with atypical pneumonia. However, some patients with asymptomatic or mildly symptomatic cases have also been reported [2,3,4]. The enormous number of confirmed COVID- 19 positive (19.2 million) and death cases (717,754) as reported in updates of August 07, 2020 suggests the failure of the entire world to cope with this critical prevailing corona situation. The fast penetration of SARS-CoV-2 across the globe had adversely affected the human health, and therefore, researchers adopted various strategies to combat and control the havoc created by this pathogen. Computational biology approaches are very fast in identifying the possible molecular targets for novel vaccine and drug designing [46]. Previous studies suggested the use of spike glycoprotein, proteases, and RNA-dependent RNA polymerases for controlling various diseases caused by diverse pathogens. The approaches such as novel vaccine designing, drug repurposing, novel drug development, and protein–protein interactions (PPIs) were applied to find a workable solution to stop this deadliest pathogen and provide a feasible protection to the masses.

To date, efforts are continuing to design small molecule inhibitor, vaccine, and many other therapeutic options are practiced. Many companies claimed success in vaccine design, but the recent re-infection and their long-term effectivity are still questionable. Therefore, designing such therapeutic small molecules or peptide inhibitors is strongly suggested which are of long-term benefit. In this regard, many studies reported thousands of small molecules and peptides, but their final therapeutic potential is still to be tested. The notion of protein–protein interactions (PPIs) is gaining popularity day by day due to its applications in deciphering many biological processes. Thus, PPIs can be used as promising targets by incorporating new modulating entities into PPIs network. Using the old drug or vaccine or peptides could aid this process to avoid such long experimental procedure. Hence, here, we have repurposed a small peptide from the previous study which reported the inhibitory effects of this peptide.

Previous study reported that peptides could be a starting point for the drug discovery against SARS. Using molecular modeling and experimental approaches, several peptides which are similar in structure and pharmacophoric features were designed and studies using molecular modeling and in vitro assays. Here, utilizing the knowledge from previous study that a peptide (ATLQAIAS) significantly halt the replication of the SARS by targeting the 3CLpro. We repurposed the previous peptide against the SARS-CoV-2 using the computational pipeline. The peptide (ATLQAIAS) was docked against the 3CLpro from SARS-CoV-2. This analysis suggests that the peptide interacts with the catalytic dyad of the 3CLpro and blocks it important residues. Considering the initial knowledge, a useful residue scan approach was utilized to improve the native peptide based on in silico mutagenesis approach. The in silico mutagenesis replaced each amino acid by 19 substitutes and search for high affinity peptides based on the dAffinity criteria. Our analysis yields only four peptides which possess stronger interactions and energy than the native wild-type peptide. The activity of these rationally designed peptides was further confirmed using molecular dynamics simulation and free energy calculations. Our analysis suggest that the designed novel peptides possess stable dynamics behavior and compactness. On the other hand, the binding free energy revealed that all our designed peptides possess stronger binding affinity than the native wild type. This was also confirmed by the per-residue energy decomposition analysis which revealed that the mutated residues changed the energy level of each residue. In addition, the clustering of the MD trajectories also revealed the variation in the dynamics of each peptide. These peptides can be directly synthesized by the experimental lab researchers and test it against the SARS-CoV-2.

In conclusion, these analyses revealed a new paradigm for the development of peptide inhibitors against the SARS-CoV-2. Using such long computational pipeline from basic molecular modeling to free energy calculations confirmed the behavior and inhibitory properties of the designed mutant peptides. Our wide-ranging analyses of binding affinities disclosed that our designed peptide owns the potential to hinder the SARS-CoV-2 and will reduce the progression of SARS-CoV-2-borne pneumonia. Our analysis strongly suggests the experimental and clinical validation of these peptides to curtail the recent corona outbreak. The only limitation of this study is the experimental validation of the final peptide candidate which is required to confirm its effectiveness in inhibiting the SARS-CoV-2 and its clinical use.