Targeting SARS-CoV-2 Receptor Binding Domain with Stapled Peptides: An In Silico Study.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has evolved into a pandemic of unprecedented scale. This coronavirus enters cells by the interaction of the receptor binding domain (RBD) with the human angiotensin-converting enzyme 2 receptor (hACE2). In this study, we employed a rational structure-based design to propose 22-mer stapled peptides using the structure of the hACE2 α1 helix as a template. These peptides were designed to retain the α-helical character of the natural structure, to enhance binding affinity, and to display a better solubility profile compared to other designed peptides available in the literature. We employed different docking strategies (PATCHDOCK and ZDOCK) followed by a double-step refinement process (FIBERDOCK) to rank our peptides, followed by stability analysis/evaluation of the interaction profile of the best docking predictions using a 500 ns molecular dynamics (MD) simulation, and a further binding affinity analysis by molecular mechanics with generalized Born and surface area (MM/GBSA) method. Our most promising stapled peptides presented a stable profile and could retain important interactions with the RBD in the presence of the E484K RBD mutation. We predict that these peptides can bind to the viral RBD with similar potency to the control NYBSP-4 (a 30-mer experimentally proven peptide inhibitor). Furthermore, our study provides valuable information for the rational design of double-stapled peptide as inhibitors of SARS-CoV-2 infection.

Introduction

During December 2019, several cases of pneumonia, resulting from an unknown virus, were reported in Wuhan, the capital of Hubei province in China. Later, it was determined the illness was caused by a novel coronavirus, which was named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).[1] Unfortunately, the disease evolved into a pandemic of monumental proportions[2] that has resulted in more than 162 million reported cases and over 3.3 million deaths around the world by May 2021 according to the World Health Organization (WHO).[3] The deleterious consequences of this pandemic have made COVID-19 an economic, social, and public health crisis.[4−6] Although some vaccines have reached the market, additional interventions are needed in the following scenarios: (1) vaccine efficacy could be affected by new viral mutations or improper storage; (2) patients not able to receive a vaccine due to health challenges or unavailability; or (3) a segment of the population averse to receiving vaccinations.

The coronaviruses (COVs) are enveloped, single-stranded RNA viruses that cause severe respiratory complications.[7] Some SARS-CoV-2 features include its rapid spread, ease of contagion, and a death toll of 3% of the diagnosed cases.[2] The SARS-CoV-2 virion is composed of four proteins: a spike protein (S), a membrane glycoprotein (M), an envelope protein (E), and a nucleocapsid protein (N).[8,9] The spike protein, which is found at the virus surface, is widely accepted as the key player in the infection process.[2,10,11] This protein can be divided into distinct areas: a receptor binding domain (RBD), a central helix/heptad repeat, and a C-terminal region that associates with the plasma membrane of human cells.[12,13] The RBD of the S protein undergoes conformational changes to maximize its association with its target, the human angiotensin-converting enzyme 2 (hACE2), to achieve human contagion.[9,14,15] Thereafter, host proteases cleave the S protein into two subdomains: an N-terminal S1 portion and a membrane-bound C-terminal S2 unit. The receptor binding event and the proteolytic cleavage work synergistically to promote viral entry into human cells.[10] It has been reported that the SARS-CoV-2 RBD presents 10–20 times higher binding affinity for hACE2 (∼15 nM) than the RBD of SARS-CoV.[13,16,17] A recent in silico study conducted by Ponga and co-workers[18] evaluated the binding affinity and bond-breaking force between the SARS-CoV-2 spike protein and hACE2 receptors. The model estimated an energy of 12.6 ± 1 kcal/mol, resulting in a dissociation constant of KD = 1.3 nM, highlighting the high affinity of the complex.

Given its role in infection, affecting the interaction between the S protein RBD and the hACE2 has been identified as a promising strategy to prevent SARS-CoV-2 contagion.[2,19−24] Crystal structures of the S protein and the hACE2 show that the RBD presents a network of H-bonds with the α1 and α2 helices and a loop linking the β3 and β4 antiparallel strands on the hACE2 protein.[17,21,25] However, most of the key interactions are between the RBD of the S protein and the α1 helix.[7,26]

The first potential peptide therapeutic against SARS-CoV-2, to our knowledge, was reported by the Pentelute group.[27] These scientists designed 23-mer peptides derived from the hACE2 α1 domain, which were refined by molecular dynamics (MD) simulations. Bio-layer interferometry on the top designed peptide showed that the structure presented a dissociation constant, KD, of 1.3 μM for the Sino Biological insect-derived SARS-CoV-2-RBD (the binding of the 23-mer sequence derived from the hACE2 was not reported), a value around 100 times higher compared to the results published by Wrapp et al. for the binding of the hACE2 full length and the spike protein, which determined a KD of ∼14.7 nM.[13] Nevertheless, this was the first precedent to demonstrate that peptides could be used as therapeutics against SARS-CoV-2. Unfortunately, the N-termini biotinylated peptide from the Pentelute group[27] did not associate with human embryonic kidney (HEK) expressing the SARS-CoV-2-RBD or other insect-derived variants. In a theoretical study, Basit and co-workers designed a truncated hACE2 peptide that was 98 amino acids long. In silico studies using ZDOCK followed by MD simulations predicted the designed peptide should bind to the S protein RBD with higher affinity than the native hACE2 protein (ΔG of −12.7 and −10.7 kcal/mol, respectively).[28] The group of Rana performed MD simulations of 136 different 23-mer peptides.[29] The researchers suggested that their top peptide should present higher affinity for the RBD domain than the native hACE2 peptide as can be deduced by the obtained binding free energies: −304.1 vs −155.8 kJ/mol, respectively. Another intriguing study was done by Sitthiyotha and Chunsrivirot.[30] These authors cleverly combined computational protein design (Rosetta) and MD simulations (AMBER) to generate a library of 25-mer peptides based on residues 21–45 of the α1 helix of hACE2. Their results predict that their peptides will display higher affinity toward the SARS-CoV-2-RBD than the peptides developed by the Pentelute group.[23] Another study that is noteworthy was reported by Chowdhury and co-workers.[31] In such study, the authors used docking to screen 50 peptides that presented known activity against SARS-CoV-1 and, after MD simulation analysis, have identified two leads (29-mer and 33-mer peptides) that should inhibit SARS-CoV-2. A theoretical study done by Han and Král[7] showed that a peptide containing elements from the hACE2 α1 and α2 helices (ca. 75–139 mer), with the addition of carefully selected linkers to stabilize the peptidic structure, presents more stability and better binding affinity to the RBD than a peptide derived only from the α1 helix (34-mer).[7] For a complete review on bioactive peptides in the perspective of their potential activity against SARS-CoV-2, we recommend to read the literature review conducted by Bhullar and co-workers.[32]

