Next Article in Journal
Concentration Sensitivity of Nucleic Acid and Protein Molecule Detection Using Nanowire Biosensors
Previous Article in Journal
Plasma-Activated Water Promotes Wound Healing by Regulating Inflammatory Responses
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

In Silico Design of Peptide-Based SARS-CoV-2 Fusion Inhibitors That Target WT and Mutant Versions of SARS-CoV-2 HR1 Domains

by Shana V. Stoddard 1,*, Felissa E. Wallace 1,2, Serena D. Stoddard 1,3, Qianyi Cheng 4, Daniel Acosta 1, Shaliz Barzani 1, Marissa Bobay 1, Jared Briant 1, Christian Cisneros 1, Samantha Feinstein 1, Kelsey Glasper 1, Munazza Hussain 1, Abigail Lidoski 1, Pranay Lingareddy 1, Grace Lovett 1, Leslie Matherne 1, Jackson McIntosh 1, Nikita Moosani 1, Lia Nagge 1, Kudzai Nyamkondiwa 1, Isaiah Pratt 1, Emma Root 1, Mary Rose Rutledge 1, Mackenzie Sawyer 1, Yash Singh 1, Kristiana Smith 1, Ubaid Tanveer 1 and Sona Vaghela 1add Show full author list remove Hide full author list
1
Department of Chemistry, Rhodes College, 2000 North Parkway, Memphis, TN 38671, USA
2
Walnut Hills High School, 3250 Victory Pkwy, Cincinnati, OH 45207, USA
3
College of Veterinary Medicine, Tuskegee University, 201 Frederick D Patterson Dr, Tuskegee, AL 36088, USA
4
Department of Chemistry, University of Memphis, 3744 Walker Ave, Memphis, TN 38152, USA
*
Author to whom correspondence should be addressed.
Submission received: 14 April 2021 / Revised: 1 July 2021 / Accepted: 2 July 2021 / Published: 8 July 2021

Abstract

:
In 2019, novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) began infecting humans, resulting in the COVID-19 pandemic. While the push for development of vaccines has yielded some positive results, the emergence of additional variants has led to concerns surrounding sustained vaccine effectiveness as the variants become the dominant strains. This work was undertaken to develop peptide-based antivirals capable of targeting both the wildtype (WT) heptad repeat 1 (HR1) domain of SARS-CoV-2 and the new HR1 variants which have developed. In silico protein mutagenesis, structural characterization, and protein–protein molecular docking were utilized to determine molecular interactions which facilitated binding of peptide-based antivirals targeting the HR1 domains. Molecular dynamics simulations were utilized to predict the final binding affinities of the top five peptide inhibitors designed. This work demonstrated the importance of hydrophobic interactions in the hydrophobic gorge and in the rim of the HR1 domain. Additionally, the placement of charged residues was shown to be essential in maximizing electrostatic interactions. The top five designed peptide inhibitors were all demonstrated to maintain good binding affinity to the WT and the variant HR1 SARS-CoV-2 domains. Therefore, the peptide inhibitors designed in this work could serve as potent antivirals which are effective in targeting both the original SARS-CoV-2 and the HR1 variants that have developed.

Graphical Abstract

1. Introduction

The 2019 coronavirus disease (COVID-19) has created an unprecedented pandemic that has taken an alarming toll on the world. The COVID-19 pandemic resulted from a novel strain of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1,2]. SARS-CoV-2 is now one of three coronaviruses, which causes severe illness in humans. The first two strains to cause severe respiratory illness in humans were severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle Eastern respiratory syndrome coronavirus (MERS-CoV) [3,4,5,6]. The current pandemic has infected 126,596,165 individuals, resulting in the loss of 2,775,970 lives in 192 countries and regions as of March 2021 [7]. Currently, 13 vaccines have been approved in a number of countries for use to decrease COVID-19 infections [8]. The current vaccines were designed for the original emergent strain of SARS-CoV-2. Unfortunately, many new strains have now emerged, with five being classified as variants of concern, namely B.1.1.7 [9,10], P.1 [11], B.1.429 [12], B.1.427 [12], and B.1.351 [13]. The emergence of variants may limit the effectiveness of these vaccines [14,15,16,17]. In light of this, antivirals must also be pursued as treatment options against COVID-19.
SARS-CoV-2 is the seventh human coronavirus and is classified as a betacoronavirus [18]. The mechanism of invasion used by SARS-CoV-2 involves a two-step process of attachment of the receptor-binding domain (RBD) to the angiotensin-converting enzyme-2 (ACE2) receptor followed by fusion to the cell membrane [19,20]. The fusion process is controlled by the binding of heptad repeat 1 (HR1) and heptad repeat 2 (HR2) on the spike (S) glycoprotein. The HR1 and HR2 domains form a fusogenic hairpin, bringing the viral capsid to the surface of the cell. Prevention of hairpin formation through the binding of a peptide-based inhibitor would ultimately prevent the virus from invading the cell. Coronaviruses have multiple copies of the S glycoprotein on their surface, making the S glycoprotein a widely available target [21]. The development of peptide-based inhibitors that target the HR1 domain and prevent the fusion process are currently being pursued by many groups [22,23]. Many of the developed peptides do not have structural information available yet from X-ray crystallography, NMR, or electron microscopy studies. Information about the key binding interactions of current peptide inhibitors to the SARS-CoV-2 HR1 domain would be essential to help further optimization efforts of developing tighter binding peptide inhibitors for SARS-CoV-2.
The S protein of SARS-CoV-2 has several known mutations in the RBD and the HR1 domains [24,25,26]. The two mutations A930V and D936Y are present in the HR1 domain. These two mutations have been shown to be present in the MT050493 and MT259281 strains which emerged in India and the USA, respectively [17]. The A930V HR1 mutation has been hypothesized to increase chances of peptide entry inhibitor resistance. While the current computational work has been used to identify potential antivirals for COVID-19 [27,28,29,30,31,32,33,34], to date studies have not addressed the impact of the A930V and/or the D936Y mutations in the design process. A peptide inhibitor that could target the wildtype (WT) and the variants would be ideal to pursue as a therapeutic intervention for SARS-CoV-2.
In this work, peptide-based inhibitors that target the WT SARS-CoV-2 HR1 and the SARS-CoV-2 HR1 variants A930V and D936Y were designed. Additionally, molecular docking was used to predict the binding of the EK1 peptide to the WT SARS-CoV-2 HR1 and characterize the binding interactions of the EK1 peptide inhibitor to the WT SARS-CoV-2 HR1 domain; a complex that has not yet been crystallized. Molecular docking which can be used to identify how two molecules interact with each other greatly assists in optimization work [27,35,36,37,38,39]. Molecular docking was used in this work to optimize the binding interactions between the WT SARS-CoV-2 HR1 domain and peptide inhibitors, followed by molecular dynamics simulations to predict binding affinities of our top antiviral designs. This study could provide guidance to further optimization efforts to develop peptide-based inhibitors that can target a broad range of SARS-CoV-2 variants.

2. Materials and Methods

2.1. Modeling of A930V and D936Y SARS-CoV-2 HR1 Domains

Homology models of the two SARS-CoV-2 HR1 variants, A930V and D936Y, were created in UCSF Chimera using the rotamer function. Additionally, a double mutant, with both A930V and D936Y, was also developed using the same process. Minimization of each model was performed to ensure that the residue mutations did not cause dramatic changes in the overall structure of the HR1 domain. Minimization was performed in UCSF Chimera. Then, 1000 steps of steepest decent followed by 10 steps of conjugate gradients were used to complete minimization of HR1 domain models.

2.2. Structural Analysis of SARS-CoV-2 HR1

Characteristics of the HR1 binding sites were evaluated to determine structural properties that would enhance binding interactions. Both electrostatic potential maps and hydrophobicity renderings were produced for the HR1 domains. All characterization was performed in UCSF Chimera [40]. The protein structures for the WT SARS-CoV-2 HR1 domain were imported into UCSF Chimera from the Protein Data Bank (PDB). PDB code for the SARS-CoV-2 used was 6LXT [22]. The homology models of the SARS-CoV-2-A930V and SARS-CoV-2-D936Y variants and the double mutant SARS-CoV-2-A930V-D936Y were used for analysis of these strains. The EK1 pan-CoV inhibitor was imported from PDB 5ZVM [41] and used as a template for understanding binding interactions between peptide inhibitors and the HR1 sites.

2.3. Conservation of HR1 and HR2 Binding Site Analysis

