Computational Design of 25-mer Peptide Binders of SARS-CoV-2.

SARS-CoV-2 is the novel coronavirus causing the COVID-19 pandemic. To enter human cells, the receptor-binding domain (RBD) of the S1 subunit of SARS-CoV-2 (SARS-CoV-2-RBD) initially binds to the peptidase domain of angiotensin-converting enzyme 2 receptor (ACE2-PD). Using peptides to inhibit SARS-CoV-2-RBD binding to ACE2 is a potential therapeutic solution for COVID-19. A previous study identified a 23-mer peptide (SBP1) that bound to SARS-CoV-2-RBD with comparable KD to ACE2. We employed computational protein design and molecular dynamics (MD) to design SARS-CoV-2-RBD 25-mer peptide binders (SPB25) with better predicted binding affinity than SBP1. Using residues 21-45 of the α1 helix of ACE2-PD as the template, our strategy is employing Rosetta to enhance SPB25 binding affinity to SARS-CoV-2-RBD and avoid disrupting existing favorable interactions by using residues that have not been reported to form favorable interactions with SARS-CoV-2-RBD as designed positions. Designed peptides with better predicted binding affinities, by Rosetta, than SPB25 were subjected to MD validation. The MD results show that five designed peptides (SPB25F8N, SPB25F8R, SPB25L25R, SPB25F8N/L25R, and SPB25F8R/L25R) have better predicted binding affinities, by the MM-GBSA method, than SPB25 and SBP1. This study developed an approach to design SARS-CoV-2-RBD peptide binders, and these peptides may be promising candidates as potential SARS-CoV-2 inhibitors.

Introduction

Coronavirus disease (COVID-19) is an ongoing pandemic that has rapidly spread worldwide. COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that is closely related to SARS-CoV in 2002–2003,[1−4] but it can transmit from humans to humans more easily than SARS-CoV.[1,5,6] Similar to other coronaviruses, SARS-CoV-2 contains four major structural proteins including the envelope (E), membrane (M), nucleocapsid (N), and spike (S) proteins.[7,8] The spike proteins of SARS-CoV-2 consist of S1 and S2 subunits that are responsible for the attachment and fusion of the virus with host-cell membranes.[4,9,10] The receptor-binding domain (RBD) of the S1 subunit of SARS-CoV-2 directly binds to the peptidase domain (PD) of the angiotensin-converting enzyme 2 (ACE2) receptor of human cells, while the S2 subunit is responsible for the membrane fusion.[6,11,12]

ACE2 is a type-I membrane metallocarboxypeptidase expressed in the lungs, intestine, heart, and kidneys.[11,13−15] The full-length structure of ACE2 consists of the N-terminal PD and the C-terminal collectrin-like domain, which ends with a single transmembrane helix and an intracellular segment.[11,16,17] The primary physiological role of ACE2 involves in the maturation of angiotensin (Ang), which is a peptide hormone that regulates vasoconstriction and blood pressure. ACE2 is involved in the SARS virus infection via the cell–surface receptor for virus entry into the lungs. SARS-CoV-2 also binds to the ACE2 receptor that facilitates virus entry into the host cells.[11,18−21] The α1 helix of the ACE2 peptidase domain (ACE2-PD) is a main recognition site of the RBD of SARS-CoV-2 (SARS-CoV-2-RBD). The α2 helix and the linker of the β3 and β4 antiparallel strands also have minor contributions to RBD binding.[11,22]

Various potential therapeutics to control viral infections including neutralizing antibodies, small-molecule drugs, and peptides have been explored.[23−27] Disrupting protein–protein interactions of ACE2 and SARS-CoV-2-RBD to prevent virus entry in human cells is a potential therapeutic solution for COVID-19. Because small molecules are often ineffective at disrupting large protein-binding interfaces,[28] inhibiting protein–protein interactions is an important challenge in traditional drug discovery. As alternative solutions, peptides can be used as inhibitors that bind at the interface to disrupt protein–protein interactions because they have larger surface area and more functional groups than small molecules that can imitate and obstruct the native protein–protein interactions[29] as in the case of enfuvirtide, which is a peptide inhibitor that can inhibit HIV entry and has been clinically approved as anti-HIV medicine.[30] Additional advantages of protein therapeutics include high specificity that leads to less potential to interfere with normal biological processes and to cause adverse effects as well as less likely to elicit immune responses in humans.[31]

A previous study showed that antiviral nanoparticles with long and flexible linkers that mimicked heparan sulfate proteoglycans had inhibitory activity against various viruses.[32] Although their applicability to SARS-CoV-2 has not been studied, multiple peptides could potentially be attached to the surfaces of these nanoparticles for multivalent binding to SARS-CoV-2 to increase the binding affinities of these peptides.[33] Furthermore, the peptide inhibitors of SARS-CoV-2 could potentially be used as inhaled therapeutics for topical lung delivery.[34]

Computational techniques have been employed to identify peptides that could potentially bind to SARS-CoV-2-RBD.[33,35−38] Moreover, a previous experimental study found that the 23-mer peptide binder (SBP1) that was derived from α1 helix of ACE2-PD bound to SARS-CoV-2-RBD with KD of 47 nM,[37] which is comparable to KD of ACE2 bound to SARS-CoV-2-RBD (14.7 nM),[39] and it could potentially be used as an inhibitor of SARS-CoV-2. However, the binding affinity of this peptide to SARS-CoV-2-RBD can be further enhanced using computational protein design and molecular dynamics (MD) to improve its effectiveness as a therapeutic agent.

Using the α1 helix (residues 21–45) of ACE2-PD as a design template, this work employed computational protein design (Rosetta) and MD (AMBER) to design SARS-CoV-2 25-mer peptide binders (SPB25) with better predicted binding affinities to SARS-CoV-2-RBD than SBP1. Our design strategy is to increase the binding affinity between SPB25 and SARS-CoV-2-RBD and avoid disrupting existing favorable interactions. The designed peptides with increased predicted binding affinities are promising candidates as potential inhibitors of SARS-CoV-2.

Methods

Structure Preparation