The work of Han and Král[7] identified the importance of α-helix stability on therapeutic peptides against COVID-19, a key feature that can be appreciated on available crystal structures (Protein Data Bank (PDB): 6M0J, 6M17). The main consideration when designing this type of therapeutics is that short peptides could lose their secondary structure,[33] leading to a disruption of the bioactive conformation in the absence of a complete protein fold.[34] As a result, the therapeutic utility of such structures would be limited. Thus, a longer peptide would be desirable to keep the α-helical character. However, it may be unpractical from an immunogenic, synthetic, and cost perspective to prepare long peptides. Shorter sequences (∼25 mer) will be attractive therapeutics due to lower production costs and higher yield. However, such molecules should be able to maintain their secondary structure, be stable to degradation, and present good water solubility.

As can be appreciated from the discussion above, there are few studies exploring peptide-based antivirals to treat COVID-19. Herein, we designed short peptide sequences that may keep their 2ry structure and display high metabolic stability. Stapled peptides are a novel therapeutic modality that allows locking a structure in its bioactive conformation through the site-specific introduction of a chemical linker.[34−36] These peptides can inhibit intracellular protein–protein interactions (PPIs), such as in the case of SARS-CoV-2/hACE2, because they are capable of covering multiple contact points.[34,37,38] Further, the addition of a staple can keep the 2ry structure in place.[36,39] An additional advantage of such bonds is their higher metabolic stability due to their improved proteolytic resistance (compared to natural amino acids).[34]

Recently, while this work was in progress, Curreli and co-workers[40] reported the design, synthesis, and biological evaluation of 30-mer double-hydrocarbon stapled peptides based on the hACE2 helix (no computer simulation was associated to their work). The reported stapled peptides showed high helical contents (50–94% helicity), no cytotoxicity at the highest dose tested, a good profile of resistance to degradation by proteolytic enzymes in human plasma and the most active peptide possessed antiviral activity (half-maximal inhibitory concentration (IC50) reported around: 1.9–4.1 mM and 2.2–2.8 mM) measured in different cell lines. Maas et al. also reported the synthesis and evaluation of hACE2-derived 35-mer peptides containing a mono-lactam bridge that are able to inhibit the RBD–hACE2 complex formation.[41] Their lead peptide presented increased affinity for the RBD (IC50: 3.6 μM, KD: 2.1 μM) compared with the control group; a 35-mer sequence extracted from the hACE2. In contrast, Morgan et al. reported that their 23-mer mono-stapled peptides effectively constrained the helical structure in solution, but none of those peptides prevented virus internalization,[42] which indicates that must be an optimal peptide length and number/position of the staples to reach antiviral activity. Moreover, these results agree and validate our theory that using stapled peptides derived from the α1 helix of hACE2 could lay the foundations for further optimization of a potential clinical candidate.

Therefore, since (to date) no SARS-CoV-2 specific drug has been described and the currently proposed pharmacological treatments are based on repurposed drugs,[43,44] aiming to develop novel therapeutics against SARS-CoV-2, we performed in silico studies of short, 22-mer stapled peptides that mimic the structure of the hACE2 receptor to identify candidates with higher affinity for the S protein RBD using different stapling chemical groups. We postulate such peptides will prevent the RBD–hACE2 association halting virus entry in human cells,[7,23] while displaying increased physicochemical and biological stabilities.

Methods

Peptide Design

We analyzed the structure of the hACE2 protein and its interactions with the RBD of the SARS-CoV-2 S protein to rationally design our inhibitors.[17] As indicated by Král[7] and Hazelhurst,[26] the α1 helix of hACE2 contains 10 out of 15 residues believed to be important for association with the RBD: Gln24, Tyr27, Asp30, Lys31, His34, Glu35, Glu37, Asp38, Tyr41, Gln42. Thus, our strategy focused on extracting a 22-mer peptide that maintains these relevant amino acids while including changes to improve α-helical stability and to optimize interactions with the RBD. As mentioned, other important residues from the hACE2 protein were not considered to avoid the design of inhibitors with a long peptide sequence, which will be unpractical from a synthetic point of view.

Others have suggested the α-helical character is key for the interaction of hACE2 with the RBD.[7] Interestingly, around residue 37 there is a kink, or bend, within the peptide that gives the appearance of two connected α-helices. This detail was considered during the design to allow the peptide to arrange in this bent conformation. Therefore, the proposed modifications to preserve the α-helical structure were located before or after residue 37. In our design, we introduced changes considering two aspects: (i) modifications to lock the secondary structure to preserve the α-helical conformation and to increase metabolic stability; and (ii) alterations to increase the affinity of the peptides for the RBD.

Regarding the first strategy, we introduced a lactam bridge by replacing residues Phe28 and Phe32 (Figure a, indicated by Xa and Xb). These two residues were chosen to achieve the optimal connection at residues i and i + 4,[45−47] between the amide-forming amino acids. We hypothesize this change will not alter the interactions because the Phe residues do not face the S protein RBD but will lock the α-helical conformation of the peptide.[45−47] Moreover, since two aromatic groups are removed, this will reduce the lipophilicity of the peptide (the calculated Log P of an FF dimer is 0.42, while the KD lactam bridge is −1.68[48]). The same lactam bridge was proposed for the Leu39 and Ser43 residues (Figure a–c, indicated by Xc and Xd). In addition, we proposed stapled peptides replacing the same residues (Phe28 and Phe32) with alkenes (Figure b, indicated by Xe and Xf). In this design, the optimal distance for connection is also at residues i and i + 4.[45,49] We postulate these new bonds will increase metabolic stability; i.e., stapled peptides will be more resistant to proteases (although they may be oxidized). Another proposed change included a double-lactam bridge between Ala25 and Phe32 using 1,4-diaminobutane (Figure c, indicated by Xg and Xh). One interesting characteristic about our designed peptides, compared to the ones proposed by Curreli et al.,[40] is the hydrophobic profile of the final structures. Our staples substitute hydrophobic residues in the original structure (i.e., Ala25, Phe28, Phe32, and Leu39). These changes should reduce the hydrophobic character of the peptides and, consequently, might display enhanced solubility in aqueous media.