The HR1 and HR2 domains of SARS-CoV-2 were analyzed to determine the degree of conservation of residues using the program ConSurf [42,43]. The 6LXT file contained both the HR1 and HR2 domains of SARS-CoV-2. Therefore, the two domains were separated in USCF Chimera to produce individual files. One file was produced containing the HR1 residues and a second containing the HR2 domain residues of SARS-CoV-2 from the 6LXT pdb file. The new pdb files were then submitted to the ConSurf server for SARS-CoV-2 and default parameters were selected for the analysis of conservation. Multiple sequence alignment for conservation analysis with 150 sequences used for conservation analysis is provided in Supplementary Table S1. UCSF Chimera was used to recolor each file to visualize the degree of conservation using the ConSurf coloring scheme.

2.4. Design of Peptide Inhibitors

The EK1 pan-CoV peptide inhibitor (PDB: 5ZVM) [41] was used as a structural template for design of novel peptide-based inhibitors, which can target the new variants that are developing. After docking of the EK1 peptide to the HR1 domain, a rational design strategy was employed to determine what in silico mutations could be introduced to promote additional molecular interactions with the HR1 domain. Existing molecular interactions which contributed to steric clashes, electronic clashes, or mismatches in hydrophobic and hydrophilic residue pairings were identified. For example, if the docked position of the peptide inhibitor resulted in the placement of a negatively charged residue to be positioned near a negatively charged residue on the HR1 domain of SARS-CoV-2, this was considered an electronic clash. If the docked position of the peptide inhibitor designed was demonstrated to position a hydrophobic residue in a hydrophilic region on the HR1 domain, this was considered to be a mismatch, resulting in an unoptimized interaction region. Mutations on the template inhibitor peptide were then created by introducing new amino acids in the positions which were hypothesized to remove the clash. Additionally, the docked positions were evaluated to determine if there were locations where residues could be introduced that would add a new interaction resulting in increased binding affinity. For example, if a polar uncharged residue on the template was near a positively charged residue on the HR1 domain, then it was hypothesized that placement of a negatively charged residue might enhance binding to the positively charged residue while maintaining any hydrogen bonding, if present. After identification of new molecular interactions to be introduced, peptide sequences were mutated in silico using the rotamer function in UCSF Chimera. In some cases (the 100 series of peptides, i.e., 101, 102, 103…), random mutations were selected by using the sequence only without inspection to produce variations in binding orientations for analysis. The most probable side chain orientation which did not produce any significant clashes upon visual inspection was selected as the rotamer orientation. Designed peptides were then evaluated using molecular docking to determine the impact of the mutation on binding. New rationally designed peptides were subsequently used as new templates for further optimization and design. Subsequent in silico mutations were introduced into these new designed peptides to produce final peptide candidates with optimized molecular interactions. Of the randomly designed peptides using the sequence, only 2 of the peptides which were predicted to improve binding were shown to do so, producing a true positive rate of 12.5%, while a total of 14 which were expected to improve binding were not shown to do so, producing a false positive rate of 87.5%. Of the rationally designed peptides, 54 of the peptides which were predicted to improve binding were shown to do so, producing a true positive rate of 78.2%, while a total of 15 which were expected to improve binding were not shown to do so, producing a false positive rate of 21.7%.

2.5. Protein–Protein Docking

The HR2 domains, the EK1 peptide template, and the designed peptide inhibitors were all docked to each of the WT SARS-CoV-2, SARS-CoV-2-A930V, SARS-CoV-2-D936Y, and SARS-CoV-2-A930V-D936Y HR1 protein structures. Protein–protein docking was performed to determine the strength of binding between the HR2 and HR1 binding site and between the peptide inhibitors and the HR1 binding sites using the program HDOCK [44,45,46,47,48]. Additionally, an RMSD calculation was performed to determine the difference in position of the HR2 docked peptide to the crystal structure in UCSF Chimera. Default parameters were used for all docking runs. HDOCK server uses an iterative knowledge-based scoring function referred to as an ITScore-PP [48] to rank the top 10 poses. Larger negative numbers indicate stronger binding interactions between the two macromolecules evaluated. This energy score has been demonstrated to correlate well to experimental binding affinities with a correlation coefficient of 0.71. The resulting HDOCK docking poses were further analyzed for patterns to understand structural implications contributing to enhanced binding scores.

2.6. Molecular Dynamics Simulations

From the binding poses of the top 5, middle 4, and bottom 3 candidates (ranked by docking score) of designed peptides, EK1, HR2, and HR1 structures were retrieved from the docking results. The Zn2+ ions were included as part of the HR1 domain in the calculation. For each of the complexes, the combined structures of HR1 and the candidate peptides/EK1/HR2 were used as the starting structures in MD simulations. The proteins, candidates, and complexes were modeled using the AMBER ff14SB force field [49], and each of the complex structures was solvated in a rectangular box with TIP3P water molecules [50]. The 12-6-4 LJ-type nonbonded model was used for Zn2+ ions [51], and counter ions (Na+/Cl) were added to neutralize the complex system when necessary. Then, the complex system was equilibrated by (1) short 500 steps of steepest descent followed by 500 steps of conjugate gradient energy minimization, (2) 50 ps heating from 0 K to 300 K in constant volume periodic boundaries (NVT), (3) 50 ps density equilibration in constant pressure periodic boundaries (NPT), (4) 500 ps NPT equilibration, and (5) 6 ns production at 300 K. In the minimization, heating, and density equilibration stage of the simulation, weak restraints of 2.0 kcal/mol Å [50] were applied on the complex, and no restraints were applied in the 500 ps equilibration and production. With hydrogen atoms constrained with SHAKE [52], simulations were run using a Langevin thermostat with a time step of 2 fs. The MD simulations were carried out with AMBER 18 [53]. The MM-PBSA method [54] was employed to calculate the binding free energies ( Ɗ G b i n d ) with a single trajectory with the MMPBSA.py program [55] in AMBER18. The entropy contribution was ignored in this study. The internal and external dielectric constant was set to 4 and 80, respectively. Other parameters were kept as the default.

3. Results

3.1. Comparision of the HR1 Domain of the WT SARS-CoV-2 HR1, A930V SARS-CoV-2, D936 SARS-CoV-2, and A930V-D936Y SARS-CoV-2

3.1.1. Evaluation of Conservation in SARS-CoV-2 HR1 Domain

Structurally, the HR1 domain is a trimer of three HR1 regions which twist around each other similarly to the collagen triple helix (Figure 1a). The surface exposed portions of each helix were demonstrated to be the most variable portions of the HR1 domains. For the purposes of this study, the surface exposed residues of the helices are called the “rim” of the HR1 domain. The buried portions of the helices that create the hydrophobic groove binding pocket for the HR2 residues are primarily conserved (Figure 1b). The conservation is, however, not continuous through the inner groove, thus creating two key regions. The first region (region A) located at the N terminal end of the HR1 domain begins with residue N914 and ends at residue S939. This region contains both the A930 residue which is slightly conserved and the D936 residue which is demonstrated to be highly variable. The second region (region B), which is longer than the first region, begins in the middle of the HR1 at the S943 domain and stops at residue E988. Supplementary Table S1 shows the residues that were identified to be variable and the corresponding residues that are found at these positions. The multiple sequence alignment for the 150 sequences used to determine residue variability is shown in Supplementary Table S2. Residues identified to be variable were typically polar uncharged or charged residues. These substitutions would be ideal for further investigations to determine any potential impacts on structural stability of future variants.

3.1.2. Evaluation of Electrostatic Charges and Hydrophobicity in SARS-CoV-2 HR1 Domain

