key: cord-0697159-oogzhetv
authors: Wang, Xiaocong; Bie, Lihua; Gao, Jun
title: Structural Insights into the Cofactor Role of Heparin/Heparan Sulfate in Binding between the SARS-CoV-2 Spike Protein and Host Angiotensin-Converting Enzyme II
date: 2022-01-21
journal: J Chem Inf Model
DOI: 10.1021/acs.jcim.1c01484
sha: a7bdb547d3de905f84c571996fd863f65099c19d
doc_id: 697159
cord_uid: oogzhetv

[Image: see text] The viral entry process of the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) requires heparin and heparan sulfates from the cell surface, functioning as a cofactor for human angiotensin-converting enzyme 2 (ACE2) for recognizing the receptor-binding domain (RBD) of the spike (S) protein on the surface of the virion. In the present study, the binding poses of an oligosaccharide with four repeating units of GlcNS6S-IdoA2S (octa) predicted by Vina-Carb in the RBD binding site were employed in molecular dynamics (MD) simulations to provide atomic details for studying the cofactor mechanism. The molecular model in the MD simulations reproduced the length- and sequence-dependent behavior observed from the microarray experiments and revealed an important planar U-turn shape for HP/HS binding to RBD. The model for octa with this shape in the ACE2–RBD complex enhanced the interactions in the binding interface. The comparisons with the ACE2–RBD complex suggested that the presence of octa in the RBD binding site blocked the movements in a loop region at the distal end of the RBD binding interface and promoted the contacts of this loop region with the ACE2 N-terminus helix. This study shed light on the atomic and dynamic details for HP/HS interacting with RBD and provided insights into their cofactor role in the ACE2–RBD interactions.

The ongoing pandemic of novel coronavirus 2019 (COVID- 19) has caused serious public health crises and enormous economic damages around the globe. 1 Significant efforts have been devoted to studying the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a new virus belonging to the β coronavirus family that has caused this pandemic. 2,3 SARS-CoV-2 is a single-stranded RNA-envelop virus, 4 which utilizes the receptor-binding domain (RBD) of its highly glycosylated trimeric spike (S) protein on the surface of the virion for the recruitment to the host receptor, human angiotensin-converting enzyme II (ACE2), to facilitate its entry via the fusion between the viral membrane and the host target membrane. 5−7 Due to the pivotal role of the S protein in viral recognition and entry, as well as their locations on the viral surface, S proteins have been employed as immunogens for generating neutralizing vaccines and targets for developing therapeutic drug molecules. 8, 9 Similarly, ACE2 has also been studied as a potential target for interfering with the ACE2targeting coronavirus infections. 10, 11 Thus, a detailed understanding of the binding interactions between SARS-CoV-2 and ACE2 is essential for elucidating the binding and entry of this new virus and its enhanced infectivity, 6,12−15 as well as rational design for effective therapeutics. 16 Previous studies have revealed that after the S protein undergoes a hinge-like conformational change that exposes the RBD, 5−7,17 two short β strands with connecting loops and a flexible loop at the distal end in RBD of the SARS-CoV-2 S protein were recognized by the N-terminal helix of ACE2. 12, 18 This recognition and binding processes engaged a network of hydrophilic interactions formed by multiple hydrogen bonding interactions and salt bridges. 7, 12 Further studies identified significant involvements of glycans in this process 19−21 after the structural analysis showed that both SARS-CoV-2 S protein and ACE2 are heavily glycosylated. For example, the glycans at N165 and N234 in RBD of the SARS-CoV-2 S protein play roles in the conformational plasticity of RBD, where they stabilize the open state of RBD; 22, 23 the glycans at N331 and N343 are shown to be critical for immune recognition. 23, 24 In the meantime, the glycans at N53 in ACE2 facilitate the stabilization of the ACE2−RBD binding surface; 21 the glycans at N90 have been theoretically shown to enhance the ACE2− RBD binding. 25−27 The exterior glycans that are not conjugated with the S protein or ACE2 have also been shown to modulate their binding interactions. A recent seminal study by Clausen et al. 28 revealed that heparin (HP)/heparan sulfate (HS) functions as a necessary cofactor to facilitate the recruitment of RBD in the S protein of SARS-CoV-2 to ACE2 by interacting with residues that are adjacent to the ACE2-binding site. As the binding sites for ACE2 and HP/HS in the RBD of the SARS-CoV-2 S protein are adjacent, 28, 29 understanding the mechanism of RBD−HP/HS interactions and their relationship with the ACE2−RBD binding interactions could exceedingly advance our knowledge of the SARS-CoV-2 infections and provide insights into a promising way for mitigating the pandemic. 30−32 HP/HS is a major class of unbranched and highly negatively charged glycocalyx polysaccharide, namely glycosaminoglycans (GAGs), which are composed of sulfated repeating [Dglucuronic acid (GlcA) or L-iduronic acid (IdoA)] and Dglucosamine (GlcN) disaccharide units and regulates a wide range of activities, including cell dynamics, inflammation, signaling pathways, etc., through their interactions with different proteins. 33−35 Liu et al. 36 performed microarray binding experiments with an extensive HP/HS oligosaccharide library and revealed that the bindings between the RBD in the SARS-CoV-2 S protein and HP/HS depend on the length and sequence of the oligosaccharides, in which HP with four GlcNS6S-IdoA2S repeating units showed stronger bindings than those with three repeating units and removing the sulfate moieties significantly reduced their binding intensities. However, the lack of the crystallographic structures of the RBD−HP/HS complexes obstructs the understanding of their interactions at the atomic level. Given that HP/HS is a linear unbranched polysaccharide that could undergo large internal motions to adopt multiple energy-favored conformations, it is also essential to study their interactions with accurate predictions of their 3D shapes and comprehensive representations for their dynamical motions in binding sites 37−39 In the present study, a model with four repeating unit of GlcNS6S-IdoA2S was built and docked to the putative binding site in the RBD of the SARS-CoV-2 S protein. Molecular dynamic (MD) simulations were performed for the complex models that were built with proper force fields. Careful validations of the RBD complex model in MD simulations were carried out by reproducing their length-and sequencedependent behavior observed in the microarray experimental results. The oligosaccharide in the validated binding pose presented a U-shape, which permits multiple residues at the turn forming stable interactions with previously identified key arginine residues simultaneously, as well as those at the reducing and nonreducing ends. This U-shape is speculated to be a general structural feature for HP/HS in the RBD binding sites, as the energy-favored glycosidic linkage conformations in HP/HS allows five consecutive monosaccharides in the same plane, which parallels with the side chains of arginine residues in the binding site and maximizes the possibilities of forming stable interactions.