The complex structure of SARS-CoV-2-RBD bound to ACE2 was obtained from the protein data bank (PDB ID: 6M0J).[40] The complex structures of 23-mer peptide (SBP1[37]) (21 IEEQAKTFLDKFNHEAEDLFYQS 43) and 25-mer peptide (SPB25) (21 IEEQAKTFLDKFNHEAEDLFYQSSL 45) bound to SARS-CoV-2-RBD were obtained from the crystal structure of α1 helix of ACE2-PD bound to SARS-CoV-2-RBD. All ionized amino acids of all systems were protonated at the physiological pH (pH 7.4). The LEaP module of AMBER18[41] was employed to construct the final structures of the complexes.

Computational Peptide Design

The structure of the SPB25/SARS-CoV-2-RBD complex was used as a design template for designing SARS-CoV-2-RBD peptide binders. Designed positions were chosen based on the following criteria: (i) they are not residues that were previously reported to form favorable interactions with SARS-CoV-2-RBD to avoid disrupting existing favorable interactions between the peptide and SARS-CoV-2-RBD and (ii) their side chains are in the orientations that could potentially form favorable interactions upon mutations with SARS-CoV-2-RBD to increase the binding affinity between the peptide and SARS-CoV-2-RBD. Each designed position was allowed to be natural amino acids except G and P. The CoupledMoves protocol[42,43] in RosettaDesign module of Rosetta3.11[44] with beta_nov16 energy function was employed to design, repack, and minimize the structure of each designed residue. The neighboring residues within 10 Å of SPB25 were also repacked and minimized. For each design, 400 independent runs were performed, resulting in the total of 400 conformations of designed sequences (some sequences might have multiple conformations). An interface analyzer[45,46] module of Rosetta3.11 was employed to calculate ΔGbind (Rosetta) of each designed conformation in REU. ΔΔGbind (Rosetta) upon mutation was calculated from the difference between the values of ΔGbind (Rosetta) of the designed sequence/conformation and SPB25. From each designed position, the designed peptides with ΔΔGbind (Rosetta) < 0 REU were selected for MD simulations to validate whether their predicted binding affinities by the molecular mechanics–generalized born surface area (MM-GBSA) method[47−49] were better than that of SPB25. The double mutants of SPB25 were also constructed from its single mutants with better predicted binding affinities, by the MM-GBSA method, than SPB25 and subjected to MD simulations to validate whether their predicted binding affinities were better than that of SPB25.

MD Simulations and Analyses

The LEaP module of AMBER18[41] was employed to construct the complex structures of ACE2/SARS-CoV-2-RBD, SBP1/SARS-CoV-2-RBD, SPB25/SARS-CoV-2-RBD, and designed peptides/SARS-CoV-2-RBD as well as peptide structures in isomeric truncated octahedral boxes of TIP3P water molecules with the buffer distance of 13 Å using protein ff14SB[50] and GLYCAM06j-1 force field parameters.[51] To remove unfavorable interactions, each system was minimized using the five-step procedure,[52−57] and all steps included 2500 steps of steepest descent and 2500 steps of conjugate gradient with different restraints on the proteins. At the beginning, the heavy atoms of protein were restrained with a force constant of 10 kcal/(mol Å2), while the hydrogen atoms and water molecules were minimized. In the second, third, and fourth steps, the force constants of 10, 5, and 1 kcal/(mol Å2) were then used to restrain the backbone of the protein, respectively, while the rest of the system was minimized. In the last step, the entire system was minimized without any restraining force.

All systems were simulated under the periodic boundary condition using the GPU (CUDA) version of the PMEMD module.[58−60] The SHAKE algorithm[61] was applied to constrain all bonds involving hydrogen atoms, allowing a simulation time step of 0.002 ps. The temperatures of all systems were controlled using the Langevin dynamics technique[62] with a collision frequency of 1.0 ps–1. All systems were heated from 0 K to the physiological temperature of 310 K for 200 ps in the NVT ensemble, while the backbones of the proteins were restrained with a force constant of 10 kcal/(mol Å2). They were subsequently equilibrated without restraint at 310 K for 300 ps in the NVT ensemble. Finally, they were simulated at 310 K and 1 atm in the NPT ensemble for 100 ns.

In terms of analyses, the root mean square deviation (rmsd) values with respect to the minimized structure were calculated to measure the stability of each system. The last 20 ns trajectories of all systems with stable rmsd values were selected for further analyses. Employing the MM-GBSA method,[47−49] the total binding free energies (ΔGbind (MM-GBSA)) of all systems were calculated to predict the binding affinity between ACE2/SBP1/SPB25/designed peptides and SARS-CoV-2-RBD. The designed peptides with better predicted binding affinity than SPB25 were further analyzed in terms of per-residue free energy decomposition and binding interactions. To identify important binding residues of ACE2/SBP1/SPB25/designed peptides in binding to SARS-CoV-2-RBD, per-residue free energy decomposition was computed. To analyze hydrogen bond interactions, hydrogen bond occupations were calculated. A hydrogen bond was considered to occur if the following criteria were met: (i) a proton donor–acceptor distance ≤ 3.5 Å and (ii) a donor-H-acceptor bond angle ≥ 120°.[52−54,63] Hydrogen bond occupations were classified into four levels: (i) strong hydrogen bonds (hydrogen bond occupations >75%), (ii) medium hydrogen bonds (75% ≥ hydrogen bond occupations >50%), (iii) weak hydrogen bonds (50% ≥ hydrogen bond occupations >25%), and (iv) very weak hydrogen bonds (25% ≥ hydrogen bond occupations >5%).[53−55] Define secondary structure of protein was used to calculate percent helicity of each peptide, and the percent helicity was computed from the summation of the percentage of α- and 310-helix structures.[64]

Results

Computational Design of SARS-CoV-2-RBD Peptide Binders