Currently, two mutations have emerged in the HR1 domain of SARS-CoV-2, A930V and D936Y [17,26]. These two mutations are present in the MT050493 and MT259281 strains which emerged in India and the USA, respectively [17]. These substitutions could contribute to changes in both the electrostatic and hydrophobic distribution of residues in the HR1 domain. A comparison of the changes in electrostatic charges and the hydrophobicity of HR1 domains of the WT, A930V, D936Y, and double mutant A930V-D936Y are shown in Figure 2 and Figure 3, respectively.
The electrostatic charges and the hydrophobicity of the WT HR1 domain of SARS-CoV-2, residues 912 to 988 (PDB: 6LXT), were evaluated to determine important characteristics that might facilitate binding capabilities. Additionally, changes in electrostatic and hydrophobicity composition which occurred as a result of the A930V, D936Y, and A930V-D936Y mutants were also evaluated to determine the impact on the binding surface of the HR1 domain. Molecular docking of the HR2 domain to the D936Y mutant produced a lower HDOCK docking score of −693.65, indicating weaker binding, compared to the WT docking score of −756.39. All HDOCK scores values are reported as ITScore-PP [48]. The D936Y mutation which decreases the negative charge of the rim was therefore demonstrated to decrease binding affinity to the HR2 domain. Residue R1185 of the HR2 domain binds electrostatically to D936 of the WT SARS-CoV-2 HR1 domain. In the D936Y and A930V-D936Y mutants, this electrostatic interaction is disrupted; however, a cation–pi interaction can be created between R1185 and the aromatic ring of Y936. Taken together, this demonstrates that the electrostatic interaction is favored over the cation–pi interaction between the HR1 and HR2 domains of SARS-CoV-2. This preference should be considered in the design of peptide-based inhibitors targeting SARS-CoV-2.
The inner grooves of regions A and B are both very hydrophobic, as expected (Figure 3). In region A, the A930V mutation contributes to an increase in hydrophobicity in the groove (Figure 3b,d). Residue L1193 of the HR2 peptide participates in a hydrophobic interaction with the HR1 domain in this pocket near A930. Molecular docking of the HR2 domain to the A930V variant produces a slightly higher docking score of −758.22 compared to the docking score of the HR2 domain to the WT SARS-CoV-2 (−756.39). The resulting increase in the hydrophobicity due to the A930V substitution subsequently produces only a moderate increase in the binding affinity of the HR2 to the HR1 domain. The D936Y mutation also creates a more hydrophobic area as the Tyr sidechain has a large aromatic ring that can also participate in hydrophobic interactions and commonly participates in pi–pi stacking hydrophobic interactions.
When considering the double mutant, the docking score of the HR2 domain to the A930V-D936Y mutant was −693.65. This is a weaker binding affinity, indicated by the smaller negative number, when compared to the WT SARS-CoV-2 HR1 domain interaction with the HR2 peptide. The result, which is similar to that of the D936Y mutation, demonstrates that when considering these mutations, the loss of the electrostatic interaction on the rim of the HR1 domain contributes to more destabilization in binding. Taken together, this demonstrates that while the core of the HR1 domain is very hydrophobic, considerations of electrostatic interactions could be more important in the design of peptide-based inhibitors targeting SARS-CoV-2 and should be incorporated into designs.

3.2. Characterization of Binding Interactions between SARS-CoV-2 HR1 and HR2 Domains and SARS-CoV-2 HR1 Domain and EK1 Peptide

The structural basis for binding of the EK1 peptide to SARS-CoV has been reported via X-ray crystallography [41]. The conservation between the SARS-CoV and SARS-CoV-2 HR1 domains is 89% homologous, suggesting a similar mode of binding. We confirmed this using protein–protein docking in this study and determined the key molecular interactions promoting binding between the WT SARS-CoV-2 HR1 domain and EK1. The EK1 peptide was removed from the 5ZVM crystal structure and docked to the WT SARS-CoV-2 HR1 domain. For comparison, the HR2 domain and the HR1 domains in the crystal structure 6LXT were separated into distinct pdb files and also subjected to protein–protein docking. The best ranked docking pose of 6LXT HR2 to the HR1 domain was identical to the crystal structure file, producing an RMSD of 0.000, computed in UCSF Chimera, demonstrating the accuracy of the docking program (Supplementary Figure S1). The EK1 peptide was determined to have a mode of binding similar to both the binding present between the HR1 and HR2 SARS-CoV-2 domains and the EK1 peptide bound to the HR1 domain of SARS-CoV (Figure 4).
HDOCK docking scores for the HR2 domain and the EK1 peptide were −756.39 and −449.70, respectively. The HR2 domain produces 22 hydrogen bonds to the HR1 domain (Figure 4a). The EK1 peptide was shown to have 15 hydrogen bonds to the WT SARS-CoV-2 HR1 domain (Figure 4b).
The HR2 and EK1 peptides can be structurally described as having three distinct regions: a helix core flanked on the N terminal and the C terminal ends by two less structured peptide tails. Residues participating in hydrophobic interactions which occur between the HR1 and HR2 domains of SARS-CoV-2 are located primarily in the gorge in each of these three regions, as expected (Figure 5). The hydrophobic residue composition of the HR2 domain consists of only the smaller hydrophobic residues Ala, Leu, Ile, and Val (Figure 5a–c). The EK1 peptide consist of Leu, Ile, and Val residues in addition to the larger hydrophobic residues Phe, Met, and Tyr (Figure 5d–f).
The HR2 domain has several charged residues which can participate in electrostatic interactions with the HR1 domain (Figure 6). Several key interactions were noted in the HR1, HR2 complex: E1202 to K921, both E1195 and E1188 to K933, and E1182 to K947 (Figure 6c,d,i,j). EK1 also has several electrostatic interactions: E41 to K921, E34 to K933, E21 to K947, K23 to D950 (Figure 6e,f,k,l). It was noticed that the EK1 also had two electrostatic clashes in region A and one in region B. Residue K31 clashes with the positively charged region created by K933; E27 clashes with D936, and E19 clashes with D950 (Figure 6k,l). Correction of these electrostatic clashes was demonstrated to enhance binding affinity which will be discussed in terms of the peptide design of inhibitors.

3.3. Design of HR1 Fusion Inhibitor Peptides

Peptide inhibitors were rationally designed through an iterative process to maximize hydrogen bonding, hydrophobic, and electrostatic interactions. Initial designs having only one or two point mutations were utilized to explore the impact of amino acid substitutions on binding to the HR1 domains using the EK1 peptide as a starting scaffold. Successive in silico mutations were then introduced based on the docking scores produced using HDOCK to maximize attractive molecular interactions between the HR1 domain and peptide inhibitors. A total of 84 peptide sequences were produced and evaluated. Protein–protein docking of inhibitors with the SARS-CoV-2 WT HR1 site was performed on all inhibitors using the HDOCK server. A full list of sequences and HDOCK binding scores are provided in Supplementary Table S3. Data show that 57 of the designed inhibitors produced docking scores better than the EK1 peptide inhibitor (Table 1), demonstrating stronger affinity for HR1 than the EK1 peptide. The 30-mer peptide inhibitor sequences and docking scores for the top five inhibitors, 130A2, 130A4, 132B, 132C, and 132D; HR2; and EK1 peptide are shown in Table 1. The top five peptide sequences demonstrated similar binding interactions and are further analyzed below to demonstrate the preferences in binding of the HR1 domain.

3.4. Analysis of Binding Interactions for Top Five Peptide Inhibitors

An analysis of the molecular interactions enhancing binding of the top five inhibitors (130A2, 130A4, 132B,132C, and 132D) to the SARS-CoV-2 WT HR1 produced similar modes of binding (Figure 7). Several new molecular interactions which had not been accessed by the HR2 domain or the EK1 peptide were identified. Each inhibitor formed 13 hydrogen bonding interactions with the WT SARS-CoV-2 HR1 domain. These inhibitors formed a hydrogen bonding interaction with residue Q935 which was not found to participate in hydrogen bonding between either the SARS-CoV-2 HR1 WT and HR2 domain or the binding of EK1 by the SARS-CoV-2 WT.
Hydrophobic interactions were optimized by the addition of a WWWEW motif in the C terminal portion of the peptides (Figure 8a). The WWWEW motif introduced new hydrophobic interactions between L922 and V915 of SARS-CoV-2 WT which were not observed with the HR2 sequence or the EK1 peptide. Residues W29 and W26 of the designed peptides participate in these hydrophobic interactions. In the helix binding region, a new hydrophobic interaction with A942 was garnered by including a Trp residue at position 15 (Figure 8b). A second new hydrophobic interaction was garnered by introducing a Met residue at position 10, which was demonstrated to interact with A944. Residues L3 and W1 found on the N terminal side of the designed peptides also produced new hydrophobic interactions not observed with either the HR2 domain or the EK1 peptide with the SARS-CoV-2 WT domain (Figure 8c). The L3 and W1 residues interact hydrophobically with residues V951 and A958 of the SARS-CoV-2 HR1 WT domain, respectively.
Therefore, the residue E19 was introduced in the Optimization of the electrostatic interactions between the designed peptides and HR1 domain was completed by removing the electrostatic clashes present with EK1 (Figure 9). In the helix region of the peptide, in silico mutations were introduced which could produce electrostatic interactions with the HR1 domain. Additionally, positions of certain residues were swapped to optimize the attractive electrostatic interactions present. For clarity, we describe any shifts of residues using the corresponding EK1 peptide positions: E19 of the EK1 peptide which clashes with D950 of the HR1 domain and E27 of EK1 which clashes with D936. Residues E19 and E27 were mutated to R8 and R17, respectively, to remove these electrostatic clashes with residues D950 and D936, respectively. Both changes produced successful electrostatic interactions and removed the aforementioned electrostatic clashes. Residue K31 on the EK1 peptide was near residue K933 of the HR1 domain. It was also hypothesized that the introduction of a negatively charged residue at position 30 on EK1 would promote an electrostatic interaction designed peptides to electrostatically interact with K933. This substitution was also demonstrated to produce an electrostatic interaction supporting stronger binding affinity.
The K30 residue of EK1 was not demonstrated to participate in any electrostatic interactions with the HR1 domain. This residue, however, if shifted to position 29 of the EK1 peptide, was hypothesized to be able to interact electrostatically with D936. Thus, residue R18 was introduced, which corresponds to a shift from position 30 to 29 on EK1 in the peptide inhibitors designed in this study. The results demonstrated that residue R18 produced an electrostatic interaction with D936. In the C terminal regions of the peptide, it was also observed that residue D38 of EK1 did not form any electrostatic interactions with the HR1 domain. It was hypothesized that a shift to position 37 would allow for the formation of an electrostatic interaction with K921. Therefore, the residue E28 was introduced in the designed peptides to strengthen the electrostatic interactions with K921. This change was demonstrated to produce the desired electrostatic interaction.

