key: cord-0971262-m78e36xx
authors: Li, Yao; Wang, Tong; Zhang, Juanrong; Shao, Bin; Gong, Haipeng; Wang, Yusong; He, Xinheng; Liu, Siyuan; Liu, Tie‐Yan
title: Exploring the Regulatory Function of the N‐terminal Domain of SARS‐CoV‐2 Spike Protein through Molecular Dynamics Simulation
date: 2021-09-02
journal: Adv Theory Simul
DOI: 10.1002/adts.202100152
sha: a2c4510c5f9d4b4f604fcd75599ee40cce352add
doc_id: 971262
cord_uid: m78e36xx

SARS‐CoV‐2 is what has caused the COVID‐19 pandemic. Early viral infection is mediated by the SARS‐CoV‐2 homo‐trimeric Spike (S) protein with its receptor binding domains (RBDs) in the receptor‐accessible state. Molecular dynamics simulation on the S protein with a focus on the function of its N‐terminal domains (NTDs) is performed. The study reveals that the NTD acts as a “wedge” and plays a crucial regulatory role in the conformational changes of the S protein. The complete RBD structural transition is allowed only when the neighboring NTD that typically prohibits the RBD's movements as a wedge detaches and swings away. Based on this NTD “wedge” model, it is proposed that the NTD–RBD interface should be a potential drug target.

The Coronavirus Disease 2019 (COVID- 19) , which has become a worldwide pandemic, is caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). [1] [2] [3] Similar to other coronaviruses, such as SARS-CoV and MERS-CoV, SARS-CoV-2 is an enveloped, positive-stranded RNA virus. [4, 5] Furthermore, SARS-CoV-2 has high homology with SARS-CoV, sharing the same human cell receptor, angiotensin-converting enzyme 2 (ACE2). [6, 7] After the completion of the genome sequencing of SARS-CoV-2, [1] a great amount of research attention has been put on elucidating the infection process and searching for potential therapies against COVID-19. The experimentally determined structures of the Spike (S) protein, [8, 9] Main protease, [10] RNA dependent polymerase, [11] as well as the complex formed by the S protein and ACE2 [12, 13] provide atomic details of SARS-CoV-2. Meanwhile, computational efforts have been exerted to improve the understanding of the infection mechanism via evolutional analysis. [14, 15] According to previous studies, the first step of the infection is the attachment of the S protein to ACE2. [16] The homo-trimeric S protein is located on the surface of the viral membrane, with each monomer consisting of two regions, S1 and S2, which are respectively responsible for receptor binding and membrane fusion. S1 is composed of a receptor binding domain (RBD) and an N-terminal domain (NTD). [17] Reported infection models for other coronaviruses suggest that the association of ACE2 and S trimers triggers a series of molecular events, including the shedding of S1, refolding of S2, and subsequent fusion of viral and host cell membranes. [18] [19] [20] [21] Multiple conformations of the SARS-CoV-2 S protein have been identified, with varying orientations of the RBDs and distinct ACE2 binding affinities. [8, 9] The closed state with all three RBDs in the downward orientation is ACE2-inaccessible and thus inactive, while the partially open state with one RBD flipped up is capable of receptor binding. [9] Recently, the semi-open state with two RBDs turned upward as well as the open state with all three RBDs turned upward have been detected by cryo-electron microscopy (Cryo-EM), [22] although the corresponding high-resolution structures have not been released yet.

With more atomic details of the S protein unveiled, the functional role of the RBD is generally understood by the research community. However, the contribution of the NTD, another component of S1, in the infection process is still elusive, despite some subtle clues. In MERS-CoV, a neutralizing monoclonal antibody 7D10 targeted at an NTD epitope precludes the RBD from binding to its cell receptor DPP4 by steric hindrance. [23] Interestingly, in SARS-CoV-2, a monoclonal antibody 4A8 that directly interacts with the NTD without blocking the RBD-ACE2 binding still exhibits neutralizing effects. [24] Recently, an obvious movement of the NTD was detected accompanying the RBD's conformational changes in the D614G variant of the S protein. [25] These studies suggest that the NTD is at least not completely independent of the functional conformational change of the RBD in the S protein and thus may be involved in the regulation of viral infection.

Molecular dynamics (MD) simulation is a computational tool for understanding the dynamics and function of macromolecules. Particularly, it can supplement the static spatial information provided by experimental methods such as X-ray crystallography and Cryo-EM. Early-stage explorations on S protein mainly focused on the interface between the RBD monomer and the receptor ACE2, aiming to identify hot spots on the interface, [26] estimate the binding free energies [27] as well as detect key residue pairs and important interactions. [28] [29] [30] Subsequently, more thorough simulation studies on the trimeric S protein at the atomistic level gradually revealed its dynamic transition among different conformational states. Gur et al. [31] used steered simulations to track the transition between closed and partially open states and identified new intermediate states.

