Posttranslational modifications optimize the ability of SARS-CoV-2 spike for effective interaction with host cell receptors

Significance SARS-CoV-2 spike protein, which forms the basis for high pathogenicity and transmissibility of the virus, is a prime target for the development of both diagnostics and vaccines for the debilitating disease caused by the virus. We present a full model of spike methodically crafted and used to study its atomic-level dynamics by multiple microsecond simulations. The results shed light on the impact of posttranslational modifications on the pathogenicity of the virus. We show how glycan–glycan and glycan–lipid interactions broaden the protein’s dynamical range and thereby, its effective interaction with the surface receptors on the host cell. Palmitoylation of the spike membrane domain, however, results in a unique deformation pattern that might prime the membrane for fusion.

The pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) continues to impact nearly all aspects of human life (1,2). SARS-CoV-2 belongs to the beta coronavirus genus, the same as SARS-CoV and Middle East respiratory syndrome coronavirus (MERS-CoV), which were responsible for the severe acute respiratory syndrome outbreak in 2003 and the Middle East respiratory syndrome in 2012, respectively (3,4). The first step in viral infection by SARS-CoV-2 and other coronaviruses in general is the binding of its extended spike glycoproteins to specific human surface receptors (5,6). Spike proteins contain all components necessary for viral entry into human cells; after binding to the surface receptors, they initiate the fusion of the viral and human membranes and subsequent release of the viral genome into the host cell (7). Characterizing the mechanism of action of the spike protein thus remains key to our understanding of the critical steps in viral infection, paving the way for development of therapeutic strategies against the virus.
Several SARS-CoV-2 spike cryo-EM structures have been resolved (8)(9)(10)(11)(12)(13)(14)(15). Although providing essential information on the structural details of the spike globular domain (spike head), which contains the receptor-binding domain (RBD) binding the Angiotensin-converting enzyme 2 (ACE2) host cell receptors (6), the structures still lack functionally important regions (Fig. 1A). These include the conserved fusion peptide involved in initiating the fusion (16,17), the stalk heptad repeat 2 (HR2) domain (18), and the transmembrane (TM) domain containing multiple palmitoylation sites that stabilize the spike in the envelope and thus, are essential for its assembly and activity in other viruses (19,20). The structures also lack critical information about the glycosylation of different residues, known to be essential in mediating protein folding and shaping viral tropism as well as shielding the virus from immune recognition (21).
Additionally, the assembly and budding of the virus are known to take place in the endoplasmic reticulum-Golgi intermediate compartment (ERGIC) (22,23), an environment that is often not taken into account in structural studies.
Significance SARS-CoV-2 spike protein, which forms the basis for high pathogenicity and transmissibility of the virus, is a prime target for the development of both diagnostics and vaccines for the debilitating disease caused by the virus. We present a full model of spike methodically crafted and used to study its atomic-level dynamics by multiple microsecond simulations. The results shed light on the impact of posttranslational modifications on the pathogenicity of the virus. We show how glycan-glycan and glycan-lipid interactions broaden the protein's dynamical range and thereby, its effective interaction with the surface receptors on the host cell. Palmitoylation of the spike membrane domain, however, results in a unique deformation pattern that might prime the membrane for fusion.
Computational simulations offer an effective technique for modeling complex systems like the spike at a detailed level. Recent classical molecular dynamics (MD) simulations of the spike head (24)(25)(26)(27)(28)(29) as well as constructed full-length spike models (30)(31)(32) have provided powerful insights into functionally relevant features [e.g., how glycans can shield the spike from the immune system (24,25,28,31,32) or the intradomain dynamics of the spike head that may facilitate the RBD to switch between the "up" and "down" (active or inactive) conformations (26)(27)(28)(29)31)]. However, we are still far from fully understanding how the global dynamics of the spike may contribute to its function, specifically in relation to effectively locating its target receptor in crowded cellular environments. In this context, characterizing how glycosylation may play a role in the protein's global conformational dynamics is as important. Furthermore, the role of palmitoylations in the TM/endodomain, which are known to participate in the modulation of the membrane curvature and in mediating cell fusion in other viruses (33)(34)(35)(36), remains lacking in SARS-CoV-2.
We constructed a full-length, membrane-embedded spike and investigated its dynamics by atomistic MD simulations. Taking a hybrid approach combining homology modeling, proteinprotein docking and MD, and the available experimental data (37)(38)(39)(40), a full-length, palmitoylated, and fully glycosylated spike is modeled. After validating the modeled parts, multimicrosecond MD simulations of the full spike were used to characterize its global dynamics, focusing particularly on the role of glycans. Control simulations of the nonglycosyslated spike highlighted the importance of the glycans in the global motion of the protein.
Palmitoylation of the protein is found to modulate the membrane curvature, potentially facilitating its fusion with the host cell. This study provides insight into the dynamics used by the spike to sample the surface of its host cell and the direct role of the glycans in regulating this process while successfully allowing the virus to evade the host immune response.

