Exploring the Role of Glycans in the Interaction of SARS-CoV-2 RBD and Human Receptor ACE2

Abstract

COVID-19 is a highly infectious respiratory disease caused by the novel coronavirus SARS-CoV-2. It has become a global pandemic and its frequent mutations may pose new challenges for vaccine design. During viral infection, the Spike RBD of SARS-CoV-2 binds the human host cell receptor ACE2, enabling the virus to enter the host cell. Both the Spike and ACE2 are densely glycosylated, and it is unclear how distinctive glycan types may modulate the interaction of RBD and ACE2. Detailed understanding of these determinants is key for the development of novel therapeutic strategies. To this end, we perform extensive all-atom simulations of the (i) RBD-ACE2 complex without glycans, (ii) RBD-ACE2 with oligomannose MAN9 glycans in ACE2, and (iii) RBD-ACE2 with complex FA2 glycans in ACE2. These simulations identify the key residues at the RBD-ACE2 interface that form contacts with higher probabilities, thus providing a quantitative evaluation that complements recent structural studies. Notably, we find that this RBD-ACE2 contact signature is not altered by the presence of different glycoforms, suggesting that RBD-ACE2 interaction is robust. Applying our simulated results, we illustrate how the recently prevalent N501Y mutation may alter specific interactions with host ACE2 that facilitate the virus-host binding. Furthermore, our simulations reveal how the glycan on Asn90 of ACE2 can play a distinct role in the binding and unbinding of RBD. Finally, an energetics analysis shows that MAN9 glycans on ACE2 decrease RBD-ACE2 affinity, while FA2 glycans lead to enhanced binding of the complex. Together, our results provide a more comprehensive picture of the detailed interplay between virus and human receptor, which is much needed for the discovery of effective treatments that aim at modulating the physical-chemical properties of this virus.

1. Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is responsible for the highly contagious coronavirus disease 2019 (COVID-19). It has led to an ongoing global pandemic where emerging mutations may require continuous development of novel therapeutic strategies. Similar to other coronaviruses, SARS-CoV-2 uses the Spike glycoprotein to interact with human host cells during viral infection [1,2]. Specifically, the receptor binding domain (RBD) of the Spike glycoprotein binds the angiotensin-converting enzyme 2 (ACE2) receptor on the host cell (Figure 1). This facilitates fusion of virus and host cell membranes, leading to entry of the virus into the cell.

Since RBD binding to ACE2 is key to viral entry, this interaction is a major target for development of antibody therapeutics and vaccine design. Mutations in the RBD can alter virus-receptor interaction [3,4], and thus, viral infectivity. Recent structural studies have revealed the residue-residue contacts between SARS-CoV-2 RBD and ACE2, suggesting possible interactions that may determine the stability of the complex [5,6,7]. However, from these static data, it is unclear which of these residues are the critical ones that form interactions most frequently. Further, both the viral Spike and host receptor ACE2 are densely glycosylated with asparagine linked N-glycans [8,9]. Since some of these glycans are spatially located near the RBD-ACE2 interface (Figure 1a–d), it is necessary to assess how they may affect the binding affinity of RBD to ACE2. Previous experiments have elucidated some roles of glycans where they, for instance, can impact antibody interactions and epitope exposure [10]. Furthermore, different glycan types with characteristic structures can be critical for pathogen-host interaction [11,12]. Other experiments have shown that disruption of ACE2 glycosylation can impair SARS-CoV-1 viral entry into cells [13]. However, due to structural complexity and heterogeneity of glycans, together with limited instrumental sensitivity [14], it is unclear how individual glycan types could distinctly modulate RBD binding. Such thorough understanding is needed to suggest more precise strategies for future experiments that seek to design effective treatments.

For this purpose, we apply extensive molecular dynamics (MD) simulations (>200 $\mathsf{\mu }$s) that elucidate the detailed interplay between SARS-CoV-2 RBD and host ACE2, and the role of different glycan types during this process. Since glycans can be extremely dynamic, numerous longer-timescale simulations starting from distinctive initial configurations are required, in order to more precisely capture glycan type-specific effects (see SI Methods for details). Here, we study RBD-ACE2 complexes that either have a) mannose-9 (MAN9) glycans attached to six ACE2 residues: i.e., Asn53, Asn90, Asn103, Asn322, Asn432, and Asn546 or b) sialated complex flucosylated 2-antennae (FA2) glycans attached to the same aforementioned ACE2 residues. In both of these glycan-included models, the Spike RBD is built with a single FA2 glycan at Asn343 (see SI Methods for details). While interprotomer interactions in the Spike trimer [16], and the effects of numerous glycans [17] are critical contributors to the conformation of the Spike RBD, the present study focuses on when the RBD of an infection-capable Spike configuration is already bound to human receptor ACE2. Specifically, we aim at elucidating the physical interactions at the interface of the RBD-ACE2 complex, and the energetic factors that impact RBD-ACE2 affinity. For this purpose, we included a single glycan at Asn343 in the RBD, since this is the only Spike glycan that is spatially located near the RBD-ACE2 interface. Moreover, by focusing on local energetic factors for the RBD-ACE2 complex, we provide a theoretical basis that may guide the design of extracellular ACE2 domain as decoy to neutralize viral infection [18,19], and aid the study of RBD immunogens that can lead to potent neutralizing antibody induction [20,21].