Through equilibrium simulations, Casalino et al. [32] found that two glycan sites in NTD modulated the conformational dynamics of RBD, implying the potential role of NTD in S protein activation. Most recently, on one of the most powerful supercomputer clusters, Zimmerman et al. [33] captured the dramatic opening process of the apo S protein using an adaptive sampling strategy.

As mentioned above, experimental and computational studies coincidentally suggest that NTD is involved in the regulation of S protein activation. Unfortunately, its detailed contribution and mechanism are still unclear. In this work, we systematically investigated this problem by running all-atom MD simulations on the glycosylated trimeric S protein. Specifically, we simulated four systems that corresponded to the closed, partially open, semi-open, and open states of the S protein, respectively. To explore the functional role of the NTD with acceptable consumptions, each system was simulated 1 s using accelerated MD (aMD) [34] with 5 repeats to thoroughly sample the relative motion between the NTD and the RBD. Our results suggest that the upward RBD is intrinsically inclined to tilt toward to the more stable downward orientation, and that the NTD either favorably interacts with the RBD as a wedge to block this motion or occasionally detaches from the RBD to allow its free transition between the upward and downward orientations. Based on this "wedge" effect of the NTD, we virtually screened potential drugs that strengthen the NTD-RBD interface, which would hypothetically prohibit the structural transition of RBD that is required for viral infection.

To study the conformational changes of the S proteins of SARS-CoV-2, we built four simulation systems of the S trimer with different initial conformations, i.e., a closed system with three "downward" RBDs, a partially open system with one "upward" RBD, a semi-open system with two "upward" RBDs, and an open system with three "upward" RBDs. The S proteins in the closed system and the partially open system were sourced from PDB structures 6VXX [9] and 6VYB [9] , respectively, where non-native mutations in 6VYB were mutated back to the wild-type counterparts, and the missing fragments were reconstructed by SWISS-MODEL. [35] [36] [37] [38] [39] Meanwhile, the initial conformations in the semiopen and open systems were built by homology modeling using SWISS-MODEL by taking the S proteins of SARS-CoV with 2 upward RBDs (PDB ID: 6NB6) [21] and 3 upward RBDs (PDB ID: 6NB7) [21] as templates, respectively. With an explicit water solvent model, each simulation system consists of nearly one million atoms. Considering the intimidating sizes of simulation systems and the simulation steps required for large-scale conformational changes, we initiated 5 independent 1 μs aMD simulations [34] following 10 ns equilibration with a timestep of 4 fs [40] to enhance the conformational sampling in each system.

As shown in Figure 1 , remarkable conformational changes are observed in the S1 region, including both NTD and RBD, while the S2 region retains structural rigidity. Hence, we used the S2 region as a basis to define the central axis of the S protein (red dashed lines shown in Figure 1 ) and used the upper part of it as a reference for evaluating the conformational changes of the S1 region (see the Experimental Section for more details). The simulation results show that NTDs may play an important role in regulating the movement of the adjacent RBDs. In the closed system, besides the trend of one NTD moving away from the central axis observed in a trajectory, all RBDs were locked in the "downward" conformation, and the other NTDs were stable in all 5 repeats (Figure 1a ), in agreement with the experimental evidence that the closed S protein is stable and forms a condensed structure with RBDs packed tightly against the adjacent NTDs. [8, [19] [20] [21] 41] In contrast, dramatic relative movements between NTDs and RBDs were observed in the other three simulation systems. In the partially open system, the upward RBD of chain B tended to tilt downward, while the adjacent NTD of chain C moved toward the central axis to prevent this downward movement; this observation was confirmed by all five repeated simulations (Figure 1b ). This blockage of RBD by NTD was also observed in the semi-open and open systems (Figure 1c,d) . Interestingly, an unusual event was detected in these two systems: when the NTD occasionally moved away from the central axis, the adjacent RBD was allowed to tilt downwards completely (Figure 1e ,f).