The validated binding pose for the oligosaccharide was merged into the model of the ACE2−RBD complex to study the mechanism of HP promoting the ACE2−RBD binding interactions. The essential stable intermolecular interactions in the concave ACE2−RBD binding interface, observed in the crystal structure of the ACE2−RBD complex, were reproduced by MD simulations with and without bound HP. But, the presence of HP would block the motion of a distal end loop region that resides between the binding sites of ACE2 and HP. This blockage reduced the freedom of the loop region and promoted its interactions with residues at the N-terminus of ACE2, which agreed with the observation that HP/HS functions as a cofactor for the recruitment of SARS-CoV-2 by the ACE2 at the host cell membranes.

The present theoretical study provided structural insights into the poses of HP in the RBD binding site and the detailed atomic mechanism of the HP promoting the ACE2−RBD interactions.

■ METHODS Structure Preparation. The initial coordinates for the S protein in SARS-CoV-2 were obtained from the Protein Data Bank (PDB entry code: 6VYB 5 ). The RBD domain was extracted from its chain B (residues 334−527). Missing residues, 455−461, 467−490, and 516−521 were built in the SWISS-MODEL webserver 40 with the template from the RBD in 6VW1. 12 The initial coordinates for the ACE2−RBD complex were also obtained from the Protein Data Bank (PDB entry code: 6VW1 12 ). For consistency, the index of residues in this study followed those in 6VW1. The initial structure for octa was built with the GAG builder on www.glycam.org. The ring conformations of GlcNS6S and IdoA2S were set to 4 C 1 and 1 C 4 prior to docking, respectively. Docking was performed with Vina-Carb 41 with the center of a docking region of 60 × 60 × 60 Å 3 at the geometry center for the side chains of R346, R355, and R466, and the top 6 ranked poses were extracted. Default Vina-Carb options (chi_coeff = 1 and chi_cutoff = 2) were used. Force field parameters for carbohydrate molecules were taken from GLYCAM06 (version j), 42 and those for proteins were taken from AMBER18 (ff99sb). 43, 44 Counterions (Na+) were added to neutralize each protein complex using the tLEaP module 44 before they were solvated in a truncated octahedral box (8 Å buffer with TIP3P water model).

Simulation Setup. A two-step energy minimization of the solvated complexes was performed under the NVT conditions: first, only the water molecules and counterions were subjected to energy minimization (500 steps SD followed by 24 500 steps CG) and the atoms in solute were restrained (100 kcal/ mol·Å 2 ); in the second step, the energy minimization circle was repeated with the restraints applied only to the Cα atoms on the protein backbone and ring atoms in carbohydrate molecules. Then, the solvated complexes were heated to 300 K within 50 ps under the NVT conditions with a restraint (10 kcal/mol·Å 2 ) applied only to the Cα atoms in the protein.

