A Comparative Analysis of SARS-CoV-2 Variants of Concern (VOC) Spike Proteins Interacting with hACE2 Enzyme

In late 2019, the emergence of a novel coronavirus led to its identification as SARS-CoV-2, precipitating the onset of the COVID-19 pandemic. Many experimental and computational studies were performed on SARS-CoV-2 to understand its behavior and patterns. In this research, Molecular Dynamic (MD) simulation is utilized to compare the behaviors of SARS-CoV-2 and its Variants of Concern (VOC)-Alpha, Beta, Gamma, Delta, and Omicron-with the hACE2 protein. Protein structures from the Protein Data Bank (PDB) were aligned and trimmed for consistency using Chimera, focusing on the receptor-binding domain (RBD) responsible for ACE2 interaction. MD simulations were performed using Visual Molecular Dynamics (VMD) and Nanoscale Molecular Dynamics (NAMD2), and salt bridges and hydrogen bond data were extracted from the results of these simulations. The data extracted from the last 5 ns of the 10 ns simulations were visualized, providing insights into the comparative stability of each variant’s interaction with ACE2. Moreover, electrostatics and hydrophobic protein surfaces were calculated, visualized, and analyzed. Our comprehensive computational results are helpful for drug discovery and future vaccine designs as they provide information regarding the vital amino acids in protein-protein interactions (PPIs). Our analysis reveals that the Original and Omicron variants are the two most structurally similar proteins. The Gamma variant forms the strongest interaction with hACE2 through hydrogen bonds, while Alpha and Delta form the most stable salt bridges; the Omicron is dominated by positive potential in the binding site, which makes it easy to attract the hACE2 receptor; meanwhile, the Original, Beta, Delta, and Omicron variants show varying levels of interaction stability through both hydrogen bonds and salt bridges, indicating that targeted therapeutic agents can disrupt these critical interactions to prevent SARS-CoV-2 infection.


Introduction
At the end of 2019, a novel coronavirus designated as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is highly transmissible, emerged in the city of Wuhan, China, and caused an outbreak called coronavirus disease (COVID-19) [1,2].
Unlike the closely related SARS-CoV in the 2003 outbreak, symptom-based screening alone could not contain it [3].It is the third human coronavirus known to co-opt the peptidase angiotensin-converting enzyme 2 (ACE2) to enter the cell [4].This SARS-CoV-2 is a positive-sense, single-stranded RNA virus that is rapidly evolving and continually accrues genomic mutations as it continues to be transmitted [5].Genetic mutations of SARS-CoV-2 will occur by two mechanisms: randomly occurring mutations followed by selection and recombination [6].More than one mutation type can be found at the same position in the spike protein sequence [7].Global Initiative Sharing Avian Influenza Data (GISAID) provides the names for the SARS-CoV-2 lineage; this scientific nomenclature is complex for ordinary people to remember or recall.As a result, people started calling these variants by the place where they were detected.To simplify this, the WHO Virus Evolution Working Group has recommended using letters of the Greek alphabet for naming variants of SARS-CoV-2 [8].These highly mutated forms of SARS-CoV-2 had enhanced transmission rates relative to previous variants and were termed Variants of Concern (VOCs)-designated Alpha, Beta, Gamma, Delta, and Omicron [2,8].
The Alpha variant of the COVID-19 virus (B.1.1.7)was first identified in the United Kingdom in September 2020 [9].Within one month, two rapidly growing lineages with large numbers of genetic changes were reported from South Africa [10] and Brazil [11].The B.1.351(Beta) variant rose in prevalence in South Africa from 11% in October to 87% by December [12].The P.1 (Gamma) variant emerged in Manaus, Brazil, a region estimated to have achieved an infection rate approaching 75% by October 2020 but experienced a surge in new cases beginning in November 2020 [13].Subsequently, the Delta variant (B.1.617.2) increased in prevalence from 2% in February 2021 to 87% in May 2021 in Maharashtra, India, as India experienced a dramatic surge in cases [14].Later on, The Omicron variant of SARS-CoV-2 (B.1.1.529)was first identified in South Africa and Botswana and was reported to the World Health Organization (WHO) on November 24, 2021 [15].The emergence of these variants highlights the ongoing evolution of SARS-CoV-2.Continued surveillance and monitoring of the virus's genetic changes will be crucial in controlling the ongoing pandemic and developing effective treatments and vaccines [16].
Experimental studies have revealed that spike proteins of both SARS-CoV-2 and SARS-CoV bind to ACE2 before entering the cell for replication [17].Traditional wetlab protein structural studies techniques involved nuclear magnetic resonance, X-ray crystallography, and cryo-electron microscopy [18].Those traditional methods were widely used before because they can investigate the dynamics and flexibility of molecules in their natural solution state and detailed views of molecular interactions.But most of the techniques are restricted by size limitations, lab setting requirements, budgets, time consumption, etc.The structural analysis finds out the different parts of mutation among those variants.The structural comparison illustrates that some variants have mutations that occurred in some locations different than the others, and some of the variants have fewer mutations while some variants have more mutations.
Many computational methods are used to study SARS-CoV-2 and its variants .For example, Khan revealed that the B.1.618variant is an antibody-escaping variant with slightly altered ACE2-RBD affinity from the perspectives of binding affinity and hydrogen bonding network [45].Kumar compared the spike proteins of Omicron and Delta variants.He observed a disorder-order transition in the Omicron variant between spike protein RBD regions 468-473, which may significantly influence disordered residues/regions on spike protein stability and binding to ACE2 [46].These analyses will help researchers to understand the structural and binding changes in the different domains and to assess the future repercussions of these changes.AlphaFold is a computational tool in the protein structural domain that can be used in protein structure prediction [47,48].The original virus and variants like Delta and Omicron sequences were predicted using I-Mutant3.0,a support vector machine-based tool for predicting protein stability changes resulting from single point mutations [46].Protein modeling, structure-based virtual screening, and comparative docking are the other methods used [49].Moreover, electrostatic intermolecular interactions are important and have been studied by different research groups [50,51], including pHdependency potential and energy of association of spike protein RBDs and hACE2 [52,53].
In this research, computational techniques, including Molecular Dynamic Simulation (MD simulation) and protein-protein interaction and analysis, are utilized to study the interaction between the spike protein of SARS-CoV-2 and its variants (Alpha, Beta, Gamma, Delta, and Omicron) and the human ACE2 protein.Tools like NAMD2, VMD, Chimera, and Prodigy aided in aligning structures, comparing mutations, performing MD simulations, and extracting salt bridges and hydrogen bonds' data from the results of MD simulations.Our analysis concludes that the Original and Omicron are most structurally alike.The Gamma variant forms the strongest interaction with hACE2 through hydrogen bonds, and Omicron easily attracts the hACE2 receptor.The Original, Beta, Delta, and Omicron variants show varying levels of interaction stability through hydrogen bonds and salt bridges, indicating that targeted therapeutic agents can disrupt these critical interactions to prevent SARS-CoV-2 infection.