We proposed novel angle and distance metrics to quantify the movements of the NTD and RBD in the simulation trajectories. Previous research used three single residues, i.e., D405-V622-V991, to define the angle for evaluating the orientation of the RBD. [42] To ensure robustness of the results, we replaced these single residues with the centroids of structurally rigid regions and defined three metrics: 1) RBD angle ( r ) that describes the upward and downward tilting poses of the RBD and any conformations in between; 2) RBD distance (d r ) that describes the distance between the RBD and the central axis; 3) NTD distance (d n ) that describes the distance between the NTD and the central axis (Figure 2a , see the Experimental Section for more details). As shown in Figure 2b-d; and Figure S1 , Supporting Information, the quantitative measurements of the relative motions between NTDs and neighboring RBDs agree with the representative conformational changes shown in Figure 1 . Specifically, the angle and the distance of the downward RBD and the distance of the corresponding NTD in the closed system exhibit no obvious changes during the whole simulation ( Figure S1 , Supporting Information). As a contrast, all motions in the simulation systems with upward RBDs are significant and can be classified into two modes. In the first mode, the NTD wedges in to block RBD movement, as detected in all three simulation systems possessing upward RBDs. As shown in the first three columns of Figure 2b -d, both angle and distance of the upward RBD decrease, but this tendency is later blocked by the quick movement of the NTD toward the central axis (i.e., a marked reduction of the NTD distance). In the second mode, NTD swings away, and RBD accomplishes its complete downward tilt, as observed in both the semi-open and open systems. As shown in the last two columns of Figure 2d , the distance of the NTD increases in the early stage of simulation, suggesting that it moves away from the central axis, releasing the steric hindrance for RBD movement. Accordingly, the highlighted upward RBDs in both semi-open system (the penultimate column) and open system (the last column) eventually approach much smaller distance and angle values compared with the blocked RBDs (the first three columns), which indicates that these highlighted RBDs tilt downward completely. To further validate these two modes, we estimated the distributions of RBD distance and angles based on the last 200 ns trajectories from all five repeats for each simulation system, while the control distributions of the downward RBDs (chain A in both partially open and semi-open systems) are calculated from the whole trajectories. Consistently, the metrics of the upward and downward RBDs are separated by a marked gap in the RBD blocked mode (see the first three columns of Figure 2e , f), while huge overlaps occur in the RBD downward tilting mode (see the last two columns of Figure 2e ,f). Moreover, one-way ANOVA analysis confirms that significant difference is present between experimental groups and control groups in the RBD blocked mode but is absent in the RBD downward tilting mode (Table S1, Supporting Information). These results imply that the NTD may act as a wedge to regulate the functional movement of the RBD.

The complete downward tilt of the RBD was not caught in the partially open system, possibly because the NTD kept wedging in during the simulations. To validate this hypothesis, we artificially forced the NTD of chain C to the swinging-out posture (using targeted MD (tMD) on a target position inferred from a reference structure with large NTD distance in the semi-open system) and then investigated the movement of the upward RBD in this system through 5 redundant 400 ns aMD simulations (Figure S2a, Supporting Information). Compared with previous simulations that started from the crystal structure, the distributions of RBD angles and distances exhibited remarkable left-shifts in the NTD-forced-out trajectories ( Figure S2b , Supporting Information), indicating that the upward RBD tilts downward to a significantly larger extent when its partner NTD is absent. We also conducted a time-structure based Independent Components Analysis (tICA) [43] to identify slow motions with high time autocorrelation on the first 200 ns of these NTD-forced-out simulations. As expected, the first tIC presents the highest correlation with the difference vector between the representative structures of the upward and downward RBDs ( Figure S2c , Supporting Information). Hence, the downward tilting movement of the RBD is the main slow motion when the NTD is forced away.

Markov state models (MSM) [44] are a powerful tool used to extract structural information from multiple simulations and has been exploited to track dynamic progress and elucidate slow motions in MD simulations. [44] In our work, considering the relative structural rigidity of the S2 region, we extracted Cartesian coordinates of C atoms of S1 to build MSM and identified intermediate conformational states in the two modes described previously. For the NTD wedging-in mode, all five repeated trajectories of the partially open system were used in the principal component analysis (PCA) to reduce the dimensionality of original coordinates. As shown in Figure S3a (Supporting Information), the top 4 principal components (PCs) were picked to represent the S1 regions of three chains since they explained 69.38% of the total variance. A lag-time of 1.0 ns was set for the time-discretization of the process and a number of 100 microstates was chosen based on the implied timescale analysis. [44, 45] These microstates were built using minibatch K-Means clustering [46] and were subsequently aggregated into a number of macrostates (or metastable states), for which pairwise transition probabilities were calculated. [44] The detailed workflow of MSM construction is illustrated in the Experimental Section.