To dissect the determinants of RBD-ACE2 binding, we performed all-atom explicit-solvent simulations of three systems: (i) the RBD-ACE2 complex without glycans (Figure 1a), (ii) RBD-ACE2 with six MAN9 glycans on ACE2 and one FA2 glycan on RBD (Figure 1b), and (iii) RBD-ACE2 with six FA2 glycans on ACE2 and one FA2 glycan on RBD (Figure 1c). These three separate sets of simulations comprise of 90 trajectories that have a combined simulated time of more than 225 $\mathsf{\mu }$s. Statistical analysis of these long simulations reveal the key residues of RBD and ACE2, which form binding contacts with higher probabilities. This provides a complementary view to structural data by quantifying the significance of identified interactions. Importantly, we find that this RBD-ACE2 contact signature is not affected by the presence of MAN9 or FA2 glycans, suggesting RBD-ACE2 contacts are inherently robust. Using our simulated interactions as foundation, we evaluate the effect of N501Y, a recent mutation that accounts for the majority of new infections in South East England. Specifically, we show that N501Y in RBD introduces additional stabilizing interactions with Y41 and K353 of ACE2, which can contribute to enhanced infectivity of the new variants that carry N501Y. Furthermore, we show that the glycan on Asn90 of ACE2 (regardless whether MAN9 or FA2) contacts the RBD with distinctly higher probabilities, compared to the remaining glycans on ACE2. This indicates that the Asn90 glycan may play a critical role in binding/unbinding events of RBD, as previously suggested by atomic force microscopy (AFM) experiments [22], biolayer interferometry [23], and other biochemical assays [24]. Finally, an analysis of binding energetics demonstrates that each glycan type is associated with distinct enthalpic contributions to binding affinity. Specifically, we show that MAN9 glycans on ACE2 lead to weaker RBD-ACE2 binding, while the FA2 glycans, with their negatively charged sialic acid tips, stabilize the complex. Together, the results provide a quantitative framework that allows future studies to more precisely modulate the infectivity of SARS-CoV-2.

2. Results

2.1. Contacts between ACE2 and RBD Are Not Altered by the Presence of Glycans

To obtain a detailed picture of how SARS-CoV-2 RBD interacts with human receptor ACE2, it is necessary to identify the residues that form close contacts at the RBD-ACE2 interface. Furthermore, we investigate how this RBD-ACE2 interaction is affected by different glycan species that are present on the ACE2 surface. To this end, we performed three separate sets of all-atom explicit-solvent simulations of (i) the RBD-ACE2 complex without glycans (Figure 1a), (ii) RBD-ACE2 with MAN9 glycans on ACE2 (Figure 1b), and (iii) RBD-ACE2 with FA2 glycans on ACE2 (Figure 1c). For each simulation set, we calculated the probability of ACE2 and RBD residues forming contacts with each other (Figure 2a–c). Here, a contact is considered formed if the smallest distance between the heavy atoms of two amino acid residues is within 4 Å. We chose this distance cutoff of 4 Å to specifically capture contacting residues that form short-range physical interactions at the RBD-ACE2 interface. Moreover, this contact definition is consistent with a recent structural study that has reported a set of RBD-ACE2 contacts, using the 4 Å cutoff [5]. Applying the same criterion for identifying contacts, we will later compare our simulated interactions to the contacts reported from experimental structural studies (see further below).

The contact maps for the non-glycan (Figure 2a), MAN9-included (Figure 2b), and FA2-included (Figure 2c) simulations are highly similar among each other, suggesting that RBD-ACE2 protein contact formation is virtually independent of the presence of glycans. To describe how glycans may affect the movement of the RBD domain relative to ACE2, we calculated the root mean square deviation (rmsd) of the C-alpha atoms of RBD after alignment of ACE2, as a function of time. From this, we calculated the distribution of rmsd values for each simulation set (Figure 2d). The distributions for the non-glycan and glycan-included simulations resemble each other closely, indicating that the glycans do not affect the orientation of RBD relative to ACE2. Specifically, in all simulations, the RBD fluctuates about a native-like basin that is located at small rmsd values ($\mathrm{rmsd}\approx 0.5\mathrm{nm}$, Figure 2d). Small rmsd values suggest that this ensemble corresponds to the RBD-ACE2 configuration used as reference for the rmsd calculations. Here, the structure by Ref. [5] (PDB ID: 6M0J, Figure 1a) was used as reference. Together, the unperturbed results for contact formation (Figure 2a–c) and rmsd distributions (Figure 2d) demonstrate that short-range interactions and relative orientations of RBD-ACE2 are robust to the presence of glycans.

2.2. Quantifying the Relative Significance of RBD-ACE2 Contacts

To further characterize the RBD-ACE2 interaction profile, we defined persistent short-range contacts as those that remain within 4 Å, in at least $60\%$ of the sampled configurations. Our simulations reveal that there are 25 RBD-ACE2 short-range contacts that form persistently, which may be the key interactions in the binding of the complex (Figure 2e and Table S1). The majority of persistent short-range contacts involve residues in the N-terminal helix of ACE2, which interact with the receptor-binding motif (RBM) of RBD (Figure 2e). Representative RBD-ACE2 interactions that occur most frequently in the simulations include: Y453-H34, F456-T27, N487-Y83, Y489-F28, N501-K353, G502-K353, and Y505-K353 (first and second labels denote RBD and ACE2 residues, respectively). These pairs form with at least $90\%$ probability, suggesting they may play a particularly prominent role in the interplay between the virus’ RBD and human receptor (see Table S1 for a complete list of all 25 persistent short-range contacts). In fact, of these seven RBD residues, two are sites where the common mutations Y453F and N501Y occur, and the remaining five are highly conserved among pandemic variants. Specifically, Y453F is a mutation that originates from minks in Denmark and has been found in over 1400 out of ∼314,000 sequences recorded globally (i.e., $0.45\%$ global frequency), as of 16 January 2021. The N501Y RBD mutation, which first became prominent in the B.1.1.7 lineage, has been found in over 18,000 out of ∼314,000 sequences (i.e., $5.7\%$ global frequency). The B.1.1.7 lineage was first detected in the UK in September 2020, and has become a dominant variant of the virus and is believed to be more transmissible. In the next section further below, we apply our contact statistics, as shown in Figure 2a, to assess how the N501Y mutation in B.1.1.7 may facilitate virus-receptor binding.