Figure 1

Template sequence and structural basis of the proposed stapled modifications. (a) Lactam bridge replacing residues Phe28 and Phe32 and Leu39 and Ser43; (b) alkene stapling brace replacing Phe28 and Phe32 and lactam bridge replacing Leu39 and Ser43; (c) double-lactam bridge between Ala25 and Phe32 and lactam bridge replacing Leu39 and Ser43. Left: three-dimensional (3D) model. Right: two-dimensional (2D) structures of the proposed stapled bonds. “X” denotes the position of the staple within each peptide.

To further increase affinity for the S protein, we included amino acids that may form additional interactions with the RBD. For example, Asp30 was replaced by Glu to maximize contacts with Lys417 on the S protein. Additionally, for some proposed structures, the Leu45 on the C termini was substituted by an acetylated Lys [K(Ac)] or a Gln to maximize contacts with Gly446. This latter change should also enhance water solubility. The sequence of all design stapled peptides is presented in Table and their 2D structures can be found in the Supporting Information (SI) (Figure S1). An interesting aspect of our peptides is the fact that they were designed based on a short domain of the hACE2 α1 helix (Figure S2). Therefore, they are expected to be active against strains possessing mutations in the RBD that are distant from the interaction points found in the 22-mer hACE2 α1 helix. Examples of such mutations are Asn439Lys[50] (N439K) and Leu452Arg (L452R). According to the Centers for Disease Control and Prevention (CDC), the last cited mutation is associated with reduced effectiveness of treatments based on monoclonal antibodies.[51,52]

Table 1

Description of the Sequence of the Proposed Stapled Peptidesa

The staples insertion positions are highlighted in blue (according to the description from Figure ), and the replacement of original residues is highlighted in red.

Preparation of Protein and Peptides

The structure of the receptor and the peptides were obtained from the crystal structure of hACE2 bound to the RBD of SARS-CoV-2 (PDB: 6M0J).[25] For the docking protocol, the receptor employed was composed of the entire structure of the glycosylated RBD from SARS-CoV-2. The file was energy minimized using UCSF Chimera software version 1.14.[53] Amber ff14SB force field was used to assign charges to the structure and the parameters for the energy optimization were: steepest descent steps were set to 1000; steepest descent step size 0.02 Å; conjugate gradient steps 1000, conjugate gradient step size 0.02 Å. Then, DockPrep tool in UCSF Chimera software was employed to: remove water molecules, repair truncated side chains, add hydrogens, and assign partial charges. Protonation states were assigned at physiological pH (∼7.4). The structures of the peptides were built from the crystal structure using the 22 residues of the hACE2 α1 helix that lie on the interface between RBD and hACE2: Gln24, Ala25, Lys26, Tyr27, Phe28, Leu29, Asp30, Lys31, Phe32, Asn33, His34, Glu35, Asp36, Glu37, Asp38, Leu39, Phe40, Tyr41, Gln42, Ser43, Ser44, Leu45. The 18 peptides were built from the hACE2 original structure minimized (to maintain its natural conformation). This sequence was altered according to the proposed modifications (see Figure S1 (SI) for more details regarding the structures of all of the proposed peptides) followed by structure energy minimization in all positions where modifications were assigned (as described before for the receptor). Additionally, for comparison purposes, the best hydrocarbon stapled peptide based on the hACE2 helix published by Curreli and co-workers,[40] peptide NYBSP-4, was also modeled using the same protocol employed with the designed peptides. The NYBSP-4 peptide was also docked and underwent MD simulations using the same approaches described below. Furthermore, after the analysis of our MD results, we decided to analyze how our best stapled peptide would behave in the presence of mutations in the RBD. Therefore, we utilized the same crystal structure (PDB: 6M0J) and modified it in UCSF Chimera using the Dunbrack backbone-dependent rotamer library[54] provided with this software (version 1.14).[53] The mutated RBD was submitted to MD simulation using the protocol described in the MD section

Generation of Docking Solutions

To predict the binding affinity of our proposed peptides to the SARS-CoV-2 RBD, we used a rigid body protein–protein docking approach employing two docking servers: ZDOCK web server[55] and PATCHDOCK web server.[56] ZDOCK uses a fast Fourier transform (FFT) algorithm to allow an efficient global docking search. The method searches for all possible binding modes in the translational and rotational space between the protein and peptide, evaluating each pose using an energy-based scoring function. The energy function used by ZDOCK is Z score, which is cumulative of pairwise shape complementarity function with desolvation and electrostatics. Finally, the ZDOCK method ranks the predicted docking poses based on the Z score generated for each solution.[57] PATCHDOCK employs an algorithm based on shape complementarity principles to obtain the docking solutions and uses PATCHDOCK score as the energy function which ranks the docked models based on desolvation energy, interface area size, and geometric score.[56,58]

Since our peptides were designed to keep the interactions of the natural sequence of hACE2 and promote new interactions with the RBD, for both server protocols, the following key interface amino acids were used as restraints, according to experimental data previously published:[17,25,59] RBD Tyr449, Tyr453, Asn487, Tyr489, Gln493, Gln498, Tyr500, Asn501 and from hACE2 Gln24, Tyr27, Lys31, His34, Glu37, Asp38, Tyr41, Gln42. Using this approach, ZDOCK and PATCHDOCK servers filter results such that only those with specified residues in the binding site are returned, focusing on the solutions that account for the known important interactions and other contacts that might arise from the modifications proposed.

Docking Solutions Refinement and Validation

Since a rigid approach was employed to generate the docking solutions and due to the fact that during protein–protein interactions, both side chains and backbone might change their conformation, we used FIBERDOCK server[60] for the refinement of the docking candidates. FIBERDOCK accounts for a high level of flexibility regarding protein interactions, providing a higher probability of finding a near-native conformation for the predicted solutions. The method used by FIBERDOCK models both side-chain/backbone flexibility followed by a rigid body optimization on the ligand orientation targeting the problem of flexibility and scoring of the solutions produced by ZDOCK and PATCHDOCK. This refinement algorithm mimics an induced fit process, in which the solution candidates generated previously are re-ranked according to the minimum global energy to identify suitable protein–peptide complexes.[60,61]

