Adaptive Evolution of Peptide Inhibitors for Mutating SARS-CoV-2.

The SARS-CoV-2 virus is currently causing a worldwide pandemic with dramatic societal consequences for the humankind. In the last decades, disease outbreaks due to such zoonotic pathogens have appeared with an accelerated rate, which calls for an urgent development of<br>adaptive (smart) therapeutics. Here, we develop a computational strategy to adaptively evolve peptides that could selectively inhibit mutating S protein receptor binding domains (RBDs) of different SARS-CoV-2 viral strains from binding to their human host receptor, angiotensin-converting enzyme 2 (ACE2). Starting from suitable peptide templates, based on selected ACE2 segments (natural RBD binder), we gradually modify the templates by random mutations, while retaining those mutations that maximize their RBD-binding free energies. In this adaptive evolution, atomistic molecular dynamics simulations of the template-RBD complexes are iteratively perturbed by the peptide mutations, which are retained under favorable Monte Carlo decisions. The computational search will provide libraries<br>of optimized therapeutics capable of reducing the SARS-CoV-2 infection on a global scale. <br>.

The very fast spreading of a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in the human population has already led to hundreds of thousands of fatalities and dire socio-economic effects worldwide. Therefore, many concepts and strategies were initiated to identify drug targets [1], develop effective therapeutics and vaccines against SARS-CoV-2 [2, 3], including the optional use of therapeutics developed for other purposes [4-7]. Some of the strategies have focused on direct targeting of the spike (S) protein [8, 9], located on the outer surface of the SARS-CoV-2 virion, which initiates the cell entry process. This process starts by binding of the receptor binding domain (RBD) on the S protein to a human host receptor angiotensin converting enzyme 2 (ACE2), followed by a proteolytic cleavage and release of its S1 subunit [10]. Then, the remaining S2 subunit undergoes a conformational change, which eventually leads to virus-host membrane fusion [11]. The S protein RBD of SARSCoV-2 binds more strongly to ACE2 than SARS-CoV [12]. This fact, together with other innovative characteristics of the S protein, including a furin preactivation of the cleavage site, and a hidden RBD, likely contribute to the large infectivity of SARS-CoV-2 [12].

SARS-CoV-2 is mutating and further adjusting to the human environment, like other novel viral pathogens. Many strains of SARS-CoV-2 have already been detected [13, 14]. Some mutations, such as D614G on the S protein, lead to rapid and enhanced viral transmission [15], causing this strain to locally dominate. The mutating SARS-CoV-2 coronavirus could also adapt to new hosts, such as domestic animals [16]. Therefore, to mitigate the large spreading and effects of SARS-CoV-2, it is important to identify classes of therapeutics that could be rapidly developed to act against multiple coronavirus strains.

Here, we address this problem by introducing libraries of adaptive peptide therapeutics that could block different S protein strains from interacting with ACE2. With the use of advanced computational methods, we show that adaptive evolution of suitable peptide templates can provide multiple inhibitors for competitive binding to different S protein variants. Using libraries of peptides adapted to different S protein variants could prevent their mutational escape, in analogy to using cocktails of S protein antibodies [9]. Our algorithm can evolve peptide inhibitors that competitively bind to (block) a set of desired targets, such as different S protein variants, S proteins with glycan shielded sites [17-19], and other related systems [20].

Results and Discussion

S proteins variants.

Over the time and geographical regions, SARS-CoV-2 virus has so far evolved into more than 104 mutated strains, shown in the mutation tree of publicly available unique genome sequences (through June 2020) in Fig. 1a [21]. Out of these mutations, so far (June 2020) 25 have been recognized in the S protein RBD, as summarized in Fig. 1b. Five of these mutations include amino acids that form a part of the ACE2-binding surface (A475V, G476S, S477I, V483A, and V503F), as highlighted in Fig. 1c.

FIG. 1:

Mutations of SARS-CoV-2. a) Time-dependent mutation tree of SARS-CoV-2 colored according to geographical regions of origin (through June 2020) [21]. b) 25 single mutations identified in RBD of the S protein. c) 5 amino acid mutations on RBD in contact with the ACE2 receptor. d) Binding free energies are evaluated with the MMGB-SA method for the ACE2-RBD complexes, including the originally reported RBD (wild type, labeled as WT) and the five mutant RBDs listed in panel c.