The structure of the 25-mer peptide, SPB25 (21 IEEQAKTFLDKFNHEAEDLFYQSSL 45) bound to SARS-CoV-2-RBD was obtained from the crystal structure of the α1 helix of ACE2-PD bound to SARS-CoV-2-RBD (PDB ID: 6M0J),[40] and it was used as a design template (Figure ). Residues I21 to L45 were chosen because residues in this region can form favorable interactions with SARS-CoV-2-RBD.[33,37] Employing Rosetta, our design strategy is to enhance the binding affinity between SPB25 and SARS-CoV-2-RBD and avoid disrupting previously reported existing favorable interactions between SPB25 residues (residues 21–45 of ACE2-PD) and SARS-CoV-2-RBD.[33,40] Designed positions were selected based on the following criteria: (i) they are not residues that were previously reported to form favorable interactions with SARS-CoV-2-RBD to avoid disrupting existing favorable interactions between the peptide and SARS-CoV-2-RBD and (ii) their side chains are in the orientations that could potentially form favorable interactions, upon mutations, with SARS-CoV-2-RBD to further enhance the binding affinity between the peptide and SARS-CoV-2-RBD. I1(21), E3(23), F8(28), L19(39), and L25(45) were selected based on these criteria. Each designed position was allowed to be natural amino acids except G and P because G and P occur infrequently in an α helix. Moreover, P can cause a destabilizing kink in a helix structure.[65] Rosetta gave the total of 67 designed peptides with single mutation (Table S1). Thirteen designed peptides have better ΔGbind (Rosetta) than SPB25 (ΔΔGbind (Rosetta) < 0 REU) and were selected for MD to validate whether their predicted binding affinities by the more accurate MM-GBSA method[47−49] (ΔGbind (MM-GBSA)) were better than that of SPB25 (ΔΔGbind (MM-GBSA) < 0 kcal/mol). These designed peptides are SPB25I1D, SPB25I1H, SPB25I1S, SPB25I1T, SPB25E3A, SPB25E3D, SPB25F8D, SPB25F8N, SPB25F8R, SPB25F8W, SPB25L25M, SPB25L25R, and SPB25L25V.

Figure 1

Structure of the SPB25/SARS-CoV-2-RBD complex that was used as a design template. SPB25 and SARS-CoV-2-RBD are colored in yellow and cyan, respectively. The designed positions (I1, E3, F8, L19, and L25) are labeled in red.

Validation by MD

To determine whether the designed peptides have better predicted binding affinity, as calculated by the MM-GBSA method, than SPB25, MD was performed on the structures of SPB25 and all 13 designed peptides in complex with SARS-CoV-2-RBD. Moreover, MD was also performed on the ACE2/SARS-CoV-2-RBD and SBP1/SARS-CoV-2-RBD complexes, and their predicted binding affinities were also used for comparison because the experimental KD values of ACE2 (14.7 nM) and SBP1 (47 nM) are available.[37,39] To analyze the stabilities of all systems, the rmsd values of all atoms and backbone atoms were calculated. Figure S1 shows that the simulations of all systems were likely to reach equilibrium around 80 ns. Therefore, the 80–100 ns trajectories of all systems were selected for further analyses.

To predict the binding affinities of all systems, the MM-GBSA method was used to calculate ΔGbind (MM-GBSA) of all systems during the 80–100 ns trajectories (Table ). The values of ΔGbind (MM-GBSA) of ACE2/SARS-CoV-2-RBD, SBP1/SARS-CoV-2-RBD, and SPB25/SARS-CoV-2-RBD complexes are −71.2 ± 0.4, −55.1 ± 0.4, and −60.3 ± 0.4 kcal/mol, respectively. Out of 13 designed peptides, 3 designed peptides such as SPB25F8N, SPB25F8R, and SPB25L25R have better ΔGbind (MM-GBSA) than SPB25 with ΔΔGbind (MM-GBSA) of −8.0 ± 0.6, −3.3 ± 0.6, and −3.4 ± 0.6 kcal/mol, respectively. Based on these three designed peptides with single mutation, the designed peptides with double mutations such as SPB25F8N/L25R and SPB25F8R/L25R were also constructed using Rosetta and subjected to MD validation. Their ΔGbind(Rosetta) and ΔGbind(MM-GBSA) are better than those of SPB25. The values of ΔΔGbind(MM-GBSA) of SPB25F8N/L25R and SPB25F8R/L25R are −6.0 ± 0.6 and −5.2 ± 0.6 kcal/mol, respectively. Moreover, the predicted binding affinity of these five designed peptides are better than that of SBP1 that is the experimentally proven peptide binder of SARS-CoV-2-RBD.[37]

Table 1

Predicted Binding Free Energies to SARS-CoV-2-RBD of ACE2, SBP1, SPB25, and Designed Peptides that were Selected for MD Simulations, as Calculated with Rosetta and the MM-GBSA Method

system	ΔΔGbind (Rosetta)a (REU)	ΔGbind (MM-GBSA) (kcal/mol)	ΔΔGbind (MM-GBSA)b (kcal/mol)
ACE2	–	–71.2 ± 0.4	–10.9 ± 0.6
SBP1	–	–55.1 ± 0.4	5.2 ± 0.6
SPB25	–	–60.3 ± 0.4	0.0 ± 0.6
SPB25I1D	–0.8	–52.0 ± 0.4	8.3 ± 0.6
SPB25I1H	–0.1	–56.7 ± 0.4	3.6 ± 0.6
SPB25I1S	–0.7	–53.4 ± 0.3	6.9 ± 0.5
SPB25I1T	–0.5	–53.0 ± 0.4	7.3 ± 0.6
SPB25E3A	–0.1	–50.6 ± 0.6	9.7 ± 0.7
SPB25E3D	–0.2	–57.3 ± 0.4	3.0 ± 0.6
SPB25F8D	–2.3	–49.1 ± 0.4	11.2 ± 0.6
SPB25F8N	–0.2	–68.3 ± 0.4	–8.0 ± 0.6
SPB25F8R	–1.2	–63.6 ± 0.4	–3.3 ± 0.6
SPB25F8W	–0.1	–34.9 ± 0.3	25.4 ± 0.5
SPB25L25M	–0.2	–53.5 ± 0.4	6.8 ± 0.6
SPB25L25R	–1.5	–63.7 ± 0.4	–3.4 ± 0.6
SPB25L25V	–0.1	–57.1 ± 0.4	3.2 ± 0.6
SPB25F8N/L25R	–2.5	–66.3 ± 0.4	–6.0 ± 0.6
SPB25F8R/L25R	–0.8	–65.5 ± 0.4	–5.2 ± 0.6

The difference between ΔGbind(Rosetta) of a system and that of SPB25.

The difference between ΔGbind(MM-GBSA) of a system and that of SPB25.