3.5. Analysis of Potential Pan-SARS-CoV-2 HR1 Variant Capability of Top Five Inhibitors

The top five inhibitors produced HDOCK scores that outperformed the EK1 peptide were further evaluated to determine if they had potential to inhibit the emerging HR1 variants of SARS-CoV-2 that might emerge with the A930V and D936Y mutations. Protein–protein docking to evaluate the binding of 130A2, 130A4, 132B, 132C, and 132D inhibitors, HR2 peptide, and EK1 peptide for comparison to the SARS-CoV-2-A930V, SARS-CoV-2-D936Y, and SARS-CoV-2-A930V-D936Y variants was also performed using HDOCK. The resulting docking scores are shown in Table 2. All peptides demonstrated the best scores in the SARS-CoV-2 WT HR1 domain except for the HR2 peptide which was demonstrated to bind slightly better to the A930V mutant. The top five peptide inhibitors outcompeted the EK1 inhibitor in each of the variant HR1 domains.

3.6. Molecular Dynamics Simulations of Top Five Binding Inhibitors

To confirm binding affinity for the best inhibitors, molecular dynamics (MD) simulations were performed to predict the binding affinity of the top five inhibitors, EK1, and HR2 peptides. The top five inhibitors were shown to have higher binding free energies in comparison to the EK1 peptide (Table 3). The data support that these compounds would have strong binding to the HR1 active sites and could serve as a peptide-based fusion inhibitor to prevent coronavirus infections. The bfactor and RMSD plots of the MD simulations are shown in Supplementary Figures S2 and S3, respectively. The RMSD fluctuates in the range of 1.6–3.5 Å, not indicating a large conformational change. Additionally, the alpha helical form was demonstrated to be preserved after MD simulations. An example Ramachandran plot of the 132B peptides is shown in Supplementary Figure S4.

4. Discussion

The top five inhibitors have two fewer hydrogen bonds to the SARS-CoV-2 WT HR1 domain compared to the EK1 peptide. These compounds, however, still outcompeted EK1. Thus, introduction of new hydrophobic contacts and the optimization of the electrostatic interactions proved to be an important feature for binding affinity and should be strongly considered in the design of new peptide inhibitors.
While many of the hydrophobic residues are located in the inner gorge, there are still some surface accessible hydrophobic contacts on the rim of the SARS-CoV-2 HR1 domain. To fully maximize hydrophobic interactions, designs should incorporate residues which can target these portions of the HR1 domains of SARS-CoV-2 variants. Including larger aromatic rings to create contacts with the rim of the gorge was demonstrated to have a positive effect on binding affinity. This is most likely because the larger but planar aromatic rings would have an increased surface area, facilitating more hydrophobic contacts. The planarity of the aromatic side chain, which also contributes to rigidity, allows the ring to “slide” in between the grooves on the outer rim. The C terminal WWWEWE motif of the designed peptides demonstrates this placement of the larger aromatic residues in the rim of the HR1 (Figure 8a). The use of Trp and Tyr is also useful in that they can still participate in hydrogen bonding to other residues on the rim of HR1. This was shown to be the case with Y9 and W25, which both sit in a hydrophobic groove, but also hydrogen bond with N960 and Q935, respectively. It is notable that the W25 to Q935 hydrogen bond has not been explored or utilized to our knowledge in the design of the EK1 inhibitor.
The inner gorge was easily accessed with flexible hydrophobic residues such as methionine. Use of larger aromatic hydrophobic residues in the gorge would be more challenging, as placing larger aromatic rings in the smaller inner pockets might contribute to steric clashes, reducing the peptide binding ability. The top five peptide inhibitors in this study have a W27 residue, which is positioned in the N terminal region of the HR1 domain. Residue W2 sits on a very hydrophobic upper wall edge between two smaller inner pockets, essentially laying over the entrance of the hydrophobic pocket, but does not fit inside the smaller inner pocket.
When considering the design of peptide-based inhibitors targeting the HR1 domains of the SARS-CoV-2 WT, electrostatic interactions play an important role in enhancing binding affinity. The strategic alignment of positively and negatively charged areas on the HR1 domain to negatively and positively charge residues, respectively, supports binding through electrostatic interactions. The removal of the electrostatic clashes from the EK1 peptide produced improvements in binding affinity even when only one electrostatic clash was corrected (SNF-M006 versus EK1 peptide; in Supplementary Table S3). The use of a positively charged residue to interact with the 936 position, which currently has two variants (D or Y), is ideal because it can electrostatically interact with D936 or participate in a cation–pi interactions with Y936. However, the cation–pi attraction may not be as favorable.
The top five designed inhibitors were demonstrated to maintain strong HDOCK scores for all SARS-CoV-2-A930V, SARS-CoV-2-D936Y, and SARS-CoV-2-A930V-D936Y HR1 domains. This suggests that they can serve as pan-SARS-CoV-2 HR1 variant inhibitors in a similar fashion to the EK1 peptide. These inhibitors, which have been designed using the post fusion state models of the HR1 domain, would be ideal to prevent the conformational change to the HR1/HR2 coiled/coil hairpin. Ideally, these peptides would bind to the HR1 domain during the pre-fusion native state or during the pre-hairpin transition state. One caveat would be that there are minor structural variations in the side chain conformations of the pre-fusion, pre-hairpin that are not fully represented in the post-hairpin state which the peptides were designed using. The inhibitors designed in this study are, however, predicted to be stronger affinity binders than EK1 through both molecular docking and MD simulations. As EK1 is known to successfully inhibit multiple coronaviruses from invading cells, the results of this study suggest that these inhibitors would also still produce inhibition of this critical step in invasion of the viral capsid. These inhibitors thus merit further investigation as pan-SARS-CoV-2 HR1 variant inhibitors. This is the first report to our knowledge of computational design of pan-SARS-CoV-2 HR1 variant inhibitors, specifically addressing the new A930V and D936Y variants that have emerged.

5. Conclusions

This work was undertaken to design peptides with potential to serve as pan-SARS-CoV-2 HR1 variant inhibitors. Additionally, this study predicted the molecular interactions between the SARS-CoV-2 HR1 domain and the EK1 peptide which facilitate binding between the two macromolecules. Using both classical molecular docking and molecular dynamics simulations, we have demonstrated the design of novel peptides that have potential to serve as pan-SARS-CoV-2 HR1 variant inhibitors. Tools to access hydrogen bonding interactions and the importance of the hydrophobic and electrostatic molecular interactions to enhance binding to the receptor were elucidated through molecular docking. The resulting 30-mer peptides demonstrated a higher affinity for the SARS-CoV-2 WT HR1 domain than the EK1 peptide. The top five inhibitors were demonstrated through molecular docking to maintain a good level of binding affinity when docked against the SARS-CoV-2-A930V, SARS-CoV-2-D936Y, and SARS-CoV-2-A930V-D936Y HR1 domains, therefore, we suggest that they can serve as pan-SARS-CoV-2 HR1 variant inhibitors. The design considerations from this work, coupled with the structural information demonstrating binding poses, will assist in the development SARS-CoV-2 HR1/HR2 fusion inhibitors and additional pan-SARS-CoV-2 HR1 variant inhibitor antivirals to combat the COVID-19 pandemic.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/biophysica1030023/s1, Figure S1: Docked HR2 domain to SARS-CoV-2 HR1 domain (tan ribbon) aligned to 6LXT crystal structure of SARS-CoV-2 HR1/HR2 complex (blue ribbon), Figure S2: MD simulation resulted B-factors for top5, middle 4, and bottom 3 peptide candidates and EK1, Figure S3: RMSD result of the MD simulation, Figure S4: Example Ramachandran plot of the 132B peptide, indicating these residues are all in the helical form, Table S1: Residue variability of WT HR1 domain of SARS-CoV-2, Table S2 Sequences selected by CONSURF for conservation analysis of SARS-CoV-2 HR1 domain, Table S3: Docking scores and sequences for HR2, EK1, and designed peptide inhibitors.