In Fig. 1d, we present the RBD-ACE2 host-receptor binding free energies, ΔG, obtained in the presence of these 5 mutations. The 5 RBD:ACE2 complexes are simulated for 30 ns, and their average binding energies are obtained in the last 15 ns. The originally reported RBD and the S477I RBD have the strongest binding to the human ACE2, ΔG ≈ −50 kcal/mol, while the other systems bind with ΔG ≈ − (40 – 50) kcal/mol. In order to block all these RBD variants, the peptide inhibitors should have comparable or lower ΔG values.

Adaptation of peptides by single mutations.

Two ACE2-based peptide structures, shown in Fig. 2a, are selected as templates for the first generation peptide inhibitors of the S protein [8]. The smaller template-1 includes single truncated α1 helix of ACE2 (amino acids 21 – 43), and the larger template-2 includes two α1α2 helices of ACE2 (amino acids 19 – 83). In the adaptive evolution search for optimized therapeutics, the selected ACE2-extracted peptide templates will be gradually mutated to increase their binding strength to RBD.

FIG. 2:

Modeling of peptide-RBD complexes. a) Complexes of S protein RBD (blue) and two peptide templates (red). Locations of five S protein mutations examined in the present work are marked by blue spheres. Amino acid residues changed in singly mutated peptides are marked by yellow spheres. b) Free energies of binding, ΔG, between the originally reported S protein RBD and the wild type or singly mutated ACE2-based peptides. Locations of mutated peptide amino acids are highlighted in panel a. c) Snapshots of initial and optimized template-1 peptides binding to the original RBD. The sequence of the optimized peptide was obtained after 100 mutation attempts, with 10 ns long MD simulation after each mutation. The final peptide with the optimized sequence was further relaxed in a 175 ns MD simulation. The initial peptide is shown as a red helix, with amino acids that are subsequently mutated shown in thin faded yellow licorice. The optimized peptide is shown as an orange helix, with mutated amino acids shown in thick yellow licorice. d) Adaptive evolution of template-1. The plot shows the binding free energies, ΔG, between the peptide and the original RBD, presented as a function of the performed mutation, where the peptide:RBD complexes are relaxed for 10 ns after each mutation attempt. e) The time evolution of ΔG between the final optimized peptide and the original RBD. The average value, obtained from the last 75 ns of the trajectory (gray), is ΔG = −57 kcal/mol. The faded green line shows the data points calculated every 0.1 ns, and the dark green line shows the running average. f) Initial and optimized sequences of template-1 peptides. The final peptides were optimized for binding to the original and mutant RBDs, with peptide-RBD complexes relaxed in 10 ns MD simulations after each attempted mutation.

In recent mutagenesis experiments, the whole ACE2 with single mutations in regions directly contacting RBD were examined for their binding to the original S protein [22]. To perform preliminary testing of our adaptive evolution search of peptide therapeutics, we first selected the most fit mutants from these experiments, and implemented their mutations in our templates-1,2. We simulated 22 peptides, i.e., the original templates and their 10 single mutants, complexed with the original S protein RBD. Their free energies of binding, ΔG, were evaluated in 100 ns simulations and presented in Fig. 2b. Template-1 binds to RBD with ΔG ≈ −19 kcal/mol, while its mutants have higher affinities giving ΔG ≈ −(24 – 35) kcal/mol. all cases, template-1 significantly In changes its conformation in the bound configuration, as the helix loses the curvature observed when within ACE2, and the hydrogen bonding between Glu35 (template-1) and Gln493 (RBD), enabled by the helix curvature, is broken.

In contrast, template-2 has more direct contacts with RBD than the shorter template-1 variants, so it binds to it more strongly, ΔG ≈ −36 kcal/mol. However, only two template-2 mutants (H34A and K31W) have higher affinities to RBD than the original template-2, having ΔG ≈ −45 kcal/mol. These simulations also revealed that peptides with a stronger binding covered larger RBD sections (Fig. S1), and reduced the RBD exposure to other potential binding partners. These results show that the experimental results obtained for mutated ACE2 [22] (Fig. S2) can provide a good guidance in the mutation of template-1, but the same mutations are less effective in the larger template-2.

Adaptive evolution of peptide inhibitors.