In terms of binding free energy components (Figure ), the main components contributing to the predicted binding affinities of ACE2, SBP1, SPB25, SPB25F8N, SPB25F8R, SPB25L25R, SPB25F8N/L25R, and SPB25F8R/L25R to SARS-CoV-2-RBD are the electrostatic interaction terms because they have the most favorable values. The van der Waals energy and non-polar solvation terms also favor their binding to SARS-CoV-2-RBD. The polar solvation terms have unfavorable contribution to the predicted binding affinity.

Figure 2

Binding free energy components of ACE2/SARS-CoV-2-RBD, SBP1/SARS-CoV-2-RBD, SPB25/SARS-CoV-2-RBD, and designed peptides/SARS-CoV-2-RBD. (A) ΔGbind (MM-GBSA), (B) van der Waals energy, (C) electrostatic interaction, (D) polar solvation, and (E) nonpolar solvation.

As shown in Figure , SPB25F8N is the designed peptide with the best predicted binding affinity, as measured by ΔGbind(MM-GBSA) (−68.3 ± 0.4 kcal/mol). Its predicted binding affinity is better than those of SBP1 and SPB25 by 13.2 ± 0.6 and 8.0 ± 0.6 kcal/mol, respectively, and worse than ACE2 by 2.9 ± 0.6 kcal/mol. The favorable binding of SPB25F8N to SARS-CoV-2-RBD is mostly caused by the increase in the favorable van der Waals energy and non-polar solvation terms, as compared to those of SBP1 and SPB25. The favorable electrostatic interaction term of SPB25F8N is better than that of SBP1 but worse than that of SPB25. The unfavorable polar solvation term of SPB25F8N is better than that of SPB25 but worse than that of SBP1. The predicted binding affinities of SPB25F8R, SPB25L25R, SPB25F8N/L25R, and SPB25F8R/L25R are also better than those of SBP1 and SPB25. The favorable binding of these peptides to SARS-CoV-2-RBD is mostly caused by the increase in the favorable van der Waals energy and non-polar solvation terms as well as the decrease in unfavorable polar solvation terms, as compared to those of SBP1 and SPB25. The favorable electrostatic interaction terms of these designed peptides are worse than those of SBP1 and SPB25.

Identification of Important Binding Residues

Per-residue free energy decomposition was calculated to identify important binding residues of ACE2, SBP1, SPB25, and designed peptides to SARS-CoV-2-RBD, as shown in Figure . In this study, an important binding residue was defined to be a residue with the total energy contribution better than −1.0 kcal/mol.[53] Focusing on residues 21–45 of the α1 helix of ACE2, Q24, T27, F28, D30, K31, H34, and Y41 were predicted to be the important binding residues to SARS-CoV-2-RBD. The important binding residues of SBP1 were predicted to be similar to those of SPB25, and they are Q4 (24), T7 (27), F8 (28), D10 (30), K11 (31), H14 (34), Y21 (41), and Q22 (42). Most of the important binding residues of these peptides were also predicted to be similar to the important binding residues of the α1 helix (residues 21–45) of ACE2 with the addition of Q22 (42) in SBP1 and SPB25. Overall, the numbers of important binding residues of SPB25F8N (9), SPB25F8R (10) SPB25L25R (10), SPB25F8N/L25R (11), and SPB25F8R/L25R (8) were predicted to be relatively similar to those of SBP1 (8) and SPB25 (9). Moreover, five residues of all designed peptides were predicted to have high binding affinity (better than −2.0 kcal/mol): Y21 has the highest binding affinity to SARS-CoV-2-RBD, followed by Q4, T7, N8/R8, and K11, respectively. Additionally, F20 of SPB25F8N, H14 of SPB25F8R, Q22 of SPB25L25R, and H14 and S24 of SPB25F8N/L25R were also predicted to have high binding affinity to SARS-CoV-2-RBD.

Figure 3

Per-residue free energy decomposition of ACE2, SBP1, SPB25, and designed peptides in binding to SARS-CoV-2-RBD. The residue number of ACE2 is in parentheses.

In terms of residue 8, the F8N and F8R mutations were predicted to substantially increase the total energy contribution of this residue from −1.4 kcal/mol in SPB25 to −4.3 and −5.7 kcal/mol in SPB25F8N and SPB25F8R, respectively. In addition, the total energy contributions of other residues such as E17, F20, and Y21 were substantially increased from 0.0, 0.0, and −3.0 kcal/mol in SPB25 to −1.4, −3.0, and −4.5 kcal/mol in SPB25F8N and to −1.3, −1.5, and −4.8 kcal/mol in SPB25F8R, respectively. The total energy contributions of Q4 was also increased from −2.5 kcal/mol in SPB25 to −3.9 kcal/mol in SPB25F8N. Additionally, the total energy contributions of D10, K11, H14, and L25 were increased from −1.4, −4.1, −1.8, and −0.4 kcal/mol in SPB25 to −2.2, −4.8, −3.3, and −1.5 kcal/mol in SPB25F8R, respectively. For residue 25, the L25R mutation was predicted to increase the total energy contribution from −0.4 kcal/mol in SPB25 to −1.6 kcal/mol in SPB25L25R. Moreover, the total energy contributions of E3, Q4, T7, and Y21 were also increased from −0.9, −2.5, −2.8, and −3.0 kcal/mol in SPB25 to −1.4, −3.7, −3.5, and −3.3 kcal/mol in SPB25L25R, respectively.

For the double mutations at residues 8 and 25, the F8N/L25R and F8R/L25R mutations were predicted to increase the total energy contributions of residues 8 and 25 from −1.4 and −0.4 kcal/mol in SPB25 to −2.6 and −1.4 kcal/mol in SPB25F8N/L25R, and to −4.4 and −1.1 kcal/mol in SPB25F8R/L25R, respectively. The total energy contributions of Q4, H14, E17, F20, Y21, and S24 were also increased from −2.5, −1.8, 0.0, 0.0, −3.0, and 0.1 kcal/mol in SPB25 to −3.1, −2.1, −1.2, −1.7, −6.8, and −2.2 kcal/mol in SPB25F8N/L25R, respectively. Furthermore, the total energy contributions of Q4, K11, E17, and Y21 were increased from −2.5, −4.1, 0.0, and −3.0 kcal/mol in SPB25 to −2.7, −4.5, −2.0, and −4.8 kcal/mol in SPB25F8R/L25R, respectively.

Hydrogen Bond and Pi Interactions

