result of MD simulation, Markov state model (MSM) is applied to catch long-time kinetic information given time-limited simulation trajectories. 16, 17 Another advantage of MSM is that it can divide a large number of protein structures from simulation into subspaces based on the extracted kinetic information, and differences among those spaces can be calculated and used for comparison.
Machine learning techniques have achieved great accomplishments in chemistry and biology including material discovery, 18 structure representation 19 and computation acceleration 20 . The great contributions from machine learning mainly come from its ability to deal with large scale data and its accurate and explainable models, 21 
The root-mean-square deviation (RMSD) is used to measure the overall conformational change of each frame with regard to a reference structure in a MD simulation trajectory.
For a molecular structure represented by Cartesian coordinate vector r i (i = 1 to N ) of N atoms, the RMSD is calculated as following:
r 0 i represents the Cartesian coordinate vector of the i th atom in the reference structure.
The transformation matrix U is defined as the best fit alignment between the protein structures along trajectories with respect to the reference structure.
The root-mean-square fluctuation (RMSF) is used to measure the fluctuation of each atom in each frame with regard to a reference structure in a MD simulation trajectory.
Where T is the total frames and r i is the average position of atom i in the given trajectory.
Distances between pairs of backbone Cα were chosen as representative features for one struc- 
Random forest is a machine learning technique that can be used for classification. 25 
In a random forest model, a quantitative evaluation of the importance for each feature used for training is calculated through training process. This feature importance is calculated using Gini impurity:
where f i is the frequency of a label at a node of being picked to split the data set and C is the total number of unique labels. A random forest model is a collection of several decision tree models. The importance of node i in decision tree is calculated as:
where w j is the weighted number of samples reading node i, C i is the impurity value of node i and m is the number of child nodes. The feature importance of feature i is calculated as f i = s l n j k∈all nodes n k (6) where s is the times of node j split on feature i. Feature importance within a decision tree is further normalized by:
The feature importance in random forest is the averaged importance of feature i in all decision tree models:
where norm f i is the normalized feature importance of one single decision tree and N is the number of decision trees. 29 Feature importance of all pairwise Cα distances were calculated using the above methods.
The feature importance of an amino acid is the summation of importances of features that are related with that amino acid. The relative accumulated feature importance of each amino acid shows the distinguishability and contribution of that amino acid among all amino acids in differentiating states. 
Transition path theory (TPT) 34 
where T ik is the transition probability from state i to state k. Sequentially, the effective flux is calculated as:
where π i is the equilibrium probability of normalized transition matrix T π and q − i is the backward-committor probability calculated as q − i = 1 − q + i . However, backward flux f ji should also be considered and subtracted when calculating net flux. Therefore, the net flux
Total flux can then be calculated as: 
The flux from initial state to target state can be decomposed to individual pathways p i . Dijkstra algorithm is implemented in MSMBuilder for pathway decomposition. A set of pathways p i can be generated along f i , which provides a relative probability by:
3 Results
Simulation trajectory analysis shows dynamical activity The plot suggested that the close state is relatively stable while the partially open state is dynamically active and undergoes significant conformational changes after 1 microsecond.
However, the simulation of open state after 6 microseconds suggested a convergence in the RMSD value and a relatively stable structure, which corresponds to the detached S1 subunit from S2 fusion machinery.
RMSF results were plotted in Figure 1 PCCA was applied to map microstates onto macrostates based on the eigenfunction structure of transition probability matrix. The resulting macrostates with transition probability are shown in Figure 3 Figure 3B to show a series of transition processes.
To better understand the shift between close and open states in S1 subunit, the pair wised Cα distances of S1 were extracted as features representing the character of protein configurations.
There are 540 residues on each chain, residue ID from 27 to 676, and total of 1620 * 1619/2 = Supervised machine learning model was applied to extract the key residues that are vital during allosteric process and study the structural differences among macrostates. For each simulation trajectory, frames were saved for every 1.2 nanoseconds (ns), resulting in 8, 334 frames. Therefore, 16, 668 samples with 1, 311, 390 features were extracted from the trajectories. Each sample was labeled based on the above macrostate result. Full dataset was further split into training set(70%) and testing set(30%). After the preparation of data, random forest model was applied to distinguish the intrinsic conformational differences among macrostates. Training scores and testing scores were plotted in Figure 4A . 7 was chosen for the depth with corresponding testing accuracy of 92.18%, which indicated that the random forest model was able to catch the conformational characteristics of each macrostate only using pair-wised Cα distances. To further investigate the relationship between chain A and two other chains, the original Cα distances related with chain A was excluded as reduced features and applied to another random forest model. Training and testing results Figure 4B . The top 500 features accounted for 74.8% percent of the overall feature importance, shown in Figure S2 . The testing accuracy with reduced features reached 88.04% at depth 8.
The top five important Cα distances were listed in Table 2 . In order to identify key residues along the transition from the closed to the partially open state, the feature importance of each Cα distance was added and accumulated to the two related residues. S1
subunit structure was plotted in Figure 5A as reference. Top 20 important residues on chain B and C, with corresponding structure and accumulated structure importance under each figure, were plotted in Figure 5 (B-C). Full results of residue importance on chain B and C are shown in S3.
The significance of the partially open state of receptor-binding domain (RBD) in SARS-CoV-2 for interacting with the host cell receptor has been extensively studied. 36, 37 The opening of S1 subunit, thus exposing RBD, is necessary for engaging with ACE2 and following cleavage of S 2 site. 38 
The spike protein is essential for SARS-CoV-2 as it destabilizes the trimer structure, causing the detachment of S1 subunit and exposing the RBD domain to host cell membrane. In this study, we used publicly available simulation trajectories of S protein and studied the asymmetric dynamics nature of the trimer structure. Markov state model was applied to divide the conformational space into 8 macrostates. The representative structures along the most probable channel in transition path theory are shown to present a clear route from the close state to the partially open state. Transition matrix was calculated to determine the probability of close-open transitions, with maximum of 9.9% from the close macrostate implied a relationship between chain A and two other chains. Yet, whether this relation originates from the mutual influence among chains or the intrinsic asymmetry in biological functions needs further investigation. Overall, our study quantitatively analyzed the motion of S1 subunit in S protein with important Cα distances and residues. Possible prevention and cure searching progress could be facilitated combining with these vital findings, which provide new opportunities for potential drug research and vaccine development. Figure S2 : Accumulated feature importance of top 500 Cα distances.