Structure Comparison
The mutations were identified by comparing the sequences of each protein, and we colored the sequence IDs where differences were found.We compared these variants chronologically (e.g.Original vs. Alpha vs. Beta).The second highest number of muta-tions is between the Original and Alpha variants (Figure 1A).Only two mutations exist between the Alpha and Beta variants (Figure 1B) or the Beta and Gamma variants (Figure 1C).The third highest is between Gamma and Delta (Figure 1D).Notably, most muta-tions occurred between the Delta variant and the Omicron variant (Figure 1E).

Omicron
35.26 26.14The PRODIGY model explains binding affinity changes in viral variants.Mutations lead to different protein-protein contacts (PPCs) and non-interacting surface (NIS) residue compositions.According to the PRODIGY model's predictions, such mutations can either synergistically improve or antagonistically reduce the bond strength.

Preparing Protein Structures
As the initial step of this study, the protein structures of the Original, Alpha, Beta, Gamma, Delta, and Omicron variants of the SARS-CoV-2 coronavirus were acquired from the Protein Data Bank, with 6VW1 [54], 7FEM [55], 7VXD [56], 7NXC [57], 7V8B [58], and 7U0N [59] being their PDB-IDs, respectively.The structural alignment process was crucial to ensure accurate comparisons using the Chimera tool [60].Chimera was used to align all the variants' structures with an original virus (6VW1), where every variant's length was trimmed to the same size as the Original and was compared (Figure 11).This structure alignment can help gain insights into these variants' evolutionary history.Structural similarities can help infer common ancestry and evolutionary relationships [61].It also helped focus on RBD, responsible for interacting with the hACE2 protein.Eliminating the part of the structure not aligned with the original virus decreased the computational time.With the same color of each variant shown above, we noticed that in most variants, the original protein structure has many mutations among these variants.The Alpha and Beta, Beta and Gamma structural comparisons show fewer mutations than other variants.Next, we visualized the protein structures among the original SARS-CoV-2 and other variants.
The RMSD values indicate how similar the two structures are.Using an online tool 2StrucCompare (http://2struccompare.cryst.bbk.ac.uk/index.php,accessed on 1 February 2023), we calculated the RMSD value between Original SARS-CoV-2 and the Alpha variant (1.44 Å), the Alpha and Beta variants (1.48 Å), the Beta and Gamma variants (1.09 Å), the Gamma and Delta variants (0.96 Å), and the Delta and Omicron variants (0.73 Å).The RMSD values of Original SARS-CoV-2 protein and its variants are as follows: Original and Alpha is 1.44 Å; Original and Beta is 0.99 Å; Original and Gamma is 0.70 Å; Original and Delta is 0.79 Å; and Original and Omicron is 0.45 Å.Therefore, the Original and Omicron variants are the two most structurally similar proteins.

Hydrogen Bonds
Based on the MD simulation, the hydrogen bonds at S protein RBDs and hACE2 were calculated, as shown in Figure 2. Figure 2 depicts the hydrogen bonds between the interaction of Original, Alpha, Beta, Gamma, Delta, and Omicron variants with hACE2.Each interaction of variants has 24, 18, 21, 23, 24, and 16 pairs of hydrogen bonds with occupancies more than 75%, with the highest occupancies being 95.8%, 94.91%, 95.5%, 97.2%, 96.4%, and 92.91% for each variant.This shows that Gamma has the highest occupancy among all the variants.The hydrogen bond pairs of Arginine amino acid (Arg) and Glutamic amino acid (Glu) with the highest occupancies are Arg559-Glu564, Arg115-Glu182, and Arg559-Glu564.If the cutoff was picked at 90%, these variants have 3, 1, 3, 2, 2, and 2 hydrogen bonds; this shows that Original and Beta form the highest number of hydrogen bond pairs with occupancy more significant than 90%, followed by the Gamma, Delta, and Omicron variants.This observation strongly suggests that hydrogen bonding significantly contributes to the overall stability of these protein-protein interactions.and Omicron variants.This observation strongly suggests that hydrogen bonding significantly contributes to the overall stability of these protein-protein interactions.Notably, among the hydrogen bond acceptors, the Arginine amino acid consistently appears in all variants, forming hydrogen bonds with high occupancy.Also, the Glutamic amino acid, especially Glu564, is the donor in most high-occupancy hydrogen bonds.Arginine and Aspartic acid are the pivotal amino acids identified in forming these high-occupancy hydrogen bonds across all variants.Their presence at the interaction interfaces signifies their critical role in defining the binding affinity of the viral spike protein to the human cell receptor hACE2.

Salt Bridges
Figure 2 depicts the salt bridges formed between the S protein of SARS-CoV-2 and its variants with hACE2.From the data obtained from the MD simulation, only the last 500 frames were considered for the analysis, as the interaction is unstable during the first 500 frames.The standard line at 3.2 Å for all the salt bridges in Figure 2 is used to determine the strength of the salt bridge.If the distance is consistently above the classic line, it is considered a weak salt bridge, whereas if the space is consistently below the standard line, it is regarded as a strong salt bridge.From Table 1, Original, Alpha, Beta, Gamma, (F) Omicron (pink).Hydrogen bonds with occupancies > 75% were considered for the analysis.The title of each graph shows the number of hydrogen bonds formed during the interaction of the variant with hACE2 and having occupancy > 75%.Each bar in these charts describes the hydrogen bond Acceptor-Donor pair with its occupancy in %.
Notably, among the hydrogen bond acceptors, the Arginine amino acid consistently appears in all variants, forming hydrogen bonds with high occupancy.Also, the Glutamic amino acid, especially Glu564, is the donor in most high-occupancy hydrogen bonds.Arginine and Aspartic acid are the pivotal amino acids identified in forming these highoccupancy hydrogen bonds across all variants.Their presence at the interaction interfaces signifies their critical role in defining the binding affinity of the viral spike protein to the human cell receptor hACE2.

Salt Bridges
Figure 3 depicts the salt bridges formed between the S protein of SARS-CoV-2 and its variants with hACE2.From the data obtained from the MD simulation, only the last 500 frames were considered for the analysis, as the interaction is unstable during the first 500 frames.The standard line at 3.2 Å for all the salt bridges in Figure 3 is used to determine the strength of the salt bridge.If the distance is consistently above the classic line, it is considered a weak salt bridge, whereas if the space is consistently below the standard line, it is regarded as a strong salt bridge.From Table 1, Original, Alpha, Beta, Gamma, Delta, and Omicron have 1, 2, 1, 1, 2, and 1 pair(s) of stronger salt bridges and 2, 0, 2, 0, 0, and 3 pair(s) of weaker salt bridges, respectively; this shows that the Alpha and Delta variants have the highest number of solid salt bridges and Gamma exhibits weaker interaction among all the variants because it has only one pair of strong salt bridge and there are no weak salt bridges.
Further deepening the complexity, our study also examined the specific amino acid pairs responsible for forming these salt bridges.Glutamic acid (Glu) and Lysine (Lys) emerged as the predominant pair, with Glutamic acid frequently functioning as the donor in these interactions and Lysine acting as the acceptor in all the stable salt bridges.Intriguingly, fewer sturdy bridges were formed when Glutamic acid was paired with Arginine, despite Arginine's role as an acceptor.This observation adds another layer to our understanding of the structural biology of these variants, indicating that not all amino acids contribute equally to the stability of salt bridges and, by extension, to the overall structural integrity of the protein.Arginine, despite Arginine's role as an acceptor.This observation adds another layer to our understanding of the structural biology of these variants, indicating that not all amino acids contribute equally to the stability of salt bridges and, by extension, to the overall structural integrity of the protein.All the salt bridges we analyzed in Figure 3 are visualized in the protein structure, as shown in Figure 4.The hACE2 receptors are colored in cyan, spike proteins: Original colored in orange, Alpha colored in gray, Beta colored in green, Gamma colored in yellow, Delta colored in purple, and Omicron colored in pink.Amino acid pairs are colored in different colors for different pairs between hACE2 receptor and spike proteins.The Original, Alpha, Beta, Gamma, Delta, and Omicron variants have 1, 2, 1, 1, 2, and 1 pair(s) of stronger salt bridges and 2, 0, 2, 0, 0, and 3 pairs of weaker salt bridges, respectively; this shows that Alpha and Delta form the more robust interaction with hACE2.When variants with less strong salt bridges are compared, like Original, Beta, Gamma, and Omicron, Gamma exhibits weaker interaction because it has no weak salt bridges, while the other variants have few.

Electrostatic Calculations
Original (Figure 5A) has a more even electrostatic potential distribution, while the positive potential is dominant in distribution; there are some neutral potential and negative potential in the middle of the front view surface and some on the sides, mainly on top and bottom.Alpha (Figure 5B) has more positive potential in the front view of the electrostatic surface and some negative and neutral potential on the side compared to the Original.The neutral potential is the second dominant potential in Alpha.Beta (Figure 5C) is mainly covered by positive potential, while some minor neutral and negative potential are on the side.Delta (Figure 5D) is dominated by positive potential.The negative potential is the second dominant, combined with minor neutral potential.Gamma (Figure 5E) is dominated by positive potential, with some neutral potential on the surface and some negative potential on the side.Omicron (Figure 5F) is dominated by positive potential, with more negative potential in the middle of the front view and top combined with some neutral potential (marked with a circle in Figure 5F).The other variants have a similar electrostatic potential distribution for the front view except for the Original and Omicron variants.The Original and Omicron both have negative and minor neutral potential on the side of the middle part of the front view, which shows high similarities between them.Note that the bottom area is the binding side.
In the top view of the Original (Figure 6A), we can see that the negative potential is more dominant than the positive and neutral potential.The top view of Alpha (Figure 6B) is dominated by neutral and positive potential, and some red possibilities are on the side and center.The top view of Beta (Figure 6C) is mainly dominated by positive and neutral potentials, with some minor weak negative potential on the side.Delta (Figure 6D) combines positive, negative, and neutral possibilities while the negative potentials are scattered.The top view of the center of the Gamma variant (Figure 6E) has more neutral potential, and the positive and negative potential are on the side.Omicron distributes the positive, neutral, and negative potentials more evenly.
neutral potential (marked with a circle in Figure 4F).The other variants have a similar electrostatic potential distribution for the front view except for the Original and Omicron variants.The Original and Omicron both have negative and minor neutral potential on the side of the middle part of the front view, which shows high similarities between them.Note that the bottom area is the binding side.In the top view of the Original (Figure 5A), we can see that the negative potential is more dominant than the positive and neutral potential.The top view of Alpha (Figure 5B) is dominated by neutral and positive potential, and some red possibilities are on the side and center.The top view of Beta (Figure 5C) is mainly dominated by positive and neutral  The top view of the electrostatic surfaces shows that Original and Omicron have similar electrostatic potential distributions, which are more likely dominated by negative potential.The Alpha, Beta, and Delta variants have similar electrostatic potential distributions because they are more likely to be dominated by the positive potential.The Gamma variant has a more even electrostatic potential distribution.
Moving over to the bottom view of the electrostatic surface of these variants, which are the binding sites of these variants.The original (Figure 6A) has a more even electrostatic potential distribution; half of the bottom (the receptor part) is dominated by negative potential along with some neutral potential, and the other half of the base is dominated by positive potential.The bottom part of the Alpha variant electrostatic potential (Figure Moving over to the bottom view of the electrostatic surface of these variants, which are the binding sites of these variants.The original (Figure 7A) has a more even electrostatic potential distribution; half of the bottom (the receptor part) is dominated by negative potential along with some neutral potential, and the other half of the base is dominated by positive potential.The bottom part of the Alpha variant electrostatic potential (Figure 7B) is dominated by positive potential.In contrast, the receptor part is covered by negative and neutral potential on the side.Beta (Figure 7C) is mainly covered by positive potential, but there are some negative potentials in the receptor part and the bottom of graph C. The bottom view of the Delta electrostatic potential (Figure 7D) differs from the previous variants; Delta is covered mainly by positive electrostatic potential, and even the receptor part is entirely covered by positive potential.Gamma (Figure 7E) is also mostly covered by the positive potential.Still, there are some weak negative potentials in the receptor part, some neutral potentials, and the bottom of Graph E. Omicron (Figure 7F) is dominated by positive potential, but the receptor part is covered by neutral and weak negative potential.
According to our previous research results [33,35,36,40], human ACE2 is overall negative, so if the electrostatic potential of these variants is generally dominated by positive potential, that means the binding affinity between the human ACE2 and the variant is high, which means the humans are more easily to be affected by the variant.
Interestingly, the pocket (leftmost outgoing part on each graph) is an active mutation site.By comparing the Original and the Delta variants, the pocket area has more negative potential in Original than in Delta.Then, from the Delta to Gamma variants, the electrostatic potential becomes negative potential again, and from the Gamma to Omicron variants, the pocket is neutral with slightly negative potential.Because the pocket is targeted in drug design, these observations may relate to drug effectiveness design.According to our previous research results [33,35,36,40], human ACE2 is overall negative, so if the electrostatic potential of these variants is generally dominated by positive potential, that means the binding affinity between the human ACE2 and the variant is high, which means the humans are more easily to be affected by the variant.
Interestingly, the pocket (leftmost outgoing part on each graph) is an active mutation site.By comparing the Original and the Delta variants, the pocket area has more negative potential in Original than in Delta.Then, from the Delta to Gamma variants, the electrostatic potential becomes negative potential again, and from the Gamma to Omicron variants, the pocket is neutral with slightly negative potential.Because the pocket is targeted in drug design, these observations may relate to drug effectiveness design.

Hydrophobic Analysis
Since the hydrophobic surface mainly measures the hydrophobicity of amino acid

Hydrophobic Analysis
Since the hydrophobic surface mainly measures the hydrophobicity of amino acid side chains.We can see that the original variant's (Figure 8A) hydrophobic distribution has more negative lipophilicity potential, but in general, the distribution looks random and does not have a specific pattern.We notice that the distribution of positive and negative lipophilicity potentials and some neutral lipophilicity potentials are like the Alpha (Figure 8B) hydrophobic surface.Notice that the shadow of the protein structure causes the black hole region in the center of graph B. The Beta variant's (Figure 8C) negative lipophilicity potential is dominant, but some positive lipophilicity potentials are near the bottom.For the Delta variant (Figure 8D), the overall distribution is still dominated by negative lipophilicity potential.For the Gamma variant (Figure 8E), the overall hydrophobic distribution is dominated by the negative potential, with some positive potentials near the bottom.For Omicron (Figure 8F), the lipophilicity distribution is dominated by negative lipophilicity potential.
By eyeballing, the lipophilicity distribution of the Beta (Figure 8C), Gamma (Figure 8D), and Omicron (Figure 8F) variants are similar, which means they also have similar hydrophilic properties.These hydrophobic surfaces are consistent with the previous electrostatic surfaces; the strong electrostatic potential tends to be more hydrophilic due to interactions with water.
In the hydrophobic top view of the original variant (Figure 9A), the hydrophobic surface is dominated by negative lipophilicity potential, and there are some positive potentials in the center of the top view.For the Alpha variant (Figure 9B), the hydrophobic surface is also dominated by negative potential, with some minor positive and neutral potentials.The negative potential dominates the Beta variant (Figure 9C), but more neutral potentials exist.The negative potential with minor neutral and positive potentials dominates the Delta variant (Figure 9D).The Gamma variant (Figure 9E) is dominated by negative potential.The Omicron variant (Figure 9F) is mainly dominated by negative potential, while some positive potentials are scattered around.Generally, the top view of the hydrophobic surface is similar among these variants.By eyeballing, the lipophilicity distribution of the Beta (Figure 7C), Gamma (Figure 7D), and Omicron (Figure 7F) variants are similar, which means they also have similar hydrophilic properties.These hydrophobic surfaces are consistent with the previous electrostatic surfaces; the strong electrostatic potential tends to be more hydrophilic due to interactions with water.
In the hydrophobic top view of the original variant (Figure 8A), the hydrophobic surface is dominated by negative lipophilicity potential, and there are some positive potentials in the center of the top view.For the Alpha variant (Figure 8B), the hydrophobic sur-  In the bottom view of the original variant (Figure 9A) hydrophobic surface, the area is covered close to evenly by negative and positive potentials along with neutral potentials, and the pocket has slightly positive potential and is close to neutral.The Alpha variant (Figure 9B) is evenly covered by positive and negative potentials and some neural potentials.We also noticed that the Alpha variant's hydrophobic structure differs slightly from the others.We noticed the left part of the receptor is not inward like the other variants (and the outward part has high lipophilicity), so the Alpha variant is more likely to bind to hydrophobic targets other than the desired target.The Beta variant (Figure 9C) has more negative than positive potentials, along with some neutral potentials.Structurally, the Beta variant's rightmost part is different from other variants.For the Delta variant (Figure 9D), the hydrophobic surface combined negative and positive potentials, mostly even with some neutral potentials.We notice that the leftmost pocket shape differs from the other variants; the shape looks like two pockets instead of one.For the Gamma variant (Figure 9E), the positive and negative potentials are almost evenly distributed, and the potentials of the pocket have a slightly positive potential.For the Omicron variant (Figure 9F), the hydrophobic surface is mainly distributed even by the negative and positive potentials along with some neutral potentials, and the pocket has slightly positive potentials and close to neutral potentials.In the bottom view of the original variant (Figure 10A) hydrophobic surface, the area is covered close to evenly by negative and positive potentials along with neutral potentials, and the pocket has slightly positive potential and is close to neutral.The Alpha variant (Figure 10B) is evenly covered by positive and negative potentials and some neural potentials.We also noticed that the Alpha variant's hydrophobic structure differs slightly from the others.We noticed the left part of the receptor is not inward like the other variants (and the outward part has high lipophilicity), so the Alpha variant is more likely to bind to hydrophobic targets other than the desired target.The Beta variant (Figure 10C) has more negative than positive potentials, along with some neutral potentials.Structurally, the Beta variant's rightmost part is different from other variants.For the Delta variant (Figure 10D), the hydrophobic surface combined negative and positive potentials, mostly even with some neutral potentials.We notice that the leftmost pocket shape differs from the other variants; the shape looks like two pockets instead of one.For the Gamma variant (Figure 10E), the positive and negative potentials are almost evenly distributed, and the potentials of the pocket have a slightly positive potential.For the Omicron variant (Figure 10F), the hydrophobic surface is mainly distributed even by the negative and positive potentials along with some neutral potentials, and the pocket has slightly positive potentials and close to neutral potentials.
Among the bottom view of these variants' hydrophobic surface, the distribution of the hydrophobic surface is primarily positive and negative potentials, and there is no apparent trend of the dominant potentials, which contributes to not making the area of binding receptor highly hydrophilic or hydrophobic and make these variants to be steady.Interestingly, the pocket area is also mostly evenly distributed by negative, positive, and neutral potentials around the pocket area.Among the bottom view of these variants' hydrophobic surface, the distribution of the hydrophobic surface is primarily positive and negative potentials, and there is no apparent trend of the dominant potentials, which contributes to not making the area of binding receptor highly hydrophilic or hydrophobic and make these variants to be steady.Interestingly, the pocket area is also mostly evenly distributed by negative, positive, and neutral potentials around the pocket area.

Binding Affinity Calculations
We deployed the PRODIGY tool to analyze the binding affinities and interaction patterns of SARS-CoV-2 variants with ACE2 receptors.Intermolecular contacts were grouped into six classes based on their polar/apolar/charged types, i.e., charged-charged, chargedpolar, charged-apolar, polar-polar, polar-apolar, and apolar-apolar contacts, as shown in Figure 10.Moreover, PRODIGY provided information about Non-Interacting Surface (NIS) properties (Table 2).The binding affinities and dissociation constants were calculated using Equation ( 2), represented in Table 3.This equation effectively quantifies the impact of various interfacial contacts and surface properties on overall binding affinity.For instance, positive coefficients observed for charged-charged and polar-polar interactions indicate their favorable contributions via such mechanisms as salt bridges and hydrogen bonds.On the other hand, negatively signed coefficients for charged-apolar and polar-apolar interactions emphasize their destabilizing effects due to misfit in complementary properties.Additionally, the equation also considers non-interacting surfaces where apolar or charged characteristics can hinder or facilitate solvation or ionic interactions, respectively, thus highlighting how complex the molecular structure-dynamics relationship is.
The Original strain (Figure 10A) has 4 charged-charged, 24 charged-apolar, and 12 apolar-apolar PPCs, resulting in substantial negative ΔG (−12.7 kcal/mol), indicating a strong binding affinity.The 33.33% apolar and 27.02% charged NIS residues, based on the PRODIGY model, create a balanced interface for binding.The Alpha variant (Figure 10B) becomes more hydrophobic with one charged-charged PPC and 41.31% apolar NIS residues, which lowers the binding energy to −11.1 kcal/mol.Beta (Figure 10C) keeps 4

Binding Affinity Calculations
We deployed the PRODIGY tool to analyze the binding affinities and interaction patterns of SARS-CoV-2 variants with ACE2 receptors.Intermolecular contacts were grouped into six classes based on their polar/apolar/charged types, i.e., charged-charged, chargedpolar, charged-apolar, polar-polar, polar-apolar, and apolar-apolar contacts, as shown in Figure 11.Moreover, PRODIGY provided information about Non-Interacting Surface (NIS) properties (Table 2).The binding affinities and dissociation constants were calculated using Equation ( 2), represented in Table 3.This equation effectively quantifies the impact of various interfacial contacts and surface properties on overall binding affinity.For instance, positive coefficients observed for charged-charged and polar-polar interactions indicate their favorable contributions via such mechanisms as salt bridges and hydrogen bonds.On the other hand, negatively signed coefficients for charged-apolar and polar-apolar interactions emphasize their destabilizing effects due to misfit in complementary properties.Additionally, the equation also considers non-interacting surfaces where apolar or charged characteristics can hinder or facilitate solvation or ionic interactions, respectively, thus highlighting how complex the molecular structure-dynamics relationship is.The Original strain (Figure 11A) has 4 charged-charged, 24 charged-apolar, and 12 apolar-apolar PPCs, resulting in substantial negative ∆G (−12.7 kcal/mol), indicating a strong binding affinity.The 33.33% apolar and 27.02% charged NIS residues, based on the PRODIGY model, create a balanced interface for binding.The Alpha variant (Figure 11B) becomes more hydrophobic with one charged-charged PPC and 41.31% apolar NIS residues, which lowers the binding energy to −11.1 kcal/mol.Beta (Figure 11C) keeps 4 charged-charged PPCs and increases apolar-apolar PPCs to 17 with 37.1% apolar NIS residues, slightly lowering binding energy to −11.6 kcal/mol.Gamma (Figure 11D) has no charged-charged PPCs but lots of charged-apolar PPCs (21) and 35.22% apolar NIS residues, maintaining high binding energy at −12.2 kcal/mol, close to the Original.Delta (Figure 11E) shifts to more charged-polar interactions (11) and fewer total PPCs (62) and shows a slight loss in binding energy at −11 kcal/mol.Omicron (Figure 11F) has the highest charged-charged PPCs at 8 but slightly lower total PPC count (64) and 35.26% apolar NIS residues, with a solid binding energy of −11 kcal/mol.
The PRODIGY model explains binding affinity changes in viral variants.Mutations lead to different protein-protein contacts (PPCs) and non-interacting surface (NIS) residue compositions.According to the PRODIGY model's predictions, such mutations can either synergistically improve or antagonistically reduce the bond strength.

Preparing Protein Structures
As the initial step of this study, the protein structures of the Original, Alpha, Beta, Gamma, Delta, and Omicron variants of the SARS-CoV-2 coronavirus were acquired from the Protein Data Bank, with 6VW1 [54], 7FEM [55], 7VXD [56], 7NXC [57], 7V8B [58], and 7U0N [59] being their PDB-IDs, respectively.The structural alignment process was crucial to ensure accurate comparisons using the Chimera tool [60].Chimera was used to align all the variants' structures with an original virus (6VW1), where every variant's length was trimmed to the same size as the Original and was compared (Figure 1).This structure alignment can help gain insights into these variants' evolutionary history.Structural similarities can help infer common ancestry and evolutionary relationships [61].It also helped focus on RBD, responsible for interacting with the hACE2 protein.Eliminating the part of the structure not aligned with the original virus decreased the computational time.

Structural Comparison
The protein structures SARS-CoV-2 and ACE2 (PDB ID 6VW1) and its variants Alpha (PDB ID 7FEM), Beta (PDB ID 7VXD), Gamma (PDB ID 7NXC), Delta (PDB ID 7V8B), and Omicron (PDB ID 7U0N) were downloaded from the Protein Data Bank (PDB) [62].Before performing the structural comparison, we resized the Alpha, Beta, and Omicron protein structures for better comparison performance.Section 2.1 compares the Original and Alpha protein structures, Alpha and Beta, Beta and Gamma, Gamma and Delta, and Delta and Omicron.We found the structural differences for each pair of the variant comparisons by comparing the DNA sequences.Then, we colorized the amino acid side chains of the DNA sequence differences in a ball stick shape.The amino acid side chains play a crucial role in proteins' folding and final 3D structure, which is also helpful in observing protein structural differences and analyzing hydrophobic interactions.

Performing Molecular Dynamic Simulations
Protein Structure Files (PSFs) were generated using the Visual Molecular Dynamics 1.9.4 (VMD) tool [63] for each variant, serving as a crucial component to define the atomistic details of the protein structures, including their composition and connectivity.PSFs contain information about the atoms, bonds, angles, and dihedrals within the protein structure.To simulate the interaction between each variant with hACE2, MD simulation was performed using the NAMD2 tool.Each simulation was performed at 300 k temperature and 2000 minimization steps, followed by 10 million degrees.It resulted in 1000 frames for each simulation of 10 ns.The first 5 ns of the simulation is unstable, So the results of the last 5 ns were considered for a more accurate result [36].

Salt Bridge and Hydrogen Bond Analysis from MD Simulation Results
Using VMD, data on salt bridges and hydrogen bonds were extracted.Hydrogen bonds with occupancies that were more significant than 75% were considered for the analysis (Figure 2).Hydrogen bond data were generated for the last 500 frames, keeping the Donor-Acceptor distance at 4.0 Å.The number of hydrogen bonds formed, their length, and the strength of the interactions were quantified for each variant.Salt bridges were removed with a cutoff at 3.2 Å, and the data from the last 500 frames of the simulation were utilized for the analysis (Figure 3).For hydrogen bonds, the cutoff was set at 4 Å.A total of 3.2 was used as a base reference line to determine the strength of the salt bridges and the salt bridges that were distanced consistently above 6 Å were excluded as they do not contribute much to the stability of the interaction.

Visualizing Salt Bridges and Hydrogen Bonds Data
The salt bridges and hydrogen bonds data, extracted from MD simulation results, were converted into separate values formatted (CSV) files using Python script.These CSV files were then utilized to visualize the results using Microsoft Excel (Figures 2 and 3).Considering the salt bridge data (Figure 2), a line chart was used to visualize the data where each graph represents a salt bridge formed during the interaction, and a standard horizontal black line at 3.2 Å was used to determine the strength of the salt bridge.Hydrogen bonds (Figure 3) were represented using bar graphs where each bar represents the hydrogen bond formed and occupancy of the respective hydrogen bond.The colors for each salt bridge and hydrogen bond were determined to be the same as those for the variants.

Electrostatic Surface
To show the electrostatic surface, we calculated the protein electrostatic potential of each variant by solving the Poisson-Boltzmann equation using Delphi and visualized it in ChimeraX.The Poisson-Boltzmann equation (PBE) is as follows: where φ(r) is the electrostatic potential, o ˛(r) is the dielectric distribution, ρ (r) is the charge density based on the atomic structures, κ is the Debye-Huckel parameter, kB is the Boltzmann constant, and T is the temperature.The calculated electrostatic potentials by Delphi were based on the electrostatic potential maps.The electrostatic potential maps were computed using the PDB2PQR server and APBS tools; we set the pH to 7.0 in the PDB2PQR configuration.After we had the electrostatic potential maps, we colored the protein surface by the electrostatic potential, with red representing the negative charges and blue representing the positive charges.We set the color scale range between −10 kT/Å and 10 kT/Å.

Hydrophobic Surface
The hydrophobicity of protein is determined by the hydrophobicity of the amino acids that compose it [64].The kd hydrophobicity is automatically assigned to amino acid residues as an attribute with the Kyte and Doolittle hydrophobicity scale values [65].Looking at the hydrophobicity distribution of a protein can help us better understand the folding of proteins.

Binding Affinity Calculation
The resulting docking data of each variant were processed and analyzed using the tools of PRODIGY software (https://rascar.science.uu.nl/prodigy/, accessed on 1 February 2023) Field [66].Within PRODIGY, a simple but robust descriptor of binding affinity is implemented based solely on the structural properties of protein-protein complexes [67].
Using the protein-protein critical affinity benchmark of Kastritis et al., 2011 [44], PRODIGY directly correlates the number of interfacial contacts (ICs) in a protein-protein complex and the binding strength.Its predictor is, therefore, based on a simple linear regression of ICs and some of the properties of the non-interacting surfaces (NISs), which have been shown to influence the binding affinity by Kastritis et al., 2014 [68].where ICs xxx/yyy is the number of Interfacial Contacts found at the interface between In-teractor1 and Interactor2 classified according to the polar/apolar/charged nature of the interacting residues (i.e., ICs charged/apolar is the number of ICs between charged and apolar residues).Two residues are defined in contact if any of their heavy atoms is within 5.5 Å.
Based on the predicted binding affinity (∆G) according to Equation (1), the dissociation constant (K d ) is calculated via the following formula: where R is the ideal gas constant (in kcal K −1 mol −1 ), T is the temperature (in K), and ∆G is the predicted free energy.By default, the temperature is set at 298.15 K (25.0 • C).

Conclusions
Our structural comparison reveals significant mutations between the Delta and Omicron variants, with fewer mutations occurring between the Alpha and Beta, Beta and Gamma variants.The RMSD calculations indicate that the Original and Omicron variants are most similar.Hydrogen bond analysis shows that hydrogen bonds significantly contribute to the overall stability of protein-protein interactions between SARS-CoV-2 variants and the hACE2 protein.From our observation, the Gamma variant has the highest hydrogen bond occupancy among all the variants with the top three hydrogen bonds, Arg559-Glu564, Arg115-Glu182, and Arg559-Glu564.This observation indicates that the Gamma variant has the strongest interaction with hACE2.
Salt bridge analysis from MD simulations revealed that the Alpha and Delta variants have the highest number of solid salt bridges, while the Gamma variant shows weaker interactions with only one pair of strong salt bridges.Glutamic acid (Glu) and Lysine (Lys) frequently form stable bridges, whereas Glutamic acid paired with Arginine forms fewer sturdy bridges.This suggests that not all amino acids equally contribute to the stability of salt bridges and the protein's structural integrity.Key residues such as Arginine and Aspartic acid are crucial in forming high-occupancy hydrogen bonds, highlighting potential targets for antiviral strategies.Therapeutic agents can be developed to disrupt these interactions, breaking the connection between SARS-CoV-2 variants and hACE2.
The electrostatic potential distribution shows high similarity between the Original and Omicron variants.Given that human ACE2 is overall negative, variants with predominantly positive potential have high binding affinity, making humans more easily affected.The pocket area, an active mutation site, shows varying potentials: Original's pocket is more negative than Delta's, turns negative again from Delta to Gamma, and becomes neutral with slightly negative potential from Gamma to Omicron.These variations suggest significant implications for drug design and effectiveness.Supported by salt bridge results from MD simulations, these insights are valuable for developing effective therapeutic drugs.
Lipophilicity distribution analysis shows similar hydrophilic properties for the Beta, Gamma, and Omicron variants, consistent with their electrostatic surfaces.The distribution of hydrophobic surfaces is balanced, contributing to the stability of the variants.These analyses are crucial for drug design, providing insights into the variants' hydrophobic features.
Non-interacting surface (NIS) analysis and binding affinity calculations show that mutations lead to different protein-protein contacts (PPCs) and NIS residue compositions.Either apolar-polar contacts (APCs) or charged-apolar contacts (CACs) are dominant in each variant.Omicron has the highest charged-charged PPCs.Throughout the whole mutation process (from Original to Omicron), PPCs, charged-polar contacts tend to have similar amounts.Apolar-apolar contacts are unstable, with APCs and CACs contributing most to protein-protein interaction.
This work enhances the collective scientific understanding of the SARS-CoV-2 virus and its variants, with potential implications for therapeutics and public health measures.The study highlights the importance of using computational tools to uncover the complex interplay of molecular interactions in biological systems, paving the way for future research.

Figure 1 .
Figure 1.Hydrogen bonds analysis between spike proteins and hACE2 RBD with occupancy over 75%.(A) Original (orange); (B) Alpha (gray); (C) Beta (green); (D) Gamma (yellow); (E) Delta (purple); (F) Omicron (pink).Hydrogen bonds with occupancies > 75% were considered for the analysis.The title of each graph shows the number of hydrogen bonds formed during the interaction of the variant with hACE2 and having occupancy > 75%.Each bar in these charts describes the hydrogen bond Acceptor-Donor pair with its occupancy in %.

Figure 2 .
Figure 2. Hydrogen bonds analysis between spike proteins and hACE2 RBD with occupancy over 75%.(A) Original (orange); (B) Alpha (gray); (C) Beta (green); (D) Gamma (yellow); (E) Delta (purple);(F) Omicron (pink).Hydrogen bonds with occupancies > 75% were considered for the analysis.The title of each graph shows the number of hydrogen bonds formed during the interaction of the variant with hACE2 and having occupancy > 75%.Each bar in these charts describes the hydrogen bond Acceptor-Donor pair with its occupancy in %.

Figure 2 .
Figure 2. Formation of salt bridges between hACE2 and spike proteins.(A) Original (orange); (B) Alpha (gray); (C) Beta (green); (D) Gamma (yellow); (E) Delta (purple); (F) Omicron (pink).Each panel depicts the distance in units of Å between donor (left) and acceptor (right) amino acids over time in units of ns.The cutoff is set to 3.2 Å and shown in the black reference line of each graph.All the salt bridges we analyzed in Figure2are visualized in the protein structure, as shown in Figure3.The hACE2 receptors are colored in cyan, spike proteins: Original colored in orange, Alpha colored in gray, Beta colored in green, Gamma colored in yellow, Delta colored in purple, and Omicron colored in pink.Amino acid pairs are colored in different colors for different pairs between hACE2 receptor and spike proteins.

Figure 3 .
Figure 3. Formation of salt bridges between hACE2 and spike proteins.(A) Original (orange); (B) Alpha (gray); (C) Beta (green); (D) Gamma (yellow); (E) Delta (purple); (F) Omicron (pink).Each panel depicts the distance in units of Å between donor (left) and acceptor (right) amino acids over time in units of ns.The cutoff is set to 3.2 Å and shown in the black reference line of each graph.

Figure 4 .
Figure 4. Electrostatic surface (front) of (A): Original, (B): Alpha, (C): Beta, (D): Delta, (E): Gamma, (F): Omicron.The electrostatic potentials are red for negative potential, white for neutral potential, and blue for positive potential, ranging from − 10 kT/e to 10 kT/e.The black circles indicate the difference in shape among all variants.

Figure 5 .
Figure 5. Electrostatic surface (front) of (A): Original, (B): Alpha, (C): Beta, (D): Delta, (E): Gamma, (F): Omicron.The electrostatic potentials are red for negative potential, white for neutral potential, and blue for positive potential, ranging from −10 kT/e to 10 kT/e.The black circles indicate the difference in shape among all variants.

Figure 6 .
Figure 6.Electrostatic surface (top) of (A): Original, (B): Alpha, (C): Beta, (D): Delta, (E): Gamma, (F): Omicron.The electrostatic potentials are red for negative potential, white for neutral potential, and blue for positive potential.The top view of the electrostatic surfaces shows that Original and Omicron have similar electrostatic potential distributions, which are more likely dominated by negative potential.The Alpha, Beta, and Delta variants have similar electrostatic potential distributions because they are more likely to be dominated by the positive potential.The Gamma variant has a more even electrostatic potential distribution.

Figure 7 .
Figure 7. Hydrophobic (front) view of spike proteins, (A): Original, (B): Alpha, (C): Beta, (D): Delta, (E): Gamma, (F): Omicron.The color changes from yellow to green to represent the amino acids lipophilicity potential, ranging from 20 to −20.The circles show the structural differences and lipophilicity potential distribution differences.

Figure 8 .
Figure 8. Hydrophobic (front) view of spike proteins, (A): Original, (B): Alpha, (C): Beta, (D): Delta, (E): Gamma, (F): Omicron.The color changes from yellow to green to represent the amino acids lipophilicity potential, ranging from 20 to −20.The circles show the structural differences and lipophilicity potential distribution differences.
charged-charged PPCs but lots of charged-apolar PPCs (21) and 35.22% apolar NIS residues, maintaining high binding energy at −12.2 kcal/mol, close to the Original.Delta (Figure10E) shifts to more charged-polar interactions(11) and fewer total PPCs (62) and shows a slight loss in binding energy at −11 kcal/mol.Omicron (Figure10F) has the highest charged-charged PPCs at 8 but slightly lower total PPC count (64) and 35.26% apolar NIS residues, with a solid binding energy of −11 kcal/mol.

Table 3 .
The predicted binding affinity (ΔG) and the dissociation constant ( ) for each.

Table 1 .
Number of strong and weak salt bridges for each variant.

Table 2 .
Percentage of non-interacting surface (NIS) residues for each variant.

Table 3 .
The predicted binding affinity (∆G) and the dissociation constant (K d ) for each.