Hydrogen bond occupations, pi–pi, cation–pi, and sigma–pi interactions of all systems were analyzed to identify important hydrogen bonds and pi interactions for the binding of ACE2, SPB1, SBP25, and designed peptides to SARS-CoV-2-RBD, as shown in Table and Table S2. Key binding interactions are shown in Figure . Y83 and K353 of ACE2 were predicted to form strong hydrogen bonds with N487 and the backbone of G502 of SARS-CoV-2-RBD, respectively. S19, Q24, and D355 of ACE2 were also predicted to form medium hydrogen bonds with the backbone of A475, N487, and T500 of SARS-CoV-2-RBD. Furthermore, weak hydrogen bonds were predicted to form between E35, E37 of ACE2 and Q493, Y505 of SARS-CoV-2-RBD, respectively. Moreover, ACE2 was predicted to additionally form 25 very weak hydrogen bonds with SARS-CoV-2-RBD. In terms of pi interactions of the ACE2/SARS-CoV-2-RBD complex, pi–pi interactions were predicted to form between Y83 of ACE2 and F486 of SARS-CoV-2-RBD. H34 and R393 of ACE2 were predicted to form cation–pi interactions with K417 and Y505 of SARS-CoV-2-RBD, respectively. Additionally, K353 of ACE2 was predicted to form sigma–pi interactions with Y505 of SARS-CoV-2-RBD.

Table 2

Numbers of Hydrogen Bonds and Pi Interactions of ACE2, SBP1, SPB25, and Designed Peptides Contributing to SARS-CoV-2-RBD Binding

	number of hydrogen bonds					pi interaction
system	strong	medium	weak	very weak	residue that forms a hydrogen bond with SARS-CoV-2-RBD using its backbone or sidechain	pi–pi	cation–pi	sigma–pi
ACE2	2	3	2	25	S19, Q24, D30, K31, H34, E35, E37, Y41, Q42, Y83, N330, K353, D355	Y83···F486	H34···K417	K353···Y505
							R393···Y505
SBP1	1	2	11	14	Q4, D10, K11, H14, E15, E17, D18, Y21, S23		K11···Y489
							H14···K417
SPB25	1	4	11	20	Q4, D10, K11, H14, E15, D18, Y21, Q22	Y21···Y505	H14···K417
SPB25F8N	0	5	9	28	E3, Q4, N8, D10, K11, H14, E15, E17, D18, F20, S24, L25		K11···F486	Y21···Y449
							K11···Y489	Y21···G496
							H14···K417
SPB25F8R	2	3	3	28	E2, Q4, R8, D10, K11, H14, E17, D18, Y21, S24	Y21···Y449	K11···Y489
SPB25L25R	2	4	12	18	Q4, D10, K11, E15, D18, Y21, Q22, R25	Y21···Y505	H14···K417
SPB25F8N/L25R	0	6	3	25	Q4, N8, D10, H14, E17, D18, Y21, S24, R25	F20···Y505	K11···Y489	K11···Y489
							H14···K417
							Y21···R403	Y21···Y505
							R25···Y449
SPB25F8R/L25R	3	3	6	16	R8, D10, K11, E15, E17, D18, Y21, Q22, R25		K11···Y489	K11···Y489
							H14···K417

Figure 4

Key binding interactions between SARS-CoV-2-RBD (cyan) and (A) ACE2, (B) SBP1, (C) SPB25, (D) SPB25F8N, (E) SPB25F8R, (F) SPB25L25R, (G) SPB25F8N/L25R, or (H) SPB25F8R/L25R. The structures of SBP1, SPB25, and designed peptides (yellow) were superimposed with ACE2 (gray). Key hydrogen bonds and salt bridges (hydrogen bond occupations >25%) are shown in blue dashed lines. These structures are the structures closest to the average structures from the 80–100 ns MD trajectories.

The total numbers of predicted hydrogen bonds and pi interactions of SPB25 are more than those of SBP1, supporting the binding free energy result that predicted SPB25 to bind better to SARS-CoV-2-RBD than SBP1. In terms of hydrogen bond interactions, the predicted numbers of strong and weak hydrogen bonds of SPB25 are similar to those of SBP1, while those of medium and very weak hydrogen bonds of SPB25 are more than those of SBP1. Y21 of SPB25 and SBP1 were predicted to form strong hydrogen bonds with the backbone of G502 of SARS-CoV-2-RBD. SPB25 and SBP1 were also predicted to form 11 weak hydrogen bonds with SARS-CoV-2-RBD. These hydrogen bonds were predicted to form between SARS-CoV-2-RBD and D10, K11, E15, D18, and Q22 of SPB25 and between SARS-CoV-2-RBD and Q4, D10, E15, E17, and S23 of SBP1. D15, D18, and Q22 were predicted to form four medium hydrogen bonds in the SPB25/SARS-CoV-2-RBD complex, while E17 and S23 were predicted to form two medium hydrogen bonds in the SBP1/SARS-CoV-2-RBD complex. Furthermore, residues of SPB25 and SBP1 including Q4, D10, H14, and D18 were predicted to form very weak hydrogen bonds with SARS-CoV-2-RBD. In terms of pi interactions, cation–pi interactions were predicted to form between H14 of SPB25/SBP1 and K417 of SARS-CoV-2-RBD. K11 of SBP1 was predicted to additionally form cation–pi interactions with Y489 of SARS-CoV-2-RBD. SPB25 was also predicted to form an additional pi–pi interaction between Y21 of SBP25 and Y505 of SARS-CoV-2-RBD.