The above results have clearly demonstrated that suitable peptide templates with appropriate mutations can acquire a strong binding to specific targets. To optimize such peptides against specific viral strains, we have developed combined mutation/selection (evolution) computational algorithms which can guide a multi-step adaptive evolution of the peptides: 1) Random mutations are introduced into random positions of the peptide templates. 2) Short MD simulation trajectory of the mutated-peptide:RBD complex are run and followed by a selection or rejection of the mutation via Monte Carlo (MC) sampling using a Metropolis criterion applied to the change of the free energy of peptide-RBD binding, ΔG (Methods). 3) The mutation/selection process is iteratively repeated until the binding affinity of peptides to the target S protein RBD is satisfactory (Methods). 4) Additional evolution of the molecules might be considered after the MC decisions to allow for a better internal relaxation of the molecules. Due to a partly stochastic nature of MD simulations, the randomness in mutations, and the MC selection, different peptides can be obtained in separate trajectories that correspond to separate local minima of the free energy surface. These peptides form a pool (ensemble) of potential therapeutics evolved for a selected viral strain, which can be further enriched by considering multiple viral strains.

In Figs. 2c–e, S3, and Table S1, coupling of template-1 to RBD was optimized in the adaptive evolution, where 10 ns MD simulation trajectories were generated after each trial mutation of template-1. Of the 100 mutations attempted, 13 mutations were accepted, and 11 amino acids were changed (individual residues can be mutated more than once). Figure 2d reveals the progression of binding free energies with the mutations of template-1, starting from ΔG = −19 kcal/mol. As detailed in Fig. 2c, during the adaptive evolution, the mutating helical peptide lost its bending (this change is independent of mutations) and multiple initial contacts with RBD. At the same time, it shifted with respect to its initial position and formed many new contacts. Peptide residue E37 formed a salt bridge with the original RBD, while residues Q24, Y26 (mutated), Q30 (mutated), S41 (mutated) and R42 (mutated) formed hydrogen bonds of varying stability with the original RBD. The resulting peptide bound to RBD with ΔG = −70 kcal/mol at the end of thirteen 10 ns-long trajectories (associated with individual accepted mutations).

As shown in Fig. 2e, additional 175 ns relaxation of this peptide resulted in a slightly less favorable ΔG ≈ −57 kcal/mol. Therefore, adaptive evolution requires sufficiently long relaxation times for a good stabilization of the whole system. Short relaxation times may result in incomplete peptide adjustments and free energies that can be misleadingly favorable. Moreover, a faster MC convergence could be achieved by considering the whole free energy changes rather than the peptide-RBD binding free energies. However, internal reorganizations of molecules inevitably take part in long trajectories, so the difference in binding energies alone might be sufficient for the MC decision, as long as additional relaxation is allowed between individual MC steps (point 4 in the method).

Next, we adaptively evolved template-1 coupled to 3 separate singly-mutated RBDs, chosen from Fig. 1d. For simplicity, 100 mutations were attempted, followed by 10 ns simulations after each attempt. The adaptive evolution gave peptides with ΔG ≈ −(45 – 70) kcal/mol, as summarized in Figs. 3a–c and Fig. 2f. Peptides targeting A475V and G476S RBDs each had 5 accepted mutations, respectively, while peptides targeting S477I RBD had 19 accepted mutations. In the A475V RBD case, one of the early accepted mutations lead to breaking of the helix secondary structure, and thus to a different peptide-RBD binding mode. This shows that individual alpha helices without additional stabilization, such as by side branching [23], might be too simplistic therapeutics.

FIG. 3:

(a-c) Adaptive evolution of template-1 coupled with singly-mutated RBDs. (d) Adaptive evolution of template-2 coupled with the original RBD. After attempted mutations, peptide:RBD complexes were relaxed in 10 ns simulation steps.

In Fig. 3d, the adaptive evolution was performed with a more stable template-2 (α1α2), but random mutations were only introduced into the α1 helix, which was in direct contact with the original RBD. After 12 accepted mutations and 10 changed amino acids (listed in Table S2), the binding strength increased from ΔG = −36 kcal/mol to −60 kcal/mol. Therefore, the adaptively evolved template-2 can compete with the whole ACE2, having ΔG ≈ −50 kcal/mol (Fig. 1c). Both the initial and optimized template-2 preserve the curvature of the α1 helix, despite the mutation of E35 (to Y35), which is observed to interact with Q493 (RBD). This feature preserves the binding pattern observed in ACE2:S protein RBD complex [1]. However, the salt bridge between D30 (initial peptide) and K417 (RBD) (shown in Fig. S4 for template-2 optimized with 1 ns simulations after mutations) is lost in the peptide optimized with 10 ns simulations after mutations.