Our simulated results provide a quantitative evaluation of the interactions reported by recent structural studies. Specifically, Lan et al. [5] have defined a list of contacting residues at the RBD-ACE2 interface for their structure (PDB ID: 6M0J), applying a contact cutoff distance of 4 Å. This is the same definition used to determine short-range interactions in the present study. Similarly, Shang et al. [7] have presented a set of RBD-ACE2 contacts for their structure (PDB ID: 6VW1), which is in line with a distance cutoff of 4 Å. Since experimental structures are based on average coordinates, it is not clear to what extent those experimentally-reported contacts statistically form or break in solution. Here, we use our explicit-solvent simulations to identify which of those contacts are more likely to persist within 4 Å. To evaluate a more comprehensive set of structurally-derived contacts, we further included the structure of Yan et al. (PDB ID: 6M17 [6]), in addition to the aforementioned structures of PDB IDs: 6M0J [5] and 6VW1 [7]. Applying the contact cutoff of 4 Å, these three configurations together implicate a total of 42 contacts at the RBD-ACE2 interface. Our simulations reveal that these interactions can have substantially diverse probabilities of forming, ranging between $p=0.04$–$0.98$ (Table S1). Notably, all 25 persistent short-range contacts captured by the simulations represent a subset of those 42 experimentally-reported interactions. Thus, our analysis demonstrates that while recent structures have revealed the possible interactions, the present calculations quantify their relative involvement, which can aid in estimating the potential impact of mutations. Together, using large sets of simulation data (∼225 $\mathsf{\mu }$s), we have provided an evaluation of host-virus contacts at a longer timescale, and have additionally shown that these short-range physical interactions are not perturbed by the presence of MAN9 or FA2 glycans. Below, we demonstrate how our contact data may be applied to evaluate potential consequences of specific mutations, such as N501Y in the recently prevalent B.1.1.7 lineage that was originally detected in the UK.

2.3. Simulated Interactions Help Assess the Effects of Mutations, Such As N501Y in RBD

The detailed contact interactions, as revealed by our simulations (Figure 2a), provide a basis that can help evaluate the effects of mutations. In recent months, the new B.1.1.7 variant of SARS-CoV-2, carrying the mutation N501Y in the RBD (Figure 3), reportedly accounts for over $60\%$ of new COVID-19 cases in South East England. This N501Y variant is associated with increased affinity to host receptor [3], but the molecular factors responsible for this is unclear. Since the N501Y mutation site is at the RBD-ACE2 binding interface (Figure 3a), we will examine how this mutation may modulate the virus-host interaction, a factor that can determine infectivity. For this, we performed in-silico mutagenesis of N501Y in the RBD-ACE2 complex, using the FoldX software [25]. Utilizing the interaction statistics from Figure 2a as basis, together with a structural and chemical-physical analysis, we argue below how N501Y facilitates the binding of RBD to ACE2, which may contribute to enhanced infectivity.

To describe the consequences of N501Y, we first identify all ACE2 residues that frequently interact with N501, the mutation site in RBD. From the simulations, N501 forms short-range contacts with Y41 and K353 of host ACE2 in $74\%$ and $93\%$ of the sampled configurations, respectively (see Figure 2a and Table S1). The significant interactions provide reason to specifically focus on these two host residues and assess how their contact formation with the N501Y mutation may alter stability (Figure 3b). For this, we note that N501Y would enable the formation of stabilizing cation-$\pi$ interaction with K353 of ACE2. In addition, the longer side chain of N501Y (relative to N501) may facilitate the intermolecular hydrogen-bonding between the OH group of N501Y (tyrosine) and ACE2 residues, including K353 (Figure 3b). It has been shown that hydrogen bonds by tyrosine OH groups can be a significant contributor to stability [26]. Hence, a tyrosine (N501Y) side chain that allows for more favorable hydrogen bonds between molecules would further facilitate binding. Moreover, the N501Y mutation introduces additional stabilizing contributions: i.e., the aromatic-aromatic interaction between N501Y and Y41 (i.e., $\mathsf{\pi }$-stacking; Figure 3b). Notably, N501Y would lead to enhanced hydrophobic effects in the inside of the binding surface. This mutation would therefore create a more protein core-like environment, which stabilizes the RBD-ACE2 coupling. Together, our analysis demonstrates how N501Y implicates numerous energetic factors whose combined effect likely facilitates the binding of the mutated virus to host cells. This interpretation is consistent with the mutational experiments by Starr et al. [3] that have demonstrated enhanced RBD-ACE2 affinity through N501Y.

2.4. RBD Forms Contacts with Asn90 Glycan of ACE2 More Frequently Than Other Glycans

Some glycans on ACE2 are located near the RBD-ACE2 interface (Figure 1b,c). While we have shown that glycans do not alter the contact signature of ACE2 and RBD, it is unclear how and to what extent glycans of the host may interact with viral RBD. Understanding precisely which glycans form close contacts with RBD will provide a more comprehensive view of the factors determining affinity and may suggest additional strategies to modulate binding.