As shown in Figure 3a ,b, all five trajectories in the partially open system were finally lumped into six macrostates and the structure closest to the centroid of each macrostate cluster in the space of the top 4 PCs was chosen as the representative structure. State 0 is the initial state with one RBD turned upward, whereas states 4 and 5 are the final states with this RBD blocked by the adjacent NTD. States 1, 2, and 3 are intermediate states between the initial and final states. With the NTD quickly moving toward the central axis in the early stage of simulation, the initial state 0 can directly transform into the final state 4 with a high probability. Meanwhile, in the intermediate state 2, the NTD initially moves We also built an MSM for the mode characterized by the NTD moving out and the RBD tilting downward from two trajectories of the semi-open system. Following the same protocol, we chose the top 5 PCs, set the microstate number as 100 and set the lagtime as 1.0 ns ( Figure S3b, Supporting Information) . As shown in Figure 3c , this mode is represented by a four-state transition network, including two transition paths between the RBD"s initial upward state (state 0) and the RBD"s final downward state (state 3), mediated by two intermediate states (states 1 and 2). In this mode, the NTD first moves away from the central axis, eliciting the upward RBD to fall (state 1) and/or moves inward to stabilize the already downward RBD (state 2). As shown in Figure 3d , the final state (state 3) exhibits the smallest values in all three metrics, implying that the upward RBD has tilted down completely to form a highly packed structure in this state. Notably, the transition probabilities here are obtained from aMD simulations with boosting potentials [47] and thus may not correlate well with the "real" probabilities estimated from an equilibrium distribution. [44] 

As demonstrated above, each state in MSM is crucial to the transition from the initial to the final state, and the whole network clearly illustrates how the overall conformational change happens step by step. In order to figure out the molecular mechanism of these conformational changes, we performed molecular mechanics generalized born surface area (MM/GBSA) [55] on each microstate with one protomer as the ligand and the remaining two as the receptor and decomposed the total binding free energy to individual residue pairs. By this means, the residue pairs that make significant contribution to the physical interactions between the NTDs and their adjacent RBD partners can be identified.

The identified residue pairs are mainly located at the NTD-RBD interface but can be clustered spatially into several clumps (Tables S2 and S3 , Supporting Information). The interaction patterns of these clumps vary substantially during the structural in (a). Final states have smaller NTD distance (d n ) but comparable RBD angle ( r ) and distance (d r ) when compared with the initial state. c) MSM in the semi-open system illustrates the mode characterized by NTD moving out and RBD tilting down. Four macrostates were extracted, among which state 0 is the initial state with an upward RBD, states 1 and 2 are intermediate states, while state 3 is the final state with the downward conformation. d) A quantitative evaluation on all macrostates in (c). Both RBD angle ( r ) and distance (d r ) in the final state are remarkably smaller than those in the initial state.

transition. Specifically, for the mode characterized by NTD wedging in and RBD being blocked in the partially open system (Figure 4a ; and Table S2 , Supporting Information), the red clump, including H519 and P230, forms in the intermediate states 2 and 3 and remains in the two final states 4 and 5. Meanwhile, the yellow clump, including P521 and Y170, forms and breaks intermittently and disappears in the final state 5. We suspect that the interactions between these residue pairs are more relevant in mediating the relative motions between the NTD and the RBD. On the contrary, the more conservative clumps that exist among most macrostates, e.g., the green clump including R357 and N165, as well as the purple clump including N360 and T167, are likely to play more basic roles such as structural stabilization. For the mode characterized by NTD moving out and RBD tilting downward in the semi-open system (Figure 4b ; and Table S3 , Supporting Information), the identified residue pairs scarcely overlap among macrostates, possibly due to the more drastic relative motions between the NTD and the RBD. However, it is noteworthy that the initial state lacks strongly interacting residue pairs if −1.5 kcal mol −1 is taken as the energy threshold. Interactions gradually form when the RBD tilts downward, reaching the final stable state with the NTD-RBD interface composed of three clumps.

We further performed evolutional analysis on key residue pairs identified in the partially open system: the residues of the RBD of chain B from 355 to 360 and from 519 to 521, and the residues of the NTD of chain C from 159 to 170 and from 230 to 233. We first compared such residue pairs with the corresponding residues in SARS-CoV and MERS-CoV. As shown in Figure 4c , most of these key residues are conserved in SARS-CoV, while only parts of them are conserved among three kinds of sequences (conservativeness reflected by the darkness of blue colors), which is consistent with the overall sequence similarity among the three coronaviruses. When evaluated on the samples of COVID-19 patients derived from Global Initiative on Sharing All Influenza Data (GISAID), [51, 52] all these residues were highly conserved regardless of region or race (Figure 4d ). These results indicate that the NTD-regulated structural transition may be generalizable for all current SARS-CoV-2 samples, which implies that potential therapy based on this mechanism may be applicable to all COVID-19 cases.