In summary, we have demonstrated that ACE2-based peptide templates can be adaptively evolved by computation using mutation/selection processes to form optimized inhibitors for a strong and competitive S protein RBD binding. The developed approach can be used to design peptide inhibitors based on templates extracted from different ACE2 polymorphs, including those from other species [24], and other proteins binding to viral pathogens. The optimized inhibitors obtained in different evolution runs can be collected to form libraries of suitable therapeutics for different RBD variants. Cocktails (ensembles) of peptide therapeutics could be delivered by different means to provide a broad-spectrum protection against different SARS-CoV-2 strains.

Methods

MD simulations.

The simulated peptide template-1 (amino acids 21 – 43 of ACE2) and template-2 (amino acids 19 – 83 of ACE2) were separately bound to the S protein RBD. All structures were directly based on the crystal structure of the human ACE2 protein bound to the SARS-CoV-2 Spike protein RBD (pdbID: 6LZG) [10]. The mutations in peptides and RBD were introduced using the psfgen plugin in VMD [25].

The systems were simulated using NAMD2.13 [26] and the CHARMM36 protein force field [27]. The simulations were conducted with the Langevin dynamics (γ = 1 ps−1) in the NpT ensemble, at temperature of T = 310 K and pressure of p = 1 bar. The particle-mesh Ewald (PME) method was used to evaluate Coulomb coupling, with periodic boundary conditions applied [28]. The time step was set to 2 fs. The long range van der Waals and Coulombic coupling were evaluated every 1 and 2 time steps, respectively. After 2, 000 steps of minimization, the solvent molecules were equilibrated for 1 ns, while the complexes were restrained using harmonic forces with a spring constant of 1 kcal/(molÅ). Next, the systems were equilibrated in 100 ns production MD runs with no restraints.

MMGB-SA calculations.

We used the Molecular Mechanics Generalized Born - Surface Area (MMGB-SA) method [29, 30] to estimate the relative binding free energies between peptides (or ACE2) and RBDs. The free energies were estimated from separate MMGB-SA calculations for three systems (peptide/ACE2, RBD, and the whole complex) in configurations extracted from the MD trajectories of the whole complex in the explicit solvent. The MMGB-SA free energies of the extracted configurations of the three systems were calculated as where E, G, G, and ΔS are the sum of bonded and Lennard-Jones energy terms, the polar contribution to the solvation energy, the nonpolar contribution, and the conformational entropy, respectively. The E, G and G terms were calculated using the NAMD 2.13 package [26] generalized Born implicit solvent model [31], with a dielectric constant of the solvent of 78.5. The G term for each system configuration was calculated in NAMD as a linear function of the solvent-accessible surface area (SASA), determined using a probe radius of 1.4 Å, as G = γSASA, where γ = 0.00542 kcal/molÅ2 is the surface tension. The ΔS term was neglected, since the entropic contribution differences nearly cancel when considering protein-protein binding of single residue mutants [32, 33]. Moreover, the entropy term, which is often calculated with a large computational cost and low prediction accuracy, is likely to be similar for the studied systems, which differ in single or several mutations. Since the G values are obtained for peptide configurations extracted from the trajectories of complexes, G doesn’t include the free energies of peptide reorganization; the correct free energies of binding should consider configurations of separately relaxed systems. The approximate binding free energies of the studied complexes were calculated as 〈ΔG〉 = 〈G(P/ACE2–RBD)–G(P/ACE2) – G(RBD)〉, where P/ACE2 represents peptides or complete ACE2, and the 〈 averaging 〉 is performed over configurations within the second half of the calculated trajectories.

Adaptive evolution algorithm.

A mutation/selection algorithm is developed and used to iteratively increase the affinity of binding between peptide templates and the S protein RBD. The algorithm involves sequences of steps combining molecular dynamics simulations and MC decisions using the Metropolis criterion [34, 35], which results in MC sampling of the peptide sequence space. Initially, the selected peptide template (directly extracted from ACE2) complexed with RBD is equilibrated for 100 ns. Then, the free energy of binding of the peptide:RBD complex, , is evaluated. Next, a random mutation is introduced at a random position in the peptide, followed by a short 1 – 10 ns equilibration in MD simulations of the complex, and evaluation of the free energy of binding of the complex. Finally, the mutation is accepted if or if the Metropolis criterion is satisfied, where r is a random number in the (0, 1) interval. If the mutation is accepted, then the new peptide becomes the new reference peptide and its ΔG becomes the reference value for the next attempted mutation. In each run of the algorithm, 100 – 200 mutations were attempted on peptide templates, as stated in the results.