Prior to data collection, the solvated complexes were equilibrated at 300 K under NPT conditions with a Berendsen thermostat 45 for 10 ns, during which all restraints were removed. MD simulations for the data collection of each solvated complex were performed for 500 ns with the GPU implementation 46 of PMEMD using the AMBER18 software package. 44 Covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm 47 in all MD simulations, which allows a simulation time step of 2 fs. A nonbonded cutoff of 8 Å was applied to van der Waals interaction energy calculations, and the particle mesh Ewald approximation was applied to the long-range electrostatic interaction energy calculations. Standard 1−4 nonbonded Journal of Chemical Information and Modeling pubs.acs.org/jcim Article scaling factors for proteins (2.0/1.2) and carbohydrate molecules (1.0/1.0) were employed. 42 Representative Conformation of Octa and Model Generations for hexa-1, hexa-2, and ACE2−RBD-octa. The conformation of octa that was most similar to its average shape in the RBD complex acquired from the MD simulations was extracted and presented as its representative conformation. The GlcNS6S-IdoA2S segment was removed from the reducing and nonreducing ends of octa in the RBD complex to generate the molecular models for RBD-hexa-1 and RBDhexa-2 complexes, respectively. The molecular model for ACE2−RBD-octa was created by adding octa to the model for the ACE2−RBD complex after superimposing the Cα atoms of the stable secondary structures (residues 353−359, 375−381, 393−403, 430−438, 507−517) in the RBD of the ACE2− RBD complex and the representative conformation of the RBD-octa complex.

Interaction Energy Calculations. Interaction energies in the complexes were calculated using the molecular mechanicsgeneralized Born solvent-accessible surface area (MM-GBSA) 48, 49 under the single trajectory methodology with the MMPBSA.py.MPI module in AMBER on 10 000 snapshots extracted evenly from each MD simulation. The GB 1 OBC model (igb = 2) 50 and internal dielectric constant (ε int ) of 4.0 were applied in all MM-GBSA calculations. 51, 52 ■ RESULTS AND DISCUSSION Predicting HP Binding Poses in RBD. HP is an unbranched oligosaccharide that could undergo large conformational changes by altering the conformations of glycosidic linkages to adopt multiple energy-favored conformations. Thus, predicting the pose of HP in the RBD binding site requires not only approximations of binding interaction contributions from individual monosaccharides but also estimations of the penalties from different glycosidic conformations that determine the overall shape of the binding pose. Thus, Vina-Carb 41 was employed for this purpose, as it was designed to include the glycosidic linkage conformations in its scoring function and has been proven to be effective in predicting the binding poses for flexible oligosaccharides in protein complexes. 53, 54 The center for the docking region was selected as the geometric center of the atoms in R346, R355, and R466, which were shown to be critical for interacting with HP in previous studies 28 and were located in a cavity next to the ACE2-binding site.

A recent microarray study on the RBD of the SARS-CoV-2 S protein using an HP/HS library suggested the lengthdependent behavior of the bindings between HP/HS oligomers and RBD. 36 An octasaccharide with four repeating units of GlcNS6S-IdoA2S displayed the strongest binding affinity, followed by a hexasaccharide with three repeating units. Further reducing the length of the oligosaccharide almost abolished the binding interactions. According to this finding, a molecular model for this octasaccharide (octa) was built and docked in the RBD binding site with Vina-Carb. The top six poses with the predicted binding energies more negative than −4.0 kcal/mol were extracted to form RBD-octa complexes for further MD simulation studies (Table S1 ). All six selected poses of octa (Figures 1 and S1) have shown to interact with at least two of the three key residues in its binding site. Obviously, octa in poses 2, 3, and 6 shared a similar Ushape in the binding site, which allows the residues in the turn form interactions with R355 and R466 simultaneously, and those at the reducing and nonreducing ends to interact with R346. Although the docking region was set to cover the entire binding site, R346, R355, and R466 appeared to be essential for the predicted octa poses, as all 20 poses showed close proximity to these residues. It is worth noting that these three selected poses were not shown to interact with K444, a key residue in RBD for interacting with HP reported in the pioneering experimental studies. Longer oligosaccharides with this U-shape are highly likely to interact with K444 through the terminal residues. The U-shape was observed in the other three poses, as well as 11 out of the remaining 14 docked poses ( Figure S2 and Table S1 ), yet their orientations and locations differed. In the MD simulations of RBD-octa complexes built with octa in poses 2, 3, and 6, the oligosaccharide displayed positional stabilities. The total RMSDs for the ring atoms in octa were less than 12 Å (Figure 1) . However, octa failed to maintain stable binding interactions in the RBD complexes built with octa in the other three poses, where the RMSDs reached up to 30 Å in the MD simulations ( Figure S1 ).

Having predicted a similar pose for octa in the RBD binding site from three different Vina-Carb results that could maintain stable interactions with key residues during MD simulations provided a solid starting geometry for further studies. Thus, the representative conformation of RBD-octa from these three MD simulations was extracted for further validating this novel binding pose and understanding the length-and sequencedependent behavior of HP binding to RBD. All MD simulations associated with this representative conformation of octa in the following study, including the RBD complexes with hexa and tetrasaccharides, and ACE2−RBD-octa complexes were determined with three replicas for establishing the reproducibility of the observed phenomena in simulations.

