Atomistic insight into the essential binding event of ACE2-derived peptides to the SARS-CoV-2 spike protein

Binding affinity of p4 and p6 towards SARS-CoV-2 RBD

We measured the K_D value for p4 and p6 to the RBD of wild-type SARS-CoV-2 using microscale thermophoresis (MST), and we obtained dissociation constants of 0.93 ± 0.21 µM and 0.45 ± 0.09 µM, respectively (Figure 2A). These results are consistent with the ones obtained from measurements with SARS-CoV-1, showing an increased antiviral activity for p6 compared to p4 indicating a stronger interaction with the RBD (Han et al. 2006). Similar values were also obtained for RBD variants alpha and delta (Supplementary Information, Table 1 and SI Figure 1), showing the broad inhibitor potential for these peptides. We thus studied p6, which is only seven ACE2 residues longer than p4 but shows a doubled binding affinity. Likewise, in a recent work we reported the interaction of recombinant ACE2 and the RBD to be in the same range (Nie et al. 2021), also similar to values found in the literature (Barton et al. 2021; Chan et al. 2020; Lan et al. 2020).

Figure 2:

Affinity and structural characterization of ACE2-derived peptides.

(A) MST binding curves from experiment using different concentrations of ACE2-derived peptides p4, p6 against constant concentration of fluorescently labeled RBD. Error bars: SEM with N = 3. (B) CD spectra for peptides p4 and p6, acquired in 20 mM phosphate buffer, and (C) in presence of 50% of trifluoroethanol as an helicity inducer.

Despite its increased affinity in comparison to p4, the affinity of p6 to the RBD is still approximately half an order of magnitude lower than that of ACE2 for the RBD. Interestingly, the helical content decreases with the C-terminal extension, as shown by comparing the helicity of p4 with p6. This effect is intensified by dissolving the peptides in TFE (Figure 2B and C). To gain further insight into the p6-RBD complex with the intention of proposing variants with enhanced binding, we performed all-atom Molecular Dynamics (MD) simulations of the complex, starting from the ACE2-RBD cryoEM structure (PDB: 7DF4) (Xu et al. 2021). During the simulations we observed that the structure of p6 fluctuates substantially from the bioactive conformation (i.e., the corresponding peptide regions within ACE2 bound to the RBD). This large conformational flexibility can be observed in the multi-modal distribution of the root mean square deviation (RMSD) of the MD-configurations from the initial structure (Figure 3A). We calculated the per-residue contribution to binding by using the molecular mechanics energies with a generalized Born and surface area continuum solvation method (MM/GBSA) (Figure 3B). The average value of the binding enthalpy for the complex was −35 ± 4 kcal/mol. We noted that the truncation of the RBD-interacting interface of ACE2 results in the disruption of hydrophobic interactions, leaving some hydrophobic residues exposed to the solvent, consequently destabilizing the overall peptide secondary structure in solution (Figure 3C).

Figure 3:

Analysis of SARS-CoV-2 RBD-p6 complex stability.

(A) Density plot of RMSD p6 residues (estimated by considering peptide backbone positions after concatenating three independent 250-ns MD simulations). (B) Per-residue contributions to the binding enthalpy and error bars show average and confidence interval (95%) from three replicas concatenated into a single trajectory. Overall binding energy of the complex is −34.8 ± 0.7 kcal/mol. Blue arrows highlight the two anchor residues T27 and K353. (C) Snapshots of p6 (black) bound to the RBD (purple), from MD simulations are shown as ribbons with the initial structure of p6 superimposed in gray. Solvent-exposed hydrophobic amino acids are shown as white sticks. The anchors are shown as blue sticks and RBD residues within 3 Å from the anchors as purple sticks. (D) K _D ^X/K _D ^p6 ratios obtained from MST experiments are shown, where X is either p4, p6(T27A) or p6(K353A) and error bars show the standard deviation. Absolute K _D values are provided in SI Table 3.

Despite the variable observed binding modes of p6, two residues, T27 and K353, stand out in respect to the number of contacts established during the simulations of the complex, which is also reflected in their favorable per-residue energy contributions (Figure 3B). To further evaluate the specific anchoring roles of T27 and K353, we synthesized p6 single point alanine mutants T27A and K353A and evaluated their affinities for the RBD using MST (SI Figure 4). The dissociation constants for each of the mutants is comparable to that of p4, indicating a 3-fold loss in affinity relative to p6 due to the deletion of the 7-residue β-hairpin segment, or one of the two anchor residues alone (Figure 3D).

Figure 4:

Characterization of p6 and variants thereof in complex to the RBD.