The interactions between the designed peptides and SARS-CoV-2-RBD are displayed in Figure . Overall, the binding positions and orientations of all designed peptides to SARS-CoV-2-RBD are relatively similar to those of ACE2. In terms of the designed peptides with single mutation, the total numbers of predicted hydrogen bonds and pi interactions of SPB25F8N are more than those of SPB25 and SBP1, supporting the binding free energy result that it has better predicted binding affinity to SARS-CoV-2-RBD than SPB25 and SBP1. The mutated residue N8 of SPB25F8N was predicted to form three medium and two very weak hydrogen bonds with N487 and Y489 of SARS-CoV-2-RBD, while no hydrogen bond was predicted to form between F8 of SPB25 and SARS-CoV-2-RBD. Other residues such as E3, Q4, D10, K11, H14, E15, E17, D18, F20, S24, and L25 of SPB25F8N were also predicted to form hydrogen bonds with SARS-CoV-2-RBD. Furthermore, there are three predicted cation–pi interactions (K11···F486, K11···Y489, and H14···K417) and two sigma–pi interactions (Y21···Y449 and Y21···G496) formed between SPB25F8N and SARS-CoV-2-RBD. The total numbers of predicted hydrogen bonds and pi interactions of SBP25F8R are more than those of SBP1 and similar to those of SPB25. However, the predicted number of strong hydrogen bonds of SBP25F8R is more than that of SPB25 and SBP1. The mutated residue R8 of SPB25F8R was predicted to form two strong hydrogen bonds with the backbone of A475 and one very weak hydrogen bond with N487 of SARS-CoV-2-RBD, respectively. Additionally, E2, Q4, D10, K11, H14, E17, D18, Y21, and S24 were predicted to form hydrogen bonds with SARS-CoV-2-RBD. Moreover, there are predicted pi–pi (Y21···Y449) and cation–pi (K11···Y489) interactions between SPB25F8R and SARS-CoV-2-RBD. The numbers of predicted hydrogen bonds and pi interactions of SPB25L25R are more than those of SBP1 and similar to those of SPB25. However, the predicted number of strong hydrogen bonds of SPB25L25R is more than that of SPB25 and SBP1. The mutated residue R25 of SPB25L25R was predicted to form one weak and five very weak hydrogen bonds with T500 of SARS-CoV-2-RBD, while L25 of SPB25 was not predicted to form any hydrogen bonds with SARS-CoV-2-RBD. Furthermore, Q4, D10, K11, E15, D18, Y21, and Q22 of SPB25L25R were predicted to form hydrogen bonds with SARS-CoV-2-RBD. Additionally, there are predicted pi–pi (Y21···Y505) and cation–pi (H14···K417) interactions between SPB25L25R and SARS-CoV-2-RBD.

In terms of the designed peptides with double mutations, the predicted number of hydrogen bonds of SPB25F8N/L25R is higher than that of SBP1 and lower than that of SPB25. However, the predicted number of pi interactions of SPB25F8N/L25R is higher than those of SBP1 and SPB25. The mutated residue N8 of SBP25F8N/L25R was predicted to form medium and very weak hydrogen bonds with N487 and Y489 of SARS-CoV-2-RBD, and the mutated residue R25 was predicted to form very weak hydrogen bonds with G446 (backbone), Y449 and Q498 of SARS-CoV-2-RBD, while F8 and L25 of SPB25 were not predicted to form any hydrogen bonds with SARS-CoV-2-RBD. Q4, D10, H14, E17, D18, Y21, and S24 of SBP25F8N/L25R were also predicted to form hydrogen bonds with SARS-CoV-2-RBD. Furthermore, there are one predicted pi–pi (F20···Y505), four cation–pi (K11···Y489, H14···K417, Y21···R403, and R25···Y449), and two sigma–pi (K11···Y489 and Y21···Y505) interactions between SBP25F8N/L25R and SARS-CoV-2-RBD. The total number of predicted hydrogen bonds of SPB25F8R/L25R is higher than that of SBP1 and lower than that of SPB25. However, the predicted number of strong hydrogen bonds and pi interactions of SPB25F8R/L25R is higher than those of SPB25 and SBP1. The mutated residue R8 of SPB25F8R/L25R was predicted to form a strong hydrogen bond with the backbone of A475 of SARS-CoV-2-RBD, while the mutated residue R25 of SPB25F8R/L25R was predicted to form four very weak hydrogen bonds with N498 and T500 of SARS-CoV-2-RBD. D10, K11, E15, E17, D18, Y21, and Q22 of SPB25F8R/L25R were also predicted to form hydrogen bonds with SARS-CoV-2-RBD. Additionally, there are two predicted cation–pi (K11···Y489 and H14···K417) and one sigma–pi (K11···Y489) interactions between SPB25F8R/L25R and SARS-CoV-2-RBD.

Peptide Helicities

Figure S2 shows the rmsd plots of peptides in water, and Figure shows the percent helicity of each peptide in water and in the peptide/SARS-CoV-2-RBD complex, as calculated from their 80–100 ns trajectories. Because of their high flexibilities, the percent helicities of the N-terminus and C-terminus of each peptide are lower than those of the residues in the middle, except for some residues in the middle of SPB25, SPB25L25R, and SPB25F8N in the peptide/SARS-CoV-2-RBD complexes. Overall, the trends of percent helicities of all peptides in water are slightly higher than those in the peptide/SARS-CoV-2-RBD complex. Moreover, the trends of percent helicities in water of SPB25F8N, SPB25F8N/L25R, SPB25F8R/L25R, SPB25F8R, SPB25L25R., and SPB25 are slightly higher than those of SBP1.

Figure 5

Percent helicities of each peptide in water (A) and in the peptide/SARS-CoV-2-RBD complex (B). SBP1, SPB25, SPB25F8N, SPB25F8R, SPB25L25R, SPB25F8N/L25R, and SPB25F8R/L25R are shown in gray, red, green, blue, purple, pink, and orange, respectively.

Discussion

SARS-CoV-2 is the novel coronavirus responsible for the COVID-19 pandemic that has caused large numbers of cases and deaths worldwide. SARS-CoV-2-RBD initially binds to ACE2-PD to enter human cells. Blocking the SARS-CoV-2-RBD binding with ACE2-PD to prevent coronavirus from entering human cells is a potential therapeutic solution for COVID-19. Using small molecules to disrupt large protein binding interfaces is often ineffective.[28] Alternatively, because peptides have a larger surface, more chemical functionalities, and more similar interactions to the native protein–protein interactions than small molecules, they can be used as inhibitors to disrupt protein–protein interactions at the protein-binding interface.[29]