Modeling the Length-and Sequence-Dependent Behavior. The predicted pose for octa in the RBD binding site needs to replicate the length-and sequence-dependent behavior observed in the microarray results using the HP/HS library before it can be used to provide insights into the RBD− HS interactions and its cofactor function for the recruitment of SARS-CoV-2 S protein by ACE2. Therefore, models for RBD complexes with three repeating units of GlcNS6S-IdoA2S, RBD-hexa-1, and RBD-hexa-2, were generated by removing the GlcNS6S-IdoA2S unit at the reducing or the nonreducing end of octa, respectively, from the representative conformation of the RBD-octa complex, and subjected to MD simulations. For a fair comparison, an MD simulation starting from the representative conformation of RBD-octa was also performed. The sequence-dependent behavior could be verified by careful comparisons of the binding contributions from different moieties in the oligosaccharide obtained by hydrogen bond analyses and MM-GBSA energy calculations to the microarray results using the HP/HS library.

All three oligosaccharides, octa, hexa-1, and hexa-2, maintained stable bindings to RBD during their MD simulations. The positional RMSDs for the ring atoms reference to the starting geometry for MD simulations in all of these oligosaccharides in their three replicas of the MD simulations were generally less than 10 Å (Figures S3−S5) . For the RBD-octa complex, the low RMSD values suggested that 55 The protein solvent-accessible surface is shown in gray and those for the key residues in RBD interacting with HP/HS are shown in blue and labeled.

Journal of Chemical Information and Modeling pubs.acs.org/jcim Article the binding pose originally predicted by Vina-Carb was maintained in the new MD simulation. For the RBD-hexa-1 and RBD-hexa-2 complexes, the low RMSD values indicated that hexa-1 and hexa-2 in MD simulations maintained the same predicted binding pose as octa, since their starting geometries were generated by removing a disaccharide from the reducing or nonreducing ends, respectively. As seen in Figure 2 , the representative conformation of hexa-1 in the RBD complex from the MD simulation displayed a similar shape as the six residues in the representative conformation for octa from the nonreducing end, while that for hexa-2 matched the six residues from the reducing end. The interaction similarities among the three oligosaccharides with their counterparts in RBD corresponded to their spatial similarities in the binding site. The stable intermolecular hydrogen bond interactions formed by hexa-1 and hexa-2 were similar to those formed by the corresponding parts in octa (Tables 1 and S2−S4). The middle sections of hexa-1 and octa both interacted with R355 and R466, while their nonreducing ends interacted with R346.

In the meantime, the nonreducing end of hexa-2, which matched the middle section of octa, interacted with R355 and R466, while its reducing end interacted with R346. Thus, the spatial and interaction similarities of these oligosaccharides in the binding site exemplified the robustness of the binding pose originally predicted by Vina-Carb and provided further confidence. The microarray results for RBD using the HP/HS library showed that the oligosaccharide with four repeating units of GlcNS6S-IdoA2S was 15% stronger in binding than that with three repeating units. As seen in Table 2 , the total MM-GBSA energies for octa were more negative than those for hexa-1 and hexa-2 (6 and 9%, respectively). The trend in the binding energy difference suggested that the sequence-dependent binding behavior was reproduced in the models of RBD complexes. The hydrogen bond analysis showed that 2-Osulfate of IdoA2S at the nonreducing end of octa, hexa-1, and hexa-2 formed stable interactions with the residues in RBD, including R346, which was previously identified as essential (Table 1) . Losing these interactions would lead to reductions of the binding strengths for the oligosaccharide, which agreed to the over 60% binding strength reduction for the hexasaccharide in the microarray experiments after replacing the nonreducing end IdoA2S with a GlcA residue. 36 The stable hydrogen bond interactions formed by the 2-O-sulfate moieties of other IdoA2S residues in hexa-1 and hexa-2 also agreed with the elimination of binding interactions of the hexasaccharide after further replacements of IdoA2S with GlcA residues. Similarly, the hydrogen bond interactions between 6-O-sulfate in GlcNS6S residues and RBD residues in the binding site corresponded to weaker binding strength for hexasaccharide after replacing GlcNS6S with GlcNS residues in the microarray experiments. Energetically, the decomposition of MM-GBSA energies showed that 2-O-sulfates in IdoA2S and 6-O-sulfate in GlcNS6S residues all contributed to the binding interactions (Table 2) , confirming their involvements in the RBD−HP interactions. Therefore, the binding pose predicted by Vina-Carb for octa reproduced the length-and sequence-dependent behavior of HP/HS interacting with RBD.