(A) Aligned sequences of p6 and p6 variants (Madeira et al. 2019). The black underline represents mutations and the gray underline represents the addition of residues. Some ACE2 residue numbers are shown as a reference. Non-consecutive regions are pointed out with an asterisk. The glycine in bold was added artificially and it is not within the ACE2. (B) ACE2-RBD interface with the RBD shown as purple cartoon and p6 shown as black cartoon. The position of the glycine that joins ACE2 helix-1 and the β-hairpin is highlighted. Residues added to p6B and to p6C, D, E are shown in gray cartoons. The gray-dashed line represents the junction made for p6C, D and E. Some ACE2 residue numbers are shown as a reference. (C) Density plot of p6 and p6 variants RMSD residues. (D) Binding enthalpies of p6 and p6 variants and error bars show average and standard error of the mean from three replicas concatenated into a single trajectory for each system. (E) Density plot of solvent surface accessible area (SASA) for key RBD residues. The gray straight line represents SASA for the key RBD residues (as defined in SI Table 2) in the CryoEM structure of the complex RBD:ACE2 (PDB:7DF4).

Analysis of p6 variants

In an attempt to design ACE2-derived peptides with antiviral activity against SARS-CoV-2 we used p6 as a lead structure peptide and evaluated five variants, denoted p6A to p6E (aligned sequences shown in Figure 4A). p6A and p6B contain mutations of several hydrophobic and non-interacting residues (i.e., residues that participate in native intramolecular ACE2 hydrophobic contacts); the former also incorporates one mutation (E22Y) to favor new interactions and the latter had two additional N-terminal residues, included to promote the stabilization of its α-helix (Petukhov et al. 2009). In variants p6C-E, we aimed at enhancing favorable peptide-RBD interactions by extending the N-terminal region with wild-type (WT) residues from ACE2 helix-2, using two residues from helix-1 to artificially link these two segments (Figure 4B); in particular p6C is the WT sequence with just this N-terminal extension, whereas p6D and p6E have the same length but with further punctual mutations intended to enhance intermolecular interactions. An ACE2 mutational analysis recently reported was taken into consideration in proposing these specific mutations (Chan et al. 2020).

We performed MD simulations of each variant bound to the RBD, so as to compare their structural stability and binding affinity. The RBD structure remained stable in all simulations and the RMSD of the backbone atoms remained below 3.5 Å during the three 250-ns trajectories (SI Figure 2A). The RMSF of the RBD residues (SI Figure 2B) presented a region with more fluctuations (T470-F490) which corresponds to a mobile loop that has been reported to enhance interactions in SARS-CoV-2 compared to SARS-CoV-1 (Wang et al. 2020a; Xu et al. 2021). The N-terminal region of ACE2 and some residues from helix-2 interact with this loop. We note that the three extended variants (p6C, p6D, and p6E) include this RBD binding region, whereas p6A and p6B do not (Figure 4A). The distributions of RMSD values for the five variants indicate that all maintained the initial bioactive conformation better than the WT p6 sequence, with the exception of p6C, which was similar to p6 in terms of its overall flexibility (Figure 4C). The estimated binding enthalpies, in contrast, showed only slight improvements in binding, and only for p6C (Figure 4D). Not only the helix formation is strongly hindered without the ACE2 context, but also the β-hairpin. During the MD runs, p6 maintains residues 351, 352, 355, 356, 357 on an antiparallel β-strand structure, as in the ACE2, 22% of the time on average, p6D and p6E about 10%, and p6A, B and C less than 5% (SI Figure 2C).

Keeping in mind that the goal of designing peptides requires they competitively inhibit ACE2-RBD binding, we considered how well the different binding modes of the variants covered residues known to be involved in the ACE2-RBD interaction. Figure 4E shows the distribution of the total solvent accessible surface area (SASA) of critical RBD residues (SI Table 2) from the combined 750-ns simulation of each variant. In the native bioactive conformation, the corresponding SASA value is approximately 700 Å² (gray reference line in Figure 4E), reflecting the tightly bound ACE2-RBD complex. As the distribution of SASA values shifts closer to this reference line, we consider the variant binding mode(s) to be more effective at “covering” the critical portions of the RBD. The analysis of the SASA distributions shows that all variants of p6 are slightly shifted to the left with respect to the WT p6 sequence, with p6E particularly standing out in this regard.

Together these analyses highlight that the mutation of solvent-exposed hydrophobic residues appears to be most effective at stabilizing the bound conformation of p6 variants (p6A and p6B, Figure 4C), whereas the extension of the N-terminal region showed more of an impact in enhancing the coverage of important RBD residues (in particular p6E, Figure 4E) and in enhancing the calculated total binding enthalpy (variant p6C, Figure 4D).