To describe RBD-glycan interaction, we calculated the probabilities of RBD residues forming contacts with any glycan that is attached to ACE2 (Figure 4). This analysis was performed for the simulation sets that either included MAN9 (cf. Figure 1b) or FA2 glycans (cf. Figure 1c) on ACE2. In both MAN9 and FA2 simulations (Figure 4a,b), ACE2 glycans on Asn53, Asn90, Asn103, and Asn322 form contacts with the RBD. Of these four glycans, the ones on Asn53, Asn103, and Asn322 have lower contact probabilities (i.e., $p\le 30\%$), suggesting that their interaction with RBD is more transient despite their proximity to the binding surface. In stark contrast, the remaining glycan on Asn90 of ACE2 contacts the RBD with distinctly higher probabilities ($p\ge 70\%$), for both MAN9 and FA2 species (Figure 4a,b and Figure S1). Specifically, this Asn90 glycan frequently forms contacts with residues of the RBD ranging approximately from 403–417 (Figure 4a–c). Notably, this RBD region does not form any contacts with ACE2 during simulation, regardless of whether glycans were present or absent (Figure 2a–c). Since RBD-ACE2 contacts are independent of the presence of glycans (as discussed earlier), this result indicates that contact formation between RBD and the Asn90 glycan of ACE2 does not compete with the protein-protein interactions between RBD and ACE2.

The prominent role of the Asn90 glycan as elucidated by our long, unrestrained simulations is consistent with previous atomic force microscopy (AFM) measurements and steered molecular dynamics (SMD) simulations (∼270 ns) [22]. In that study, AFM and SMD analyses suggest that the Asn90 glycan can affect the association and disassociation of RBD and ACE2. As shown in Figure 4, the frequent contacts between Asn90 glycan and RBD implies that there are significant steric effects associated with this glycan. Since RBD-ACE2 interface interactions are not perturbed by glycans (as discussed above), the sterics of Asn90 glycan would hinder the unbinding of RBD, hence stabilizing the RBD-ACE2 complex. This result helps explain why, in AFM unbinding experiments [22], separating RBD and ACE2 requires higher forces when ACE2 glycans are not removed. For the case when ACE2 is not bound to RBD, the excluded volume of the Asn90 glycan is likely to interfere with the binding interface, thereby impeding the association of RBD. Consistent with this interpretation, mutation studies have shown that removal of the Asn90 glycan leads to increased RBD binding events [23,24]. Together, our simulations provide a mechanistic basis for how the sterics of Asn90 glycan can impede the unbinding and binding of RBD.

In previous SMD simulations [22], which reported on glycan-RBD interactions, all glycans of ACE2 were modeled as N-glycan core pentasaccharide, which is a minimum structure for all N-glycans. Thus, the results from Ref. [22] may not apply to specific glycan types that have distinct structures or chemical-physical properties. Here, we specifically included MAN9 or FA2 glycans in ACE2 in our simulations and have shown that they are associated with comparable glycan-RBD contact interactions (Figure 4a,b). While the contacts made with RBD may be similar for MAN9 and FA2, different chemical-physical properties of the glycans can lead to distinct binding energetics of the RBD-ACE2 complex, as will be discussed in the following section.

2.5. RBD-ACE2 Binding Energetics Depend on Glycan Type

To elucidate how the energetics of RBD-ACE2 binding may be affected by MAN9 or FA2 glycans, we applied the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach [27,28]. This technique may be applied to distinguish energetic stabilities of biomolecular states. Moreira et al. [29] have performed these calculations to provide an approximation of the relative stabilities between different Spike protein conformations. Here, we employ this end-state free-energy calculation to estimate the relative changes in binding of RBD and ACE2, with and without glycans. In the present MM-PBSA calculations, the binding energy of a ligand-receptor complex (RBD and ACE2) is approximated by changes in molecular mechanics and solvation energies. Here, the sum of the two contributions is referred to as MM-PBSA energy, which is composed of enthalpic terms and an approximation for the solvent entropy (see section “MM-PBSA calculations” of the SI Methods for details). To assess the free energy of binding (i.e., affinity), one needs to further determine changes in conformational entropies. Separately from the MM-PBSA energy calculation, the configurational entropy was estimated using a quasi-harmonic approach. We found that the relative changes in entropy are similar for the MAN9 and FA2 simulations, and therefore, did not include them in the binding energy calculations (see SI Methods and Figure S2a). Hence, the MM-PBSA energies discussed below can be used to approximate relative changes in binding energetics, allowing us to infer how each glycan type can influence RBD-ACE2 affinity.

To dissect the effect of different glycan types on RBD-ACE2 stability, we calculated MM-PBSA energies (i.e., stability) for the simulations where glycans were not included, and where MAN9 or FA2 glycans were bound to ACE2 (Figure 5a). For each simulation set, we used large numbers of representative snapshots for MM-PBSA analysis to ensure statistical convergence (see SI Methods for details). We find that the simulations with MAN9 glycans on ACE2 result in a $14.7\%$ decrease in RBD-ACE2 binding stability, relative to the simulations without glycans (Figure 5a). In contrast, the trajectories with FA2 glycans on ACE2 lead to a $9.1\%$ increase. While this MM-PBSA analysis does not report free energies, the distinct combinations of enthalpy and solvent entropy demonstrate that RBD-ACE2 affinity can be dependent on glycan species. Specifically, since the entropies associated with MAN9 and FA2 glycans are comparable (Figure S2a), our calculations suggest that FA2 glycans in ACE2 would have a stabilizing effect on the RBD-ACE2 complex, while the MAN9 type would lead to weaker binding affinity.