Hypothetical General Shape of HP in RBD. In addition to the theoretical reproductions of the length-and sequencedependent interaction behavior of RBD and HS, understanding the U-shape of octa in the RBD binding site could provide further insights into the RBD−HS interactions and advanced knowledge of rational drug designs.

It is of interest to determine the energetic preferences toward the planar U-turn shape in the oligosaccharide composed of repeating units of GlcNS6S-IdoA2S. First, an oligosaccharide model with 8 repeating units of GlcNS6S-IdoA2S, whose glycosidic linkages were in the energy-favored ϕ/ψ conformations suggested in the CHI energy study, 41, 56 appeared to be a helical structure ( Figure 3A,B) . In the helices, five consecutive residues in the oligosaccharide could generate a planar U-turn geometry. As shown in Figures S6−S8 , the glycosidic linkages in the oligosaccharide appeared to maintain the initial conformations suggested in the CHI energy study, 41 suggesting that such a helical structure is energetically stable in explicit solvents. Finally, the glycosidic linkage conformations were altered manually to make the conformation of this oligosaccharide linear ( Figures S9−S11) , which was later subjected to MD simulations. Interestingly, the oligosaccharide spontaneously formed a helical conformation, although not as ordered as seen in the previous MD simulations. Therefore, theoretical models of an oligosaccharide with repeating units of GlcNS6S-IdoA2S confirmed that a helical conformation could be energy-favored in solution and potentially facilitate its binding to the RBD. In the planar U-turn shape, the sulfate moieties extended vertically to the plane ( Figure 3B ). Both the 

GlcNS6S

No stable interactions observed.

Journal of Chemical Information and Modeling pubs.acs.org/jcim Article planar shape of the oligosaccharide and extended exo-cyclic moieties favor the formation of intermolecular hydrogen bonds and salt bridges between HP/HS and the receptor. The planar shape was observed between 2nd and 4th GlcNS6S of octa in the representative conformation of the RBD-octa complex ( Figure 3C ). The sulfate moieties also appeared to be vertical to the oligosaccharide plane and formed multiple intermolecular interactions with residues in RBD simultaneously (Tables S2−S4) . It is worth noting that octa in the binding site of RBD was not in the helical form as shown in the ideal model in Figure 3A ,B. The glycosidic linkage distributions for octa during MD simulations showed that the terminal residues in octa could adopt multiple stable positions, as their glycosidic linkages adopted multiple conformations. On the other hand, the glycosidic linkages for most of the monosaccharides in the same plane adopted only the energy-favored conformations. As the shapes of hexa-1 and hexa-2 related to the corresponding sections in octa, both of their poses contained a planar U-turn. The absence of a repeating unit from each end of octa caused only a slight reduction in the binding affinity instead of a significant reduction or elimination. These observations of oligosaccharide conformations and changes of binding affinities suggested that this planar U-turn in five consecutive residues is important for the binding of HP/ HS to RBD. Furthermore, the microarray study showed that the binding strength of a tetrasaccharide with only two repeating units of GlcNS6S-IdoA2S was under detection limits. For comparisons, three tetrasaccharide models were generated by removing two repeating units from the reducing end of octa (tetra-1) or two repeating units from the nonreducing end (tetra-2), or one repeating unit from both ends (tetra-3). MD simulations of their complexes with the RBD all displayed Table 2 . Per-Residue MM-GBSA Interaction Energies (kcal/mol) of octa, hexa-1, and hexa-2 in the Binding Site of RBD Residue not present. Figure S12 ) compared to those observed in the RBD-octa complexes. Both the experimental and modeling studies suggest that two repeating units of GlcNS6S-IdoA2S with an incomplete planar U-turn would significantly harm the binding strength. Furthermore, the binding strength for a tetrasaccharide with only two repeating units of GlcNS6S-IdoA2S with RBD was under the detection limit in the microarray study, suggesting that an incomplete planar U-turn would significantly harm the binding strength.

The U-turn shape was also observed in other studies, both free form in the solution phase and the protein complex 57, 58 ( Figure S13 ). Therefore, one can speculate that this planar Ushape with five residues may be a general structural form for the HP/HS binding in the RBD binding site. Collectively, the model of the octasaccharide was validated for representing the dynamic and energetic features of the RBD−HS complex through the reproduction of the lengthand sequence-dependent binding behavior by MD simulations starting from the predicted pose in the RBD binding site. It would be further employed in generating the molecular model for the ACE2−RBD−HS complex, as the binding sites for HS and ACE2 in RBD are separated.

Mechanism of the Cofactor Role for HP in ACE2−RBD Binding Interactions. Thus, having confirmed the validity of the molecular model for the RBD-octa complex, the representative conformation for octa in the RBD binding site was utilized to study the cofactor mechanism of HP/HS for the recognition of the SARS-CoV-2 S protein by ACE2. Prior to studying the mechanism, the molecular model of ACE2−RBDocta was examined by its structural stabilities and ability to replicate the experimentally observed key interactions at the ACE2−RBD binding interface in MD simulations.