Given all the consistent observations and comprehensive analysis demonstrated above, we propose the "wedge" model to illustrate the crucial regulatory role of the NTD in the functional conformational change of the S protein. As shown in Figure 5 , the NTD acts as a wedge in two aspects: wedging in for RBD blockage and moving out to allow RBD tilting. In general, the receptor-accessible upward RBD is inclined to tilt downward to form the more compact structure. The free reorientation of the RBD, however, will be blocked when the neighboring NTD moves toward the central axis and forms favorable interactions with the RBD. Occasionally, the NTD may swing away and detach from the RBD, sterically allowing the free downward movement of the RBD. Notably, albeit not directly observed in our simulation, the upward tilting of RBD Table S2 , Supporting Information, for more details). b) Key residue pairs of each macrostate in MSM for the mode characterized by NTD moving out and RBD tilting downward (see Table S3 , Supporting Information, for more details). c) Evolutional analysis of key residue pairs in (a) among SARS-CoV-2 (PDB ID:6VYB), [9] SARS-CoV (GenBank: AY278488.2), [48] and MERS-CoV (GenBank: NC_01 9843.3). [49] Results were visualized by Discovery Studio Visualizer. [50] d) Evolutional analysis of key residue pairs in (a) among all samples of COVID-19 patients from GISAID. [51, 52] Data were collected on September 28th, 2020. Results were visualized by WebLogo. [53, 54] Hydrophilic, neutral and hydrophobic residues are colored in blue, green, and black, respectively. has been visualized in previous works. [33] Therefore, we speculate that NTD may play as a wedge in this opposite motion of RBD. According to this model, the NTD may be a potential drug target for COVID-19 therapy. Particularly, drugs targeting the NTD-RBD interface is supposed to augment the NTD "wedge" effect and lock the S protein into an individual conformational state (complete "up" or complete "down" states), which would inhibit the functional structural transition required for viral infection.

Given the "wedge" model of the NTD, the partially open state was employed to recruit drug compounds that enhance interactions on the NTD-RBD interface. The NTD of chain C from the last frame in the equilibrium simulation of the partially open system was chosen as the receptor. We used AutoDock Vina [56, 57] to perform virtual screening against the sub dataset of "in-trials" in ZINC15, [58] [59] [60] including all approved drugs and investigational compounds in clinical trials. 9800 chemical compounds were evaluated, and top compounds with high binding affinities are shown in Table S4 (Supporting Information). The compound ZINC000261494659, a Paclitaxel analogue, has the highest binding affinity of −13.8 kcal mol −1 and binds at the interface between NTD of chain C and the "up" RBD of chain B. The drug compound strengthens the NTD-RBD binding by forming new interactions at the interface, including hydrogen bonds (e.g., with N360 of the RBD of chain B and V171 of the NTD of chain C), hydrophobic interactions (e.g., with P330 of the RBD and V159 of the NTD) and van der Waals (Figure 6a) . ZINC000256109538 can also bind at the NTD-RBD interface with high affinity but in a different pose (Figure 6b ), forming interactions with T108, Q115, and T167 of the NTD as well as I233, R357, and N394 of the RBD. In addition, we also performed virtual screening against the database of TCMSP, [61] which consists of 14 249 chemical compounds extracted from Chinese herbs. Compounds with high binding affinities to the NTD-RBD interface were selected using a threshold of −9.5 kcal mol −1 , which were also found to coincide with the ingredients of the prescriptions recommended by the National Health Commission of China for the treatment Figure 5 . The schematic diagram for the NTD "wedge" effect in controlling conformational changes of the S protein. The leftmost structure is the simplified initial conformation, where only one NTD (orange) and its neighboring RBD (cyan) are showed in this sketch. When the NTD wedges in and favorably interacts with the RBD, the free tilting of the RBD is blocked and the S protein is locked into individual conformation states. Occasionally, when the NTD swings out and detaches from the RBD, the RBD's reorientation is restored, thus enabling the interstate transition of the S protein.

Hypothetically, the opposite motion of the RBD tilting upward may occur during the activation of the S protein, which is illustrated by dotted arrows.

of COVID-19 (Tables S5-S13, Supporting Information). For example, in the prescription of "Qing Fei Pai Du Tang," 18 constituents from 8 herbs were detected with high binding affinities (Figure 6c ). These results imply that the NTD-RBD interface may present the potential binding target for anti-SARS-CoV-2 drugs, which awaits computational and experimental validation in the future.

