Evolution of Stronger SARS-CoV-2 Variants as Revealed Through the Lens of Molecular Dynamics Simulations.

Using molecular dynamics simulations, the protein-protein interactions of the receptor-binding domain of the wild-type and seven variants of the severe acute respiratory syndrome coronavirus 2 spike protein and the peptidase domain of human angiotensin-converting enzyme 2 were investigated. These variants are alpha, beta, gamma, delta, eta, kappa, and omicron. Using 100 ns simulation data, the residue interaction networks at the protein-protein interface were identified. Also, the impact of mutations on essential protein dynamics, backbone flexibility, and interaction energy of the simulated protein-protein complexes were studied. The protein-protein interface for the wild-type, delta, and omicron variants contained several stronger interactions, while the alpha, beta, gamma, eta, and kappa variants exhibited an opposite scenario as evident from the analysis of the inter-residue interaction distances and pair-wise interaction energies. The study reveals that two distinct residue networks at the central and right contact regions forge stronger binding affinity between the protein partners. The study provides a molecular-level insight into how enhanced transmissibility and infectivity by delta and omicron variants are most likely tied to a handful of interacting residues at the binding interface, which could potentially be utilized for future antibody constructs and structure-based antiviral drug design.

Introduction

The coronavirus disease 2019 (COVID-19) is a highly infectious respiratory disease caused by severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2), a novel coronavirus [1]. SARS-CoV-2 is the causal agent of the worldwide COVID-19 pandemic, which has brought about major health, economic, and social problems with unprecedented consequences. This virus has claimed over 5.5 million lives globally and continues to spread amongst human populations as mutations arise in the SARS-CoV-2 spike (S) protein, one of the four structural proteins that comprise the virus. Since the start of the pandemic, biochemical research has led to the development of various COVID-19 vaccines and communities have established transmission mitigation strategies to prevent disease spread [2]. However, the emergence of SARS-CoV-2 variants is causing major concern in global public health. Variants have shown a higher transmissibility, virulence, re-infection frequency, and resistance to the action of monoclonal and polyclonal antibodies. They have been observed to contain a high number of mutations in the S protein, particularly in the receptor-binding domain (RBD) [3, 4]. Protruding out of the surface of virions, the spike protein undergoes extensive glycosylation [5], a type of post-translational modification that helps the virus particle to avoid being detected by our cellular defense system. Once attached, the primary function of the S protein is to stay latched on to the host cell surface and initiate viral entry. In particular, the RBD of the S protein binds specifically to the protease domain (PD) of the angiotensin-converting enzyme 2 (ACE2), the human cell surface receptor, which is also heavily glycosylated. Recent atomistic simulations, supported by experimental studies demonstrated that the dynamics of the RBD attachment to the host receptor is impacted by these glycans [6-9]. This is a critical initial step that allows SARS-CoV-2 to penetrate human host cells and cause infection [10]. ACE2 is a membrane protein found on the surface of brain, lungs, kidneys, heart, and intestine cells and normally functions to help regulate the cardiovascular system through the hydrolysis of angiotensin II [11].