The overall RMSDs for Cα atoms in ACE2 and RBD were less than 10 Å, suggesting that proteins in the ACE2−RBDocta complex were structurally stable during the MD simulation. In addition, similar to the RBD-octa complex, octa in the ACE2−RBD-octa complex also displayed low structural variations along the course of MD simulations ( Figure S14 ). The RBD in the SARS-CoV-2 S protein binds to ACE2 through a network of hydrophilic interactions of multiple hydrogen bond interactions and salt bridges 7 in their binding interface that is formed between two short β strands with connecting loops in RBD and an N-terminal helix in ACE2. Two hotspots in the binding interface of the ACE2− RBD complex, hotspots 31 and 353, centered with K31 and K353 on the ACE2-binding ridge were previously identified by crystallographic measurements. 7 The hydrogen bond interactions between K31 and Q493 in RBD in the crystal structure of the ACE2−RBD complex 12 were observed in the MD simulation of the ACE2−RBD-octa complex (Table 3) . Although the key interactions between K353 and the backbone oxygen atom of G496 in the crystal structure were not found to be stable in the MD simulation, K353 formed stable hydrogen bond interactions with Y495 and Q498 that are adjacent to G496 (Table 3 and Figure 4 ). In addition to the essential interactions formed by the two key hotspot residues, more hydrogen bond interactions were observed within the binding interface (Table 3 and Figure 4 ). Both Y41 and D355 in ACE2 formed stable interactions with the side chain of T500 in RBD in the MD simulation; in the crystal structure of the ACE2− RBD complex, Y41 and D355 were found to be within 4 Å of the side chain of T500. Similarly, stable interactions were seen between the side chain of N330 in ACE2 and the backbone of T500 in RBD, in which the distances were less than 4 Å in the crystal structure. Additional hydrogen bond interactions were observed between Q24 in the ACE2-binding ridge and N481 in the loop region of the RBD binding interface, which could potentially enhance the contacts between two regions. In the crystal structure, although the loop region presented a different conformation that separated two residues, the interactions between these two regions were still observed as F486 formed hydrophobic contacts with the ACE2 N-terminus.

For comparison purposes, the model of the ACE2−RBD complex (without octa) was also built and subjected to MD simulations. Similar to the ACE2−RBD-octa complex, the overall RMSDs for Cα atoms in the ACE2−RBD complex were less than 10 Å ( Figure S14) , showing similar structural stability to those observed with the presence of octa in the RBD binding site. Furthermore, the intermolecular hydrogen bond interactions in the ACE2−RBD binding interface without octa in the complex were also similar to those previously observed in the ACE2−RBD-octa complex (Tables 3 and S5 and Figure 4) . The key interaction pairs involving hotspot residues, K31−Q493 and K353−Q498, were observed. Instead of the backbone of Y495, K353 formed interactions with the backbone of G496, as seen in the crystal structure. In addition to the hotspot residues, D355 was also observed to interact with T500. The pairing of N330−T500 at the distal end of the binding interface was not observed as stable interactions in the MD simulation of the ACE2−RBD complex, yet, the interactions between Q42 and Q498 were seen as stable. Collectively, the models for the ACE2−RBD complex with and without the presence of octa in the binding site of RBD maintained structural stabilities and reproduced essential intermolecular hydrogen bond interactions observed in the ACE2−RBD binding interface. Therefore, the dynamic and thermodynamic differences of these two complexes were employed for understanding the cofactor mechanism of HP/ HS in the recognition of RBD in the SARS-CoV-2 S protein by ACE2.

The stable hydrogen bond pair of Q24−N481 was observed when octa was bound to RBD, but not when octa was absent. The loop region at the distal end of the binding interface that contained the N756 residue showed distinct conformations in the representative conformations for both complexes. Without the presence of octa, this loop region was seen to be further away from the ACE2 N-terminus helix in the binding interface compared to that in the complex with the presence of octa ( Figure 4) . The positions of this loop region in the static representative conformation in the ACE2−RBD complex extracted from the MD simulation appeared to be different from that observed in the crystal structure. 7 Thus, to confirm the dynamic behavior of the loop region beyond the static representations, it is necessary to examine the positional variations of the loop region during the course of the MD simulations in ACE2−RBD complexes with and without the presence of octa in the RBD binding site. The loop region showed larger positional variations in the ACE2−RBD complex (as high as 30 Å) than those in the ACE2−RBDocta complex (under 10 Å) as seen from their positional RMSD values ( Figure S15 ). Other theoretical studies also showed that this loop region displayed higher dynamical Journal of Chemical Information and Modeling pubs.acs.org/jcim Article motions without the presence of the oligosaccharide. 59 It can be speculated that the dynamic motions in the loop region were significantly reduced in the ACE2−RBD complex with the presence of octa in the RBD binding site. The principal component analysis (PCA) for the motions of the loop region was performed to further study the correlations between the dynamic motion changes for the loop region and the presence of octa in the RBD binding site. The PCA results agreed with the positional variations observed through RMSD calculations ( Figure 5 ). The motions of the loop were stronger as indicated by the longer arrows when octa was absent than those when octa was present. More importantly, the top motion component of the loop in RBD without octa was toward the binding site of octa, suggesting that the vacancy of the binding site for HP/HS in RBD offered more freedom for the loop region. Conversely, the spatial freedom of the loop would be reduced when octa is present in the binding site, and so that the motion. Given that the position of the loop region in RBD is between the binding sites for ACE2 and HP/HS, the presence of octa in the RBD binding site would block the movements of the loop and reduce its motions, which in turn would promote its contact and binding interactions with the N-terminal helix of ACE2, as observed in the MD simulation of the ACE2−RBD-octa complex.

