The outbreak of Coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) presents a tremendous threat to global public health. It caused over 35 million confirmed cases and more than 1 million deaths as of October, 2020. Due to unavailability of antiviral medicines or approved vaccines, the current treatment strategy is supportive care to relieve symptoms and isolation of infected individuals to reduce transmission, which has placed a huge burden on the public healthcare system and led to massive social and economic distress.
SARS-CoV-2 is an enveloped virus with a positive-sense single-stranded RNA genome 1 . The spike (S) protein anchored in the viral envelope is a class I fusion protein that mediates receptor binding and host cell entry by interacting with human angiotensin converting enzyme-2 (ACE2) [2] [3] [4] , and it is also the target of a variety of neutralizing antibodies [5] [6] [7] [8] . S protein is a homo-trimer and each monomer has two subunits (S1 and S2) separated by a cleavage site that is recognized by host proteases 9 . A number of recently published structural studies using cryogenic electron microscopy (cryo-EM) have provided a good understanding of the S protein structure at near-atomic resolution 4, 10-12 . The S1 subunit responsible for receptor binding is composed of the signal peptide (SP), N terminal domain (NTD), and receptor binding domain (RBD), and the S2 subunit responsible for membrane fusion is composed of the fusion peptide (FP), two heptad repeats (HR1 and HR2), transmembrane domain (TM), and cytoplasmic domain (CP). The three RBDs on the top of S protein head are conformationally variable. In closed conformations, all three RBDs lay flat with the receptor binding motif occluded by the RBDs on the neighboring monomers. In open conformations, one or multiple RBDs lift up and expose the receptor binding motif(s).
Although the cryo-EM structures of S protein have provided crucial information about its overall structure, highly flexible protein regions such as loops and stalk still remain unresolved. Molecular dynamics (MD) simulation provides molecular-level insight into the underlying mechanisms of biological functions that are difficult to elucidate only with experiments. Recently, cryo-electron tomography (cryo-ET) and MD simulation have been used to explore the conformational variability of S protein stalk that gives the head orientational freedom and allows S protein to scan the host cell surface 13, 14 . However, it still remains unclear what structural freedoms in the stalk portion are determinant to the overall shape of S protein and its orientation, and how they affect the binding to ACE2. In addition, MD simulation along with accessible surface area (ASA) calculations have been used to estimate the impact of glycan shielding on the exposure of antibody epitopes 15 . Mutations of two glycosylation sites have been performed to study the role of two N-linked glycans in stabilizing an RBD open conformation 16 . Further investigation is still required to evaluate whether the ASA difference between glycosylated and non-glycosylated structures truly reflect the impact of glycan shielding on antibody binding, and whether glycans have more functional roles than stabilizing the open-state RBD.
In this work, we present all-atom MD simulations of fully-glycosylated full-length S protein in a viral bilayer, and multiple µs-long trajectories were generated for RBD in open and closed states and S stalk built from different models. We also performed multiple µs-long simulations of non-glycosylated S head-only systems. Our results provide deeper insight into functional roles of glycans that provide not only shielding for immune evasion, but also contribute to the trimer stability and transition of RBD open and close states. Moreover, our simulations give insight into essential structural roles of highly flexible stalk conformations in S protein binding to ACE2.
An illustrative snapshot of a fully-glycosylated full-length S protein structure in a viral membrane is shown in Figure 1 . When there are many missing residues and domains, the initial models for MD simulations need to be carefully built and validated against available experimental data. We have built the models using GALAXY protein modeling suite [17] [18] [19] for missing residues and domains, ISOLDE 20 for refining initial models against experimental density maps, and CHARMM-GUI Glycan Modeler [21] [22] [23] and Membrane Builder 24-26 for glycan and membrane building. The head of S trimer was built based on cryo-EM structures (PDB ids: 6VSB 4 and 6VXX 12 ). All three chains of 6VXX have RBD in a closed conformation. One chain of 6VSB (A chain) has RBD in an open conformation and the other two chains have RBD in a closed conformation. Two models were selected for each of HR2 linker, HR2-TM, and CP, resulting in a total of 16 structures after the domain by domain assembly. The glycan sequences selected for 22 N-linked and 1 O-linked glycosylation sites of each monomer were based on the mass spectrometry data 27, 28 . The detailed model generation is described in reference 29 . The model name follows the model numbers used for HR2 linker, HR2-TM, and CP structures. For example, "6VSB_1_2_1" represents a model based on 6VSB with HR2 linker model 1 (M1), HR2-TM model 2 (M2), and CP model 1 (M1). All 16 S protein simulation systems and trajectories are available in CHARMM-GUI COVID-19 archive (http://www.charmm-gui.org/docs/archive/covid19). Two models for RBD/NTD, HR2 linker, HR2-TM, and CP are enlarged on the right panel. The three individual chains of S protein are colored in yellow, gray, and white, respectively, while glycans are represented as red sticks. The palmitoylation sites of S protein are highlighted in cyan. The phosphate, carbon, and hydrogen atoms of the viral membrane are colored in green, gray, and white, respectively. For clarity, water molecules and ions are omitted. All illustrations were created using Visual Molecular Dynamics (VMD) 30 .
We have performed 1.25-µs all-atom MD simulation of each of 16 models (i.e., a total of 20 µs) each containing about 2.3 million atoms (see Methods). Conformational Analysis Tools (CAT, http://www.md-simulations.de/CAT/) were used for high-throughput analysis of all simulation trajectories. The root-mean-square-deviation (RMSD) time series in Figure S1 show that the head region of S protein (residue 1-1140) remains stable during the simulations with RMSD of 4 to 5 Å in most systems. The stalk region, however, exhibits highly flexible motions at the HR2 linker and HR2-TM (see Movies S1-2), which is consistent with S protein structures observed in cryo-ET 13 . To further understand the flexible stalk motion, bending characteristics of the two linker regions, defined as θ1 and θ2 (Figure 2A) , were quantified. Figure 2B shows the distributions of θ1 and θ2 for each model. Both M1 and M2 of θ1 show similar angle distributions centered at 150º (±15º) and 155º (±12º), respectively. The HR2-TM region, however, exhibited different bending motions. The M1 of HR2-TM shows a narrow distribution centered at 172º (±4º), whereas the M2 of HR2-TM shows a wide distribution centered at 155º (±14º). Twisting motions were also dependent on the HR2-TM model ( Figure  S5A ). While both M1 and M2 of HR2-linker show similar twist angle distributions ( ) centered at 66º (±46º) and 68º (±49º), respectively, the M1 of HR2-TM shows a narrow distribution centered at 99º (±18º) and the M2 of HR2-TM shows a wide distribution centered at 98º (±71º).
These bending and twisting characteristics are consistent with the secondary structure analysis. The secondary structures of HR2 linker M1 and M2 models are mostly in coil conformations during the simulation although local folding and unfolding occur in both models ( Figure S2 ). The secondary structure of initial HR2-TM M1 model mainly consists of helical structures that are mostly retained during the simulation time. On the other hand, the secondary structure of M2 initially modeled with turn and bend shows low helicity in the range of L1200-K1215 ( Figure S3 ). This indicates that the flexible motions of HR2-TM linker are strongly influenced by the secondary structure and initial model (within the current simulation time). Although the secondary structures of CP domains are different in between two models, they have no significant effect on the motions of stalk ( Figure S4 ). To further characterize the bending motions, the Pearson correlation coefficient (r) was calculated for all combinations of HR2 linker and HR2-TM models ( Figure S5B ). The r values of all cases range from -0.16 to 0.13, indicating that there is no correlation between bending motions of HR2 linker and HR2-TM and thus each linker acts as an independent hinge. Although 16 x 1.25-μs MD simulations were performed, it does not cover all possible configurations of S protein especially with such flexible two linkers. To increase sampling, utilizing the independent θ1 and θ2 characteristics, S protein orientation was resampled based on three regions: head-HR1, HR1-HR2, and HR2-TM. First, 30 HR2-TM conformations were randomly extracted from each trajectory (excluding HR2-TM M1 models), and their TM domain was superimposed to the TM domain of the initial model to resample HR2 domain motion. Second, 30 HR1-HR2 conformations were randomly extracted and superimposed to each of the resampled HR2 domains. Finally, 30 head-HR1 conformations were randomly extracted and superimposed to the resampled HR1 domains. In summary, 27,000 configurations of S protein on a viral membrane were generated. Figure 2C shows the tilt angle of resampled configurations of S protein using M1 and M2 for HR2 linker and only M2 for HR-TM. The S protein can tilt by up to 90º towards the membrane and tilt angles around 48º are most probable. This tilt angle distribution agrees well with the experimental observation 11, 13 . However, if M1 was used for HR2-TM, the tilt angle distribution of S protein becomes narrow ( Figure S6 ). This indicates that both M1 and M2 of HR2 linker are reliable models, but for the HR2-TM, M2 is more appropriate to represent S protein configurations. To further understand the contribution of each independent hinge motion on the tilt angle, S protein was resampled separately with HR2-TM only and with HR1-HR2 only. In both cases, the resampled S protein shows a narrow angle distribution compared to the experimental observation ( Figure 2C) , indicating that both linkers are necessary for full tilting motions of S protein observed in experiment.
To explore the effect of flexible stalk motion on ACE2 binding, we performed structural alignment of S protein to ACE2. The RBD in complex with full-length human ACE2 in the presence of neutral amino acid transporter B 0 AT1 (PDB: 6M17 10 ) was used for alignment. Fully independent bending and twisting motions of two stalk linkers allow us to increase the number of S protein samples. 125 head-HR1, HR1-HR2, and HR2-TM-CP conformations were separately extracted from each trajectory with 10-ns interval. Each RBD of head-HR1 conformations was first superimposed to the RBD-ACE2-B 0 AT1 complex. Then, the HR1-HR2 conformations were superimposed to each of HR1 from the previous step. Finally, the HR2-TM-CP conformations were superimposed to each of HR2 from the previous step. Figure 3A shows one of the most probable configurations of S protein-ACE2 complex. The tilting angle (θ) is defined in Figure 2 , and the distance (d) is defined by an arc length between the centers of mass (COMs) of two TM domains. As shown in Figure 3B , d ranges from 240 Å to 350 Å and θ ranges from 30º to 60º. At the most probable configuration, d and θ are about 290 Å and 46º, respectively. Note that there is approximately one S protein per 1,000 nm 2 (316 Å × 316 Å ) on the viral surface 14 . This sparse distribution of S protein suggests that receptor binding can be promoted by enough space to have orientational degrees of freedom for RBD. Moreover, it is reported that the most probable tilting angle of prefusion state is about 40º ~ 50º 11, 13 (also see Figure 2 ). This tilting angle appears to maximize the accessibility of the receptor binding motif to ACE2 (when the RBD is in an open conformation), which could account for the high infection rate of SARS-CoV-2. 
To explore RBD and NTD motions, we measured two structural features ( Figure 4A) : RBD-NTD distance (d) defined by the minimum distance between RBD and NTD, and RBD orientation angle (θ) defined by two points at each end of RBD and the third point on the center axis of S trimer. One RBD (in both open and closed states) forms a U-shaped pocket with the NTD in the same monomer, which is occupied by the neighboring RBD in a closed state. Therefore, d estimates the RBD-NTD pocket size and θ estimates the extent of RBD opening. The time series of d and θ in 6VSB and 6VXX are shown in Figures 4B and S7 , respectively. In cryo-EM structures available in the PDB, θ in the open-state RBD ranges from 134° to 153°, and θ in the close-state RBD ranges from 88° to 93°. The trajectories of fully-glycosylated models completely cover the RBD orientation angles observed in cryo-EM structures, and explore a wider range of conformational space. In particular, θ of 6VSB_A ranges from 120° to 170°, indicating that the open-state RBD is much more flexible than the closed-state RBD, and it is consistent with the RMSD and root-mean-square-fluctuation (RMSF) results shown in Figure S1 . For both 6VSB_A and 6VSB_C, the simulations cover the NTD-RBD distances observed in the PDB cryo-EM structures. In 6VSB_B, the pocket formed by NTD and RBD (chain B) is empty due to the opening of neighboring RBD (chain A), and consequently, the NTD moved close to RBD.
To investigate the impact of glycans on the transition between RBD open and closed states, we built and simulated non-glycosylated head-only systems (3 replicates for 6VSB and 6VXX) by removing all N-/O-linked glycans and truncating the stalk. It is worth noting that the RBD in non-glycosylated 6VSB_A started to close at the very beginning of trajectories in all three replicas ( Figure 4B ). In the trajectory (colored in brown), θ decreases to 110°, which is about in the middle of open and closed states, and in the trajectory (colored in purple), the RBD reverted to opening at around 0.25 μs. Since the transition between RBD open and closed states is a complicated process involving the motions of multiple domains and attached glycans, it may require much longer simulation time to capture the conformation in which the geometries of both protein and glycan satisfy the requirement for the transition to occur. Nonetheless, dramatic transitions from RBD open to more closed states in non-glycosylated 6VSB_A indicates that the RBD open state is unstable when the glycans are removed. This implies that glycosylation has critical roles in the viral entry as the RBD needs to be open in order to interact with ACE2. 16 . Such a difference could attribute to the initial models and simulation lengths. When RBD is closed (all except 6VSB_A), the RBD forms a sandwich-like arrangement with two glycans. N165 glycan is located above the RBD and N234 glycan is located below the RBD (Figure 5B ). Both glycans frequently interact with the closed-state RBD (Figure 5C ), which makes transition to an open conformation hard. The glycan attached to N343 on RBD orients toward the solvent and hardly interacts with other domains when RBD is open (6VSB_A) (Figure 5D ). When an RBD is closed (6VSB_B) but the neighboring RBD (6VSB_A) is open, N343 glycan orients toward the pocket between NTD and RBD and interacts with N165 glycan, which makes open to closed state transition of neighboring RBD difficult. When RBD is closed (all expect 6VSB_A and 6VSB_B) and the neighboring RBD is also closed, N343 glycan makes extensive interactions with the neighboring RBD (Figure 5E) , which contributes to the stability of both RBDs in closed conformation The atom contacts between RBD N343 glycan and neighboring RBD exist in more than 95% of snapshots when both RBDs are closed (Figure 5F ). Given that RBD is shown to constantly transit between open and closed states in experiment, we propose that glycans serve as a clutch that temporarily holds the RBD in an open or closed conformation, which modulates the lifetime of both open and closed states. In addition, we calculated the accessible surface area (ASA) reduced due to formation of S trimer. For example, the ASA reduction for chain A was calculated by SA + SBC -SABC, where SA, SBC, and SABC are the ASA of chain A, chains BC-only complex, and chains ABC complex (i.e., S trimer), respectively. The ASA reduction due to trimer formation was split into the portion from protein only and the portion involving glycans. The trimer interface interactions involving glycans is about 30% when the entire S1 subunit is considered, and it increases to about 40% when only RBD and NTD are considered (Figure S8 ). This suggests that glycans makes significant contributions to the stability of S trimer, which is different to the common belief that protein-protein interactions are the only dominating factor to the stability of a protein multimer.
Viruses evolve to minimize the immunogenicity by coating the exposed viral proteins with nonimmunogenic or weakly-immunogenic glycans. It is commonly believed that the glycans on viral envelope shield viruses from host immune system. To get an impression of such glycan shielding, we aligned the S head in each trajectory snapshot and the glycan distributions are shown in Figure 6A . Most glycans are very flexible and they move around in a wide range of space, which covers most of the trimer surface. However, the role of glycan is not limited to shielding. During the past decade, many glycan-dependent HIV neutralizing antibodies have been discovered and extensively studied, which target both envelope protein and glycans [31] [32] [33] .
In the cryo-EM structure of S trimer in complex with S309 antibody (PDB ids: 6WPS and 6WPT 34 ), S309 interacts with the glycan attached to N343.
To explore the role of glycans in antibody binding, we used TM-align 35 to superimpose the RBD in RBD-antibody complex structures in the PDB onto the RBD in each simulation snapshot, and calculated the number of glycan heavy atoms that have clashes with the antibodies (NCLASH, Figure S9) . In this work, we discuss four RBD-targeting antibodies, namely B38, CR3022, H11-D4, and S309 ( Figure 6A) . The epitopes of B38 and CR3022 are irrelevant to glycans. B38 binds to the same region of RBD as ACE2 does. This epitope is fully exposed in the open-state RBD, but when the RBD is closed, B38 epitope is masked by the neighboring RBD in either open or closed state (Figure S10A) . The epitope of CR3022 is in the inner surface of RBD, and it is only accessible when all three RBD are open (Figure S10B) .
The epitope of H11-D4 is next to the epitope of B38, and it is also fully exposed when RBD is open. When RBD is closed, the glycans attached to N165 and N343 on the neighboring chain are located near this epitope (Figures 6B and S10C) . As shown in the distribution of NCLASH with H11-D4 (Figure 6C) , N343 glycan rarely interferes with the antibody, but N165 glycan has high probabilities to make severe clashes with the antibody when RBD is closed. For comparison, we also aligned a nanobody to the RBD in each simulation snapshot. Though N165 glycan still frequently makes clashes with the nanobody, the frequency of severe clashes is much lower ( Figure S11) and thus there is a chance a nanobody can bind to this epitope as shown in PDB id 6Z43.
The epitope of S309 is surrounded by four glycans. Two of them are attached to N331 and N343 on the same RBD, and the other two are attached to N122 and N165 on the neighboring NTD (Figures 6D and S10D) . The distributions of NCLASH with S309 are shown in Figure 6E . N165 and N331 glycans rarely interfere with S309 antibody in both open-and closed-state RBD. N343 glycan has minor clashes with the antibody in most snapshots, and such minor clashes are not sufficient to block antibody binding, as these clashes can be easily removed with small changes in glycan conformation and orientation. The antibody-glycan interactions can also contribute to the antibody binding, which is observed in the cryo-EM structures. In more than half of all snapshots, N122 glycan has severe clashes with the antibody when RBD is open, but it moves away from the superposition of antibody in the remaining snapshots. This suggests that S309 epitope in the open state RBD is blocked by N122 glycan in more than half of the simulation time, but it is still accessible to the antibody. As a comparison, we calculated the ASA of S309 and H11-D4 epitopes using a probe radius of 7.2 Å that is commonly used to approximate the size of hypervariable loops of antibody, and compared the epitope ASA with and without glycans ( Figure S12) . For H11-D4, we observed significant decreases of the epitope ASA in all chains except 6VSB_A and 6VSB_C, which is generally consistent with the frequency of clashes between glycans and antibodies. However, for S309, the epitope ASA decreases significantly in all chains when glycans are present. This is contradictory to the result that only N343 glycans occasionally have only slight clashes with the superimposed antibody when the RBD is close (all except 6VSB_A). In addition, the PDB structures of S trimer in complex with S309 shows that N343 glycan interacts with the antibody. In the calculation of ASA, a point on the surface is considered inaccessible even if the probe sphere has a very tiny clash with the molecule. However, the shape of antibody is not a sphere and it can have narrow shaped regions that extends deep into the pocket in the epitope. A glycan like the one attached to N343 can reduce the epitope ASA even though it may contribute to antibody binding. Therefore, in some cases, simple comparison of protein ASA with and without glycans is likely to overestimate the impact of glycan shielding.
In this work, we present multiple µs-long all-atom MD simulations of fully-glycosylated fulllength S protein in a viral membrane. Our MD simulations reveal the overall shape of S protein and its orientation on the membrane surface are determined by highly flexible stalk composed of two independent joints. Importantly, S protein models from our simulation allow us to predict possible configurations of S protein-ACE2-B 0 AT1 complex with allowable orientations and distances between two S proteins on the membrane surface. The simulation here also provides insight into how glycans influence the open/closed state change of RBD and the antibody binding to RBD epitopes. We identify glycans attached to multiple glycosylation sites that stabilize the open and/or closed states of RBD by making a high energetic barrier between the open-closed transition. The simulation of non-glycosylated systems shows that the open-state RBD becomes unstable when glycans are removed and the transition to close state occurred at the early stage of simulation. By aligning RBD-antibody complex structures to the simulation trajectories, we reveal that the impact of glycan shielding is overestimated by a simple ASA analysis. More importantly, the glycan does not only serve as shielding for immune evasion, but it can also contribute to antibody binding. Our work sheds light on the full structure and dynamics of S protein and we hope our work to be useful for development of vaccines and antiviral drugs.
In this study, the CHARMM36(m) force field was used for proteins 36 , lipids 37, 38 , and carbohydrates [39] [40] [41] . TIP3P water model 42 was utilized along with 0.15 M KCl solution. The total number of atoms is 2,343,394 (6VSB_1_1_1: 668,899 water molecules, 2,128 K + , and 1,857 Cl -); see CHARMM-GUI COVID-19 Archive (http://www.charmm-gui.org/docs/archive/covid19) for other system information. The van der Waals interactions were smoothly switched off over 10-12 Å by a force-based switching function 43 and the long-range electrostatic interactions were calculated by the particle-mesh Ewald method 44 with a mesh size of ~1 Å.
All simulations were performed using the inputs generated by CHARMM-GUI 45 and GROMACS 2018.6 46 for both equilibration and production with LINCS algorithm 47 . Temperature was maintained using a Nosé-Hoover temperature coupling method 48, 49 with a t of 1 ps, for pressure coupling (1 bar), semi-isotropic Parrinello-Rahman method 50, 51 with a p of 5 ps and compressibility of 4.5 x 10 -5 bar -1 was used. During equilibration run, NVT (constant particle number, volume, and temperature) dynamics was first applied with a 1-fs time step for 250 ps. Subsequently, the NPT (constant particle number, pressure, and temperature) ensemble was applied with a 1-fs time step (for 2 ns) and with a 2-fs time step (for 18 ns). During the equilibration, positional and dihedral restraint potentials were applied, and their force constants were gradually reduced. Production run was performed with a 4-fs time-step using the hydrogen mass repartitioning technique 52 without any restraint potential. Each system ran about 20 ns / day with 1,024 CPU cores on NURION in the Korea Institute of Science and Technology Information.