In our work, we found that the NTD acts as a wedge to regulate the movement of the adjacent upward RBD: holding it to stand up (the "wedging in" mode) or allowing it to tilt down (the "moving out" mode). Both modes were observed in the simulation systems with upward RBDs, and several metastable states were identified by MSM. This NTD "wedge" model should be further validated by biochemistry and molecular biology experiments. Notably, previous research has proven that the closed and the partially open structures coexist in vivo [62] and the full-atom MD simulations on trimeric S protein [32, 33] including ours revealed the dynamic transitions between the two structure states. Among the three simulation studies, Zimmerman's simulation reached an unprecedented time scale of 0.1 s with the powerful distributed computational system Folding@home and thus directly observed bidirectional conformational changes of RBDs. Although their results shed light on the infection mechanism and the design of antivirals, the relative motions and potential regulatory functions of NTDs were not systematically investigated. In the works of Casalino et al. and ours, MD simulations were only performed to microseconds due to the magnitude of the system size (≈1 000 000 atoms). Casalino et al. laterally confirmed the involvement of NTD in S protein activation by identifying the two glycosylation sites N165 and N234 in this domain. We moved one step forward and proposed a full mechanistic model to show the regulatory role of NTD in the RBD conformational change. Therefore, the "wedge" model proposed in this work complements previous computational and experimental studies on S protein and deepens our understanding of NTD-RBD interactions during viral infection.

Furthermore, the semi-open and open states are more flexible with more than one upward RBD, and thus their simulations exhibit high diversities. As shown in Figures S4 and S5 (Supporting Information), besides the two kinds of relative movements between the NTD and the RBD, we also observed that the upward RBDs formed interactions with each other. In the semiopen system, two upward RBDs tilted down, and the three downward RBDs in the final structure formed stable mutual interactions ( Figure S4 , Supporting Information). Similar behavior was also observed in the open system ( Figure S5 , Supporting Information), where all three upward RBDs tilted downward completely to form a compact packing. Occasionally, two RBDs tilted downward, while the remaining upward RBD only formed close interactions with one of the other two RBDs ( Figure S5b , Supporting Information). Nevertheless, the complicated interactions with multiple upward RBDs and their potential influence on viral infection still await further investigation.

Our main observations were obtained from aMD simulations, which disrupted the Boltzmann distribution by applying energy boosts. Hence, unperturbed distributions of conformational states need to be recovered from the microstate transition analysis in MSM, which are essential for the evaluation of kinetic properties, such as transition probabilities, and thermodynamic properties, such as the potential of mean force. In addition, although we have screened out dozens of promising compounds with high binding affinities to the NTD-RBD interface from the ZINC15 and TCMSP databases, more stringent calculations and experimental validations should be conducted to evaluate the identified drug candidates in the future.