During the refinement, the following parameters were used: an antibody–antigen complex type was chosen; a restricted side-chain optimization in both receptor and ligand was performed; the Monte Carlo cycles were set to 100; an atomic radii scale of 0.95 was used. Following this first step, a second refinement was performed using the FIBERDOCK candidate with the best alignment at the interface according to the orientation of the natural hACE2 binding sequence. In this second FIBERDOCK run, we used a full side-chain optimization and the best solution from each peptide was ranked according to minimum global energy (Tables and 3). The visualization of all complexes was performed using UCSF Chimera version 1.14.[53] The final selection of docked complexes was based on the global energy of the bound predicted complexes after the two steps of refinement and also on visual inspection to assess their spatial orientation in the interface guided by the crystal structure of the SARS-CoV-2 RBD–hACE2 complex. From that point, we selected the candidates that presented a right orientation in the interface, the lowest root-mean-square deviation (RMSD) compared to the original 22-mer peptide derived from the α1 helix of hACE2, and the lowest global energy (higher stability) calculated by the scoring function of the FIBERDOCK method. The efficiency of this approach was assessed by a re-docking protocol of the original peptide sequence extracted from the crystalized hACE2 (PDB: 6M0J). Then, we analyzed the orientation in the binding site and measured the RMSD of the docked predicted 22-mer peptide derived from the α1 helix of hACE2 (QAKTFLDKFNHEAEDLFYQSSL) and compared to the same structure extracted from the crystal.

Table 2

Docking Energy Information Resulting from PATCHDOCK and Refined with FIBERDOCKa

peptide	global energy	attractive eVdW	repulsive eVdW	ACE	HB
NYBSP-4	–54.30	–30.03	14.27	4.01	–1.94
modification 15	–51.52	–24.10	10.76	3.39	–2.40
modification 10	–49.90	–22.66	9.04	3.32	–2.35
modification 9	–49.75	–23.38	13.91	3.39	–2.06
modification 4	–49.65	–23.20	12.87	3.36	–1.60
modification 14	–49.35	–24.96	18.36	3.52	–1.83
modification 11	–49.10	–23.89	7.31	5.28	–1.06
modification 5	–48.19	–23.53	13.92	3.36	–1.66
modification 1	–46.74	–22.88	8.11	5.66	–1.52
original 3	–46.09	–23.73	12.62	6.11	–1.68
modification 6	–46.07	–23.38	9.43	5.48	–0.76
modification 2	–45.84	–22.98	10.73	5.85	–1.46
modification 3	–45.12	–24.19	9.08	6.80	–0.59
modification 13	–45.00	–25.14	9.28	6.89	–0.47
original 1	–44.65	–21.68	5.65	5.80	–1.51
original 2	–44.41	–23.27	12.57	5.97	–0.98
modification 8	–44.34	–24.48	10.67	6.90	–0.49
hACE2	–43.35	–23.34	9.09	5.63	–0.99
modification 7	–37.26	–24.58	13.72	7.67	–1.16
modification 12	–37.26	–24.58	13.72	7.67	–1.16

Energy presented in kcal/mol; VdW: van der Waals interaction; ACE: atomic contact energy (desolvation energy); HB: hydrogen bond.

Table 3

Docking Energy Information Resulting from ZDOCK and Refined with FIBERDOCKa

peptide	global energy	attractive eVdW	repulsive eVdW	ACE	HB
modification 3	–53.81	–26.31	10.06	5.55	–0.20
NYBSP-4	–53.42	–25.70	12.20	2.62	–0.24
original 3	–47.75	–24.67	7.33	5.64	–0.83
modification 13	–46.60	–27.81	14.30	7.13	–2.33
hACE2	–46.49	–24.03	11.84	5.73	–1.58
modification 12	–45.84	–27.40	11.83	6.63	–0.48
modification 2	–45.20	–26.04	12.04	5.20	–0.36
modification 11	–44.65	–25.88	16.10	6.38	–0.85
modification 15	–43.90	–26.24	14.97	6.24	–0.85
modification 8	–43.15	–28.55	18.28	7.45	–0.81
modification 7	–42.91	–29.47	18.32	8.46	0.00
modification 4	–42.74	–21.60	9.25	1.99	–0.51
original 1	–42.71	–23.05	14.70	5.25	–0.97
modification 10	–41.98	–25.10	13.35	2.56	–1.60
modification 14	–41.35	–20.16	9.19	1.95	–1.00
modification 5	–39.77	–21.46	12.62	2.47	–1.30
modification 8	–39.32	–24.07	16.13	7.38	–2.22
modification 6	–38.72	–23.63	10.62	6.25	–1.00
modification 1	–38.54	–24.29	9.50	6.31	–0.95
original 2	–37.83	–20.38	5.82	5.03	–0.09
modification 9	–32.11	–19.31	7.95	1.63	–0.36

Energy presented in kcal/mol; VdW: van der Waals interaction; ACE: atomic contact energy (desolvation energy); HB: hydrogen bond.

Molecular Dynamics Simulations

MD simulations were performed to analyze the stability and dynamic properties of the RBD of the S protein bound to the most promising stapled peptides. As controls, the crystal structure of the original 22 residues from the hACE2 α1 helix, which lies on the interface between RBD and hACE2 (residues 24–45), and the experimentally proven antiviral NYBSP-4 double-stapled peptide were employed.

The MD simulations were conducted using Desmond version 3.0 provided with the Schrödinger package 2020–4,[62,63] according to the following protocol: the RBD of SARS-CoV-2 bound to the designed peptide complex was placed in a cubic box with periodic boundaries and a minimum of 1 nm from any box edge to the solute. The box was filled with TIP4P water and Na+ or Cl– ions were added, as needed by the software, to neutralize the system and to achieve a physiological concentration (150 mM). Then, the system was relaxed using the standard protocol provided in Desmond, which consists of a mixture of predefined minimizations and previous molecular dynamics executions using a constant number of particles, volume, and temperature (NVT) and a constant number of particles, pressure, and temperature (NPT) ensembles designed to slowly relax the system, while not deviating substantially from the initial protein coordinates.[63] After, the simulations for all systems were performed under NPT ensemble with the temperature set to 300 K and the pressure set to 1.013 bar using the Nosé–Hoover[64] and the Martyna–Tobias–Klein[65] algorithms. Production of MD simulations was set to 500 ns, and they were performed using Desmond GPU with the OPLS3e force field[66] through the Holland Computing Center.[67] The OPLS3e force field[66] is based on the optimized potentials for liquid simulations (OPLS) first developed by Jorgensen and co-workers.[68] OPLS3e displays significant improvements in the representation of secondary structure elements in simulated peptides and native structure stability over a number of proteins.[66]