As shown earlier, of the six glycans that are on ACE2, only the one at Asn90 forms significant contacts with RBD, regardless of ACE2-glycan type (cf. Figure 4a,b and Figure S1). Since only glycan Asn90 forms close contacts with RBD, one could ask if the binding energy perturbations (as presented in Figure 5a) resulted solely from the effects of glycan Asn90. To answer this, we removed glycan Asn90 from both trajectories where either MAN9 or FA2 glycans were on ACE2, and repeated the MM-PBSA analysis. When glycan Asn90 is excluded (Figure S2b), the MAN9 trajectory leads to a $9.6\%$ decrease in RBD-ACE2 binding stability, relative to the simulations without glycans (contrast this to the $14.7\%$ decrease for when glycan Asn90 is present, cf. Figure 5a). For the FA2 simulations, excluding glycan Asn90 results in an increase of $3.6\%$ in stabilizing energy (versus the $9.1\%$ increase for when glycan Asn90 is included, cf. Figure 5a). This comparison demonstrates that while glycan Asn90 is a dominant contributor to the observed affinity perturbations, it is not the sole contributor out of the six ACE2 glycans. Since the remaining glycans on ACE2 lead to appreciable stability changes (Figure S2b) while not forming close contacts with RBD (Figure 4a,b and Figure S1), it suggests that glycans affect RBD-ACE2 affinity through long-range electrostatic interactions, as discussed in more detail below.

2.6. Electrostatic Effects of Different Glycan Types Lead to Distinct RBD-ACE2 Binding Energetics

To identify which energetic contributions are responsible for the distinct affinities as shown in Figure 5a, we dissected individual components of the MM-PBSA energy for the simulations without glycans, and those with MAN9 or FA2 glycans in ACE2. Relative to the non-glycan simulations, the ones with MAN9 or FA2 glycans in ACE2 each lead to a similar increase in stabilizing van der Waals energetics of roughly 30% (Figure S2c). In terms of solvation energy (summation of polar and nonpolar), the MAN9 and FA2 simulations are associated with greater stabilizing contributions by 35% and 31%, compared to the non-glycan case. Regarding electrostatic energy, the trajectories with MAN9 or FA2 glycans in ACE2 result in decrease of stabilizing contributions by 42% or 32%, relative to the value for non-glycan simulations. The summation of these energetic components constitute the MM-PBSA results shown in Figure 5a. Our energy partitioning demonstrates that electrostatic contributions, specific to each glycan type, are mainly responsible for the observed variation in binding affinity of RBD-ACE2.

As just discussed above, electrostatic effects are the dominant factor for the observed difference in MM-PBSA properties associated with MAN9 and FA2 glycans. Distinct electrostatic features can stem from the charged sialic acids in FA2 and the characteristic solubility of each glycan type. It has been suggested that the ACE2 binding surface is overall negatively charged, while the corresponding RBD interface is positive [30]. Here, we further determined that ACE2 and RBD surface-potential distributions are also mainly negative and positive, respectively, for regions well beyond the binding interface (Figure 5b). Hence, interaction of the positive RBD with the negative sialated tips of the FA2 glycans on ACE2 would lead to stabilizing MM-PBSA contributions, as compared to MAN9 glycans that do not have these sialated tips. Since the RBD-ACE2 contact map is independent of the presence of glycans, the observed MM-PBSA effects should be mainly the result of long-range electrostatic interactions of the sugars that are flexible and sample wide regions around the RBD. These findings are in agreement with numerous studies that have identified sialic-acid interaction sites away from the binding interface, which may be critical for virus-host binding [31,32].

2.7. The Role of the RBD-Asn343 Glycan in Contributing to Infectivity

There is growing evidence that the glycan at Asn343 in RBD plays a critical role during viral infection [10,33]. In the present study, we modeled the RBD glycan at Asn343 as a complex FA2 glycoform, in accordance with previous experimental studies [8,34]. In our simulations, the Asn343 glycan of RBD forms contacts with the Asn53 and Asn322 glycans of ACE2 in $20\%$ and $40\%$ of the sampled configurations (regardless of ACE2 glycan type, see Figure S3a). Interestingly, the RBD-Asn343 glycan forms more frequent contacts only with the two aforementioned glycans of ACE2, but does not form significant short-range interactions with any protein residues of ACE2 (less than $10\%$ contact probabilities, see Figure S3a).

To elucidate the energetic contribution of the RBD-Asn343 glycan to virus-host affinity, we removed this glycan and repeated the MM-PBSA analysis (Figure S3b). For the simulations with MAN9 glycans in ACE2, deleting the RBD-Asn343 glycan leads to a slight increase in RBD–ACE2 stability of $10.6\%$ (i.e., compare pink with magenta bar in Figure S3b). To rationalize this result, we note that since the RBD glycan is of FA2 type, its deletion would remove repulsive electrostatics between the negatively charged FA2-sialated tips and the overall negative ACE2 surface (cf. Figure 5b), thereby increasing stability. Similarly, for the simulations with FA2 glycans in ACE2, we observe enhanced RBD–ACE2 stability (by $6.8\%$) when the RBD-Asn343 glycan is excluded (i.e., compare cyan with blue bar in Figure S3b). In this case, RBD-glycan deletion would eliminate repulsive interactions of the negative RBD-FA2 tips with (a) the negative FA2 glycan tips of ACE2, and (b) the overall negative surface of ACE2. Relative to the simulations with MAN9 in ACE2, the ones with FA2 in ACE2 exhibit greater stability after removal of RBD-Asn343 glycan (i.e., cyan vs. pink bar in Figure S3b), which may be the consequence of stabilizing electrostatics between the FA2 glycans of ACE2 and the positive RBD surface, as discussed above. In summary, the MM-PBSA calculations suggest that the FA2 glycan at RBD Asn343 has a destabilizing effect on RBD-ACE2 binding, regardless of which glycan type is present in ACE2.