The total binding energy contribution for RBD residues in the ACE2−RBD binding interface (Table S6) , ones within 10 Å of ACE2, were −49.2 ± 3.2 kcal/mol with the presence of octa, which was 15% more negative than that without the presence of octa (−42.9 ± 3.1 kcal/mol). This enhanced binding interaction agreed energetically with the cofactor role of octa in the recognition of SARS-CoV-2 S protein by ACE2. There were 11 residues (T478 to C488) in the loop region, which were included in the binding energy calculations. The difference of the total contributions from these 11 residues in the ACE2−RBD-octa (−4.6 ± 1.8 kcal/mol) and ACE2−RBD (0.7 ± 0.7 kcal/mol) complexes were −5.3 kcal/mol, while the difference of the contributions from the other binding residues, including those interacting with the hotspot residues, was only −1.0 kcal/mol. The disproportion of the binding contribution differences caused by the presence/absence of octa between the loop region at the distal end and the main binding interface suggested that the loop region is a key variable for the ACE2− RBD binding strength.

In summary, the presence of octa in the RBD binding site blocked the movements of the loop region at the distal end of the ACE2−RBD binding interface, which in turn promoted the contacts of the loop region with the ACE2 N-terminus helix and enhanced the binding interactions of the ACE2−RBD complex.

Glycosaminoglycans involve in many aspects of different viral activities. In the infections of SARS-CoV-2, HP/HS has shown to be a cofactor for the recognition of SARS-CoV-2 S protein by ACE2. Thus, understanding the mechanisms for HP/HS interacting with RBD and promoting ACE2−RBD interactions could provide useful theoretical insights into the therapeutics against COVID-19. The microarray study for the RBD of SARS-CoV-2 S protein using the HP/HS library has revealed that the binding of HP/HS to RBD is length-and sequencedependent; nevertheless, the lack of crystallographic data for the RBD−HP/HS complexes has been an obstacle for further studying the mechanisms at the atomic level.

In the present study, putative poses for octa, an oligosaccharide with four repeating units of GlcNS6S-IdoA2S, in the binding site of RBD of the SARS-CoV-2 S protein were generated by Vina-Carb that includes penalties for energetic unfavored glycosidic linkage conformations in the docking scoring functions. Only three of the top six ranked poses maintained their initial docked poses and bound stably to RBD along the course of their MD simulations. Coincidently, they shared similar conformational and spatial characters. The representative conformation of octa by MD simulations starting from these three docked poses was employed in building models for RBD complexing with hexasaccharides, RBD-hexa-1 and RBD-hexa-2, to reproduce the previously observed length-and sequence-dependent behavior before further mechanism studies. Hexa-1 and hexa-2 displayed weaker binding strengths than octa, which agreed with the length-dependent behavior. The 2-O-sulfate and 6-Osulfate groups of all three IdoA2S and GlcNS6S residues, respectively, in hexa-1 and hexa-2 formed stable intermolecular hydrogen bond interactions and contributed to the binding interaction energies, which agreed with the sequence-dependent behavior of the IdoA2S residues. A planar U-turn shape was observed in the five consecutive residues in octa, and such a shape was in favor of the HP/HS binding to RBD. The oligosaccharides with four and three repeating units of GlcNS6S-IdoA2S showed binding to RBD experimentally, but not the tetrasaccharide with two repeating units, which suggested that this planar U-turn may be a general structural form for HP/HS binding to RBD. These agreements suggested that the docked pose for octa in the RBD possessed the binding characters of HP/HS in RBD in the natural cells, if the actual binding pose is not reproduced.