Author Contributions

Conceptualization, S.V.S.; methodology, S.V.S. and Q.C.; formal analysis, S.V.S., F.E.W., S.D.S., Q.C., D.A., S.B., M.B., J.B., C.C., S.F., K.G., M.H., A.L., P.L., G.L., L.M., J.M., N.M., L.N., K.N., I.P., E.R., M.R.R., M.S., Y.S., K.S., U.T. and S.V.; investigation, S.V.S., F.E.W., S.D.S., S.V.S., F.E.W., S.D.S. and Q.C.; writing—original draft preparation and review and editing, S.V.S. and Q.C.; visualization, S.V.S.; supervision, S.V.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study are available in the supplementary materials or on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhu, N.; Zhang, D.; Wang, W.; Li, X.; Yang, B.; Song, J.; Zhao, X.; Huang, B.; Shi, W.; Lu, R.; et al. A Novel Coronavirus from Patients with Pneumonia in China, 2019. N. Engl. J. Med. 2020, 382, 727–733. [Google Scholar] [CrossRef] [PubMed]
  2. Zhou, P.; Yang, X.-L.; Wang, X.-G.; Hu, B.; Zhang, L.; Zhang, W.; Si, H.-R.; Zhu, Y.; Li, B.; Huang, C.-L.; et al. A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 2020, 579, 270–273. [Google Scholar] [CrossRef] [Green Version]
  3. Ksiazek, T.G.; Erdman, D.; Goldsmith, C.S.; Zaki, S.R.; Peret, T.; Emery, S.; Tong, S.; Urbani, C.; Comer, J.A.; Lim, W.; et al. A Novel Coronavirus Associated with Severe Acute Respiratory Syndrome. N. Engl. J. Med. 2003, 348, 1953–1966. [Google Scholar] [CrossRef] [PubMed]
  4. Drosten, C.; Günther, S.; Preiser, W.; Van Der Werf, S.; Brodt, H.R.; Becker, S.; Rabenau, H.; Panning, M.; Kolesnikova, L.; Fouchier, R.A.M.; et al. Identification of a Novel Coronavirus in Patients with Severe Acute Respiratory Syndrome. N. Engl. J. Med. 2003, 348, 1967–1976. [Google Scholar] [CrossRef] [PubMed]
  5. Kuiken, T.; Fouchier, R.A.; Schutten, M.; Rimmelzwaan, G.F.; Van Amerongen, G.; Van Riel, D.; Laman, J.D.; De Jong, T.; Van Doornum, G.; Lim, W.; et al. Newly discovered coronavirus as the primary cause of severe acute respiratory syn-drome. Lancet 2003, 362, 263–270. [Google Scholar] [CrossRef] [Green Version]
  6. Zaki, A.; Van Boheemen, S.; Bestebroer, T.; Osterhaus, A.; Fouchier, R. Isolation of a Novel Coronavirus from a Man with Pneumonia in Saudi Arabia. N. Engl. J. Med. 2012, 367, 1814–1820. [Google Scholar] [CrossRef]
  7. Johns Hopkins University & Medicine. Coronavirus Resource Center. Available online: https://0-coronavirus-jhu-edu.brum.beds.ac.uk/map.html (accessed on 27 March 2021).
  8. Regulatory Affairs Professionals Society. Available online: https://www.raps.org/news-and-articles/news-articles/2020/3/covid-19-vaccine-tracker (accessed on 27 March 2021).
  9. Rambaut, A.; Loman, N.; Pybus, O.; Barclay, W.; Barrett, J.; Carabelli, A.; Connor, T.; Peacock, T.; Robertson, D.L.; Volz, E. Preliminary Genomic Characterisation of an Emergent Sars-Cov-2 Lineage in the UK Defined by a Novel Set of Spike Mutations. Virological Database. Available online: https://virological.org/t/preliminary-genomic-characterisation-of-an-emergent-sars-cov-2-lineage-in-the-uk-defined-by-a-novel-set-of-spike-mutations/563 (accessed on 27 March 2021).
  10. Volz, E.; Mishra, S.; Chand, M.; Barrett, J.C.; Johnson, R.; Geidelberg, L.; Hinsley, W.R.; Laydon, D.J.; Dabrera, G.; O’Toole, Á.; et al. Assessing transmissibility of SARS-CoV-2 lineage B.1.1.7 in England. Nature 2021, 593, 266–269. [Google Scholar] [CrossRef]
  11. Faria, N.R.; Claro, I.M.; Candido, D.; Moyses Franco, L.A.; Andrade, P.S.; Coletti, T.M.; Silva, C.A.; Sales, F.C.; Manuli, E.R.; Aguiar, R.S. Genomic Characterisation of an Emergent SARS-CoV-2 lineage in Manaus: Preliminary Findings. Virological da-Tabase. Available online: https://virological.org/t/genomic-characterisation-of-an-emergent-sars-cov-2-lineage-in-manaus-preliminary-findings/586 (accessed on 27 March 2021).
  12. Deng, X.; Garcia-Knight, M.A.; Khalid, M.M.; Servellita, V.; Wang, C.; Morris, M.K.; Sotomayor-González, A.; Glasner, D.R.; Reyes, K.R.; Gliwa, A.S.; et al. Transmission, infectivity, and antibody neutralization of an emerging SARS-CoV-2 variant in California carrying a L452R spike protein mutation. medRxiv 2021. [Google Scholar] [CrossRef]
  13. Tegally, H.; Wilkinson, E.; Giovanetti, M.; Iranzadeh, A.; Fonseca, V.; Giandhari, J.; Doolabh, D.; Pillay, S.; San, E.J.; Msomi, N.; et al. Emergence and rapid spread of a new severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) line-age with multiple spike mutations in South Africa. medRxiv 2020. [Google Scholar] [CrossRef]
  14. Wang, P.; Nair, M.S.; Liu, L.; Iketani, S.; Luo, Y.; Guo, Y.; Wang, M.; Yu, J.; Zhang, B.; Kwong, P.D.; et al. Antibody resistance of SARS-CoV-2 variants B.1.351 and B.1.1.7. Nature 2021, 593, 130–135. [Google Scholar] [CrossRef]
  15. Madhi, S.A.; Baillie, V.; Cutland, C.L.; Voysey, M.; Koen, A.L.; Fairlie, L.; Padayachee, S.D.; Dheda, K.; Barnabas, S.L.; Bhorat, Q.E.; et al. Efficacy of the ChAdOx1 nCoV-19 Covid-19 Vaccine against the B.1.351 Variant. N. Engl. J. Med. 2021, 384, 1885–1898. [Google Scholar] [CrossRef] [PubMed]
  16. Weisblum, Y.; Schmidt, F.; Zhang, F.; DaSilva, J.; Poston, D.; Lorenzi, J.C.; Muecksch, F.; Rutkowska, M.; Hoffmann, H.-H.; Michailidis, E.; et al. Escape from neutralizing antibodies by SARS-CoV-2 spike protein variants. eLife 2020, 9, 1–31. [Google Scholar] [CrossRef] [PubMed]
  17. Singh, P.K.; Kulsum, U.; Rufai, S.B.; Mudliar, S.R.; Singh, S. Mutations in SARS-CoV-2 Leading to Antigenic Variations in Spike Protein: A Challenge in Vaccine Development. J. Lab. Physicians 2020, 12, 154–160. [Google Scholar] [CrossRef] [PubMed]
  18. Hasöksüz, M.; Kiliç, S.; Saraç, F. Coronaviruses and SARS-COV-2. Turk. J. Med. Sci. 2020, 50, 549–556. [Google Scholar] [CrossRef]
  19. Shang, J.; Ye, G.; Shi, K.; Wan, Y.; Luo, C.; Aihara, H.; Geng, Q.; Auerbach, A.; Li, F. Structural basis of receptor recognition by SARS-CoV-2. Nature 2020, 581, 221–224. [Google Scholar] [CrossRef] [Green Version]
  20. Brielle, E.S.; Schneidman-Duhovny, D.; Linial, M. The SARS-CoV-2 Exerts a Distinctive Strategy for Interacting with the ACE2 Human Receptor. Viruses 2020, 12, 497. [Google Scholar] [CrossRef] [PubMed]
  21. Wrapp, D.; Wang, N.; Corbett, K.S.; Goldsmith, J.A.; Hsieh, C.-L.; Abiona, O.; Graham, B.S.; McLellan, J.S. Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation. Science 2020, 367, 1260–1263. [Google Scholar] [CrossRef] [Green Version]
  22. Xia, S.; Liu, M.; Wang, C.; Xu, W.; Lan, Q.; Feng, S.; Qi, F.; Bao, L.; Du, L.; Liu, S.; et al. Inhibition of SARS-CoV-2 (previously 2019-nCoV) infection by a highly potent pan-coronavirus fusion inhibitor targeting its spike protein that harbors a high capacity to mediate membrane fusion. Cell Res. 2020, 30, 343–355. [Google Scholar] [CrossRef] [Green Version]
  23. Vanpatten, S.; He, M.; Altiti, A.; Cheng, K.F.; Ghanem, M.H.; Al-Abed, Y. Evidence supporting the use of peptides and peptidomimetics as potential SARS-CoV-2 (COVID-19) therapeutics. Future Med. Chem. 2020, 12, 1647–1656. [Google Scholar] [CrossRef]
  24. Jia, Y.; Shen, G.; Zhang, Y.; Huang, K.S.; Ho, H.Y.; Hor, W.S.; Yang, C.-H.; Li, C.; Wang, W.L. Analysis of the mutation dynamics of SARS-CoV-2 reveals the spread history and emergence of RBD mutant with lower ACE2 binding affinity. bioRxiv 2020. [Google Scholar] [CrossRef] [Green Version]
  25. Lokman, S.M.; Rasheduzzaman; Salauddin, A.; Barua, R.; Tanzina, A.Y.; Rumi, M.H.; Hossain, I.; Siddiki, A.Z.; Mannan, A.; Hasan, M. Exploring the genomic and proteomic variations of SARS-CoV-2 spike glycoprotein: A computational biology approach. Infect. Genet. Evol. 2020, 84, 104389. [Google Scholar] [CrossRef]
  26. Ahamad, S.; Kanipakam, H.; Gupta, D. Insights into the structural and dynamical changes of spike glycoprotein mutations associated with SARS-CoV-2 host receptor binding. J. Biomol. Struct. Dyn. 2020, 1–13. [Google Scholar] [CrossRef]
  27. Stoddard, S.V.; Stoddard, S.D.; Oelkers, B.K.; Fitts, K.; Whalum, K.; Whalum, K.; Hemphill, A.D.; Manikonda, J.; Martinez, L.M.; Riley, E.G.; et al. Optimization Rules for SARS-CoV-2 Mpro Antivirals: Ensemble Docking and Exploration of the Coronavirus Protease Active Site. Viruses 2020, 12, 942. [Google Scholar] [CrossRef] [PubMed]
  28. Han, Y.; Král, P. Computational design of ACE2-based peptide inhibitors of SARS-CoV-2. ACS Nano 2020, 14, 5143–5147. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Arun, K.G.; Sharanya, C.S.; Abhithaj, J.; Francis, D.; Sadasivan, C. Drug repurposing against SARS-CoV-2 using E-pharmacophore based virtual screening, molecular docking and molecular dynamics with main protease as the target. J. Biomol. Struct. Dyn. 2020, 20, 1–12. [Google Scholar] [CrossRef]
  30. Shah, B.; Modi, P.; Sagar, S.R. In silico studies on therapeutic agents for COVID-19: Drug repurposing approach. Life Sci. 2020, 252, 117652. [Google Scholar] [CrossRef]
  31. Guy, R.K.; DiPaola, R.S.; Romanelli, F.; Dutch, R.E. Rapid repurposing of drugs for COVID-19. Science 2020, 368, 829–830. [Google Scholar] [CrossRef]
  32. Pawar, A.Y. Combating devastating COVID-19 by drug repurposing. Int. J. Antimicrob. Agents 2020, 56, 105984. [Google Scholar] [CrossRef]
  33. Linsky, T.W.; Vergara, R.; Codina, N.; Nelson, J.W.; Walker, M.J.; Su, W.; Barnes, C.O.; Hsiang, T.-Y.; Esser-Nobis, K.; Yu, K.; et al. De novo design of potent and resilient hACE2 decoys to neutralize SARS-CoV-2. Science 2020, 370, eabe0075. [Google Scholar] [CrossRef]
  34. Dai, W.; Zhang, B.; Jiang, X.-M.; Su, H.; Li, J.; Zhao, Y.; Xie, X.; Jin, Z.; Peng, J.; Liu, F.; et al. Structure-based design of antiviral drug candidates targeting the SARS-CoV-2 main protease. Science 2020, 368, 1331–1335. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Miceli, L.A.; Teixeira, V.L.; Castro, H.C.; Rodrigues, C.R.; Mello, J.F.R.; Albuquerque, M.G.; Cabral, L.M.; De Brito, M.A.; De Souza, A.M.T. Molecular Docking Studies of Marine Diterpenes as Inhibitors of Wild-Type and Mutants HIV-1 Reverse Transcriptase. Mar. Drugs 2013, 11, 4127–4143. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Liu, G.; Wang, W.; Wan, Y.; Ju, X.; Gu, S. Application of 3D-QSAR, Pharmacophore, and Molecular Docking in the Molecular Design of Diarylpyrimidine Derivatives as HIV-1 Nonnucleoside Reverse Transcriptase Inhibitors. Int. J. Mol. Sci. 2018, 19, 1436. [Google Scholar] [CrossRef] [Green Version]
  37. Stoddard, S.V.; May, X.A.; Rivas, F.; Dodson, K.; Vijayan, S.; Adhika, S.; Parker, K.; Watkins, D.L. Design of Potent Panobinostat Histone Deacetylase Inhibitor Derivatives: Molecular Considerations for Enhanced Isozyme Selectivity between HDAC2 and HDAC8. Mol. Inform. 2018, 38, e1800080. [Google Scholar] [CrossRef] [PubMed]
  38. Balasubramaniam, S.; Vijayan, S.; Goldman, L.V.; May, X.A.; Dodson, K.; Adhikari, S.; Rivas, F.; Watkins, D.L.; Stoddard, S.V. Design and synthesis of diazine-based panobinostat analogues for HDAC8 inhibition. Beilstein J. Org. Chem. 2020, 16, 628–637. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Stoddard, S.V.; Dodson, K.; Adams, K.; Watkins, D.L. In silico Design of Novel Histone Deacetylase 4 Inhibitors: Design Guidelines for Improved Binding Affinity. Int. J. Mol. Sci. 2019, 21, 219. [Google Scholar] [CrossRef] [Green Version]
  40. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Xia, S.; Yan, L.; Xu, W.; Agrawal, A.S.; Algaissi, A.; Tseng, C.-T.K.; Wang, Q.; Du, L.; Tan, W.; Wilson, I.A.; et al. A pan-coronavirus fusion inhibitor targeting the HR1 domain of human coronavirus spike. Sci. Adv. 2019, 5, eaav4580. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Glaser, F.; Pupko, T.; Paz, I.; Bell, R.E.; Bechor-Shental, D.; Martz, E.; Ben-Tal, N. ConSurf: Identification of Functional Regions in Proteins by Surface-Mapping of Phylogenetic Information. Bioinformatics 2003, 19, 163–164. [Google Scholar] [CrossRef] [Green Version]
  43. Landau, M.; Mayrose, I.; Rosenberg, Y.; Glaser, F.; Martz, E.; Pupko, T.; Ben-Tal, N. ConSurf 2005: The projection of evolutionary conservation scores of residues on protein structures. Nucleic Acids Res. 2005, 33, W299–W302. [Google Scholar] [CrossRef]
  44. Yan, Y.; Tao, H.; He, J.; Huang, S.-Y. The HDOCK server for integrated protein–protein docking. Nat. Protoc. 2020, 15, 1829–1852. [Google Scholar] [CrossRef]
  45. Yan, Y.; Zhang, D.; Zhou, P.; Li, B.; Huang, S.-Y. HDOCK: A web server for protein–protein and protein–DNA/RNA docking based on a hybrid strategy. Nucleic Acids Res. 2017, 45, W365–W373. [Google Scholar] [CrossRef]
  46. Yan, Y.; Wen, Z.; Wang, X.; Huang, S.-Y. Addressing recent docking challenges: A hybrid strategy to integrate template-based and free protein-protein docking. Proteins 2017, 85, 497–512. [Google Scholar] [CrossRef]
  47. Huang, S.-Y.; Zou, X. A knowledge-based scoring function for protein-RNA interactions derived from a statistical mechanics-based iterative method. Nucleic Acids Res. 2014, 42, e55. [Google Scholar] [CrossRef]
  48. Huang, S.-Y.; Zou, X. An iterative knowledge-based scoring function for protein-protein recognition. Proteins 2008, 72, 557–579. [Google Scholar] [CrossRef]
  49. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. Ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [Google Scholar] [CrossRef] [Green Version]
  50. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. [Google Scholar] [CrossRef]
  51. Li, P.; Merz, K.M. Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10, 289–297. [Google Scholar] [CrossRef] [Green Version]
  52. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H.J.C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. [Google Scholar] [CrossRef] [Green Version]
  53. Case, D.A.; Ben-Shalom, I.Y.; Brozell, S.R.; Cerutti, D.S.; Cheatham, T.E.I.; Cruzeiro, V.W.D.; Darden, T.A.; Duke, R.E.; Ghoreishi, D.; Gilson, M.K.; et al. Amber 2018; University of California: San Francisco, CA, USA, 2018. [Google Scholar]
  54. Srinivasan, J.; Cheatham, I.T.E.; Cieplak, P.; Kollman, A.P.A.; Case, D.A. Continuum Solvent Studies of the Stability of DNA, RNA, and Phosphoramidate−DNA Helices. J. Am. Chem. Soc. 1998, 120, 9401–9409. [Google Scholar] [CrossRef]
  55. Miller, I.B.R.; McGee, J.T.D.; Swails, J.M.; Homeyer, N.; Gohlke, H.; Roitberg, A.E. MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321. [Google Scholar] [CrossRef] [PubMed]
  56. Mills, J.E.J.; Dean, P.M. Three-dimensional hydrogen-bond geometry and probability information from a crystal survey. J. Comput. Mol. Des. 1996, 10, 607–622. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Ribbon structure and conservation of SARS-CoV-2 HR1 domain. (a) Ribbon structure of WT SARS-CoV-2 HR1 domain is shown from N terminus (blue) to C terminus (green). (b) Conservation plot of WT HR1 domain of SARS-CoV-2. (c) Conservation mapped onto ribbon diagram of WT SARS-CoV-2 HR1 domain. Magenta and pink residues indicate strongly or moderately conserved, respectively. Cyan and light blue represent highly or slightly variable, respectively. White indicates regions that are not conserved or variable.