Based on the above computational analyses of the peptide-RBD simulations, p6D was excluded for the experimental characterization because neither the binding energy of the complex nor the coverage of key RBD residues was enhanced. Even though the binding enthalpies of p6A and p6B to the RBD did not show improvement with respect to p6, both showed promise in terms of conformational stability in the bioactive conformation (Figure 4C) and coverage of key RBD residues (Figure 3E).

We synthesized and characterized p6A, p6B, p6C and p6E by CD spectroscopy and evaluated affinity for the SARS-CoV-2 spike RBD using MST measurements. All peptides were soluble until a 0.25 mM concentration in physiological phosphate buffer at room temperature and the peptide solutions were translucent along the concentrations evaluated in the MST experiments. The helicity in trifluoroethanol (TFE) was lower than that of p6 for all the variants and the binding affinity did not have an improvement compared to p6 (Figure 5, see also SI Table 3). Nevertheless, we observed that variants with higher helical content had higher affinities to the RBD, with the logical exception of the single-point alanine substitutions of the anchor residues.

Figure 5:

K _D ^p6/K _D ^X ratios obtained from MST experiments with respect to the % helical content measured from CD spectra of peptides in PBS:TFE (1:1) solution. P6-cyc is highlighted in red.

To further explore the hypothesis that a stabilization of the α-helical structure could enhance peptide-RBD binding, we designed a cyclized p6 variant, denoted p6-cyc. Starting from p6 sequence, we replaced A36 for Fmoc-(R)-2-(7-octenyl)Ala-OH, and S43 for Fmoc-(S)-2-(4-pentenyl)Ala-OH, and subsequently performed Ring Closing Metathesis (RCM) as described in supplementary information. The positions selected for RCM implied substitution of aliphatic amino acids that were not important for the interaction with the RBD with an adequate distance in the peptide sequence to reach a (i, i+7) type of stapling. The cyclized peptide effectively not only showed a higher helical content in TFE, but also a 25% helical content in buffer, and in accordance with our hypothesis this translated with a slightly better affinity towards the spike protein RBD, as demonstrated by MST experiments. This indicates that peptide cyclization is a promising strategy to incorporate into antiviral peptide design efforts, as also described in other approaches for other systems (Curreli et al. 2020; Maas et al. 2021; Morgan et al. 2021).

Systems preparation for MD simulations

Six ACE2-derived peptides (p6 and five variants named as p6A-E) were simulated and analyzed in complex with wild-type SARS-CoV-2 RBD (residues 319 to 541). The initial configuration of all systems were taken from 7DF4 PDB [7DF4.pfb]. p6 was constructed by truncating ACE2 from residues 22 to 44 and 351 to 357. A glycine was added between the two non-consecutive parts of p6. p6A is a p6-mutated peptide as well as p6B, but the latter also has two more residues from the α-helix, 20 and 21. p6C is an extended version compared to p6 and it was designed based on p6 structure with the addition of ACE2 residues 79 to 83 at the N-terminal. In order to fill the gap between residues 83 and 22 (Figure 4B), residues 20 and 21 were artificially moved between residues 83 and 22, allowing the correct joint of these other non-consecutive parts. Therefore, p6 and p6C peptides present no mutations from the wild-type structure, but spatial rearrangements and the addition of the above-mentioned Gly. As p6D and p6E are also extended peptides, they were prepared with the same protocol as p6C. The RBD residues taken from the PDB structure were 319–541 from the monomer in the up-conformation (chain B).

The protonation state of ionizable residues were fixed to be those most probable at neutral pH, after estimation of side chain pKa values with propka 3.0, an empirical pKa predictor for buried and surface protein residues (Olsson et al. 2011). All of the His residues were determined to be protonated solely at the Nε (HIE).

All systems were modeled with the Amber ff14SB force field (Maier et al. 2015). All systems were solvated with the TIP3P water model (Jorgensen et al. 1983) using the tLeAP module of the Amber simulation suite, in octahedral boxes using a distance of 10 Å from the protein surface, resulting in the addition of approximately 25,000 TIP3P residues for the systems, respectively. A physiological salt concentration of 0.15 M was achieved through the addition of approximately 70 Na⁺ and Cl⁻ ions, respectively, after neutralizing the system charge in each case (p6: –3e; p6A: –4e; p6B: –2e; p6C: –3e; p6D: –3e; p6E: –3e; RBD has +7e additional charge for the complex). Hydrogens and other atoms not resolved in the cryo EM structures, as well as the disulfide bonds (UniProt Consortium 2021), were added using the tLeAP module of the Amber simulation suite (Case et al. 2021).

MD simulations production and analysis protocols