Finally, the root-mean-square deviation (RMSD) for both protein and peptide ligands, after MD simulations, were evaluated to understand the relative stabilities of the complexes. In addition, the trajectories of the most promising complexes were analyzed and recorded in video format using the Maestro interface from the Schrödinger package 2020–4.[62,63] The protein–stapled peptide interactions were monitored as well during the entire period of the simulation. Then, the interactions between the RBD of SARS-CoV-2 and the stapled peptides were analyzed using the analyze trajectory script[69] provided by Schrödinger, which reads a trajectory file and identifies interactions occurring between the defined sets.

Binding Free Energy Analysis

The estimation of the free energies for the binding between the RBD of SARS-CoV-2 and the most promising stapled peptide complexes (ΔGbind in kcal/mol), according to the analysis of the MD results, were computed using the molecular mechanics with generalized Born and surface area (MM/GBSA)[70] continuum solvation method[71] implemented in the Schrödinger package. The MM-GBSA binding free energy was estimated as follows: ΔGbind = Gcomplex – Greceptor – Gligand, where ΔGbind is the binding free energy and Gcomplex, Greceptor, and Gligand are the free energies of complex, receptor, and ligand, respectively. The thermal MM-GBSA script provided by Schrödinger[69,71] was used to calculate the ΔGbind for the studied complexes. This script takes in a Desmond MD trajectory, splits it into individual frame snapshots, and runs each one through MM-GBSA analysis. During the MM-GBSA calculation, 702 snapshots from the 500 ns MD simulation were used as input to compute the average binding free energy. The predicted free energies of binding are presented as average values (ΔGbind, Table ) along with the energy components used in the calculation.

Table 4

Predicted Free Energies of Binding (ΔGbind) from the MM-GBSA Analysis and the Energy Components Used in the Calculationa

peptide	ΔGbind	ΔGCoul	ΔGCov	ΔGHbond	ΔGLipo	ΔGVdW	ΔGPack	ΔGSolGB	ΔGSC
Mod11	–90.96 ± 15.05b	–140.02	4.55	–6.09	–22.95	–79.72	–5.45	158.62	0.09
NYBSP-4	–86.04 ± 18.81b	–127.81	4.97	–5.72	–23.88	–78.30	–3.61	150.14	0.13
Mod15	–81.47 ± 9.23b	–82.91	5.68	–4.18	–25.95	–75.06	–5.55	106.53	–0.03
hACE2	–75.48 ± 10.79b	–123.95	0.84	–4.88	–14.90	–62.58	0.57	129.53	–0.10

All energy terms are presented as average of all values calculated for each snapshot generated from the MD trajectory. All values are presented in kcal/mol. ΔGCoul: Coulomb energy contribution; ΔGCov: covalent energy contribution; ΔGHbond: hydrogen-bonding contribution; ΔGLipo: lipophilic energy contribution; ΔGVdW: van der Waals energy contribution; ΔGPack: π–π packing energy contribution; ΔGSolGB: generalized Born electrostatic solvation energy contribution; ΔGSC: self-contact energy contribution. For more information on the energy components, see Li et al.[71]

Standard deviation.

Results and Discussion

Docking Results

The results for the PATCHDOCK/FIBERDOCK and ZDOCK/FIBERDOCK studies are presented in Tables and 3, respectively. As it can be seen, for the first approach several proposed modifications presented a lower global energy (higher stability) compared to the sequence derived from the α1 helix of hACE2 (QAKTFLDKFNHEAEDLFYQSSL). Only two proposed modifications showed higher energies than the original sequence (modification 7 and 12). Regarding the second approach, only three proposed stapled peptides showed better energy results compared to the hACE2 original sequence. Additionally, for both docking servers, our best ranked stapled peptides present docking scores similar to the scores of the control NYBSP-4 (experimentally proven SARS-CoV-2-RBD peptide binder[40]). For the PATCHDOCK study, the best proposed peptide (modification 15) shows a docking score of −51.52 kcal/mol, while NYBSP-4 presents −54.30 kcal/mol. For the ZDOCK approach, the best docking solution (modification 3) presents a docking score slightly better (−53.81 kcal/mol) than the NYBSP-4 peptide (−53.42 kcal/mol). These findings suggest that the highest-ranked docking solutions from our study are comparable to the control in affinity toward the RBD of SARS-CoV-2.

Comparing the results of the two docking approaches indicates that even though the overall energies are similar (−51.52 kcal/mol for the best solution using the PATCHDOCK protocol and −53.81 kcal/mol for ZDOCK best’s solution), the top-ranked candidates are distinct for each method. This divergence reflects the distinct methods employed to generate the docking poses in each server.

Regarding the top predicted candidates from the PATCHDOCK protocol, the replacement of Leu45 by K(Ac) may increase the affinity for the SARS-CoV-2 RBD since all top 3 ranked predictions (modifications 15, 10, and 9) present this modification (Figures S1 and 2B). This change allows an additional hydrogen bond with Gly446 in the RBD increasing the affinity to the S protein. Additionally, the substitution of Asp30 by a Glu seems to favor the affinity between the systems since the top 2 ranked solutions present this designed replacement (Figure A). These observations help to explain the superior hydrogen-bond energy contribution of these structures compared to the hACE2 peptide (−2.40, −2.35, and −2.06 kcal/mol of the top 3 predictions vs −0.99 kcal/mol in hACE2).

Figure 2

Spike protein (RBD) of SARS-CoV-2 (yellow) and modification 15 (blue) complex from PATCHDOCK prediction. (A) Detailing of the polar interactions between the substituted Glu residue and Lys417. (B) Detailing of the polar interactions determined by acetylated Lys in modification 15. The nonpolar hydrogens were hidden for visualization purpose.