The functional relevance of the RBD-Asn343 glycan may be reflected by the fact that this glycan is highly conserved in current GISAID SARS-CoV-2 sequences. Specifically, it is lost in only 3 out of $313,826$ sequences, according to cov.lanl.gov Spike alignment (as of 16 January 2021). SARS-CoV-2 pseudovirus essays by Li et al. [10] have shown that removal of the RBD-Asn343 glycan leads to a 20-fold decrease in infectivity. Infectivity is determined by factors such as RBD-ACE2 binding stability [3,4], and/or the relative population of the all-down-RBD Spike state and infection-capable (RBD-up) Spike states [16]. Regarding the latter, Sztain et al. [33] have demonstrated that the RBD N-glycan at Asn343 facilitates the transition toward the RBD-up Spike conformation. Since the present MM-PBSA analysis indicates that the Asn343 glycan decreases RBD-ACE2 affinity, it suggests that RBD opening, which is facilitated through Asn343 glycan (as shown by Sztain et al. [33]), could be a dominant mechanism for the infectivity differences observed by Li et al. [10]. Hence, the combined findings from these studies, together with our energetics analysis, have revealed the disparate roles of the RBD-Asn343 glycan during viral infection. That is, while this glycan may reduce RBD affinity to ACE2 host receptor, it nevertheless enhances infectivity by shifting the Spike toward the RBD-up ensemble, which facilitates binding events to the host receptor.

3. Discussion

To develop effective strategies for treating COVID-19, one requires thorough insights into the factors that determine the binding between viral RBD and human ACE2. To this end, recent structural analyses [5,6,7] and mutational experiments [3,4] have alluded to possible RBD-ACE2 interactions, as well as the role that glycans may play during this process. To extend these initial findings, one must obtain a detailed mechanistic understanding of the dynamics associated with numerous contributors to binding. For this purpose, we have performed long all-atom simulations of the RBD-ACE2 complex in explicit solvent. To dissect each contributor, we simulated separate models where different glycan types are bound or excluded (Figure 1). Interestingly, the simulations show that MAN9 or FA2 glycans have virtually no effect on the interface contacts between RBD and ACE2 (Figure 2a–c), nor the fluctuations of RBD relative to ACE2 (Figure 2d). Our statistical evaluation of contacts has identified which RBD-ACE2 interface residues are most likely to form interactions, thereby pinpointing the critical sites for binding. This analysis is based on the most extensive RBD-ACE2 simulations to date, comprising over $225\mathsf{\mu }$s, which is one to two orders of magnitude longer than previous simulations. Hence, our results represent more accurate statistics that may be used as reference for interpreting experimental measurements, or gauging the impact of specific mutations. As an example of application, we have used the simulations as guide to evaluate how the novel N501Y mutation may modulate the affinity between virus and host (Figure 3). Specifically, our analysis suggests that N501Y leads to additional stabilizing interactions (i.e., with Y41 and K353 of host ACE2, see Figure 3b) that can facilitate virus binding.

The present simulations suggest that excluded volume effects of the Asn90 glycan in receptor ACE2 may determine the binding and unbinding of viral RBD (Figure 4). Precise characterization of glycan structure and function in experiments has not been possible so far, because glycan structures are highly complex, heterogeneous, and flexible [35,36]. All currently available SARS-CoV-2 Spike and human ACE2 structures only include glycans that are partially resolved (i.e., up to the minimal stem conformation). Furthermore, these incomplete glycan structures represent only a very limited conformational ensemble, since the flexible glycans can sample extremely large conformational spaces. Therefore, to accurately describe the effect of glycans, the simulations must capture sufficiently large number of configurations. To account for this, we have generated 40 long simulations with MAN9 or FA2 glycans, where each run was initiated using a distinctive starting configuration (see Section 4 for details). The simulated trajectories have revealed that the glycan on Asn90 of receptor ACE2 forms more significant contacts with viral RBD than any other ACE2 glycan. This suggests that steric effects from Asn90 glycan may impede the unbinding of RBD, thus stabilizing the RBD–ACE2 complex. These results are consistent with recent atomic force microscopy (AFM) experiments, where separating RBD and ACE2 requires higher forces when the Asn90 glycan is present [22]. On the other hand, when ACE2 is not in complex with RBD, the significant steric presence of Asn90 glycan indicates that it would obstruct the binding site, thereby hindering the RBD from associating with ACE2. This notion agrees with mutational experiments demonstrating that removal of the Asn90 glycan leads to increased RBD binding events [23,24]. Here, by uncovering the prominent role of Asn90 glycan in forming contacts with RBD, our simulations indicate how Asn90-glycan sterics may impede the binding, as well as the unbinding of RBD. In addition to these findings, a recent AFM study has elucidated the thermodynamics and kinetics of virus-host binding, highlighting the critical role of the RBD-ACE2 interface during this process [37]. Thus, the combined perspectives from experiments and computational efforts are continually dissecting the numerous energetic factors, allowing future studies to more precisely regulate the landscape that governs this interaction.

With increasing focus on the use of soluble extracellular domains of ACE2 as decoy inhibitors, it is critical to understand how specific glycan types on ACE2 can affect the binding energetics of RBD–ACE2. Our calculations suggest that distinct electrostatic features of different glycan types can be decisive for binding affinity (Figure 5). Specifically, we have shown that stabilizing energetics arise from interaction between the overall positively charged RBD and the negative sialated tips of complex FA2 glycans on ACE2. These observations are in line with experiments that reported stabilizing binding sites for sialic acids on Spike proteins in various SARS and MERS viruses [31,38]. Sialation has also been known to improve thermal stability and solubility [39,40], as well as better recognition of “self” versus “non-self” by Siglecs [41], which are all favored features for effective immunogens. Together, our simulations, along with these supporting experiments, help establish how glycosylation variability on ACE2 can contribute to the slight discrepancies in binding affinities of SARS-CoV-2 for host receptor [1,5,7,42,43], possibly explaining the broad range of host-immune responses in the human population [44].

