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

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.


Introduction
The 2019 coronavirus disease  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.

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.

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.

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.

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%.

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.

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 Zn 2+ 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 Zn 2+ 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 (DG bind ) 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. 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.
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 (∆ ) 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.

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.

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 Figures 2 and 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.

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  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.

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. Proteinprotein 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. 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.

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. Suc-

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. 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.

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. 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. with K933. 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.

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 in- 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.

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. Proteinprotein 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.

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.

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.

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://www.mdpi.com/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. 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.