As for the ZDOCK top predictions, the Asp30 substitution by Glu also seems to play an important role in the predicted binding affinity since modifications 3 and 13 show this feature (Figure S1, SI). Another interesting point from the ZDOCK solutions is the fact that two out of the three top-ranked stapled peptides present a double-lactam bridge between Ala25 and Phe32 (modification 13 and original 3, Figures and S1). The same type of staple is also present in the top candidate from PATCHDOCK protocol (modification 15), which indicates the positive effect of this modification in the top predictions.

Re-docking Evaluation

Regarding the re-docking protocol employing the original 22 amino acids from hACE2, the overall procedure was considered satisfactory after the RMSD analysis. According to the Critical Assessment of PRediction of Interactions (CAPRI)-inspired threshold for success,[72] a docking protocol finds a near-native conformation if the peptide lies within 4.0 Å backbone RMSD of the native peptide bound to the receptor (i.e., the CAPRI criterion for an acceptable peptide–protein docking prediction). The prediction made by ZDOCK server produced a docking solution for the original structure with an RMSD calculated of 2.782 Å, compared to the crystal structure. After the two refinement steps using FIBERDOCK, the best obtained prediction presented an RMSD of 0.556 Å (Figure S3, SI). For the PATCHDOCK server, the best predicted solution presented an RMSD of 3.144 Å, and after the use of FIBERDOCK, the difference was 0.373 Å (Figure S3, SI). Therefore, we concluded the docking protocol is suitable for producing acceptable peptide–protein docking predictions regarding this specific complex, and the refinement process can improve the quality of the prediction. The low RMSD obtained with both methods reinforce the likelihood of the approach in producing predictions close to a native conformation.[58]

Molecular Dynamics Simulation Analysis

MD determines, in silico, how motion can affect protein–peptide interactions, the stability of those complexes, and their conformation variations due to binding.[73,74] The first parameter analyzed from the MD trajectories was the RMSD, which measures the average distance between a group of atoms.[75] Analysis of this parameter during the simulation time allows us to determine the level of protein or peptide conformational changes. Monitoring the RMSD of the protein (left Y-axis in Figure , black lines) can provide an understanding of the structural conformation of the studied protein, giving insights into the stability. Ligand RMSD (right Y-axis in Figure , blue and red lines) is an indication of ligand stability with respect to the protein, as well as the evolution of its internal conformation. Promising inhibitors should have low RMSDs during the simulation.[76]

Figure 3

RMSD fluctuations of the original 22-mer hACE2 α1 helix (A), NYBSP-4 stapled control (B), modification 15 (C), and modification 11 (D) spike protein of SARS-CoV-2 complexes. The 500 ns MD simulations were monitored with the first frame as a reference. RMSD based on carbon α of the protein (black) (left Y-axis) and ligand RMSD (right Y-axis). Lig fit Prot in blue, and Lig fit Lig in red.

Based on the RMSDs values and RMSD standard deviation for the complexes composed of the spike protein and top-ranked docking candidates (Table S1, SI), the complex S protein–modification 15 (best docking prediction from PATCHDOCK web server) showed the best results (Figure ). As it can be seen in Figure and Table S1 (SI), when analyzed together, the ligand RMSD (“Lig fit Prot” and “Lig fit Lig”) and the standard deviation values for modification 15 complex presented the lowest values. As presented in the plots in Figure , this complex displays a very stable profile along the 500 ns trajectory. Modification 11 also displays a stable profile with low values for ligand RMSD and ligand RMSD standard deviation. The protein RMSD (Table S1 and Figure ) for these complexes reveals a stable profile, even though considering the average value for this parameter, the complexes with modification 15 and modification 11 do not present the lowest values (see Table S1 (SI), Prot CA column). Furthermore, these cited complexes perform better in the RMSD overall analysis than the original 22 residues from the hACE2 α1 helix, NYBSP-4, and the modification 3 complex (best docking prediction from ZDOCK web server, Figure S4, SI).

The Lig fit Prot (blue lines in Figure ) represents the RMSD of a ligand when the protein–ligand complex is first aligned on the protein backbone and then the RMSD of the ligand heavy atoms is measured. If the values observed are significantly larger than the RMSD of the protein, then it is likely that the ligand has diffused away from its initial binding,[76] which can be observed for modification 14, original 3, and modification 9 (Figure S4, SI). In the case of modification 15, this parameter suffers fluctuations during the simulation, especially at the beginning of the trajectory until the equilibration of the system; then, it maintains stability throughout the end of the simulation.

Furthermore, the Lig fit Lig (red lines in Figure ) shows the RMSD of a ligand that is aligned and measured just on its reference conformation (frame 1 in this case). Regarding this parameter, considering both the average Lig fit Lig RMSD value and the standard deviation, modification 15 and modification 11 complexes showed the best results of all studied complexes. Notably, for these cited complexes, the Lig fit Lig RMSD reveals a stable profile of its ligand atoms internal fluctuations (3.11 ± 0.23 Å for modification 15 and 3.20 ± 0.48 Å for modification 11 against 4.10 ± 0.51, 3.87 ± 1.00, and 3.76 ± 0.42 Å for hACE2 control, NYBSP-4, and modification 3, respectively). The RMSD and standard deviation values for the others studied peptides are presented in the SI (Table S1 and Figure S4).

The superior stability, based on ligand RMSDs values, especially for modification 15 but also observed for the other stapled peptides tested (Table S1, SI), compared to the hACE2 control, can be attributed to the staples present in these structures. The presence of the two staples resulted in the maintenance of the α-helical character (key for the interaction of hACE2 with the RBD[7]) along the trajectory. The ability to maintain the secondary structure during most of the simulation time, despite some occasional fluctuations, can be observed in the movie generated from the MD simulation trajectory of modification 15 (Movie S1, SI). The same feature is not present in the control hACE2 (Movie S2, SI) or even in modification 3 (Movie S3, SI). Additionally, modification 15 also performs better in the RMSD stability analysis compared to modification 10 and original 3 (second best docking candidates from PATCHDOCK and ZDCOK server, respectively, Figure S4 and Table S1, SI). In summary, the RBD-modification 15 complex was the best solution, showing a stable RMSD profile for both protein and ligand evaluations, and reaching the equilibrium state at the end of the MD simulation despite some occasional oscillations along the process. Modification 11 also performed well in the stability analysis, which supports both structures as potential RBD inhibitors.