The ACE2−RBD complex with octa in its validated representative conformation and without octa in the RBD binding site both reproduced the experimentally observed intermolecular interactions between the residues in the ACE2binding ridge and those in the corresponding RBD binding surface. In addition to the key hydrogen bond interaction pairs from hotspot 31 and 353, several other stable hydrogen bond interactions were also observed in its MD simulation. The total binding energy contributions from the RBD residues in the ACE2−RBD binding surface in the ACE2−RBD-octa complex were greater than those in the ACE2−RBD complex, agreeing to the experimental observations that the presence of HS would promote the binding between ACE2 and RBD. The interactions between N481 from the loop in RBD and Q24 in the ACE2 N-terminus were observed in the ACE2−RBD-octa complex, but not in the ACE2−RBD complex in the MD simulations, which possibly caused the binding energy differences between these two complexes. Further analyses showed that the loop containing N481 in RBD displayed more positional variations in the ACE2−RBD complex than the ACE2−RBD-octa complex. The decompositions of the motions showed that the loop carried movements toward the HS binding site in the ACE2−RBD complex. Therefore, given that the loop is adjacent to both octa and ACE2, it is reasonable to speculate that octa in the RBD binding site blocked such motion and reduced the freedom of the loop, which promoted the loop to interact with ACE2 on the other side. 

The proximal origin of SARS-CoV-2

A Novel Coronavirus from Patients with Pneumonia in China

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

Veesler, D. Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

Exploitation of glycosylation in enveloped virus pathobiology

Structure, Function, and Evolution of Coronavirus Spike Proteins

An update on ACE2 amplification and its therapeutic potential

Targeting transcriptional regulation of SARS-CoV-2 entry factors ACE2 and TMPRSS2

Structural basis of receptor recognition by SARS-CoV-2

Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved alpha-ketoamide inhibitors

Molecular Architecture of the SARS-CoV-2 Virus

Structure-based design of prefusion-stabilized SARS-CoV-2 spikes

Distinct conformational states of SARS-CoV-2 spike protein

Structural and simulation analysis of hotspot residues interactions of SARS-CoV 2 with human ACE2 receptor

Analysis of the SARS-CoV-2 spike protein glycan shield reveals implications for immune recognition

Virus-Receptor Interactions of Glycosylated SARS-CoV-2 Spike and Human ACE2 Receptor

Glycosylation of SARS-CoV-2: structural and functional insights

Vulnerabilities in coronavirus glycan shields despite extensive glycosylation

Deducing the N-and O-glycosylation profile of the spike protein of novel coronavirus SARS-CoV-2

Site-specific glycan analysis of the SARS-CoV-2 spike

How glycobiology can help us treat and beat the COVID-19 pandemic

Biomechanical characterization of SARS-CoV-2 spike RBD and human ACE2 protein-protein interaction

Polysulfates Block SARS-CoV-2 Uptake through Electrostatic Interactions

Characterization of heparin and severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) spike glycoprotein binding interactions

Heparin Inhibits Cellular Invasion by SARS-CoV-2: Structural Dependence of the Interaction of the Spike S1 Receptor-Binding Domain with Heparin

The 2019 coronavirus (SARS-CoV-2) surface protein (Spike) S1 Receptor Binding Domain undergoes conformational change upon heparin binding

The Structure of Glycosaminoglycans and their Interactions with Proteins

Heparin -Protein interactions

Heparin and heparan sulfate: structure and function

Heparan Sulfate Proteoglycans as Attachment Factor for SARS-CoV-2

Extension and validation of the GLYCAM force field parameters for modeling glycosaminoglycans

Perspective on computational simulations of glycosaminoglycans

No. e1388. (39) Pomin, V. H. Solution NMR conformation of glycosaminoglycans

SWISS-MODEL: homology modelling of protein structures and complexes

Improving Glycosidic Angles during Carbohydrate Docking

GLYCAM06: A generalizable Biomolecular force field. Carbohydrates

Simmerling, C. ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution

Molecular-Dynamics with Coupling to An External Bath

Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born

Algorithms for Macromolecular Dynamics and Constraint Dynamics

Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models

py: An Efficient Program for End-State Free Energy Calculations

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Assessing the Performance of the Molecular Mechanics/Poisson Boltzmann Surface Area and Molecular Mechanics/Generalized Born Surface Area Methods. II. The Accuracy of Ranking Poses Generated from Docking

Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations

Downstream Products are Potent Inhibitors of the Heparan Sulfate 2-O-Sulfotransferase

Engineering the affinity of a family 11 carbohydrate binding module to improve binding of branched over unbranched polysaccharides

3D implementation of the symbol nomenclature for graphical representation of glycans

Importance of Ligand Conformational Energies in Carbohydrate Docking: Sorting the Wheat from the Chaff

Semi-Rigid Solution Structures of Heparin by Constrained X-ray Scattering Modelling: New Insight into Heparin-Protein Complexes

Structural Snapshots of Heparin Depolymerization by Heparin Lyase I

Regulation Mechanism for the Binding between the SARS-CoV-2 Spike Protein and Host Angiotensin-Converting Enzyme II

The authors declare no competing financial interest. Vina-Carb is free for academic use and available at www. glycam.org. AMBER18 is available at www.ambermd.org. The MD simulation trajectory files are available upon request.