The RBD of the S protein consists of five antiparallel beta sheets (β1, β2, β3, β4 and β7) with short helices and loops that form the core. There is an insertion between sheets β4 and β7 that contains β5 and β6 strands, α4 and α5 helices, and loops [12]. This inserted sequence consisting of 70 residues (S438-Y508) is known as the receptor-binding motif (RBM) and contains the majority of the contacting residues that bind to ACE2 [13]. There are four pairs of disulfide bonds made from eight cysteines (three pairs in the core, which stabilize the beta-sheet (Cys336–Cys361, Cys379–Cys432, and Cys391–Cys525), and the other pair connects the loops in the RBM (Cys480–Cys488) [12, 14]. The nominal (titratable) charge distribution analysis suggests that at physiological pH, the RBD is likely to remain positively charged, relative to the rest of the S protein [15]. The RBM forms a concave outer surface that interacts with the N-terminal helices of ACE2 (Fig. 1) [12].

Fig. 1

Site of missense mutations observed for the SARS-CoV-2 variants shown on the RBD backbone (shown in blue-colored cartoon representation). Mutation sites: green spheres represent the mutation sites found in alpha, beta, gamma, delta, eta, and kappa variants. Sites color-coded in pink are found in the omicron variant. The ACE2 is shown in red-colored cartoon representation, with a catalytically important zinc ion at its peptidase active site (Color figure online)

ACE2 is a metallopeptidase, whose claw-like extracellular head domain occurs as a homodimer with a deep central channel containing a zinc-bound active site (Fig. 1) [8, 16]. The channel is surrounded by ridges that contain loops, helices, and a fraction of a β-sheet. At the surface of the ACE2, ridges are consisting of two large α-helices containing Ser19 through Tyr83 [17] that serves as a binding motif for the concave RBD of SARS-CoV-2. The two α-helices near the N-terminal (residues 19–54) act as shaft allowing the concave surface of the SARS-CoV-2 to anchor on it. [17, 18]

A large volume of structural studies depicts key aspects of the binding of RBD to the ACE2 receptor leading to the fusion of the viral protein to the human cell [12, 18–22]. Computational studies have also produced interaction details at the binding interface [14, 23–29]. ACE2 consists of a homodimer and structural elucidation suggests that a dimeric ACE2 can bind to two S protein trimers [8, 20]. Upon activation of the spike protein [6, 7] the RBD of the spike protein undergoes a glycan-shielded ‘down’ to ‘up’ (exposed) conformational change that favors binding to the ACE2 [19]. A steric clash between ACE2 and the S-protein trimer constraints the binding of ACE2 homodimer with only the “up” conformation of the RBD. [8, 20] The binding of the S protein to the ACE2 further triggers a cascade of events. The S protein has two subunits—S1 and S2, and binding of the S1 to the ACE2 receptor facilitates dissociation of the S1 subunit from the S2 subunit by utilizing the host’s cellular proteases [30]. The cleavage forces the S2 subunit to go into a more stable form that is critical for membrane fusion [12, 21]. Therefore, the entry of coronavirus is a complex process that involves a concerted action of receptor-binding and proteolytic processing. [30]

The binding of the S protein to the ACE2 receptor plays a central role in the COVID-19 infection. It is known that mutations in the S protein create new variants of SARS-CoV-2 with varying degrees of transmissibility and virulence. Therefore, a systematic study of the changes in the structural and dynamic properties of complexes formed between RBD of SARS-CoV-2 variants and ACE2 could hold the key information about the viral infectivity. In particular, the identification of the inter-residue interactions between human S protein are critical for assessing the impact of mutations on viral transmissibility. Herein, we report a molecular dynamics (MD) simulation study to probe the RBD-ACE2 complexes of wild-type SARS-CoV-2 and its seven variants. The studied mutants are alpha, beta, gamma, delta, eta, kappa, and omicron [31]. The simulated systems were analyzed to investigate the impact of mutations on the interactions, energetics, and dynamics of the RBD-ACE2 complex. Using the simulation data, an attempt was also made to explore the thermodynamic basis of unusual transmissibility of the wild-type, delta and omicron variants that could shed lights on the infectivity of newer variants to emerge in the future. [32]

Methods and Procedures

Computational Setup

MD simulations were carried out using the NAMD program [33, 34] on the hybrid GPU-CPU Cluster with 61 nodes and 3904 cores at the Blugold Center for High-Performance Computing, UW-Eau Claire. Each node contains two CPU, equipped with AMD EPYC 7452, 2.3 GHz/32-core and connected through Hewlett Packard Slingshot internode connection. Additionally, GPU nodes have NVIDIA Tesla V100S 32 GB GPU cards that enabled NAMD to have hardware acceleration through the use of GPUs or CPU vectorization. Molecular mechanical calculations were carried out using CHARMM36 force field [35, 36]. Non-bonding interactions were modeled using a switching function with a ‘switchdist’ of 10 Å, a cutoff of 14 Å, and a ‘pairlistdist’ of 16 Å. Electrostatic interactions were evaluated using the particle mesh Ewald method [37]. The leapfrog Verlet algorithm [38] was employed for integration and a time step of 2 fs was used to compute atomic velocities and displacements. For conformational sampling, a modified Nosé-Hoover method [39, 40] was employed, where fluctuations in the barostat was controlled using Langevin dynamics [41, 42]. A periodic boundary condition was used, which controls the pressure by dynamically adjusting the unit cell volume and rescaling atomic coordinates. The sampled conformations constitute an isothermal-isobaric (NPT) ensemble, which yielded enthalpic changes for the study.

All mutations and structural editing were performed using a homemade script and were visualized using Visual Molecular Dynamics (VMD) [43]. The structure of the native SARS-CoV-2 S protein in complex with human ACE2 (6M0J [12]) was obtained from the Protein Data Bank (PDB) [44]. Using VMD [43], the complex was separated into two protein segments, which were saved as two separate PDB files. The list of mutations on the RBD was obtained from the Stanford University Coronavirus Antiviral Research Database [31]. Next the RBD of the wild-type S protein was mutated to form a specific mutant variant using a homemade script. This allowed the mutated segment to be built and combined with the ACE2 receptor forming the new assemblies with RBD bound to the ACE2. While setting up these assemblies, a total of 80 water molecules, which were detected during crystal structure determination, were preserved. The zinc ion (Zn2+) bound to the core of the ACE2 was also retained in all models; its non-bonding (Lennard-Jones [45]) interaction parameters included in the CHARMM36 forcefield [35, 36] were used during MD simulations. The RBD-ACE2 complex for all variants were then solvated using TIP3P model of water molecules with a 2.4 Å cutoff between protein and solvent atoms [46]. In order to minimize the solvent volume, the protein–protein complex was rotated in 10-degree increment and a 25 Å padding was maintained in all three dimensions. The negatively charged solvated RBD-ACE2 complex was made neutral by the addition of sodium ions. The neutral solvated RBD-ACE2 complex systems had about 200,800 atoms of which about 12,500 were protein atoms.

Following the solvation and ionization of the system, 50,000 steps of minimization was run using NAMD implementation of Langevin dynamics under periodic boundary conditions [47] and a constant temperature of 298 K. [33, 37] Once the system was optimized and strain on the system was released, a 100 ns MD simulation was carried out for each mutant system.

Interaction Energy Calculation

The intermolecular interaction energies between RBD and ACE2, a major contributor to the binding energy, were calculated for the conformations extracted from the MD trajectories using standard protocol of NAMD [33]. The enthalpies were obtained as a combined electrostatic and van der Waals’ interactions. Uncertainties were calculated using block-averages as discussed in Allen and Tildesley [48] using 25 ns block of data. The entropies were calculated using the standardized procedure as described by Sun et al. [49] The pair-wise interaction free energy, was calculated from Eq. 1 by combining the entropic and the enthalpic contributions:

The difference in the interaction free energy for a particular variant was computed from the difference of the interaction free energies between the variant and that of the wild-type.

Essential Dynamics Analysis

The MD simulation trajectories were further analyzed by studying the principal component of motions [50], which is also known as essential dynamics analysis (EDA). In addition, the root-mean square deviation per frame, the amplitude of principal components, and per-residue fluctuations were determined using the program CARMA [51]. The root-mean square deviation (RMSD) for ith frame is calculated from the square root of the mean square of the deviations, averaged the overall Cα atoms for that specific frame using Eq. 3where N is the number of Cα atoms, and and are the position vectors for the jth Cα atom observed in the ith frame and the 0th (i.e., starting frame), respectively. The per-residue backbone fluctuation was recorded from the average thermal fluctuation of individual residues. The amplitude of fluctuation along a certain principal component was calculated from its maximum and minimum values

The relative amplitude was computed using Eq. 5

Results and Discussion

The list of residues mutated in the RBD of alpha, beta, gamma, delta, eta, and kappa variants are illustrated in Table 1. The mutations on the RBD sequence for the omicron variant are tabulated in Table 2. The simulated systems were studied in detail by performing several analyses of their MD trajectories.

Table 1

List of residues mutated in alpha, beta, gamma, delta, eta, and kappa variants. [31]

Protein systems	Variant name (pango-lineage)	Missense mutations	Location
Alpha	B.1.1.7	N501Y
Beta	B.1.351	K417N, E484K, N501Y
Gamma	P.1	K417T, E484K, N501Y
Delta	B.1.617.2	L452R, T478K
Eta	B.1.525	E484K
Kappa	B.1.617.1	L452R, E484Q

Table 2

Missense mutations present in the receptor-binding domain of the omicron variant. [31]

Protein systems		Variant Name (pango-lineage)B.1.1.529
Omicron	Missense Mutations	G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H
	Location

Quality of the Simulated Data

The quality of the simulation for each mutant was evaluated by following their conformational evolution during the 100 ns of simulation time. For each protein system, conformations were stored after 10 ps and for the duration of 100 ns. The overall conformational change was monitored by measuring how each conformation deviated from the starting conformation during the entire period of 100 ns simulation. In particular, to quantify the conformational change, the root-mean-square deviation (RMSD) of each conformation, averaged over all Cα atoms, with respect to the starting conformation (Figure S1) was plotted for the entire duration of MD simulation i.e., 100 ns. The resulting RMSD plots of these proteins suggest that equilibrations were reached for each system after 20–30 ns (Figure S1). An RMSD change of < 2 Å was observed for the last 50 ns (Figure S1) indicating the stability of these protein complexes.

Analysis of Interactions at the Protein–Protein Binding Interface

Several variants contain more than one mutation in the RBD of the S protein (Tables 1–2) and most of them are located on or near the RBM. The omicron variant has 15 mutations (Table 2). These mutations are expected to change the inter-residue interactions at the protein–protein interface. Therefore, the simulation data was analyzed to examine any alteration in the inter-residue interaction at the interface of the RBD and PD domain. Earlier MD simulation studies [17, 18] reported several key residues at the interface. In particular, a recent study by Fossum et al. [17] identified three distinct zones in WT (Table 3), where interacting residue clusters were observed to be essential in providing a tighter binding between the RBM to the ACE2. As shown in Table 3, the left zone consists of two aromatic–aromatic interactions, namely, Y489(OH)···F28(N) and Y489(OH)···Y83(OH). The central zone is dominated with polar (Q493(OE1)···K31(NZ), Q493(NE2)···E35(OE2)) and ion-pair (K417(NZ)···D30(OD1) interactions. In contrast, the right zone involves one cation–π ((Y505(OH)···R393(NH2)), one anion–π ((Y505(OH)···E37(OE2)), and an ion-pair (R403(NH2) ···E37(OE1)) interactions. These residues were chosen for the subsequent analysis of the SARS-CoV-2 variants.

Table 3

Various interaction distances between residues (in Å) observed in the three zones described in the analysis of protein–protein binding interface

Protein systems/ zones	Inter-residue interactions	WT	Alpha	Beta	Gamma	Delta	Eta	Kappa
Left	Y489(OH)···F28(N)	4.4 (0.7)	4.6 (0.5)	5.4 (1.2)	4.3 (0.7)	4.6 (0.8)	5.3 (0.7)	4.2 (0.5)
	N487(ND2)···Q24(OE1)	4.9 (1.7)	6.5 (1.8)	12.5 (3.8)	4.3 (1.0)	4.2 (1.1)	4.4 (1.2)	4.4 (1.2)
	N487(OD1)···Y83(OH)	3.6 (1.3)	5.1 (1.7)	10.0 (3.3)	2.9 (0.4)	3.0 (0.6)	2.8 (0.4)	3.3(1.0)
	Y489(OH)···Y83(OH)	4.2 (0.8)	5.0 (0.6)	5.5 (1.2)	4.0 (0.7)	4.0 (0.8)	3.1 (0.6)	4.2 (0.6)
Central	K417(NZ)···D30(OD1)	4.9 (2.1)	4.9 (2.5)	6.1 (1.2)b	8.7 (1.3)c	3.9 (1.0)	8.2 (3.1)	4.8 (1.7)
	Y453(OH)···H34(NE2)	3.7 (0.6)	4.2 (0.5)	4.5 (0.7)	4.4 (0.9)	4.3 (0.8)	5.1 (0.7)	4.0 (0.7)
	Q493(OE1)···K31(NZ)	3.5 (1.0)	4.3 (1.4)	4.9 (1.7)	4.1 (1.8)	4.4 (1.7)	5.5 (1.8)	4.2 (1.6)
	Q493(NE2)···E35(OE2)	3.9 (1.0)	5.7 (1.3)	4.5 (1.4)	4.4 (1.4)	4.2 (1.2)	4.2 (1.0)	4.0 (1.1)
Right	Y449(OH)···Q42(NE2)	6.4 (1.7)	6.2 (1.7)	9.4 (5.3)	4.4 (1.4)	5.7 (1.6)	4.7 (1.5)	5.1 (1.4)
	Y449(OH)···D38(OD1)	4.6 (1.6)	3.6 (1.3)	6.6 (5.1)	4.6 (1.2)	4.3 (1.4)	4.5 (1.4)	3.8 (1.2)
	Y449(OH)···K353(NZ)	5.8 (1.4)	9.7 (1.8)	7.5 (4.8)	7.5 (2.1)	5.8 (1.1)	5.5 (1.1)	5.2 (0.9)
	Q498(OE1)···K353(NZ)	5.3 (2.4)	10.4 (1.7)	8.5 (1.4)	8.3 (2.1)	5.4 (2.2)	5.0 (1.9)	4.6 (1.6)
	T500(OG1)···Y41(OH)	3.4 (0.6)	6.7 (0.9)	4.2 (0.8)	3.2 (0.4)	3.2 (0.5)	3.4 (1.0)	3.0 (0.3)
	N501(ND2)···Y41(OH)	4.9 (1.8)	3.4 (0.5)a	5.5 (1.2)a	6.5 (1.0)a	4.0 (0.6)	4.2 (1.2)	3.9 (0.7)
	Y505(OH)···R393(NH2)	3.6 (0.5)	10.7 (3.1)	8.8 (1.2)	7.6 (1.4)	3.6 (0.4)	5.5 (1.9)	4.5 (1.2)
	Y505(OH)···E37(OE2)	3.3 (0.8)	8.8 (2.1)	7.1 (1.4)	7.4 (1.1)	3.2 (0.7)	5.1 (1.7)	4.3 (1.3)
	R403(NH2)···E37(OE1)	6.1 (0.6)	8.0 (1.2)	7.1 (1.0)	9.1 (1.1)	5.9 (0.6)	7.8 (1.5)	8.0 (1.5)

These distances were averaged over an ensemble of conformation collected in 100 ns MD simulation for different variants. The quantity in parenthesis represents the uncertainty in the distances

aFor these mutants, the oxygen atom of the hydroxy group of the tyrosine at 501 position has been used

bFor this mutant, the nitrogen atom of the amide group of the asparagine at 417 position has been used

cFor this mutant, the oxygen atom of the hydroxy group of the threonine at 417 position has been used

Alpha Variant

The alpha variant has only one mutation in the RBD: N501Y (Table 1). Based on the variation of the interaction distances, it appears that the average distances are slightly larger in the left zone and remained unchanged in the central zone as compared to the wild-type (Table 3, columns 1 and 2). The main ion-pair interactions K417(NZ)···D30(OD1) remained intact in the central zone as well as Y453(OH)···H34(NE2) and Q493(OE1)···K31(NZ). In contrast, the interaction distances at the right zone were significantly larger in this variant compared to the wild-type. Although, the interaction between Y501(OH)···Y41(OH) was stronger because of the N501Y mutation, other interactions Y449(OH)···K353(NZ) and Y505(OH)···R393(NH2) became much weaker as indicated by the increase in their interaction distances (Table 3, column 2) as well as the uncertainties indicated in the parenthesis. The steric hindrance due to phenyl rings of Y501, Y449, and Y505 is responsible for this change. This demonstrates that the single missense mutation will decrease the overall interaction between the RBD and ACE2.

Beta Variant

In the beta variant, there are three significant mutations within the RBD: N501Y, K417N, and E484K (Table 1). Residue N501 changes from a polar to a nonpolar side chain, K417 changes from a positively charged side chain to a polar side chain, and E484 mutates from a negatively charged side chain to a positively charged one. The left and right zones exhibited significant loss of interactions, and the central zone appeared to be weakened in this mutant as compared to the wild-type (Table 3, columns 1 and 3). In the left zone, N487 moved away from ACE2 as evident from longer N487(ND2)···Q24(OE1) and N487(OD1)···Y83(OH) interactions distances along with significantly larger uncertainties (Table 3, column 3). Three inter-residue interactions located in the central zone of the RBD, Y453(OH)···H34(NE2), Q493(OE1)···K31(NZ), and Q493(NE2)···E35(OE2), were found to be moderately strong despite the K417N mutation that destroys the ion-pair interactions. In the right zone, many of the inter-residue interactions were seen to have much weaker interactions compared to those for the wild-type SARS-CoV-2 (Table2, columns 1 and 3). In summary, the beta mutants appear to have significantly weaker interactions in all three regions at the binding interface.

Gamma Variant

The gamma variant has three mutations in the RBD: K417T, E484K, and N501Y (Table 1). The K417T is expected to interrupt the native salt bridge interactions, and the E484K mutation changes the charge of the residue from negative to positive, likely altering the electrostatic environment and thereby the local interactions. As previously observed in the case of alpha, N501Y mutation increases steric clashes. Compared to the wild-type, the impact of mutations in the gamma RBD on the left region was minimal, however significantly weakened interactions in the central and right zones (Table 2, columns 1 and 4) were observed. As expected, the T417(OH)···D30(OD1) interaction distances at the central zone were much longer indicating a loss of this interaction. Additionally, many of the inter-residue contacts in the right zone of the RBD for the gamma variant also became much weaker (Table 3, column 4) as compared to the wild-type. In summary, an overall loss of binding strength was observed for the gamma variant as compared to the wild-type. The gamma variant is 2.5–2.8 times more transmissible than the wild-type, therefore the observed loss of binding does not correlate with the enhanced transmissibility.

Delta Variant

The delta variant contains mutations in the RBD: L452R and T478K (Table 1). Both mutations introduce a positively charged residue, likely changing the electrostatic environment within the binding interface. Strong inter-residue interactions were retained throughout the left, central, and right zones of the RBD as seen from the 100 ns simulation (Table 3, column 5). The interactions at the left and central zone remained almost unchanged when compared to the wild-type (Table 3, columns 1 and 5). The strongest interactions observed in the left and central zone of the RBD were N487(OD1)···Y83(OH) and K417(NZ)···D30(OD1), respectively (Table 3, Column 5). In the right zone, the ion-pair interactions R403(NH2)···E37(OE1) was preserved as evident from its smaller (0.6 Å) uncertainties. All other interactions in the right zone were quite comparable to the wild-type (Table 3, columns 1 and 5). This observation is aligned with the hypothesis that the binding of the delta RBD to ACE2 receptor will be as strong as the wild-type.

Eta Variant

A variant of interest, Eta, involves a single missense mutation in the RBD: E484K (Table 1). This mutation introduces a positive charge to the RBD, likely changing the electrostatic interactions at the binding interface. The inter-residue interactions measured during the 100 ns simulation showed strong interactions between N487(OD1)···Y83(OH) in the left zone, Q493(NE2)···E35(OE2) in the central zone, and T500(OG1)···Y41(OH) in the right zone (Table 3, column 6). The key ion-pair interaction between K417 and D30 at the central zone was completely gone. This is because of the charge-charge repulsion caused by the mutated residue K484. The other ion-pair interaction R403(NH2)···E37(OE1), which is preserved in the wild-type and delta variant was also weakened in the right zone (Table 3, column 6). Overall, this analysis predicts a reduction of binding strength for this mutant.

Kappa Variant

The kappa variant contains two mutations: L452R and E484Q (Table 1). The results of the 100 ns simulation (Table 2) showed that the inter-residue interactions in the left and central zones were found to remain as strong as observed for the wild-type (Table 3, columns 1 and 7). Interaction N487(OD1)···Y83(OH) in the left zone and interactions Y453(OH)···H34(NE2) and Q493(NE2)···E35(OE2) in the central zone showed the strongest binding in their respective zone. Most of the inter-residue interactions were found to have strong binding affinities in the right zone as well except the ion-pair interaction R403(NH2)···E37(OE1),.which is weakened. This mutant RBD is therefore expected to exhibit moderately strong binding with ACE2.

Omicron Variant

As illustrated in Table 2, there are 15 missense mutations on the RBD of the omicron variant. The location of these mutations on the RBD backbone is shown in Fig. 1. The main difference in the intermolecular interactions in the omicron variant occurs in the form of a new strong network of electrostatic interactions occurring almost at the center of the binding interface. The strongest among these interactions are located about the middle of the long helix (residues 19–53) of ACE2 (Fig. 2). The residue at site 493 in RBD is mutated from glutamine to arginine and both guanidinium nitrogen interacts strongly with D38 of ACE2. This is evident from the average distances of R493(NH1)···D38(OD2) and R493(NH2)···D38(OD2) of 3.0 Å (Table 4). The central part is also stabilized by new interactions due to the two other mutated residues—S496 (mutated from glycine) and R498 (mutated from glutamine). As indicated in Table 4, both S496 (OG1)···D38(OD1) and R498(NH1)···D38(OD1) interactions are fairly strong with an average of 3.1 Å with standard deviations of 0.5 and 0.7 Å, respectively. Taken together, the mutations in the omicron variant are expected to increase interactions with the ACE2 receptor.

Fig. 2

The key interaction networks in the delta and omicron variants that exhibited strong interactions. The ACE2 and RBD are shown in red and blue ribbons, respectively. The omicron variant has a strong interaction network (top inset panel) at the center of the long α-helix (residues 19–53), while the delta variant exhibits two stronger interaction zones (bottom inset panel), one in the central zone and the other in the right zone separated by the broken black line (Color figure online)

Table 4

Interaction distances between residues (in Å) observed in the three zones of omicron variant, described in the analysis of interactions at the protein–protein binding interface

Zones	Inter-residue interactions	Omicron
Left	Y489(OH)···F28(N)	4.2 (0.5)
	N487(ND2)···Q24(OE1)	5.6 (1.7)
	N487(OD1)···Y83(OH)	5.9 (2.5)
	Y489(OH)···Y83(OH)	4.1 (0.7)
Central	N417(ND2)···D30(OD1)	9.1 (1.8)
	Y453(OH)···H34(NE2)	4.2 (0.8)
	R493(NH1)···D38(OD2)	3.0 (0.7)
	R493(NH2)···D38(OD2)	3.3 (0.5)
	S496 (OG1)···D38(OD1)	3.1 (0.5)
	R498(NH1)···D38(OD1)	3.1 (0.7)
Right	Y449(OH)···Q42(NE2)	6.6 (2.4)
	Y449(OH)···D38(OD1)	8.8 (2.6)
	T500(OG1)···Y41(OH)	3.6 (0.3)
	Y501(OH)···Y41(OH)	5.7 (0.6)
	H505(ND1)···K353(NZ)	5.6 (0.9)
	R403(NH1)···E37(OE1)	7.8 (1.0)

These distances were averaged over an ensemble of conformation collected in 100 ns MD simulation. The quantity in parenthesis represents the uncertainty in the distances

Although, the transmissibility data of all variants [52, 53] does not exactly correlate with the observed pair-wise interaction energetics, the current results demonstrate a tighter binding for delta and omicron, which is consistent with the significantly enhanced infectivity caused by these variants [53]. The trend observed in this computational study correlates well with the binding study conducted using biolayer interferometry method [54], where the delta and omicron was observed to have tighter binding.

Analysis of Protein Dynamics

The essential dynamics analysis enabled to examine how the change in inter-residue interactions at the binding interface impacted the overall dynamics of protein complexes. As reported earlier [17, 18], the principal motion involves an oscillatory protein dynamics—a hinge bending motion arising out of the concave binding surface of the RBD anchored on to the long N-terminal α-helix of ACE2 consisting of residues 19–54 [16, 17]. The motion corresponding to the first principal component of all variants as well as the wild-type has been illustrated in Fig. 3. The red to blue color of the protein backbone represent the starting and ending conformations. As shown in Fig. 3, the spread of the red to blue ribbons suggests that the overall motion is reduced in the wild-type as well as in delta and omicron variants. This is also supported by the calculated changes in the amplitude of PC1 motion relative to that of the wild-type, reported in Table S1. The relative amplitude of PC1 for all variants is positive except delta and omicron variants (Table S1, column1) indicating that the extent of oscillatory motion has reduced for both delta and omicron variants. The relative amplitude for PC2 and PC3 for all mutants did not follow any particular trend and varied in both positive and negative directions indicating that these components of motion either increased or decreased in various mutants as compared to the wild-type system. Since the backbone flexibility of a particular system is the measure of the net dynamics contributed by all PCs, the backbone fluctuation of the RBD in different variants were computed relative to that of the wild-type system and were plotted against the residue number of the RBD sequence (Figure S2). Closer scrutiny of the 70 residue-containing RBM (residues 438 to 508) indicates that the backbone fluctuations decreased substantially for delta, kappa, and omicron variants confirming a more stable complex formation in these mutants. The raw data for the backbone fluctuations of the RBM averaged over 70 Cα atoms are given in Table S1 (right-most column).

Fig. 3

Essential dynamics analysis of the wild-type and various variants, using the 100 ns simulation trajectories. The red to blue color of the backbone represents the starting and ending conformations, respectively. The extent of the backbone fluctuation is a measure of the amplitude of the motion. The higher amplitude of the dynamics correlates with the reduced binding affinity between RBD and ACE2 in these variants (Color figure online)

Pair-Wise Interaction Energy Decomposition

The computation of the RBD-ACE2 interaction energy for the wild-type and mutant variants reveals that the interaction enthalpy is relatively more negative for the wild-type, delta, and omicron variants (Table 5). The interaction entropy was positive and varied significantly among all the variants. Compared to the wild-type, the free energy of interaction for the delta variant is more favorable;  = –28 kcal/mol. Also, the interaction free energy is –147.5 kcal/mol more favorable for the omicron variant with respect to the wild-type S protein (Table 5). The gamma variant exhibits the least favorable interactions with of 115 kcal/mol. These interaction energy calculations demonstrated that some variants may be more contagious compared to others because of stronger inter-residue interactions between the human ACE2 receptor and the RBD of the S protein.

Table 5

Decomposition of the pair-wise interaction energy between the receptor-binding domain (RBD) and the angiotensin-converting enzyme 2 (ACE2) receptor for SARS-CoV-2 variants

Protein systems	Interaction enthalpy Δ_interacH^o(kcal/mol)	Interaction entropy Δ_interacS^o(kcal/(mol.K)	Interaction gibbs free energy Δ_interacG^o(kcal/mol)	ΔΔ_interacG^o(kcal/mol)
Wt	– 233 (4)	0.460 (0.013)	– 371 (4)	0
Alpha	– 220 (13)	0.486 (0.043)	– 366 (16)	6
Beta	– 172 (9)	0.578 (0.030)	– 346 (9)	26
Gamma	– 145 (7)	0.372 (0.023)	– 257 (7)	115
Delta	– 235 (11)	0.546 (0.037)	– 399 (11)	–28
Eta	– 207 (13)	0.524 (0.043)	– 365 (13)	7
Kappa	– 217 (16)	0.429 (0.053)	– 346 (16)	26
Omicron	– 380 (14)	0.461 (0.047)	– 519 (14)	–148

energies were calculated using Eq. 2 for 27 °C. The uncertainties energies and entropies computed using block average method described in Sect. 2.2 are given in parenthesis

The interaction energy for three strong interactions is juxtaposed with their distances (Fig. 4), which demonstrate that the two distinct residue clusters shape the binding between the protein partners, RBD and ACE2. The first of these interaction clusters is located at the central domain and includes primarily of K417···D30 ion-pair (Table 2, Fig. 2) for WT, delta, alpha, and kappa variants. The strength of this network is evident from the interaction energy (Fig. 4a), which ranged from – 60 to – 80 kcal/mol with shorter interacting distances (~ 4 Å). The absence of the lysine residue at position 417 in beta, gamma, and omicron variants resulted in significant destabilization (Fig. 4a). For all but omicron variants, the Q493(OE1)···K31(NZ) hydrogen bonding has only a very small amount of contribution in terms of interaction energy (Fig. 4b). However, the Q493R mutation in omicron led to an alternate set of ion-pair: R493···D38 formation (Fig. 2), which preserved the strong in the central zone interactions (Table 4). This is evident from the strong interaction energy (– 70 kcal/mol) and reduced interacting distances (Fig. 4b). The second interaction cluster occurs across the right zone, where the WT and the delta variant have stronger ion-pair interactions R403(NH2)···E37(OE1) (Fig. 4c,Table 3), compared to rest of the variants.

Fig. 4

Interaction energies and distances for the major electrostatic interactions observed wild-type and its variants during 100 ns simulations. The * in b indicates R493-D38 interaction for the omicron variant. Uncertainties in both measurements are indicated by the error bars

Conclusions

The study of the RBD-ACE2 complexes in wild-type and its seven mutants reveals that a handful of inter-residue interactions at the binding interface is responsible for tighter binding of the S protein to the human receptor. The characterization of these key residue interaction networks was made through a 100 ns MD simulation for each of these protein–protein complexes. Subsequently, a detailed study was carried out to probe inter-residue interactions at the binding interface, essential protein dynamics, backbone flexibility, and interaction energy of these simulated systems.

Analysis of the interactions at the binding interface, essential dynamics, and pair-wise interaction energies indicate that the mutations in the RBD of delta and omicron variants increase favorable interactions between the S protein and human ACE2 receptor, compared to the wild-type. In contrast, the interaction between the receptor binding domain and the ACE2 receptor was much weaker in alpha, beta, gamma, eta, and kappa mutants, which is evident from larger interaction distances between interface residues and elevated pair-wise interaction free energies (Tables 3, 4, and 5). A decrease in the amplitude of the first principal component of the motion was noted for a strongly bound complexes (wild-type, delta, and omicron), which exhibits a classic hinge oscillatory motion originated from the concave surface of the RBD latched on to the long helix of the ACE2 (Figs. 1 and 2). Alpha, beta, gamma, eta, and kappa variants revealed a weakening in RBD-ACE2 interactions, where a large amplitude hinge bending motion was observed (Fig. 3).

Although a direct comparison of the energetics with experimental energies is not possible, a review of the literature provides conflicting results. An experimental study by Burrato et al. using biolayer interferometry method demonstrates that the binding affinity of the RBD for ACE2 in the alpha, beta, and gamma variants decreased, while the binding affinity in eta remained similar to the WT. [54] Thus, computed results of the present study are consistent with the experimental findings. However, Koehler et al. has observed a different trend using atomic force microscopy [55], where variants like delta and omicron exhibited weaker binding affinity compared to the WT, while beta, gamma, and kappa exhibited tighter binding. The discrepancy in the two experimental results warrants more studies on these variants.

One of the key findings of this study is the existence of two distinct residue networks that shape the binding between the protein partners, RBD and ACE2. The first of these interaction clusters is located at the central domain. The cluster consists of K417···D30 ion-pair, which is the strongest in WT and delta (Fig. 4a) but missing in other variants including omicron. However, the Q493R mutation in omicron led to an alternate set of ion-pair: R493···D38 formation (Fig. 2), which preserved the strong interactions (Fig. 4b) throughout the simulation. The second interaction cluster occurs across the right zone, where the WT and the delta variant have moderately strong ion-pair R403(NH2)···E37(OE1) (Fig. 4c, Table 3). For omicron, the Q498R mutation creates a stronger interaction at the right zone and a new ion-pair R498(NH1)···D38(OD1) (Table 4).

The favorable pair-wise interaction energy for WT, delta and omicron demonstrate that these interaction networks play a critical role in the tighter binding and hence could be responsible in significantly greater transmissibility of the WT, delta, and omicron variants. As far as one can predict the thermodynamic consequence of these mutations, a viral RBD, where K417 is preserved like in delta variant and Q493R, G496S and Q498R mutations like in omicron variant, could pose a significant risk in terms of enhanced transmissibility and infectivity for COVID-19 infection. The findings in the present study are of significance as it demonstrates the unique importance of the specific residue network in forging enhanced binding interactions between the protein partners. The information can be utilized for future antibody design and structure-based antiviral drug design.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDB 978 kb)

Supplementary file2 (PDB 977 kb)

Supplementary file3 (PDB 978 kb)

Supplementary file4 (PDB 978 kb)

Supplementary file5 (PDB 977 kb)

Supplementary file6 (PDB 978 kb)

Supplementary file7 (PDB 981 kb)

Supplementary file8 (PDB 977 kb)

Supplementary Materials consist of figures of RMSD of conformations, backbone fluctuations, and a table of relative amplitudes of principal components observed during the 100 ns simulations. Supplementary file9 (DOCX 436 kb)