At the time of this study and to the best of our knowledge, three reports by Curreli et al.,[40] Maas et al.,[41] and Morgan et al.[42] described the design and synthesis of stapled peptides also based on the native hACE2 α1 helix. In the first report, a pseudoviral assay showed that a 30-mer double-hydrocarbon-stapled peptide effectively inhibited viral entry. Maas et al.[41] predicted that 35-mer lactam-stapled peptides can inhibit the S protein RBD–hACE2 complex formation. However, the last study pointed out that mono-stapled peptides can successfully constrained α-helical structure in solution but do not prevent virus internalization. Thus, according to the suggestion by Morgan et al.,[42] the double stapling is a viable approach for inducing α1-helicity, which may prevent virus internalization as described by Curreli et al.[40] when smaller peptides (∼30 mer) are evaluated. Therefore, it is possible to assume that a double-stapled approach, as the one we describe here, compass a superior strategy to develop shorter stapled peptides based on the hACE2 as a new therapeutic to prevent viral infection.

The antiviral activity and superior α-helicity observed in the experimental study of Curreli et al.[40] (which possess a double-stapled) might be partly understood by the results of our MD simulation study. According to Figure , in which the initial and final structures from the MD studies of the control peptide and modification 15 are superimposed, the second stapled at the right side of the designed peptide helps to sustain this part of the structure in place in the complex with the S protein. This prevents the designed peptide from moving around and also to maintain the α-helicity of this side of the structure in place during the simulation. Thus, it can be theorized that the second staple prevents slight rotations of the α helix segments, which can lower affinity. These observations can also be seen in Movies S1, S4, and S5, SI. It is noteworthy to mention that most of our proposed modifications contain lactam bridges, which should be less lipophilic that a staple peptide containing an alkene group.

Figure 4

Superimposition of the initial structure (red) and final structure (blue) after 500 ns of simulation of the SARS-CoV-2 spike protein bound to (A) control α1 helix and (B) designed peptide modification 15.

Protein–Stapled Peptide Interactions Obtained from MD Simulation

The important interactions (more than 10% persistence) between the RBD of SARS-CoV-2 and modification 15 or 11 from the MD trajectory are shown in Figure . For modification 15/RBD complex, the amino acids Lys417, Glu484, and Tyr505 (RBD) and Glu30, Lys31, and Asp38 (peptide) exhibited the most remarkable interactions over the 500 ns simulation period. For modification 11/RBD complex, the residues Lys417, Gln474, Arg403, Gln493, and Tyr489 (RBD) and Glu30, Gln24, Asp38, and Lys31 (peptide) are particularly important for the interaction between the S protein and this stapled peptide.

Figure 5

Most persistent interactions (>10% of simulation time) from modification 11 (A) and modification 15 (B) with SARS-CoV-2 RBD. Interacting residues from RBD of SARS-CoV-2 (left Y-axis) and from the studied stapled peptides (right Y-axis).

One interesting feature of modification 15 and modification 11 interaction profiles is the persistence of the most important interaction during the simulation time (Figure ). This interaction is represented by a salt–bridge between Lys417 in the RBD and the substituted Glu residue from our stapled peptides (89.72 and 79.40% of persistence in the simulation time for modifications 11 and 15, respectively). This fact is particularly relevant when analyzed according to the report of Ghorbani et al.[77] In such a study, the binding free energy for the SARS-CoV-2 and RDB interaction were decomposed into a per-residue-based binding energy to find the residues that contribute strongly to the binding and are responsible for the higher binding affinity of SARS-CoV-2 for the hACE2 compared to SARS-CoV. Among their findings, the authors point out that from all of the interface residues, Lys417 had the highest contribution to the total binding energy (−12.34 ± 0.23 kcal/mol) by interacting with Asp30 from hACE2 (also via a salt–bridge). The same salt–bridge interaction between Lys417 and Asp30 from our control hACE2 shows less than half of the persistence value (35.22%) of our studied stapled peptides. Regarding the control NYBSP-4, the persistence for the interaction Lys417 and Asp30 is 74.75%. This observation directed us to conclude that: (i) the presence of two staples in modification 15, modification 11, and the control NYBSP-4 helps to maintain the bioactive conformation allowing the interacting residues to stay in the right position to determine interactions with their partner residues in the viral RBD. This ability is lost in the nonstapled control hACE2 used here; (ii) substitution of Asp30 in modifications 15 and 11 by a Glu30 is effective in placing this residue closer to its partner interacting residue (Lys417). This substitution leads to an increased interaction persistence along the simulation compared to the controls, which possess an Asp30 in the same position (Figure ). The interaction profile of the controls is available in Table S2, SI.

Figure 6

Detailing of the aspartic acid 30 substitution by glutamic acid in modification 15. Spike protein (RBD) of SARS-CoV-2 (yellow). Modification 15 (blue) and hACE2 control (cyan) overlaid. The nonpolar hydrogens were hidden for visualization purpose.

Stability and Profile Interaction in the Presence of E484K RBD Mutation

The second most important interaction of modification 15 is a salt–bridge between Glu484 (RBD) and Lys31. According to the CDC, laboratory studies suggest that specific monoclonal antibody treatments may be less effective for treating cases of COVID-19 caused by variants with the E484K substitution in the spike protein.[51,78−80] Additionally, there are three circulating variants of concern, which present such mutation: B.1.351 lineage (also known as (a.k.a.) 20H/501Y.V2) first described in South Africa;[81] P.1 lineage (a.k.a. 20J/501Y.V3) first reported in Japan in four travelers from Brazil;[82,83] and finally, the B.1.1.7 lineage (a.k.a. 20I/501Y.V1) first detected in the United Kingdom[84] (this mutation is found in some but not all sequences[51]). Therefore, we decided to analyze the interaction profile and stability of modification 15 in the presence of the spike protein RBD presenting the E484K mutation. Even though this specific interaction does not show a relevant persistence (>10%) in the interaction profile of modification 11 and the viral RBD, we decided to analyze the behavior of this peptide in the presence of the cited mutation. This was done with the purpose to assess if the stable profile demonstrated by this stapled peptide in the simulation would be somehow affected by the presence of SARS-CoV-2 variants presenting this mutation.