Figure 1. Ribbon structure and conservation of SARS-CoV-2 HR1 domain. (a) Ribbon structure of WT SARS-CoV-2 HR1 domain is shown from N terminus (blue) to C terminus (green). (b) Conservation plot of WT HR1 domain of SARS-CoV-2. (c) Conservation mapped onto ribbon diagram of WT SARS-CoV-2 HR1 domain. Magenta and pink residues indicate strongly or moderately conserved, respectively. Cyan and light blue represent highly or slightly variable, respectively. White indicates regions that are not conserved or variable.
Biophysica 01 00023 g001
Figure 2. Electrostatic surface of HR1 domains of SARS-CoV-2: (a) WT, (b) D936Y, (c) A930V, and (d) A930V-D936Y variants. The red color indicates regions of negative charge, the blue indicates regions of positive charge. The white color indicates regions of neutral charge.
Figure 2. Electrostatic surface of HR1 domains of SARS-CoV-2: (a) WT, (b) D936Y, (c) A930V, and (d) A930V-D936Y variants. The red color indicates regions of negative charge, the blue indicates regions of positive charge. The white color indicates regions of neutral charge.
Biophysica 01 00023 g002
Figure 3. Hydrophobicity plots shown on surface of HR1 domains of SARS-CoV-2: (a) full structure WT. Expanded surface of the regions found to contain mutated residues of the SARS-CoV-2. (b) WT, (c) D936Y, (d) A930V, and (e) A930V-D936Y variants. The blue color indicates hydrophilic regions, the brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. The darker the shade, the greater the degree of either hydrophobicity or hydrophilicity.
Figure 3. Hydrophobicity plots shown on surface of HR1 domains of SARS-CoV-2: (a) full structure WT. Expanded surface of the regions found to contain mutated residues of the SARS-CoV-2. (b) WT, (c) D936Y, (d) A930V, and (e) A930V-D936Y variants. The blue color indicates hydrophilic regions, the brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. The darker the shade, the greater the degree of either hydrophobicity or hydrophilicity.
Biophysica 01 00023 g003
Figure 4. Hydrogen bonding interactions between: (a) the SARS-CoV-2 WT HR1 domain and HR2 domain and (b) SARS-CoV-2 WT HR1 domain and the EK1 peptide. Hydrogen bonds are shown with blue or orange lines. Blue hydrogen bonds fall within the precise geometric criteria of 2.5 angstroms and theta is greater than 90 degrees [56]. Orange hydrogen bonds fall within the tolerance criteria for H-bond relaxed constraints of 0.4 angstroms and 20.0 degrees.
Figure 4. Hydrogen bonding interactions between: (a) the SARS-CoV-2 WT HR1 domain and HR2 domain and (b) SARS-CoV-2 WT HR1 domain and the EK1 peptide. Hydrogen bonds are shown with blue or orange lines. Blue hydrogen bonds fall within the precise geometric criteria of 2.5 angstroms and theta is greater than 90 degrees [56]. Orange hydrogen bonds fall within the tolerance criteria for H-bond relaxed constraints of 0.4 angstroms and 20.0 degrees.
Biophysica 01 00023 g004
Figure 5. Hydrophobic interactions occurring between SARS-CoV-2 HR1 WT domain and HR2 and SARS-CoV-2 WT HR1 and EK1 peptide. HR2 and EK1 residue side chains are shown as either ribbon or side chains. The HR1 domain of the SARS-CoV-2 is shown as a hydrophobic surface where blue color indicates hydrophilic regions and brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. (a) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of HR2 of SARS-CoV-2 WT domain (green) which flanks the core alpha helix. (b) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of HR2 domain (side chains only) of SARS-CoV-2 WT. (c) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of WT SARS-CoV-2 HR2 domain (green) which flanks the core alpha helix. (d) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of EK1 peptide (salmon) which flanks the core alpha helix. (e) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of EK1 peptide (only side chains shown) of SARS-CoV-2 WT. (f) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of EK1 peptide (salmon) which flanks the core alpha helix.
Figure 5. Hydrophobic interactions occurring between SARS-CoV-2 HR1 WT domain and HR2 and SARS-CoV-2 WT HR1 and EK1 peptide. HR2 and EK1 residue side chains are shown as either ribbon or side chains. The HR1 domain of the SARS-CoV-2 is shown as a hydrophobic surface where blue color indicates hydrophilic regions and brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. (a) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of HR2 of SARS-CoV-2 WT domain (green) which flanks the core alpha helix. (b) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of HR2 domain (side chains only) of SARS-CoV-2 WT. (c) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of WT SARS-CoV-2 HR2 domain (green) which flanks the core alpha helix. (d) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of EK1 peptide (salmon) which flanks the core alpha helix. (e) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of EK1 peptide (only side chains shown) of SARS-CoV-2 WT. (f) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of EK1 peptide (salmon) which flanks the core alpha helix.
Biophysica 01 00023 g005
Figure 6. Electrostatic interactions between SARS-CoV-2 WT HR1 and HR2 domains and SARS-CoV-2 WT HR1 domain and EK1 peptide shown as either ribbon diagrams or electrostatic potential plot formats. In electrostatic plots, red color indicates regions of negative charge, blue indicates regions of positive charge, and white color indicates regions of neutral charge. (a,b) SARS-CoV-2 WT HR1 binding region that interacts with the C terminal portion of the peptides. (c,d) C terminal region of the HR2 and SARS-CoV-2 HR1 WT domain. (e,f) C terminal region of the SARS-CoV-2 WT HR1 binding site and EK1. (g,h) Helix binding region of SARS-CoV-2 WT HR1 binding site and EK1. (i,j) SARS-CoV-2 WT HR1 and helix region of HR2 domain. (k,l) SARS-CoV-2 WT HR1 and helix region of EK1. Region corresponding to the binding site for the N terminal portion of SARS-CoV-2 not shown.
Figure 6. Electrostatic interactions between SARS-CoV-2 WT HR1 and HR2 domains and SARS-CoV-2 WT HR1 domain and EK1 peptide shown as either ribbon diagrams or electrostatic potential plot formats. In electrostatic plots, red color indicates regions of negative charge, blue indicates regions of positive charge, and white color indicates regions of neutral charge. (a,b) SARS-CoV-2 WT HR1 binding region that interacts with the C terminal portion of the peptides. (c,d) C terminal region of the HR2 and SARS-CoV-2 HR1 WT domain. (e,f) C terminal region of the SARS-CoV-2 WT HR1 binding site and EK1. (g,h) Helix binding region of SARS-CoV-2 WT HR1 binding site and EK1. (i,j) SARS-CoV-2 WT HR1 and helix region of HR2 domain. (k,l) SARS-CoV-2 WT HR1 and helix region of EK1. Region corresponding to the binding site for the N terminal portion of SARS-CoV-2 not shown.
Biophysica 01 00023 g006
Figure 7. Hydrogen bonding interactions between the SARS-CoV-2 WT HR1 domain and peptide inhibitor 132D. Hydrogen bonds are shown with blue or orange lines. Blue hydrogen bonds fall within the precise geometric criteria of 2.5 angstroms and theta is greater than 90 degrees [56]. Orange hydrogen bonds fall within the tolerance criteria for H-bond relaxed constraints of 0.4 angstroms and 20.0 degrees.
Figure 7. Hydrogen bonding interactions between the SARS-CoV-2 WT HR1 domain and peptide inhibitor 132D. Hydrogen bonds are shown with blue or orange lines. Blue hydrogen bonds fall within the precise geometric criteria of 2.5 angstroms and theta is greater than 90 degrees [56]. Orange hydrogen bonds fall within the tolerance criteria for H-bond relaxed constraints of 0.4 angstroms and 20.0 degrees.
Biophysica 01 00023 g007
Figure 8. Hydrophobic interactions occurring between SARS-CoV-2 HR1 WT domain and designed peptide inhibitors. Side chains of the designed peptide inhibitors are shown as ribbon structures and/or atoms and bonds. The HR1 domain of SARS-CoV-2 is shown as a hydrophobic surface where blue color indicates hydrophilic regions and brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. (a) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of peptide inhibitor (pink) which flanks the core alpha helix. (b) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of peptide inhibitor (atoms and bonds). (c) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of designed peptide inhibitors (purple) which flanks the core alpha helix.
Figure 8. Hydrophobic interactions occurring between SARS-CoV-2 HR1 WT domain and designed peptide inhibitors. Side chains of the designed peptide inhibitors are shown as ribbon structures and/or atoms and bonds. The HR1 domain of SARS-CoV-2 is shown as a hydrophobic surface where blue color indicates hydrophilic regions and brown color indicates hydrophobic regions. The white color indicates regions neither strongly hydrophilic nor hydrophobic. (a) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the C terminal region of peptide inhibitor (pink) which flanks the core alpha helix. (b) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues in the helix region of peptide inhibitor (atoms and bonds). (c) HR1 of SARS-CoV-2 WT domain (surface) interactions with hydrophobic residues on the N terminal region of designed peptide inhibitors (purple) which flanks the core alpha helix.
Biophysica 01 00023 g008
Figure 9. Electrostatic interactions between SARS-CoV-2 HR1 domain (electrostatic surface) and (a) C terminal region of 130A2, (b) 130A2 helix region, and (c) 132B helix region. The red color indicates regions of negative charge, the blue indicates regions of positive charge. The white color indicates regions of neutral charge. Electrostatic interactions between SARS-CoV-2 HR1 domain, (d) C terminal region of 130A2, (e) 130A2 helix region, and (f) 132B helix region in ribbon and atom and bond format.
Figure 9. Electrostatic interactions between SARS-CoV-2 HR1 domain (electrostatic surface) and (a) C terminal region of 130A2, (b) 130A2 helix region, and (c) 132B helix region. The red color indicates regions of negative charge, the blue indicates regions of positive charge. The white color indicates regions of neutral charge. Electrostatic interactions between SARS-CoV-2 HR1 domain, (d) C terminal region of 130A2, (e) 130A2 helix region, and (f) 132B helix region in ribbon and atom and bond format.
Biophysica 01 00023 g009
Table 1. HDOCK score and sequences of HR2, EK1, and top five peptide inhibitors to SARS-CoV-2-WT HR1 domain.
Table 1. HDOCK score and sequences of HR2, EK1, and top five peptide inhibitors to SARS-CoV-2-WT HR1 domain.
PeptideHDOCK Docking Score (ITScore-PP)Sequence
HR2−756.39DVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQE
132D−702.46WMLYLWLREMMKKMWRMREGMMETWWWEWE
132C−698.67WILYLWLREMMKKMWRMREGMMETWWWEWE
132B−693.35WMLYLWLREMMKKMWRMREGMMETWWWEWE
130A2−693.34WVLYMWLREMMKRMWQRREGMMETWWWEWE
130A4−684.38WVLYLWLREMMKRMWQRREGMMETWWWEWE
EK1−449.70NVTFLDLEYEMKKLEEAIKKLEESYIDLKE *
* Sequence with corresponding 3D structure from 5ZVM.
Table 2. HDOCK score of peptides for HR2, EK1, and top five peptide inhibitors to SARS-CoV-2 HR1 A930V, D936Y, and A930V-D936Y mutants.
Table 2. HDOCK score of peptides for HR2, EK1, and top five peptide inhibitors to SARS-CoV-2 HR1 A930V, D936Y, and A930V-D936Y mutants.
PeptideWT *
(ITScore-PP)
A930V
(ITScore-PP)
D936Y
(ITScore-PP)
A930V-D936Y
(ITScore-PP)
HR2−756.39−758.22−693.65−693.65
132D−702.46−693.10−618.74−612.06
132C−698.67−688.92−616.69−608.98
132B−693.35−688.92−610.76−602.74
130A2−693.34−684.00−587.71−580.86
130A4−684.38−675.04−579.20−575.68
EK1−449.70−449.91−448.01−455.22
* Scores for WT shown for comparison.
Table 3. Predicted binding affinity for top five inhibitors, HR2, and EK1 peptides to WT SARS-CoV-2 HR1 domain.
Table 3. Predicted binding affinity for top five inhibitors, HR2, and EK1 peptides to WT SARS-CoV-2 HR1 domain.
Peptide or InhibitorBinding EnergyPBSA (INDI = 4)
HR2−93.2
132B−75.1
130A2−73.4
132D−73.0
130A4−70.8
132C−66.5
EK1−63.1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Stoddard, S.V.; Wallace, F.E.; Stoddard, S.D.; Cheng, Q.; Acosta, D.; Barzani, S.; Bobay, M.; Briant, J.; Cisneros, C.; Feinstein, S.; et al. In Silico Design of Peptide-Based SARS-CoV-2 Fusion Inhibitors That Target WT and Mutant Versions of SARS-CoV-2 HR1 Domains. Biophysica 2021, 1, 311-327. https://0-doi-org.brum.beds.ac.uk/10.3390/biophysica1030023

AMA Style

Stoddard SV, Wallace FE, Stoddard SD, Cheng Q, Acosta D, Barzani S, Bobay M, Briant J, Cisneros C, Feinstein S, et al. In Silico Design of Peptide-Based SARS-CoV-2 Fusion Inhibitors That Target WT and Mutant Versions of SARS-CoV-2 HR1 Domains. Biophysica. 2021; 1(3):311-327. https://0-doi-org.brum.beds.ac.uk/10.3390/biophysica1030023

Chicago/Turabian Style

Stoddard, Shana V., Felissa E. Wallace, Serena D. Stoddard, Qianyi Cheng, Daniel Acosta, Shaliz Barzani, Marissa Bobay, Jared Briant, Christian Cisneros, Samantha Feinstein, and et al. 2021. "In Silico Design of Peptide-Based SARS-CoV-2 Fusion Inhibitors That Target WT and Mutant Versions of SARS-CoV-2 HR1 Domains" Biophysica 1, no. 3: 311-327. https://0-doi-org.brum.beds.ac.uk/10.3390/biophysica1030023

Article Metrics

Back to TopTop