Simulation System Design and Configuration: The initial structures of the closed system and the partially open system were chosen from the PDB structures of 6VXX [9] and 6VYB [9] , respectively. These structures are trimer glycoproteins with 1281 residues with N-glycosylation in each chain. In the research work, non-native proline mutations were mutated back to the wild type lysine and valine, respectively. Since several fragments were not resolved in these Cryo-EM structures, SWISS-MODEL [35] [36] [37] [38] [39] was employed to replenish the missing atoms. The GLYCAM-Web online server (http://glycam.org) was used to add glycans (biantennary LacNAc Nglycans) on 13 aspartic acids and changed the residue name of "ASN" to "NLN" for residues modified by glycans in the PDB files. In the semi-open and the open simulation systems, the initial structures of the S proteins were generated by SWISS-MODEL with SARS-CoV [21] structures of 2 upward RBDs (PDB ID: 6NB6) and 3 upward RBDs (PDB ID: 6NB7) taken as templates. Glycans were also added at the same positions.

Each simulation system was solvated by a cubic TIP3P water box. Considering that the protein mainly stretched in the x and y directions, the boundaries of the water box in each simulation system were extended more in the x and y dimensions than in the z dimension (Table S14 , Supporting Information) to ensure that the protein stayed in the periodic cell during the simulation process. 3 Na + ions were added to neutralize the simulation systems. The force fields "leaprc.GLYCAM_06j-1," [63] "leaprc.protein.ff14SB," [64] and "leaprc.water.tip3p" [65] in Amber20 were employed to parameterize glycans, the protein and water molecules (as well as ions), respectively. Topology files suffixed with ".prmtop" and coordinate files suffixed with ".inpcrd" were generated by running the program "tleap" in Amber20. The program "cpptraj" was executed to fix the orders of all atoms.

In order to accelerate simulation speed, a timestep of 4 fs via hydrogen mass repartitioning, which repartitioned the mass of hydrogen atoms and their attached heavy atoms was used. [40] SHAKE was used to remove the stretching degrees of freedom of hydrogen-involved bonds.

Simulation Process: We first performed 5000 steps of energy minimization to relax steric clashes, including 2500 steps of steepest descent minimization and 2500 steps of conjugate gradient minimization. The temperature of each simulation system was set to 300 K, and the preequilibration was run in the NVT ensemble for 300 ps. A radius cutoff of 10 Å was used for nonbonded interactions. Harmonic positional constraints of 10 kcal mol −1 Å −2 were applied during the pre-equilibration to gradually release improper forces in the system. Before the equilibrium simulation, the solvent was first optimized, followed by optimization of the protein side chains, and then restraints were gradually removed to allow the system to approach a low-energy state. Then 10 ns equilibrium simulation at 300 K for each system was performed. The minimization, pre-equilibration and equilibrium simulation were performed with the binary file "pmemd.cuda" in Amber20.

By adding bias potential to the original potential, aMD [34] can promote the protein to overcome local energy barriers, which would enhance the efficiency of conformational sampling. Energy boosting was applied to both the whole potential and the torsion term (iamd = 3). The aMD parameters adopted in the simulations included the average energy thresholds and inverse strength boost factors for the total potential and the dihedral potential. When estimating these parameters, the total number of atoms is scaled by 0. 16 in the partially open, semi-open, and open systems but by 0.20 in the closed system, since the last system was more stable. 1000 ns aMD was run for each simulation system with 5 repeats.

tMD was performed to detect the mode characterized by NTD moving out and RBD tilting downward in the partially open system. The distance was calculated between the NTD and the central axis of the S protein in the trajectory of the semi-open system and chose the frame with the largest distance as the target structure. The initial structure was chosen as the last frame of 100 ns equilibrium simulation on the partially open system. tMD [66] applies an additional term to the energy function based on the mass-weighted root mean square deviation (RMSD) of a set of atoms in the current structure compared to the reference structure 2 (1)

where TGTMDFRC stands for the force constant, NATTGTRMS stands for the number of C atoms of NTD, CURRMSD stands for the RMSD of the current structure compared with the target structure and TGTRMSD stands for expected RMSD value, which is equal to 0 in the work.

To conduct a quantitative analysis of relative motions between the NTD and the RBD, four stable regions were first defined by calculating the residue-wise positional fluctuations (also referred to as root-mean-square fluctuations or RMSF)

where x i stands for the Cartesian coordinates of the C atom of a residue in each frame of the whole trajectory and <x i > stands for the mean value along the simulation trajectory. By selecting residues with small RMSFs, the four stable regions were defined as follows: 1) the RBD core region including resides 334-527, with residues 454-491 excluded; 2) the rotation axis of RBD including residues 541-586; 3) the upper part of the central axis including residues 976-996; and 4) the NTD core region including resides 27-303. The centroids of each region were calculated and three evaluation metrics were defined: 1) RBD angle ( r ) is the scalar angle between centroids of the RBD core region, the rotation axis and the central axis; 2) RBD distance (d r ) describes the distance between the RBD core region and the central axis; 3) NTD distance (d n ) describes the distance between the NTD core region and the central axis. tICA [43] was used to evaluate motions of the upward RBD in the partially open system when the NTD was forced out by tMD. Instead of using structural similarity as a proxy for kinetic similarity, tICA sought to encode kinetic similarity in the states explicitly by choosing the states along the slowest degrees of freedom. Given a multidimensional dataset, tICA is able to find the linear combinations (tIC) of the input coordinates that maximize the autocorrelation function of the corresponding projection, while restraining this linear combination to be uncorrelated to the previous ones. In practice, tICA was performed on the first 200 ns of the aMD trajectory with a lag time of 2 ns.