Each system had an initial energy minimization, heating from 0 to 200 K over 50 ps with a 1 fs time step with strong restraints (100 kcal/mol/Ų) on the protein, and then heating from 200 to 300 K under the same conditions with weaker restraints on the protein (50 kcal/mol/Ų). The system was then equilibrated for an additional 50 ps at constant temperature and volume, without restraints, and finally at constant temperature and pressure without restraints for 50 ps. For constant pressure simulations, isotropic position scaling was performed with the Berendsen barostat and a pressure relaxation time of 2.0 ps. The Langevin thermostat was used with a collision frequency of 1.0 ps⁻¹ (Sindhikara et al. 2009; Uberuaga et al. 2004). The SHAKE algorithm was used to constrain bonds to nonpolar hydrogens (Ryckaert et al. 1977). An 8.0-Å cutoff was used for nonbonded interactions.

For each complex, three replicas of 250 ns, under the same NVT conditions described above, with a 2.0 fs timestep, using the PMEMD cuda module of the Amber simulation suite (Case et al. 2021).

Trajectory analysis was performed using the CPPTRAJ module of the Ambertools 19 package of programs. The binding energy was computed from the first 200 ns of the molecular dynamics trajectories extracting snapshops with a 0.2 ns timestep, using the MM-GBSA (molecular mechanics, general-Born surface area) approach implemented in AMBER. The binding energy, ΔEbind,solv, was calculated as the sum of the terms:

where ΔEbind,vacuum corresponds to the substrate binding energy in vacuum and ΔEsolv is the corresponding energy for the solvation process of each species. We performed an energy decomposition for the protein residues involved in interactions with the peptide. For the Solvent-Accessible Surface Area (SASA), we first made a detailed analysis of the interactions observed between the ACE2 and the RBD (Figure 1, SI Table 2). We compiled data already reported in the literature and calculated the SASA, using CPPTRAJ with a solvent-probe radius of 1.4 Å (Weiser et al. 1999), for the key RBD residues involved: 417, 446, 449, 453, 455, 456, 473–476, 474, 475, 476, 484, 486–490, 493–496, 498, 500–502, 505 (Lan et al. 2020; Wang et al. 2020a, 2020b; Xu et al. 2021; Yan et al. 2020). We estimated the per-residue average of antiparallel β-strand structure with the CPPTRAJ module of the AmberTools package.

Peptide synthesis, purification and characterization

All the peptides were synthesized on Rink Amide using automated SPPS standard protocols. For the cyclized version of p6, stapling was performed on a resin following known methodologies for ring closing metathesis: the resin was incubated with a 10 mM solution of bis(tricyclohexylphosphine)-benzylidene ruthenium (IV) (1st generation Grubb’s catalyst) in 1,2-dichloroethane for 1 h twice. In all the cases, the peptide cleavage from the resin was performed using a 95:2.5:2.5 mixture of TFA:water:TIPS. Subsequent purification was performed on prepHPLC systems using C18 chromatographic columns, with different mixtures of water:acetonitrile, and after lyophilization peptides were obtained as pure white powders, their purity being determined by 15 min chromatography runs on a C18 UPLC-MS instrument. Characterization was mainly done by determining helical content by circular dichroism spectra, and high resolution mass spectrometry of the pure samples (SI Figures 6–23).

Peptide binding affinity measurements

For affinity measurements His-tagged RBD of the wild-type SARS-CoV-2 (2019-nCoV) S-protein from recombinant expression in mammalian cells was kindly provided by Dr. Coskun (Technical University Dresden, Germany). S-protein RBD variants alpha (B.1.1.7) and delta (B.1.617.2) were provided from SinoBiological.

RBD labeling and purification was performed using a 2nd generation NHS-Red labeling kit (NanoTemper). RBD solution was diluted with water to 10 µM, rebuffered in sodium carbonate aqueous solution at pH 8.0, and incubated in the dark at 300 µM of the red dye. Size exclusion chromatography allowed to obtain the labeled RBD free of unreacted dye. The labeling efficiency was analyzed from spectroscopy measurements to be approximately 1:1 (protein:dye). Affinity measurements were conducted in DPBS (without Ca²⁺ and Mg²⁺) supplemented with 0.05% (v/v) Tween20, and premium capillaries (NanoTemper).

Microscale thermophoresis (MST) was then measured by making 16 sequential 1:1 dilutions of each peptide, using PBS + 0.05% Tween 20 as diluent, each with a final volume of 10 µL. Each diluted sample was then incubated with 10 µL of the 10 nM labeled RBD stock solution, therefore always keeping a 5 nM concentration of the labeled target protein in every sample.

MST signal was then acquired in a NanoTemper Technologies Monolith NT.115 Pico instrument, at an excitation power of 20% and a MST power of 40%, the signal was analyzed 1.5 s after the start of the IR-laser, and the obtained data fitted as shown previously (Bhatia et al. 2017). These conditions were kept the same for all the measured samples. All measurements were done at least by triplicate (SI Figures 3–5).