The ongoing outbreak of the novel coronavirus [2] [3] [4] , which causes fever, severe respiratory illness, and pneumonia, poses a major public health and governance challenges. The emerging pathogen has been characterized as a new member of the betacoronavirus genus (SARS-CoV-2) 5,6 , closely related to several bat coronaviruses and to severe acute respiratory syndrome coronavirus (SARS-CoV) 7, 8 . Compared with SARS-COV, SARS-CoV-2 appears to be more readily transmitted from human to human, spreading to multiple continents and leading to the World Health Organization (WHO)'s declaration of a pandemic on March 11, 2020 9 . According to the WHO, as of June 9, 2020, there had been >7,000,000 confirmed cases globally, leading to >400,000 deaths.
In the initial step of the infection, the coronavirus enters the host cell by binding to a cellular receptor and fusing the viral membrane with the target cell membrane 10, 11 . For SARS-CoV-2, this process is mediated by spike glycoprotein which is a homo-trimer 12 . In the process of viruses fusing to host cells, the spike protein undergoes structural rearrangement and transits from a metastable prefusion conformational state to a highly stable post-fusion conformational state 13, 14 . The spike protein comprises of two functional units, S1 and S2 subunits; when fused to the host cell, the two subunits are cleaved. The S1 subunit is responsible for binding to the angiotensin-converting enzyme 2 (ACE2) [15] [16] [17] receptor on the host cell membrane and it contains the N-terminal domain (NTD), the receptor-binding domain (RBD) and the C-terminal domain (CTD). NTD in the S1 subunit assists recognize sugar receptors. RBD in the S1 subunit is critical for the binding of coronavirus to the ACE-2 receptor [18] [19] [20] [21] . CTD in the S1 subunit could recognize other receptors 22 . The binding of RBD to ACE2 facilitates the cleavage of the spike protein and promotes the dissociation of the S1 subunit from the S2 subunit 23 . S2 contains two heptad repeats (HR1 and HR2), a fusion peptide, and a protease cleavage site (S2'). The dissociation of S1 induces S2 to undergo a dramatic structural change to fuse the host and viral membranes. Thus, the spike protein serves as a target for development of antibodies, entry inhibitors and vaccines 24 . Coronavirus transits from a metastable prefusion state to a highly stable post-fusion state as part of the spike protein's role in membrane fusion. The instability of the prefusion state presents a significant challenge for the production of protein antigens for antigenic presentation of the prefusion antibody epitopes that are most likely to lead to neutralizing responses. Thus, since the prefusion spike protein exists in a thermodynamically metastable state 1 , a stabilized mutant conformation is critical for the development of vaccines and drugs.
Computational mutagenesis is an effective approach to finding mutations that are able to stabilize proteins. We have previously developed a protein design platform, Eris 25, 26 , which utilizes a physical force field 27 for modeling inter-atomic interactions, as well as fast side-chain packing and backbone relaxation algorithms to enable efficient and transferrable protein molecular design. Originally, Eris has been validated on 595 mutants from five proteins, corroborating the unbiased force field, side-chain packing and backbone relaxation algorithms. In many later studies, Eris has been validated through prediction of thermodynamically stabilizing or destabilizing mutations [28] [29] [30] [31] [32] [33] , and direct protein design efforts [34] [35] [36] [37] .
In this work, we propose a pipeline to automatically stabilize spike proteins through computational mutagenesis. Within the pipeline, we first analyze the conservation score and solvent accessible surface area (SASA) of residues in the protein. We then perform discrete molecular dynamics (DMD) [38] [39] [40] [41] simulations to calculate the root mean square fluctuation (RMSF) of residues to analyze their flexibility. Based on this information, we select appropriate residues (2 < Conservation Score < 5; SASA > 0.4; RMSF > 3.5 Å; see Discussion) as mutation sites. We subject the selected residues to computational redesign using Eris to find the stabilizing mutations by calculating the change in free energy ∆∆ = ∆ − ∆ , where ∆ and ∆ are the free energies of the mutant protein and wild type proteins correspondingly. We utilize this pipeline to identify stabilization mutants of the spike protein. Next, we describe our methods in detail and provide a list of stabilizing mutations for spike protein. 
We propose a pipeline to automatically find stabilized mutants of proteins ( Figure   1 ). The pipeline can be divided into two stages. In the first stage, users designate the protein of interest, and then the pipeline will analyze the 3D structure of the protein by using different metrics. In the second stage, users designate the mutation sites according to the analysis of the 3D structure of the protein, and finally the pipeline will determine the stabilizing mutations of the protein. Users can either upload the 3D structure of the protein or input the PDB ID of the protein to designate the protein of interest. In the first stage, the first step is to remodel the 3D structure of the protein of interest to complete the missing atoms and residues.
We integrate MODELLER into our pipeline to remodel proteins. Next, the pipeline utilizes ConSurf 48 to calculate the conservation score of each residue in the protein of interest. The conservation score indicates the importance of the residue in maintaining protein structure and/or function. Subsequently, the pipeline utilizes DMD to analyze the flexibility of each residue in the spike protein through RMSF.
The technique has already been used to efficiently study the protein folding thermodynamics and protein oligomerization and allows for a good equilibration of the structures. Then, the pipeline will calculate SASA of residues in the protein.
In the second stage, users designate the mutation sites according to the conservation score, RMSF, and SASA. A high conservation score (≥ 7) indicates the residue may play important roles in the function or the stability of the structure of the protein; residues of high RMSF (> 3.5 Å) are likely the culprit to undermine the stability of the structure of the protein, hence we select residues that have a low conservation score or high RMSF; residues with SASA < 0.5 are considered buried and residues with SASA ≥ 0.5 are considered exposed to solvent. After the designation of the mutation sites, the pipeline utilizes Eris to determine the changes in free energies of the mutants. For each residue in the mutation sites, 
We utilize MODELLER 42 to re-model the 3D structure of the spike protein by using a structure deposited to Protein DataBank (PDB), PDBID:6VSB 12 , the cryo-EM structure of the prefusion state of the spike protein, as the template structure to complete the missing atoms and residues ( Figure 2A&B ). Next, we use Chiron 44 All following computational mutagenesis study are performed using the structure with the 3 RBDs in down conformation.
To validate the ability of the pipeline to identify stabilization mutants, we use the pipeline to calculate the free energy changes of several known prefusion spike protein stabilization mutants. The 2P mutation strategy (K986P and V987P) has been proved effective for the stabilization of spike protein of SARS-COV-2 and other betacoronavirus 12,20,50 . Hsieh and coworkers 51 
The spike protein is mainly composed of S1 and S2 subunits ( Figure S1 ). We select residues for mutation from the NTD and RBD domains in S1, and we also select residues from the HR1 (heptad repeat 1) and CH (central helix) domains in S2. We don't select residues from HR2 domain in S2 because the structure of the HR2 domain has not been solved in the cryo-EM structure 6VSB. At the outset, we calculate the conservation score ( Figure 3A&D ) of all residues by using ConSurf.
Based on the conservation score, most residues in HR1/CH are conservative, while residues in NTD and RBD are prone to mutation in evolution. Next, we use Pymol 47 to calculate SASA of all residues in the spike protein ( Figure 3B&E) . SASA indicates the level of residues exposed to the solvent in a protein and usually most of the functional residues are located on the protein structure's surface 52 . All four domains have both low SASA residues and high SASA residues. Then, we perform 1,000,000 steps DMD simulation for the spike protein to calculate the RMSF of each residue ( Figure 3C&F ). The residues in HR1/CH have extremely low RMSF, while the residues in NTD and RBD domains have moderate to high RMSF.
We select residues that have different conservation scores in these four domains for mutagenesis. In the NTD, we select 5 residues (Table S1 ) with the conservation score ranging from 1 (highly variable) to 9 (highly conservative). The SASA of these residues range from 0.01 (buried) to 0.75 (exposed). The RMSF range from 1.61 Å (frozen) to 4.88 Å (flexible). Likewise, in the other three domains (Table S2- 3), we select 5 to 7 residues, respectively. These residues also have diverse conservation scores, SASA, and RMSF. Of note, to avoid affecting the function of the spike protein, these residues are all not chosen from the functional sites of the spike protein, such as the ACE2 binding site in RBD.
We utilize Eris to calculate the free energy changes of mutants relative to the wild type ( Figure 4A , Figure S2&3 , and Table S4-6). In the NTD, 5 residues, E169, K113, I203, R246, and L270 are selected for mutagenesis. Among them, the free energy changes of nearly all mutations of residues I203, R246, and L270 are positive, indicating that they are destabilizing the structure. In contrast, most mutations on residues E169 and K113 have negative free energy changes, suggesting that they are stabilizing the structure. The mutant that has the most negative free energy change is R246C. However, cysteine is prone to forming disulfide bond with other cysteine, which may affect the correct folding of the protein structure, so we recommend K113M to be a better choice as the stabilization mutant.
In the RBD, 5 residues (A411, T415, Y505, N439, and D428) (Table S2) are selected for mutagenesis. The free energy changes calculated by Eris ( Figure S2 , In the HR1/CH domain, 7 residues (L948, I1018, A1026, Y1007, S1003, T961, and V976) are selected for mutagenesis as shown in Table S3 . In stark contrast to NTD and RBD, most residues have extremely high free energy changes (> 30 kcal/mol),
suggesting that these residues are not very good choices for mutagenesis. The high free energy changes also implicate that they may play important roles in stabilizing the structure so that they are irreplaceable to some extent. This finding is also in concert with the high conservation scores of these residues. That said, we can still find stabilization mutants for these residues, such as T961A and S1003M ( Figure S3 ).
Compared to experimental mutagenesis, such as random mutagenesis 53, 54 and site-directed mutagenesis 55, 56 , computational mutagenesis 25 is an efficient alternative that lays the foundation of large-scale mutation screening. However, performing computational mutation screening for all residues in the spike protein trimeric structure, which consists of 1288 residues in each monomer, is still timeconsuming and inefficient, so we seek to select critical residues to perform mutagenesis. In this work, to interrogate how to select residues as mutation sites, we select residues that have different conservation scores, SASA, and RMSF. We find that the probability of stabilizing mutations for specific residues is correlated with the conservation score, SASA, and RMSF of these residues ( Figure 4B -D).
We calculate the average and the minimum free energy change of all 19 mutations of each residue. We find that the average free energy change of mutations of residues with either high (>5) or low (<2) conservation score are typically larger than 0. Only mutations in residues that have moderate (2~5) conservation score have negative average free energy change, thus indicating a possibility to find stabilizing mutations. We posit that conservative residues are typically playing critical roles in structural stability or protein functioning 57 , making them a likely target for finding stabilizing mutations. However, if the conservation score of the residue is too high, the residue may be irreplaceable and any mutation will destabilize the structure. On the other hand, residues with low conservation scores may be less critical to the structural stability, reducing the chances for finding stabilizing mutations at these positions. Thus, we select residues that have moderate conservation score (2~5) for mutagenesis. Similarly, we find that residues with large SASA or large RMSF have more stabilization mutants than residues with low SASA or low RMSF, respectively. Overall, we select residues that have moderate conservation score (2~5), high SASA (>0.4), and high RMSF (>3.5 Å) for mutagenesis. In addition, although we only select residues in NTD, RBD, and HR1/CH domains to perform mutagenesis, residues in other regions can also be used as mutation sites. For example, the known 2P mutation strategy (K986P and V987P) has been proved effective for the stabilization of spike protein of SARS-COV-2 and other betacoronavirus 12, 20, 50 .
In this work, we propose a pipeline to automatically stabilize proteins through computational mutagenesis. We analyze the conservation score, RMSF, and SASA of residues in the spike protein through the pipeline. We propose criteria based on the conservation score, RMSF, and SASA to identify residues for mutation. Finally, we utilize Eris to calculate the free energy change and find stabilizing mutants.
All source codes are deposited in: https://bitbucket.org/dokhlab/proteinstabilization. The 3D structure of the spike protein colored SASA. The red/blue colors indicate exposed/buried residues. (F) The 3D structure of the spike protein colored by RMSF. Red means flexible and blue means frozen.