The previous experimental study designed SBP1 (23-mer peptide) based on residues 21–43 of the α1 helix of ACE2-PD bound to SARS-CoV-2-RBD and found that SBP1 bound to SARS-CoV-2-RBD with KD comparable to ACE2 and could potentially be used as an inhibitor to disrupt the binding of ACE2 to SARS-CoV-2-RBD.[37] To design novel peptides that could potentially bind to SARS-CoV-2-RBD with better predicted binding affinity than SBP1, computational protein design (Rosetta) and MD (AMBER) were employed in this study. Residues 21–45 of the crystal structure of the α1 helix of ACE2-PD bound to SARS-CoV-2-RBD were used as a template (SBP25) because a large number of residues in this region form favorable interactions with SARS-CoV-2-RBD, and some residues could potentially be mutated to increase the binding affinity to SARS-CoV-2-RBD. Our design strategy is to enhance the binding affinity between SPB25 and SARS-CoV-2-RBD and avoid disrupting the previously reported existing favorable interactions between residues 21–45 of the α1 helix of ACE2-PD and SARS-CoV-2-RBD. I1(21), E3(23), F8(28), L19(39), and L25(45) were selected as designed positions because they are not residues that were previously reported to form favorable interactions with SARS-CoV-2-RBD, and their side chains are in the orientations that could potentially form favorable interactions, upon mutations, with SARS-CoV-2-RBD to further enhance their binding affinities. Each designed position was allowed to be natural amino acids except G and P. After computational protein design, Rosetta proposed the total of 67 designed peptides with single mutation. Thirteen designed peptides have better Rosetta binding energies than SPB25 (ΔΔGbind (Rosetta) < 0 REU) and were selected for MD to validate whether their predicted binding affinities, by the more accurate MM-GBSA method, were better than that of SPB25 and SBP1.

The predicted binding affinities of SPB25 (the design template), SBP1, and ACE2 were also calculated by the MM-GBSA method. MD results show that the predicted binding affinity to SARS-CoV-2-RBD of ACE2 is better than that of SBP1 (−55.1 ± 0.4 kcal/mol), supporting the experimental results that ACE2 bound to SARS-CoV-2-RBD (KD = 14.7 nM)[39] with better affinity than SBP1 (KD = 47 nM).[37] The predicted binding affinity to SARS-CoV-2-RBD of SPB25 (−60.3 ± 0.4 kcal/mol) is better than that of SBP1 but worse than that of ACE2, suggesting that SPB25 could potentially bind to SARS-CoV-2-RBD better than SBP1 but worse than ACE2. ACE2 was predicted to bind better than both SPB25 and SBP1 because it is larger than SPB25 and SBP1 and has more residues interacting with SARS-CoV-2-RBD including residues in the α2 helix and the linker of the β3 and β4 antiparallel strands[11,22] in addition to residues 21–45 and 21–43 that were used to derived SPB25 and SBP1, respectively. SPB25 was predicted to bind to SARS-CoV-2-RBD better than SBP1 probably because SPB25 was predicted to have more residues interacting with SARS-CoV-2-RBD than SBP1 as supported by the MD results that the total numbers of predicted hydrogen bonds (involving Q4, D10, K11, H14, E15, D18, Y21, and Q22) and pi interactions (involving H14 and Y21) of SPB25 are more than those of SBP1. Moreover, the results from per-residue free energy decomposition suggest Q4, T7, F8, D10, K11, H14, Y21, and Q22 as important binding residues of SPB25 and SBP1. D18 was also predicted to be an important binding residue of SPB25.

Thirteen designed peptides with single mutations that were predicted by Rosetta to bind to SARS-CoV-2-RBD better than SPB25 (ΔΔGbind (Rosetta) < 0 REU) were subjected to MD, and their predicted binding affinities were further validated by the MM-GBSA method, which gives more accurate predicted binding affinities than Rosetta. Three designed peptides such as SPB25F8N, SPB25F8R, and SPB25L25R were predicted to bind to SARS-CoV-2-RBD better than SPB25 with ΔΔGbind (MM-GBSA) of −8.0 ± 0.6, −3.3 ± 0.6, and −3.4 ± 0.6 kcal/mol. These peptides with single mutation were employed to create the designed peptides with double mutation such as SPB25F8N/L25R and SPB25F8R/L25R using Rosetta, and they were subjected to validation by MD. SPB25F8N/L25R and SPB25F8R/L25R were predicted to bind to SARS-CoV-2-RBD better than SPB25 with ΔΔGbind (MM-GBSA) of −6.0 ± 0.6 and −5.2 ± 0.6 kcal/mol, respectively. Their predicted binding affinities were also better than SPB25F8R and SPB25L25R but worse than SPB25F8N. Most importantly, these five designed peptides (SPB25F8N, SPB25F8R, SPB25L25R, SPB25F8N/L25R, and SPB25F8R/L25R) were predicted to bind to SARS-CoV-2-RBD better than SBP1, suggesting that they should be able to bind to SARS-CoV-2-RBD better than SBP1, experimentally. The ranking of the predicted binding affinities of the designed peptides, SBP1 and SPB25 (best to worst) is SPB25F8N > SPB25F8N/L25R ≈ SPB25F8R/L25R > SPB25F8R ≈ SPB25L25R > SPB25 > SBP1. Additionally, the binding positions and orientations of SPB25 and all designed peptides to SARS-CoV-2-RBD are relatively similar to that of ACE2, suggesting that they could potentially disrupt the binding interactions between SARS-CoV-2-RBD.