Results and Discussion
Here, we first develop a full spike structure in its membranebound state and use it for microsecond-scale MD simulations to characterize its conformational dynamics. We describe these conformational changes in terms of bending and twisting motions of the spike head with respect to the stalk and examine the potential role of glycan-lipid and glycan-glycan interactions in these motions. Sequence comparison with other human coronaviruses allows us to delineate multiple conserved hinge regions in the stalk that likely promote the conformational motions of the spike head. Furthermore, we quantify the effects of palmitoylation in the TM domain, offering support for the role of this posttranslational modification in modulating the local membrane curvature.

Construction of the Membrane-Embedded Full Spike Protein.
We combined homology modeling utilizing multiple template structures from homologous coronaviruses, like SARS-CoV and MERS-CoV, with fragment-based, ab initio modeling and secondary structure prediction to construct the missing regions in the spike head and the stalk ( Fig. 1 A and B and SI Appendix, Fig. S1). Additionally, we carried out an expansive procedure for deriving a stable model of the TM domain in a native membrane (SI Appendix, Table S1), where we combined sequence-based secondary structure prediction for the TM region (SI Appendix, Fig. S2) with protein-protein docking for generating the initial configurations of TM trimers (SI Appendix, Table S3) followed by MD relaxation of suitable TM trimeric configurations in membrane (SI Appendix, Figs. S3 and S4). The above-described modeling procedure was combined with the biochemical data on the palmitoylation sites (37) and the recently reported glycomics data to determine the glycosylation compositions of specific Asn residues (38, 39) ( Fig. 1C and SI Appendix, Table S2) to develop the full-length, palmitoylated, and fully glycosylated spike structure (Fig. 1D).
We have also carried out structural comparison of our fulllength spike model with previously published models (30,31,41). The individual monomers of the spike head display relatively stable conformations after the initial equilibration (SI Appendix, Fig. S5). The HR2 domain in the previous studies was modeled using a coiled-coil template from Salmonella autotransporter adhesin (31), using a de novo coiled-coil modeling program (41), or using the HR2 from SARS-CoV (30). Here, we also make use of the HR2 template from SARS-CoV (42), resulting in a stable model throughout the simulations (SI Appendix, Fig. S5). As SARS-CoV and SARS-CoV-2 display complete sequence identity in the HR2 domain, our stable model is expected to provide a closer physiological representation of this domain. The TM domain in the previous spike models (30,31) was constructed using a template structure from HIV (43). It is noteworthy that the TM helices in HIV each contain a "GxxxG" motif, known to be involved in inter-TM helix interactions (44). Consequently, the TM domain assembly in HIV displays interhelical contacts in this region with the C-terminal half only held together by polar contacts (43). As the GxxxG motif is present in the central region of the TM helices in the SARS-CoV-2 spike (Fig. 1A), its contacts between the TM helices can be expected to lie in this region instead. The TM domain in our spike model is a trimeric assembly, with the majority of the interhelical contacts formed between the central regions of the helices containing the GxxxG motifs ( Fig. 2 A and B). Additionally, the TM domain maintains a relatively stable and symmetric structure during the simulations (Fig. 2 C-F and SI Appendix, Fig. S5). The correct arrangement of the different stalk domain structures, including the trimeric HR2 and TM domains, can be critical to the correct description of the spike's global dynamics.
Bending and Twisting Motions of the Spike Head around the Stalk. To further describe the dynamics employed by the spike to sample the crowded cellular surface, we characterized its global motions in terms of the movements of the spike head with respect to the stalk region during the simulations (Movie S1). These motions were quantified in terms of head-twist and head-bend angles as well as head distances calculated separately with respect to the spike neck, the HR2 domain, or the TM domain (  fully glycosylated, native form, the spike can sample large headtwist angles (in the xy plane) with respect to both the HR2 and TM domains, highlighting the large conformational flexibility available along this degree of freedom ( Fig. 3 B, Upper, F, and G).
A wide normal-like distribution of head-twist angles is observed with respect to the HR2 domain, with the highest sampled configurations turned by almost 90 • relative to the starting structure.
With respect to the TM domain, although lowly populated, the spike head can sample conformations that are almost oppositely facing (turned by as much as 160 • ) compared with the starting orientation. A relatively smaller range of motion is allowed around the spike neck, with the most populous regions lying close to the starting structure ( Fig. 3 B, Upper and E). These results can be directly related to the length of the linkers connecting the different domains of spike. The domains around which the spike head displays larger motions (HR2 and TM domains) are connected by longer linkers (linker 2: 13 residues; linker 3: 11 residues) compared with the domain around which the spike head displays relatively more restricted motions (spike neck), which are connected to the spike head by a smaller linker (linker 1: 7 residues). Furthermore, direct comparison of the head-twist motions in the glycosylated form with the control nonglycosylated system shows surprisingly that these motions become restricted in the latter (Fig. 3 B, Lower), especially with respect to the headtwist angles calculated around the HR2 and TM domains. These results point to an unexpected enhancing role of glycans in the overall dynamics of the spike. Similarly, the spike head displays a broad range of head-bend angles with respect to the HR2 domain, sampling configurations almost orthogonal to the starting structure ( Fig. 3 C, Upper and I ). Comparatively, smaller head-bend motions are observed around the spike neck and the TM domain, with the respective populations showing distributions over a range of 0 • to 60 • angles ( Fig. 3 C, Upper, H, and J ). The extent of these bending angles is also directly related to the changes in the head distances with respect to these domains, with the largest distance changes observed in conformations with the largest head-bend angles (Fig.  3 D, Upper). Thus, the spike head can move as much as 125Å with respect to the HR2 domain compared with its starting configuration. Again, comparison of the glycosylated and nonglycosylated forms shows that the spike domains rigidify (display a relatively smaller range of motion) in the absence of glycans (Fig. 3 C and D, Lower). The potential role of the glycans in modulating global motions of the spike is further discussed in later sections.
Mechanistically, these motions (bend/twist/translation) of the spike head represent major degrees of freedom for the RBD (located in the spike head) to locate and bind to the target ACE2 receptors in the host cell membranes. As spike proteins have been found to be sparsely distributed on the virion envelope [recent cryo-EM and tomography studies (12,45)], the global flexibility of the spike head and its allowed range of motions may play an important role in expanding the search space covered by the spikes. Cryo-EM and previous computational studies have similarly reported the conformational diversity of the spike head around its stalk (41,46). Compared with our results, a smaller conformational variability of the spike head was observed by Turoňová et al. (41); as stated by the authors, the crowded conditions generated by inclusion of four spike proteins in the simulation system might have prevented ample sampling of the available space. The crowded distribution used in that study does not represent the physiological density reported for the full virion in later cryo-EM studies (approximately one spike per 1,000 nm 2 of the viral surface area) (45). While previous simulation studies on the spike protein have concentrated on the interdomain motions in the protein (41,46) or on the conformational dynamics of the RBD (open/closed transitions) (27)(28)(29) and its response to glycosylation (28,31), in the present study we focus more on how glycans may modulate the global motions of the entire spike, particularly with respect to the positioning of the head domain with respect to the viral envelope (discussed in later sections).
Sequence Conservation in Multiple Hinge Regions. In order to examine potential generality of the observed global motions in the SARS-CoV-2 spike, we analyzed the sequence similarity of the spike protein, specifically focusing on the flexible regions connecting the different domains. These include the small connecting linkers between the spike head and the neck (linker 1: NTVYDPL), between the spike neck and the HR2 domain (linker 2: HTSPDVDLGDISG), and between the HR2 domain and the TM domain (linker 3: GKYEQYIKWP). The three linkers are found to be rich in Pro and Gly, residues known to contribute to the conformational flexibility in proteins (47,48) (Fig. 4A). Additionally, the Pro residues in these regions may be important for the hinge or recoil motion suggested to take place in the refolding of the spike during the fusion process, which may result in an elongated postfusion state of the protein (15).
To further elucidate the possible functional importance of these residues in the mobility of the spike, multiple sequence alignment was carried out with other coronaviruses infecting humans. Interestingly, all Pro and Gly residues present in the flexible regions are completely conserved in the homologous bat-CoV and human SARS-CoV and to a large extent, in human MERS-CoV, where only Pro-1163 is not conserved (Fig. 4B). Most of these residues are not conserved in other human coronaviruses except Pro-1213, which is completely conserved in all coronaviruses. We propose that the presence of these conserved flexible residues in the regions connecting the different domains not only promotes the observed conformational flexibility of the spike head around the stalk domains but compared with other coronaviruses, also provides an evolutionary advantage to SARS-CoV-2, SARS-CoV, and MERS-CoV, coronaviruses known to be highly infective in humans.
High glycan densities are also observed around both these flexible regions and the stalk overall in addition to the high densities observed around the spike head, providing a global shielding mechanism for the spike against the immune response and neutralizing antibodies (SI Appendix, Fig. S6A). The accessibility of the spike protein, measured independently for the glycosylated and nonglycosylated systems, is found to increase monotonically with increasing probe radius (SI Appendix, Fig. S6B), consistent with the previous MD results (31). In addition to the local, residue-wise profile of glycan shielding (28,32), our calculation of the global shielding effects examined here provides a more global view of how glycosylation may prevent antibodies targeting the spike.
The spike protein is composed of two subunits: S1, which contains the RBD that recognizes and binds to the human ACE2 (hACE2) receptor, and S2, which includes the stalk domain along with the fusion machinery that mediates the viral membrane fusion after cleavage and separation from the S1 subunit. Recent studies have shown that antibodies targeting the S2 subunit exhibit specific neutralizing activity against SARS-CoV-2 (49). Analyzing clinical samples with virus-reactive memory B cells suggests the presence of circulating immunoglobulin G antibodies reactive to the S2 subunit (50). Furthermore, in silico docking studies have also predicted binding of naturally occurring compounds like flavonoids to the S2 subunit (51). Even though direct evidence for the involvement of the stalk domain and the Pro-and Gly-rich regions of the S2 subunit in interacting with above-mentioned antibodies/small molecules is still lacking, the dynamics of the hinge regions seem to be a highly relevant feature of the protein B A interacting with other molecules and may provide an additional strategy in developing therapeutics against the virus.

Stalk Bending Motions Modulated by Glycan-Lipid and
Glycan-Glycan Interactions. The extensive glycosylation of the spike surface is key to the viral escaping from the immune system (28,31,32,41,46). Moreover, the RBD opening/closing transition, which has been proposed to equip the virus with a "conformational masking" mechanism to hide active RBD from the immune system (28), may be affected by the glycans. Sztain et al. (29) performed over 130 μs of weighted ensemble simulations to reveal a gating role of the glycans on the RBD, and Zimmerman et al. (28) identified cryptic epitopes by capturing dramatic RBD openings calculated for both glycosylated and unglycosylated spikes. However, the underlying mechanism for the dynamics of its stalk region is still undercharacterized. Using our sampling, we have explored how glycans near the hinge regions may modulate their flexibility and affect the global conformation of the spike. We also find that the mobility of the spike stalk would be altered in the absence of glycosylation in our control simulation.
We noted that the glycans also form contacts with each other, as well as with the lipid bilayer, during the simulations (examples of these phenomena are in Movies S2 and S3). To obtain deeper insights into the impact of these glycan interactions, especially with respect to the spike global motions, we evaluated the relationship between the inclination of the different stalk domains and the glycan-lipid or glycan-glycan contacts formed both in the native glycosylated spike simulations and in the control nonglycosylated spike. In the glycosylated spike when the inclination angle of the HR2 domain exceeded a threshold of around 20 • , we observed formation of more glycan-lipid or glycan-glycan contacts, indicating a correlation between the bending of the HR2 domain and the glycan interactions (Fig. 5). On examining the distribution of the HR2 tilt angles toward the membrane in the glycosylated spike simulation, only ∼20% of the trajectory was found to lie between 0 • and 20 • (vertical orientation) compared with ∼80% of the nonglycosylated trajectory that was below the 20 • inclination. It can be inferred that in the fully glycosylated spike when no glycan-lipid interactions are formed, the HR2 domain can move freely within the 20 • inclination (Fig. 5 A and B); the glycan-lipid interactions appear when the inclination angle is between 20 • and 50 • , further restraining the flexibility of the HR2 hinge motion in this range. A Spearman's coefficient of 0.59 ± 0.08 between the glycan-lipid contacts and HR2 bending further supports the observed correlation.
The HR2 neck inclination shows similar properties as the HR2 membrane inclination, where larger inclination angles, ranging between 20 • and 90 • , mostly correspond with high glycanglycan contacts (Fig. 5 C and D). around the flexible linkers can further enhance the relative motion between the stalk domains, in turn increasing the probed space by the spike head. N -glycans have been reported to be involved in mediating receptor recognition and facilitating cell adhesion (52). The selfrecognition and conjugation have been previously demonstrated for N -linked complex glycans with multiple N -acetylglucosamine (GlcNAc) termini (53,54). The glycans located at the bottom end of the neck region (N1158) as well as the top (N1173) and bottom (N1194) ends of the HR2 domain are all complex glycans containing multivalent GlcNAc (SI Appendix, Table S2). Additionally, these glycans are completely conserved in human coronaviruses (Fig. 4B). Overall, we propose that these complex glycosylation sites not only shield the stalk from epitopes by occupying the space around the flexible hinges connecting the neck, HR2, and TM domains but also, contribute to the flexibility of the stalk, which may, in turn, increase the range of motion of the spike head, allowing for more effective sampling of the human cell surface.
On further examination of the glycan-lipid interactions, we observed that the majority of these contacts are formed with POPC or POPE, followed by PSM (SI Appendix, Fig. S7). The majority of the interactions between glycans and lipids consist of hydrogen bonds formed between O or N atoms of the N -acetyl-Dneuraminic acid or D-galactose sugars present in the glycans and N atoms from PE head group or O atoms in the POPE/POPC phosphate groups. These interactions may further stabilize the inclined conformations of the HR2 domain with respect to the membrane. The contacts with cholesterol are minimal, which is expected as cholesterol is known to insert more deeply in the membrane core compared with the phospholipids (55,56).
Although because of the unavailability of structural information for unglycosylated spike proteins, no direct evidence for modulation of spike dynamics by glycosylation currently exists, previous studies on other glycoproteins indicate that glycans may play a role in their structure and dynamics. These include mass spectrometry results showing that oligosaccharides in immunoglobulins affect the conformational behavior of their hinge regions (57), where mutations of the glycosylation sites resulted in a change in the rigidity of the interdomain linkers (58).
Furthermore, large-scale bending motion of a membrane glycoprotein ectodomain observed by electron microscopy was found to be modulated by the conserved N -glycosylations proximal to the joint between its domains (59). It can, therefore, be expected that similar conformational modulatory effects of glycans may be also present in the SARS-CoV-2 spike protein.
Correlated Interdomain Motions in the Spike. The global conformational dynamics of the spike were further examined for the presence of correlated motions between its different domains. A correlation between the motions of the spike head and the TM domain about the HR2 domain is evident in the fully glycosylated spike (Fig. 6 A and B and Movie S4). A strong Pearson's coefficient of 0.84 ± 0.07 between the Head1 and TM vectors (Fig. 6  A and B), representing the spike head's coplanar motion with respect to the HR2-TM domain bending, further confirms the observed correlation in the motions of these domains. This is an interesting observation in that the orientational motions of the different domains of the spike appear to remain interrelated, even though they are connected by flexible hinge regions. Without glycosylations, not only is the magnitude of interdomain motions dampened (SI Appendix, Fig. S8), but at the same time, the correlation between them is also reduced, with the nonglycosylated spike showing a relatively lower Pearson's coefficient, 0.76 ± 0.07, for the in-plane motion of the spike head with respect to the HR2-TM domain bending. In general, the correlated motions between the different spike domains can be attributed to its unique tertiary structure and can be considered an intrinsic property of the protein, but at the same time, glycosylation may further amplify these motions through the extensive glycan-lipid and glycan-glycan interactions in the system (Fig. 5). For example, interactions of the glycans located at the bottom end of the HR2 domain with the lipids seem to promote its bending toward the membrane that may, in turn, facilitate the formation of more contacts between the glycans located at the top/opposite end of the HR2 domain and glycans in the neck region, promoting bending of the spike neck and head in the opposite direction.
The the spike head and the TM domain are oriented in the same direction about the HR2 domain for 88% of the simulation in the fully glycosylated spike (Fig. 6C). The nonglycosylated spike displays a similar behavior, with the orthonormal vectors oriented in the same direction for 89% of the time, further highlighting that these motions are largely a property of the spike structure itself. Topologically, this reflects a configuration where the spike head and the TM domain orient diametrically opposing each other around the HR2 domain, a configuration that may further extend the reach of the spike head and the RBD for the target cellular receptors (Fig. 6D). Additionally, considering the fact that the inclination of TM domain in the membrane can be substantial, ranging between 0 • and 30 • (SI Appendix, Fig. S9), the concerted, opposing bending of the spike head and the TM domain may pull the head back to an upward position with respect to the membrane, preventing it from severe inclination and falling into the membrane and thus, maintaining a proper configuration for binding with the ACE2 receptor.

Role of Palmitoylations in Modulating Membrane Curvature.
Protein S-palmitoylation, the covalent lipid modification of the Cys residues, widely exists in viral membrane proteins and is reported to contribute to membrane fusion and other processes (33)(34)(35). To characterize the effect of spike palmitoylations on the membrane shape, we compared the membrane curvature for the control simulation of nonpalmitoylated spike with the palmoytilated (and fully glycosylated) one, taking into account the periodicity in fitting the lipid head groups. The two simulations share a similar membrane size, but the palmitoylated system seems to deform the membrane more than the nonpalmitoylated system ( Fig. 7 A and B). Keeping in mind that the periodic boundary conditions used in the simulations can suppress longrange curvature effects, comparing the local curvature (in the range of 20 to 50Å from the center of the TM domain) shows that the palmitoylated TM domain is found within a positively curved membrane, whereas the nonpalmitoylated TM prefers a saddle point. Thus, the palmitoyl tails may induce a positive curvature locally in the vicinity of the TM domain. The observed effect of the palmitoyl tails on the local membrane properties in the vicinity of the protein is complex, primarily and largely dominated by the protein segment they are connected to (Fig.  7C). Such tethering effects may promote any pattern of membrane curvature (positive or negative) depending on the protein context.
The lipidation of other viral proteins is relatively well studied. In influenza virus, cryogenic electron tomography (cryo-ET) measurements have shown membrane curvature induced by hemagglutinin (HA) palmitoylation (34). Apart from this, nonacylated HA mutants of influenza virus show a severe negative effect on the membrane fusion step of the infection process (60). As shown by the results above, SARS-CoV-2 spike palmitoylation may similarly contribute to the bending of the membrane,  especially locally, and thus, may play a role in the viral fusion with the host cell membrane.

Conclusions
The spike protein of SARS-CoV-2 constitutes a major component in the viral infection and therefore, the main target for both diagnostics and vaccines developed for the disease caused by the virus. To better understand the many mechanistic steps the spike is involved in, it is imperative to characterize its structure and dynamics in the most detailed manner and under the most realistic conditions possible. We report here a full model for the full-length spike in its native, glycosylated, palmitoylated, and membranebound form, which we use for several microseconds of atomistic simulations highlighting some of the dynamics employed by the spike to search for and locate the host cell receptors effectively (Fig. 8). The hinge regions identified in the stalk domain directly modulate the conformational landscape of the spike head, where the RBDs reside, and may thus represent important structural elements regulating the effective role of the spike in binding its receptor. These highly conserved hinges may offer targets for developing alternative antibodies and therapeutics that bind them and modulate their structure/dynamics and therefore, the effectiveness of the spike protein. In addition, we propose possible roles for posttranslational modifications of spike, specifically glycosylation and palmitoylation, in modulating its dynamics. The results provide a deeper understanding of the intricate roles and functions of the different elements of the spike, each evolved in order to maximize the chance of successful infection and transmission of the genetic material to the host cell. Extended simulations studying the spike behavior in a full viral envelope, as recently done for the SARS-CoV-2 Delta variant (61), will allow further characterization of the impact of a more crowded environment of the envelope and its proteins in modulating the structure and dynamics of the spike and will provide additional insight into the inner working of the virus.

Materials and Methods
An extended version of the methods is available as SI Appendix.
As a first step, a full-length spike was constructed using recent cryo-EM structures of the spike head from SARS-CoV-2 [Protein Data Bank (PDB) ID codes 6VYB and 6VXX (9)], SARS-CoV [PDB ID codes 5X58 (62), 5XLR (63), and 6CRW (64)], and MERS-CoV [PDB ID code 6Q04 (65)] along with the SARS-CoV-2 RBD structure in complex with human receptor ACE2 [PDB ID code 6M17 (66)], to which missing protein segments (e.g., the fusion peptide) (16,17) were added in MOE (67) using either template-or fragment-based modeling. Missing regions without a suitable template were constructed using the Robetta protein structure prediction server (68,69). The spike neck was also modeled using Robetta, and the HR2 domain was homology modeled using the structure of SARS-CoV protein [PDB ID code 2FXP (42)]. The modeled regions were tested against secondary structure predictions by JPred4 (70).
The TM domain, which is important for spike stabilization (19,71), shows low homology to the available HIV virus (30.4% identity), making it unsuitable for template-based modeling (43). After prediction of the TM region from the sequence in Transmembrane Helices; Hidden Markov Model (TMHMM) (72), a TM monomer was constructed and embedded in a membrane with the lipid composition of ERGIC (22,73), where coronaviruses assemble (23) (SI Appendix, Table S1). After 100 ns of relaxation of the monomer in the membrane, multiple trimeric configurations were generated using Multimer Docking in ClusPro (74). The two best trimeric TM configurations were embedded in the membrane, solvated, and equilibrated for 200 ns each. The stability of these models was evaluated by calculating rmsd, TM tilt/inclination in the membrane, and intermonomeric coordination number [coordNum in the collective variables (COLVARS) module (75,76) of Visual Molecular Dynamics (VMD) (77)]. The most stable trimeric conformation was used to construct the full-length spike structure. The structure for the luminal C-terminal region was modeled with Robetta. The structure of the full-length spike was then constructed by assembling the experimental spike head and the individually constructed domains. Missing loops were then filled with MODELER (78). Based on mass spectrometry data (38,39), we then added a total of 22 N-glycans and 1 O-glycan (39) to each spike monomer (SI Appendix, Table S2) using Glycan Reader & Modeler (30,79) in CHARMM-GUI (80,81). Palmitoylations at the endodomain (37) were modeled based on SARS-CoV data (37) (residues 1,236, 1,240, and 1,241).
MD simulations of the full-length, membrane-embedded spike structure in both glycosylated and nonglycosylated forms (5 μs each) were then performed in NAMD (82,83) using CHARMM36m (84,85). The simulations were performed as an NPT ensemble (310 K and 1 bar) employing Langevin thermostat and barostat (86,87) with the particle mesh Ewald method for electrostatic forces (88) and applying SHAKE to constrain bonds with hydrogen atoms (89).
The global dynamics (SI Appendix, Fig. S10) of the protein were characterized in terms of the positional and orientational behavior of the spike head with respect to the different stalk domains. The potential effect of glycosylation on these motions was analyzed in relation to their direct interactions with the lipid bilayer as well as contacts between the glycans at the interfaces of different domains. Since the palmitoylation tails contribute to membrane deformation (34), we also quantified the potential effect of palmitoylation on the membrane curvature. Multiple sequence alignment was carried out using the Multiple Alignment using Fast Fourier Transform (MAFFT) program (90) and visualized using Jalview (91). The solvent accessible surface area was calculated with the Shrake-Rupley algorithm (92).