The present study elucidates the numerous contributors to the interplay between viral RBD and human ACE2, including protein-protein interactions, and the effect of specific glycan types on the binding of the RBD-ACE2 complex. The simulations have provided evidence that the stability of RBD-ACE2 is dependent on which glycan type is bound to the host receptor. Remarkably, ACE2 glycans can affect virus binding affinity through electrostatic effects, while not perturbing the physical contacts that are formed between virus and host. Of the numerous glycans on ACE2, the calculations have revealed that glycan Asn90 is a dominant contributor to affinity perturbations, and may play a critical role in the binding and unbinding events of RBD. Together, our results provide a theoretical framework that may be used to design more precise experiments that aim at regulating the infectivity of SARS-CoV-2.

4. Methods

To elucidate the interaction between RBD and ACE2, we performed all-atom explicit-solvent molecular dynamics (MD) simulations of the RBD-ACE2 complex (PDB ID: 6M0J [5]; Figure 1a). In all models used for simulations, the SARS-CoV-2 RBD domain is defined by residues T333–G526, and the ACE2 receptor by S19–D615. To avoid artificial charges at the protein ends, we introduced N-terminal acetylated and C-terminal N-methylamide capping groups. The ACE2 structure contains a zinc ion that is coordinated by H374, H378, E402, and one water molecule [5,45]. Zinc coordination plays a critical role in maintaining the structural integrity and stability of a protein [46]. To properly account for this in the simulations, we introduced bonded terms between the zinc ion and H374, H378, E402 in ACE2. Equilibration values for distances and angles of these bonded terms were defined by the values found in the RBD-ACE2 configuration of PDB ID: 6M0J [5]. The force constants for bonds and angles of zinc interactions were set as $k_{\mathrm{b}}=5\times {10}^5$ kJ·mol${}^{-1}·$ nm${}^{-2}$ and $k_{\theta }=1\times {10}^3$ kJ·mol${}^{-1}·$ rad${}^{-2}$, respectively. The strength of these force constants are equivalent to covalent bond interactions, as described in the CHARMM36m forcefield [47], which was applied for the present simulations. We determined these force parameters using preliminary simulations of C-Raf CRD. In those C-Raf simulations, the parameters used here led to structural fluctuations of the zinc site, whose scale at $310$K is consistent with the NMR ensemble from Ref. [48].

To dissect the role of distinctive glycan species in the interaction of RBD and ACE2, we performed additional separate sets of simulations where (a) MAN9 and FA2 glycans are attached to ACE2 and RBD, respectively (Figure 1b) and (b) FA2 glycans are on both ACE2 and RBD (Figure 1c). Generally, glycans have complex structures and are highly flexible, allowing them to sample very broad conformational spaces. Accordingly, simulations may not accurately capture the effect of glycans if these calculations depend on the choice of initial configuration. To prevent this artifact, we prepared 40 complementary initial configurations for the MAN9 and FA2 simulations (i.e., 20 for each simulation set). Each glycan initial structure was modeled based on a simulated RBD-ACE2 configuration from preliminary trajectories (see section “Choice of glycosylation and glycan modeling” of the SI Methods for details). By initiating the glycan simulations from many distinctive configurations, we capture larger conformational ensembles that can more accurately partition the contribution of each glycan type.

Simulation Details

All-atom explicit-solvent simulations were performed with the AMBER 16 software package [49]. The CHARMM36m protein [47] and CHARMM36 carbohydrate [50,51] forcefields were used, with TIP3P water model [52]. Each configuration was solvated and centered in a cubic box. The size of the cubic box was chosen to create at least 15 Å padding on each side along the largest atom-atom distance of the molecule. Each system was neutralized with an excess of $150$mM KCL. Energy minimization was performed using the steepest descent algorithm. Equilibration simulations were first carried out under the constant number-volume-temperature (NVT) ensemble for $2$ ns, and then under the constant number-pressure-temperature (NPT) ensemble for $10$ ns. During both equilibration stages, harmonic position restraints were imposed on all non-hydrogen atoms of the molecule. Constant temperature was maintained at $310$K using velocity Langevin dynamics [53], with a relaxation time of $1$ps. Constant isotropic pressure of $1$bar was achieved by employing the Berendsen barostat [54], with a relaxation time of $4$ps and compressibility of $4.5\times {10}^{-5}$bar${}^{-1}$. Covalent bond lengths were constrained with the SHAKE algorithm [55]. Van der Waals interactions were evaluated using a cutoff where forces smoothly decay to zero between $1.0$–$1.2$ nm. Coulomb interactions were computed using the particle-mesh Ewald (PME) method [56], with Fourier grid spacing of $0.08$–$0.10$ nm and fourth order interpolation. Unrestrained production simulations were performed in the NPT ensemble, with an integration time step of $4$ fs, which was enabled through hydrogen mass repartitioning [57].

For the non-glycan RBD-ACE2 model (Figure 1a), 50 simulations were performed, with a total simulated time of $165\mathsf{\mu }$s (i.e., each trajectory includes roughly $3.3\mathsf{\mu }$s). For the complex with MAN9 glycans in ACE2 (Figure 1b), 20 simulations were performed for an aggregated time of over $30\mathsf{\mu }$s (each replica of this set is about $1.5\mathsf{\mu }$s long). Finally, the simulation set for the model with FA2 glycans in ACE2 (Figure 1c) contains 20 simulations that exceed $30\mathsf{\mu }$s of accumulated time (each trajectory of this set is roughly $1.5\mathsf{\mu }$s long).

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/v13050927/s1, Figure S1: Total probability of glycan-protein contact formation for each glycan, Figure S2: Analysis of RBD-ACE2 binding energy, Figure S3: Interactions and energetic effects associated with RBD-Asn343 glycan, Table S1: Quantitative evaluation of contacts previously implicated by experimental structures.