Out of the five designed peptides, SPB25F8N is the most promising designed peptide because its predicted binding affinity is better than SPB25, SBP1, and all designed peptides. The results from per-residue free energy decomposition suggest Q4, T7, N8, D10, K11, H14, E17, F20, and Y21 as important binding residues. The enhanced binding affinity of SPB25F8N is probably caused by the increase in the total numbers of predicted hydrogen bonds (involving E3, Q4, N8, D10, K11, H14, E15, E17, D18, F20, S24, and L25) and pi interactions (involving K11, H14, and Y21) of SPB25F8N. Because Q4 (24), D10 (30), E15 (35), E17 (37), D18 (38), and Y21 (41) were also reported to form favorable interactions in the crystal structure of SARS-CoV-2-RBD bound to ACE2,[40] binding of SPB25F8N to SARS-CoV-2-RBD could potentially disrupt the binding of ACE2 and SARS-CoV-2-RBD, supporting the potential use of SPB25F8N as an inhibitor. Moreover, the F8N mutation caused a substantial increase in the total energy contribution of this residue and significant increase in the total energy contributions of other residues such as Q4, E17, F20, and Y21, as compared to those of SPB25. For SPB25F8R, the results from per-residue free energy decomposition suggest Q4, T7, R8, D10, K11, H14, E17, F20, Y21, and L25 as important binding residues, and these results are similar to those of SPB25F8N. The binding affinity of SPB25F8R was predicted to be better than those of SBP1 and SPB25. This result is supported by the fact that its total numbers of predicted hydrogen bonds (involving E2, Q4, R8, D10, K11, H14, E17, D18, Y21, and S24) and pi interactions (involving Y21 and K11) are more than those of SBP1, and the number of predicted strong hydrogen bonds of SPB25F8R (involving R8) is more than that of SPB25 and SBP1. Moreover, the F8R mutation was predicted to cause a substantial increase in the total energy contribution of this residue and increase in the total energy contributions of other residues such as D10, K11, H14, E17, F20, Y21 and L25, as compared to those of SPB25. For SPB25L25R, the results from per-residue free energy decomposition suggest E3, Q4, T7, F8, D10, K11, H14, Y21, Q22, and R25 as important binding residues. The binding affinity of SPB25L25R was predicted to be better than those of SBP1 and SPB25. This finding is supported by the fact that the numbers of predicted hydrogen bonds (involving Q4, D10, K11, E15, D18, Y21, Q22, and R25) and pi interactions (involving H14 and Y21) of SPB25L25R are more than those of SBP1, and the predicted number of strong hydrogen bonds (involving Q4 and Y21) of SPB25L25R is more than that of SPB25 and SBP1. Additionally, the L25R mutation was predicted to cause an increase in the total energy contribution of this residue and other residues such as E3, Q4, T7, and Y21, as compared to those of SPB25.

In terms of designed peptides with double mutation, the predicted binding affinities of SPB25F8N/L25R and SPB25F8R/L25R are comparable. The results from per-residue free energy decomposition suggest Q4, T7, N8/R8, D10, K11, E17, Y21, and R25 as important binding residues for both peptides. Additionally, H14, F20, and S24 were suggested to be important binding residues of SPB25F8N/L25R. Their binding affinities were predicted to be better than those of SBP1 and SPB25 probably because they were predicted to have better interactions with SARS-CoV-2-RBD than SBP1 and SPB25. Specifically, their total numbers of predicted hydrogen bonds and pi interactions are more than that of SBP1, and their predicted number of pi interactions is more than those of SBP1 and SPB25. The predicted binding affinity of SPB25F8R/L25R is better than those of their respective single mutants (SPB25F8R and SPB25L25R) probably because the F8R/L25R mutations caused the increase in the total numbers of hydrogen bonds and pi interactions of SPB25F8R/L25R, as compared to those of SPB25F8R and SPB25L25R. The predicted binding affinity of SPB25F8N/L25R is worse than those of their respective single mutants (SPB25F8N and SPB25L25R) probably because the F8N/L25R mutations caused the decrease in the total number of hydrogen bonds of SPB25F8N/L25R, as compared to those of SPB25F8N and SPB25L25R. Although the double mutation caused the slight increase in the total number of pi interactions, these favorable effects were overcome by the unfavorable effects of the decreased hydrogen bond interactions. Residues that were predicted to form hydrogen bonds with SARS-CoV-2-RBD for both peptides are N8/R8, D10, E17, D18, Y21, and R25. Additionally, Q4, H14, and S24 of SPB25F8N/L25R were predicted to form hydrogen bonds with SARS-CoV-2-RBD, while K11, E15, and Q22 of SPB25F8R/L25R were predicted to form hydrogen bonds with SARS-CoV-2-RBD. K11 and H14 of both peptides were predicted to form pi interactions with SARS-CoV-2-RBD. Moreover, F20, Y21, and R25 of SPB25F8N/L25R were predicted to form pi interactions with SARS-CoV-2-RBD. Furthermore, the F8N/L25R and F8R/L25R mutations were predicted to increase the total energy contributions of residues 8 and 25. Additionally, the total energy contributions of Q4, H14, E17, F20, Y21, and S24 of SPB25F8N/L25R were also enhanced from those of SPB25. The total energy contributions of Q4, K11, E17, and Y21 of SPB25F8R/L25R were also increased from those of SPB25.

In terms of peptide helicities, our results show that the trends of percent helicities in water of SPB25F8N, SPB25F8N/L25R, SPB25F8R/L25R, SPB25F8R, SPB25L25R., and SPB25 are slightly higher than those of SBP1, suggesting that their stabilities in water may be slightly higher than that of SBP1 that is the experimentally proven binder of SARS-CoV-2-RBD. Therefore, these results suggest that our designed peptides should be stable enough to be used as peptide binders of SARS-CoV-2. The trends of percent helicities of these peptides in water are slightly lower than those in the peptide/SARS-CoV-2-RBD complexes probably because some residues in the peptides also interacted with the residues of SARS-CoV-2-RBD and could not participate fully in forming hydrogen bonds important for helix formation.

Using computational protein design and MD techniques, this study developed an approach that can be used to design peptide binders of SARS-CoV-2-RBD. Five designed peptides are promising candidates that could potentially be used as inhibitors of SARS-CoV-2 because their predicted binding affinities are better than that of SBP1 (the experimentally proven peptide binder of SARS-CoV-2). These designed peptides could potentially be employed as inhaled therapeutics for topical lung delivery[34] to prevent the binding of SARS-CoV-2-RBD and ACE2 in the lung. Moreover, one possible way to further increase the binding affinities of peptide binders to SARS-CoV-2 is attaching multiple peptides to the surfaces of nanoparticles for multivalent binding to SARS-CoV-2.[32,33]

Conclusions

Employing computational protein design and MD, we developed an approach to design promising 25-mer peptides, which were derived from residues 21–45 of the α1 helix of ACE2-PD, with better predicted binding affinity to SARS-CoV-2-RBD than SBP1 (the experimentally proven SARS-CoV-2-RBD peptide binder). We employed the design strategy aiming to enhance the binding affinity between SPB25 and SARS-CoV-2-RBD and avoid disrupting existing favorable interactions. This study designed SPB25F8N, SPB25F8N/25R, SPB25F8R/L25R, SPB25F8R, and SPB25L25R that were predicted to bind to SARS-CoV-2-RBD better than SBP1 and SPB25, achieving the design goal. These five designed peptides are promising candidates as potential inhibitors of SARS-CoV-2.