Markov State Model: By using MSMBuilder v3.8.0, [46, 67] the MSM was built for the partially open system with the following procedures: 1) To remove the degrees of freedom of rotation and translation, all frames were aligned of all trajectories to the first frame of trajectory 1. 2) As the study focused on relative motions of the NTD and its partner RBD of the S1 region, the Cartesian coordinates of C atoms of S1 regions were extracted from five trajectories as input data for MSM. 3) PCA was performed to reduce the dimensions of the input data and projected the input data on the top 4 PCs. 4) The choice of lag-time is important since it allows the discrete microstates to be Markovian and ensures kinetic accessibility. [46, 67] According to the implied timescale analysis, [45] if the implied timescale no longer increases at any lag-time longer than , then is chosen as a proper lag-time to resolve the process. The lag time was chosen based on MSMs built with a series of microstate numbers including 50, 100, 250, 500, 1000, and 2000. At the beginning, the implied timescale rose gradually, with the 0.2 ns stepwise increment of lag-time increasing. When the lag time reached 1.0 ns, the implied timescales tended to be relatively constant, and no significant increase was detected with any longer lag times in all MSMs with different numbers of microstates. Therefore, 1.0 ns as the lag time was chosen. 5) The generalized matrix Rayleigh quotient (GMRQ) as a metric was used to choose a proper number of microstates. GMRQ is a criterion that evaluates how well the MSM eigenvectors generated on the training dataset explains kinetics variance in the test dataset. [68] MSM was built with a lag time of 1.0 ns and with varying numbers of microstates (50, 100, 250, 500, 1000, and 2000), and finally chose 100 microstates because of the highest value of GMRQ in fivefold cross validation. With all hyperparameters fixed, all subsequent MSMs were constructed with a lag time of 1.0 ns and microstate number of 100, using minibatch K-Means. [46] 6) 100 microstates were lumped into 6 macrostates by perron-cluster cluster analysis. [69] 7) For the partially open system, MSM was built with 6 macrostates and transition probabilities were calculated between each pair of macrostates. 8) The structures closest to the center of each macrostate as the representative structures were selected.

MSM was also built for the mode characterized by the NTD moving out and the RBD tilting downward in the semi-open system with the same procedures. 100 microstates were clustered in a 5D PC space and a lag time of 1.0 ns was chosen for this model. 4 macrostates were selected in the final model.

Free Energy Decomposition: MM/GBSA [55] was used for free energy decomposition and detected key residue pairs on the NTD-RBD interface. For each macrostate, 100 structures that had the smallest RMSD with the representative structure and discarded structures with RMSD larger than 4 Å were first picked up. MM/GBSA [55] was performed in Amber20 on each microstate for the selected structures. The program "Parmed" was used to split the topology file of the original system into two parts, i.e., a chain with the selected NTD acting as the ligand and the remaining two chains acting as the receptor. The total free energy on pairs of residues was decomposed and residue pairs were identified in the NTD of the ligand and in the corresponding RBD of the receptor according to decomposed binding free energies. After ranking the pairwise decomposed binding free energy, it is found that most of the residue pairs with high bind affinities were located on the NTD-RBD interface.

Virtual Screening by Docking: Virtual screening was performed by AutoDock Vina [56, 57] with the following procedures:

1) The receptor file for docking was generated. The structure of the last frame of the equilibrium simulation in the partially open system was selected as the receptor. Charges were added for each residue. 2) The virtual screening was performed against two compound databases. 9800 chemical compounds were culled from the category of "in-trails" of ZINC15 [58] [59] [60] and 14249 compounds from the database of TCMSP. [61] The compounds culled from ZINC15 were approved drugs and investigational compounds in clinical trials with the "sdf" format. The compounds culled from TCMSP were chemical compounds extracted from 502 kinds of Chinese herbs with the "mol2" format. Chemical compounds were converted into the PDB format by OpenBabel [70] and were further optimized using the Generalized Amber Force Field (gaff). Charges were added and rotatable bonds were assigned for each compound.

3) The grid box was then designed and the configuration file was set. The NTD-RBD interface of the partially open system as the target for docking was chosen. The Cartesian coordinates of the centroid of the NTD were calculated and this NTD was placed into a rectangular box with each dimension extended by 5 Å. 4) AutoDock Vina [56, 57] was run to screen compounds in two databases with the same configuration. The top 20 compounds with the highest binding affinities in ZINC15 were selected, and the detailed interactions were visualized using Discovery Studio Visualizer. For compounds in TCMSP, all compounds were selected with binding free energy lower than −9.5 kal mol −1 . Then the selected compounds that coexisted in 10 prescriptions recommended by National Health Commission of China for the treatment of COVID-19 were analyzed.

Supporting Information is available from the Wiley Online Library or from the author.

Proc. Natl. Acad. Sci

Proc. Natl. Acad. Sci

Discovery Studio Modeling Environment, Release

Global Challenges

Y.L., T.W., and J.Z. are co-first authors and contributed equally to this work. The authors thank Zhuoxuan Liu at the Department of Electrical Engineering of Tsinghua University for his support for drawing the schematic diagram for the NTD "wedge" effect. This work was partially founded by the National Natural Science Foundation of China (Nos. 31670723, 81861138009, and 91746119) to H.G. and by the Beijing Advanced Innovation Center for Structural Biology to H.G.