Author Contributions

K.N., S.C., R.A.M., B.K., and S.G. designed research. K.N. and S.C. performed simulations. S.C., K.N., and R.A.M. prepared structural models. S.C., K.N., R.A.M., B.K., and S.G. analyzed data. K.N. and S.C. created figures and wrote the paper. K.N., S.C., R.A.M., B.K., and S.G. edited the paper. B.K. and S.G. obtained funding. All authors have read and agreed to the published version of the manuscript.

Funding

K.N. and S.G. are partially supported by the DOE Office of Science through the National Virtual Biotechnology Laboratory, a consortium of DOE national laboratories focused on response to COVID-19, with funding provided by the Coronavirus CARES Act. S.C. is supported by the Center of Nonlinear Studies Postdoctoral Program. R.A.M. was supported by a Los Alamos National Laboratory (LANL) Director’s Postdoctoral Fellowship. B.K. and S.G. are partially supported by LANL LDRD project 20200706ER.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data will be available upon reasonable request.

Acknowledgments

This research used computational resources provided by the LANL Institutional Computing Program.

Conflicts of Interest

The authors declare no competing interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Simulations of the RBD-ACE2 complex in the absence and presence of different glycan types. (a) Structural representation of human receptor ACE2 (silver) bound to RBD of SARS-CoV-2 (pink). The yellow bead in ACE2 depicts the zinc ion that is coordinated by His374, His378, Glu402, and one water molecule. (b) Same protein complex as in panel (a) except that six MAN9 glycans (magenta) are bound to ACE2 (at Asn53, Asn90, Asn103, Asn322, Asn432, Asn546) and one FA2 glycan (blue) is bound to RBD (at Asn343). (c) Same complex as in (b) except that the six ACE2 glycans are FA2. Panels (a–c) represent three systems for which separate sets of simulations were performed. All molecular graphics were created using VMD [15]. (d) Schematic representation of N-glycans. Oligomannose (MAN9) and complex (FA2) glycans are depicted using the Symbol Nomenclature for Glycans (SNFG).

Figure 2. Contact signature of ACE2 and RBD is robust to the presence and types of glycans. (a) Probability of contact formation between ACE2 and RBD residues. This contact map was calculated from the simulations where glycans are not present (cf. Figure 1a; see Table S1 for list of contacts). To elucidate the role of glycans, we repeated this analysis for simulations that either include (b) MAN9 (cf. Figure 1b) or (c) FA2 glycans (cf. Figure 1c) on ACE2. The contact probabilities for both glycan systems (b,c) are virtually the same as panel (a), demonstrating that the RBD-ACE2 contact signature is robust. (d) Distributions of Cα-rmsd of RBD after least square fitting of ACE2, calculated for non-glycan and glycan simulations. (e) Structural representation of the RBD-ACE2 binding interface highlighting the residues that form persistent short-range contacts (cyan in ACE2 and ice blue in RBD). Persistent short-range contacts are those that form within a distance cutoff of 4 Å, with at least 60% probability (see Table S1).

Figure 3. N501Y mutation introduces additional stabilizing interactions between RBD and ACE2. (a) Dashed circle indicates the location of the N501Y mutation site in RBD, which is at the binding interface with ACE2. (b) Zoomed-in view of N501Y location (i.e., dashed circle region of panel (a)). The wild type residue N501 of RBD is shown in green, and the mutation N501Y in red. The simulations show that N501 forms persistent short-range contacts with Y41 and K353 (cyan) of ACE2 (cf. Figure 2a). Based on this, we assess how the N501Y mutation may alter the interactions with these ACE2 residues. As shown by dashed lines and their labels, N501Y would introduce additional stabilizing interactions with Y41 and K353 of ACE2, which will increase the binding affinity between RBD and ACE2.

Figure 4. Only the glycan on Asn90 of ACE2 forms significant contacts with RBD. (a) Probability of contact formation between RBD and MAN9 glycans of ACE2. (b) Contact probability between RBD and FA2 glycans of ACE2. Both (a,b) show that only the contacts involving Asn90 glycan occur with higher probabilities (see region indicated by arrow). RBD residues that form these protein-glycan interactions are roughly between 403–417. This RBD region does not form contacts with ACE2 in the non-glycan or glycan-included simulations (cf. Figure 2a–c). This suggests that the glycans do not compete with protein-protein interactions at the RBD-ACE2 interface. (c) Structural description of the ACE2 Asn90 glycan (magenta) and RDB residues 403–417 (green), forming glycan-protein interactions.

Figure 5. RBD-ACE2 binding affinity is glycan dependent.(a) To evaluate the binding energy between RBD-ACE2, the MM-PBSA approach was applied. Binding energy was calculated for simulations without glycans (white bar), and with MAN9 glycans (magenta bar) or FA2 glycans (blue bar) on ACE2. Simulations with MAN9 glycans on ACE2 are associated with a 14.7% decrease in stability, relative to the non-glycan simulations. In contrast, simulations with FA2 glycans on ACE2 result in a 9.1% increase in stability. (b) Electrostatic surface potential calculated for ACE2 and RBD. Two complementary views show that the ACE2 surface is overall negatively charged, while the surface of RBD is overall positive.