The ligand RMSD values calculated from the MD trajectory of E484K-mutated RBD/modification 15 and 11 reveal stable profiles of the peptides with respect to the protein and to their internal fluctuations (Figure S5). Specifically, modification 15 presented the following values for the Lig fit Prot and Lig fit Lig RMSD: 4.26 ± 0.43 and 3.46 ± 0.25 Å. For modification 11, the values for the same parameters are: 4.76 ± 0.59 and 3.19 ± 0.51 Å. These ligand RMSD and standard deviation values for the mutated complexes are comparable to the values seen in the original RBD and the same stapled peptides. As presented in Table S1 (SI), the original RBD/modification 15 shows 4.12 ± 0.49 and 3.11 ± 0.23 Å and RBD/modification 11 shows 4.57 ± 0.72 and 3.20 ± 0.48 Å for the same parameters. Overall, both mutated RBD and stapled peptide complexes present better stability profiles than the hACE2 α1 helix control and the stapled NYBSP-4 control (RMSD plot for the tested stapled peptides and mutated RBD trajectories are presented in Figure S5, and RMSD for the controls is shown in Figure and Table S1).

Important interactions between the mutated RBD of SARS-CoV-2 and the tested peptides are shown in Figure . The interaction profile of modification 15 and mutated RBD shows a shift in the character of the major interaction types compared to the trajectory of the natural RBD/modification 15 complex (see Figures B and 7B). In the trajectory resulting from the MD simulation of modification 15 and the natural RBD, the most important interaction type is represented by salt–bridge interactions with high persistence (Lys417/substituted Glu: 79.40%; Glu484/Lys31: 66.80%). However, in the E484K-mutated RBD and modification 15 complex, the interactions of this specific type have reduced persistence (Figure B). To compensate this reduction, increased number of hydrogen-bond (H-bond) interactions are displayed presenting lower persistence. This increase in the H-bond interactions must contribute to the good stability of modification 15/E484K-mutated RBD complex, as demonstrated by the stable profile seen in the RMSD parameters. Note that even though the salt–bridge interactions have a lower influence in the interaction profile of modification 15 and mutated RBD, the interaction between Lys417 and substituted Glu30 suffers a low reduction in the persistence displaying a value of 72.85%.

Figure 7

Most persistent interactions (>10% of simulation time) from modification 11 (A) and modification 15 (B) with viral RBD containing the mutation E484K. Interacting residues from RBD of SARS-CoV-2 (left Y-axis) and from modification 15 (right Y-axis).

Furthermore, for modification 11-mutated RBD complex, some interactions were lost, and some were created (see Figures A and 7A), leading to an interaction profile with the same number of contacts presented by the natural complex. The most remarkable difference is the appearance of a new salt–bridge between Arg403 and Glu37 in the mutated complex with a persistence of 86.20% (Figure A). Similarly to modification 15, the analysis of persistence for the salt–bridge interaction between Lys417 and substituted Glu in modification 11 shows a maintenance of this important interaction. Interestingly, the mutated complex shows a slightly increasing in persistence going from 89.72% in the original RBD/peptide interaction to 94.01% in the mutated complex. Together these results lead us to hypothesize that the E484K mutation present in different SARS-CoV-2 lineages will not hamper the ability of modification 15 or modification 11 to bind the viral RBD.

Binding Free Energy Analysis

Molecular docking provides a fast/efficient way to predict protein–protein interactions (PPIs) and to rank them accordingly to the docking score.[85] However, most scoring functions used by docking programs are developed to enhance computational efficiency and, consequently, present a reduced accuracy of prediction.[86] Methods that combine molecular mechanics energy and implicit solvent models such as MM/GBSA are theoretically more rigorous than docking scoring functions and are powerful tools to predict the binding affinities for protein–peptide systems.[86,87] To predict the binding affinity of our proposed peptides to the RBD of SARS-CoV-2 and to compare this predicted binding affinity with the experimental validated NYBSP-4, the MM/GBSA calculation was employed. The results of such analysis are presented in Table for the controls, modification 15, and modification 11, and in Table S3 (SI) for the other best docking ranked stapled peptides.

According to the results presented in Table , modification 11 is predicted to have better binding affinity to the viral RBD. Interestingly, this modification exhibits superior predicted affinity compared to the experimentally validated NYBSP-4. In addition, NYBSP-4 has a slightly better binding affinity compared to modification 15. However, NYBSP-4 is a longer peptide (30-mer stapled peptide) than both modification 15 and modification 11 (22-mer stapled peptides). Moreover, modifications 15 and 11 being smaller structures can be easier and cheaper to prepare. Further, both of our proposed stapled peptides should provide a better water solubility profile.

Overall, modification 15 and modification 11 showed good stability results for the parameters analyzed in this study. Modification 11 showed better results in the binding free energy prediction calculated using the MM-GBSA method. Furthermore, modification 15 displays a predicted binding affinity comparable to the control NYBSP-4 when analyzed accordingly to their size. Additional experimental studies are required to determine the activity of these peptides against SARS-CoV-2. Nevertheless, our in silico study provide initial evidence to confirm that double-stapled peptides derived from the hACE2 α1 helix present a therapeutic strategy worthy of further investigation.

Conclusions

Through the use of structure-based design, it was possible to propose 18 (22-mer) stapled peptides derived from the hACE2 α1 helix. These peptides were designed to retain the α-helical character of the natural structure, to have enhanced binding affinity between the peptides and SARS-CoV-2-RBD, to avoid disrupting existing favorable interactions, and to display a better solubility profile compared to bigger stapled peptides available in the literature. Furthermore, using docking techniques and a refinement protocol, we selected the most promising binders to perform further analysis using MD simulation and MM-GBSA free energy of binding prediction. According to our study, we identified modifications 11 and 15 as our best candidates. We predict that these peptides can bind to SARS-CoV-2-RBD with potency higher than or similar to the control NYBSP-4 (experimentally proven SARS-CoV-2-RBD 35-mer peptide binder) showing the advantages of being smaller peptides. Our most promising stapled peptides showed stable profiles in the MD simulation and could retain important interactions with the RBD even in the presence of the E484K RBD mutation. Moreover, our study provides valuable information for the rational design and development of stapled peptide inhibitors against SARS-CoV-2 infection.