key: cord-0924977-2crqj8so
authors: Rezaei, Shokouh; Sefidbakht, Yahya; Uskoković, Vuk
title: Comparative molecular dynamics study of the receptor-binding domains in SARS-CoV-2 and SARS-CoV and the effects of mutations on the binding affinity
date: 2020-12-17
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1860829
sha: f5760d7332fdf9250f751d7392768e8323f79a6b
doc_id: 924977
cord_uid: 2crqj8so

Here, we report on a computational comparison of the receptor-binding domains (RBDs) on the spike proteins of severe respiratory syndrome coronavirus‐2 (SARS-CoV-2) and SARS-CoV in free forms and as complexes with angiotensin-converting enzyme 2 (ACE2) as their receptor in humans. The impact of 42 mutations discovered so far on the structure and thermodynamics of SARS-CoV-2 RBD was also assessed. The binding affinity of SARS-CoV-2 RBD for ACE2 is higher than that of SARS-CoV RBD. The binding of COVA2-04 antibody to SARS-CoV-2 RBD is more energetically favorable than the binding of COVA2-39, but also less favorable than the formation of SARS-CoV-2 RBD-ACE2 complex. The net charge, the dipole moment and hydrophilicity of SARS-CoV-2 RBD are higher than those of SARS-CoV RBD, producing lower solvation and surface free energies and thus lower stability. The structure of SARS-CoV-2 RBD is also more flexible and more open, with a larger solvent-accessible surface area than that of SARS-CoV RBD. Single-point mutations have a dramatic effect on distribution of charges, most prominently at the site of substitution and its immediate vicinity. These charge alterations alter the free energy landscape, while X→F mutations exhibit a stabilizing effect on the RBD structure through π stacking. F456 and W436 emerge as two key residues governing the stability and affinity of the spike protein for its ACE2 receptor. These analyses of the structural differences and the impact of mutations on different viral strains and members of the coronavirus genera are an essential aid in the development of effective therapeutic strategies. Communicated by Ramaswamy H. Sarma

Betacoronaviruses (b-CoVs) are the genera of coronaviruses (CoVs) that are of particular clinical significance due to their virulence and the ability to jump across species. Prior b-CoVs that proved infectious to humans included severe acute respiratory syndrome coronavirus (SARS-CoV) and middle east respiratory syndrome coronavirus (MERS-CoV) (Yin & Wunderink, 2018) , which were both highly pathogenic. The symptoms of infection with any of these two b-CoVs were severe and the corresponding mortalities relatively high, which allowed the medical teams to stop the spread of the virus before it reached pandemic proportions. In December 2019, however, in China, severe respiratory syndrome coronavirus-2 (SARS-CoV-2) was identified as the third pathogenic virus of the b-CoV lineage .

Only a month after it had been reported in China, the spread of SARS-CoV-2 reached global proportions and in March 2020 the World Health Organization declared the pandemic, which is still unfolding. Unlike SARS-CoV and MERS-CoV, SARS-CoV-2 is typified by a broad range of symptoms produced in humans, spanning from total asymptomatic conditions to a bedridden and life-threatening state. This has allowed the virus to rapidly spread through the population, using seemingly healthy individuals as carriers. As of November 2020, the pandemic has not wound down and the virus continues to spread across the population at an accelerating rate. This state of affairs has led to the recruitment of the scientific community in the effort to provide scientific findings that would assist in curbing this ongoing pandemic.

One of the effects of the rapid spread and exchange of the viral material between people is the equally rapid mutability of the virus, which SARS-CoV-2 itself is particularly prone to in the first place owing to its ribonucleic genetic material. RNA viruses tend to be less stable than their DNA counterparts, but their mutability is significantly higher (Uskokovi c, 2020) . This has led to concerns that the virus may mutate into even less treatable forms than those that have spread through the community until now. Recent epidemiological analyses do confirm these expectations, indicating that mutated SARS-CoV-2 appearing in the most recent waves of infection in urban areas may be more contagious than the earlier versions (Long et al., 2020) . Detailed analyses of the effects of mutations on the molecular characteristics of SARS-CoV-2 are therefore of benefit, like the one we undertake here. Also, because any rational design of drugs targeting and neutralizing SARS-CoV-2 would be based on achieving a precise molecular recognition effect between the drug molecule and a receptor on the viral surface (or the interior under special circumstances), understanding these molecular interactions is of paramount importance.

CoVs in general contain 16 non-structural (nsp1-16) proteins and four structural ones. The structural ones include spike (S), membrane (M), nucleocapsid (N), and envelope (E) proteins . The envelope-anchored spike protein is an essential determinant of the specificity of the viral interaction with the host because it protrudes the viral particle most prominently and is the first of its components to come into contact with the host cell. The spike protein can be cleaved by host proteases into the N-terminal of S1 subunit and the C-terminal of S2 subunit (Huang et al., 2020) . The molecular recognition and attachment of the viral b-CoV particle to the host cell owes itself to the binding of the receptor-binding domain (RBD) on the S1 subunit of the spike protein to human angiotensin-converting enzyme 2 (ACE2) (Gui et al., 2017) . The RBD possesses two conformations: 'up' and 'down'. While the 'up' conformation corresponds to the receptor-accessible state, the 'down' conformation corresponds to the receptor-inaccessible state. Out of the two conformations, the 'up' conformation is thought to be less stable (Roy et al., 2020; Wrapp et al., 2020) . In addition, the RBD contains a receptor-binding motif (RBM), which is its most functional segment because of its direct involvement in the attachment to ACE2 (Li et al., 2005) and whose physiochemical properties are the direct determinants of this binding affinity. Based on previous studies, it is clear that both SARS-CoV-2 and SARS-CoV infect host cells by the binding of the RBD to ACE2 (Lan et al., 2020; Lukassen et al., 2020) .

In the last few decades, computational approaches have developed into a multitude of tools for forecasting the antigenic evolution of viruses. Accordingly, researchers investigate amino acid substitutions in protein regions corresponding to the antigenic change, after which conclusions of relevance for the efficient method of vaccine development are drawn (Klingen et al., 2018) . With this in mind, this study is based on one such computational assessment of the structural comparison of the RBDs of SARS-CoV-2 and SARS-CoV, with and without the interface with ACE2. In addition, over 40 mutations detected so far in the RBD structure of SARS-CoV-2 were investigated in terms of their effect on selected physicochemical properties of these molecules. With regards to the method, the RBD structures of SARS-CoV-2 and SARS-CoV were subjected to 50 ns molecular dynamics (MD) minimizations. In addition, the residue interaction networks (RINs) approach was applied to investigate proteinprotein interactions and mutation effects on the affinity of the RBD toward ACE2. In the RIN model, a network is computed and analyzed, such that its nodes and edges represent residues and bonds between the residue side chains, respectively (Piovesan et al., 2016) . This approach can give insight into the mechanisms of protein folding, conformational alterations, and the effects of amino acid mutations on the protein structure and function (Dehury et al., 2020; Fonseca et al., 2020) , with the major implications on the binding of antibodies such as COVA2-04 and COVA2-39 to the given mutated RBDs . As a result, these results are of significance for understanding the immune response to the infection, but also for developing novel therapies based on precisely targeted pharmacodynamics.

The sequences of wild type spike proteins isolated from the SARS-CoV-2 and SARS-CoV viruses were retrieved from UniProtKB (http://www.uniprot.org/), under the identifier numbers of P0DTC2 and P59594, respectively. SARS-CoV-2 and SARS-CoV RBDs consist of 223 (319-541) and 222 (306-527) residues, respectively . The automated sequence alignment was performed using the EMBOSS Needleman-Wunsch method (Needleman & Wunsch, 1970; Rice et al., 2000) . This tool indicates the percentage of identity and similarity between RBDs. I-TASSER server (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) was utilized to build the 3D structure of RBDs. I-TASSER, an online server for automated protein structure prediction, was used earlier for several biological studies (Yang & Zhang, 2015) . The predicted structures (5 models) were saved in the PDB format and sorted according to the C-scores. The best model was selected based on its C-score. The quality of the predicted structure was confirmed by PROCHECK (Laskowski et al., 1993 (Laskowski et al., , 1996 and ProSA, which calculate an overall quality score for 3D structures (S1 (A-B)) (Wiederstein & Sippl, 2007) . The local environment of amino acids was checked through the WHAT IF coarse packing quality control, which should stay above À5. Mutations were also predicted with the use of the WHAT IF mutation prediction tool (Vriend & Sander, 1993) .

A complete list of SARS-CoV-2 and SARS-CoV RBD mutations in the GISAID database (https://www.gisaid.org/) was utilized. Then, multiple sequence alignments between the wild type and the variants were performed to obtain sequence variations in the SARS-CoV-2 RBD region. In addition, a list of SARS-CoV RBD mutations in UniProt was considered and one of the mutations reported therein and not reported elsewhere (R426G) was included in the analysis.

2.2.1. Effect of mutations on thermodynamic parameters 2.2.1.1. Free energy binding and DDG (affinity and destabilizing). The structures of SARS-CoV-2 and SARS-CoV RBD-ACE2 complexes (PDB codes 6M0J and 2AJF) (Lan et al., 2020) and the structure of RBD-COVA2-04 and RBD-COVA2-39 complexes (PDB codes 7JMO and 7JMP) were downloaded from the Protein Data Bank . Then, the PRODIGY (PROtein binDIng enerGY) web server was applied to predict the binding affinity of protein-protein complexes (Vangone & Bonvin, 2015) . DynaMut (Rodrigues et al., 2018) and SAAMBE (Pahari et al., 2020) were used to predict DDG (destabilization and affinity) for mutations with respect to the given protein-protein complexes.

2.2.1.2. Solvation free energy. ProWaVE (https://www.prowave.org/) was used to calculate the solvation free energy (SFE) of the RBD and its variants. ProWaVE is an automated server for calculating the SFE based on the molecular theory of solvation (Chong et al., 2011; Chong & Ham, 2014) . ProWaVE investigates the structural properties and solvation thermodynamics for the wild and mutant types using fully atomistic, explicit-water MD simulations as well as threedimensional reference interaction site model (3D-RISM) theory. Thereupon, the following aspects were tried to be elucidated with the use of the 3D-RISM theory: 1) the protein structure and the folding kinetics alteration after inducing a single-point mutation; 2) the impact of mutations on the hydrophobicity of the protein.

2.2.1.3. Free energy surface. Protein processes such as folding or aggregation can be explained in terms of the molecular free energy, DG:

where k B is the Boltzmann constant, P is the probability distribution of the molecular system along the coordinate R, and P max denotes its maximum, which is subtracted from P(R) to ensure that DG ¼ 0 for the lowest free energy minimum. Common choices for R are the root-mean-square deviation of atomic positions (RMSD), radius of gyration (R g ), number of hydrogen bonds (H-bonds) or native contacts, or principal components. Here, the free energy was typically plotted at 310 K along axes defined by two parameters, e.g. R g and RMSD.

The degrees of protonation of titrable side chains at pH 7.4 were predicted through the Hþþ server (http://biophysics.cs. vt.edu/Hþþ), using known pKa values (Anandakrishnan et al., 2012) . The electrostatic representation of the molecule was calculated using UCSF Chimera 1.14 (Pettersen et al., 2004 ).

Energy minimization and MD simulations were carried out by GROMACS (2020.2) utilizing the Optimized Potential for Liquid Simulations (OPLS) force field (Jorgensen & Tiradorives, 1988 ). The structure produced by I-TASSER was used as an input. Water box was created with at least 1 nm (10 Å) distances from the protein using the SPC water model and applying boundary conditions. The system neutralization was done by adding Na þ and Clions at the concentration of 0.1 M. The MD simulation was carried out to examine the quality of the model structures by checking their stability via performing 50 ns simulations at a constant temperature 310 K (NVT). Energy minimization was performed in 50,000 steps to avoid any bad contacts generated while solvating the system. Then the NPT optimization was done for 100 ps.

To increase the likelihood of achieving the appropriate structure, the MD simulation was performed for 50 ns using the OPLS force field. Finally, GROMACS tools were applied for the trajectory analysis, including the parameters such as RMSD, RMSF, R g , SASA, and the number of H-bonds.

Eigenvectors of the covariance matrices and the projections of the first two principal components were calculated. The principal component analysis (PCA) was performed using GROMACS to reduce the dimensionality of the MD simulations data and identify the configuration space of a harmonic motion with only a few degrees of freedom.

RING 2.0 (http://old.protein.bio.unipd.it/ring/) (Martin et al., 2011; Piovesan et al., 2016) and SARS-CoV-2 RBD-ACE2 complex (PDB ID: 6M0J) were used to generate the SARS-CoV-2 RBD-ACE2 RIN. Then, the RING output was visualized in the Cytoscape. MCODE (Bader & Hogue, 2003) plugin in Cytoscape was used to find clusters of the network. In general, clusters in the RINs for RBD-ACE2 and free RBD (including five or more nodes), stress and betweenness were used to identify the key residues. In addition, interactions at the contact surface of the RBD and ACE2 were investigated.

Key residues were defined using stress and betweenness as the two elementary parameters of the local metrics. High betweenness values here indicate that the node is a mediator of interactions with other nodes and that it may be a key structural residue (Hu et al., 2014) . In contrast, stress is defined as the number of shortest paths passing through a node.

CHARMM-GUI Membrane Builder was applied to simulate the complex biological membrane and produce the representation of its topology (Allouche, 2011) . Also, APBS-1.5 was used to calculate the dipole moment (Jurrus et al., 2018) .

Investigating the effects of sequence variations in SARS-CoV-2 RBD is essential for the understanding of pathogenesis and the development of safe and effective prevention and treatment strategies for COVID-19, as through viral detection, vaccine design, and development of drugs . One such study, gaining an insight into the basic aspects of thermodynamics and pharmacodynamics of the interaction between wild-type and mutated SARS-CoV-2 RBDs and its receptor in the host cell has been undertaken here.

The results of the amino acid sequence alignment revealed that the RBD of SARS-CoV-2 shares 73.1% identity and 82.1% similarity with the RBD of SARS-CoV ( Figure 1 ). The secondary structure of the SARS-CoV-2 RBD contains five antiparallel b sheets (b1, b2, b3, b4 and b7) with short connecting helices and loops. The extended region between the b4 and b7 strands is the receptor-binding motif (RBM), stretching from R403 to Y508 in SARS-CoV-2 and from K403 to Y508 in SARS-CoV. The SARS-CoV-2 RBD has two b sheets less than the SARS-CoV RBD (Figure 2 (A)), which is the first indicator of the higher structural flexibility of SARS-CoV-2 RBD compared to that of SARS-CoV RBD. Also, as shown in Figure 1 , interestingly, the a-helix (a3) is formed at two different positions in the RBDs of SARS-CoV-2 and SARS-CoV. It appears that the presence of P384 in the (-VSPTKLN-) region prevents the formation of a-helix in SARS-CoV-2 RBD. Further, a total of nine cysteine residues were detected in the RBD of both viruses, eight of which involved four pairs of disulfide bonds. Three of those pairs aid in the stabilization of the b sheet structure, including Cys336-Cys361, Cys379-Cys432 and Cys391-Cys525; the remaining pair (Cys480-Cys488) is located in the distal end of the RBM loop. Also, an additional insertion (V483) was found in the (-AGSTPCNGVEGF-) loop of the SARS-CoV-2 RBD that is not found in the same region of the SARS-CoV RBD (-PDGKPCTPPAL-) ( Figure 2 (B)). Sequence comparisons showed that some of the loop residues, such as V483, E484, G485, F486 and N487, bind to antibodies. Therefore, it is likely that V483 provides an additional site for the binding of antibodies, such as COVA2-04. It lies in-between the epitope residues A475 and F486, which were previously delineated as the binding sites for neutralizing antibodies (Yi et al., 2020) . In general, investigating the impact of mutations, particularly on the conformation of the interacting proteins, is essential because mutations in proteins can have a dramatic effect on the protein folding and stability. Also, mutations alter the kinetics and thermodynamics of protein-protein interactions , and these mutations can be either selectively beneficial to the organism via evolution or straightforwardly detrimental (Jubb et al., 2017) . Therefore, the study of the effect of mutations is essential for different biomedical applications, including personalized medicine stemming from disease-associated mutation analyses or highly specific drug design and remedial interventions.

3.2.1. Effect of mutations on thermodynamic parameters 3.2.1.1. Binding free energy and DDG (affinity and destabilizing). Charting the thermodynamics of protein-protein interactions is an important precondition for revealing their mechanisms of function and understanding the effects of disease-related conditions. The total binding free energies (DG Binding ) for SARS-CoV-2 and SARS-CoV RBD-ACE2 complexes were calculated. As it is shown in Table 1 , DG Binding for SARS-CoV-2 RBD-ACE2 was À11.9 kcal mol À1 , 1.1 kcal mol À1 lower than that of SARS-CoV RBD-ACE2 (À10.8 kcal mol À1 ), clearly indicating that the binding affinity of the SARS-CoV-2 RBD to ACE2 is higher than that of the SARS-CoV RBD.

In addition, DG Binding for two SARS-CoV-2 RBD-antibody complexes was computed. The results showed that DG Binding for RBD-COVA2-04 (À11.6 kcal mol À1 ) was higher than that for RBD-COVA2-39 (À7.3 kcal mol À1 ) (Table 2) . Interestingly, DG Binding of the RBD-COVA2-04 complex was very close to the free energy of the SARS-CoV-2 RBD-ACE2 complex. Still, based on these thermodynamic considerations, the binding of the SARS-CoV-2 RBD to ACE2 is more energetically favorable than its binding to COVA2-04 or COVA2-39, which is only one aspect where the virus has an advantage over the immune protection. In addition to this, the dissociation constant of the SARS-CoV-2-ACE2 complex is relatively low at 4E-09 (Table 1) , especially in comparison with that of the RBD-COVA2-39 complex, which is almost three orders of magnitude higher (Table 2) .

A difference in the binding affinity to the epitope between the wild protein domain (DG W ) and the domain containing mutant residues (DG M ) is defined as DDG Affinity (DDG ¼ DG m -DG w ). This free energy difference can be a descriptor of how mutations influence the protein stability. Overall, DDG values below zero indicate that the mutation has caused destabilization of the protein; otherwise, it has caused its stabilization. In this analysis, the change in the binding affinity upon the mutation (DDG Affinity ) and in the folding free energy (DDG Destabilizing ) at the RBD-ACE2 (Figure 4(A,B) ) and RBD-antibody ( Figure 6(A,B) ) interfaces were considered.

Among 25 analyzed mutations of the SARS-CoV-2 RBD-ACE2 complex, 14 had DDG < 0 (Table S1) . Interestingly, all of the X!F mutations (L445F, V503F, and Y453F) had DDG Destabilizing > 0 (Table S1) (Figures 4(A) and 5). This positive value may be due to the p-electrons of the phenyl ring, which can undergo stacking with the neighboring aromatic rings of amino acids such as F, W or Y. Therefore, promoting hydrophobic interactions, as in this case, can lead to protein folding and structural stability. In general, Figure 4 (A) shows mutations leading to changes in protein interactions, which affect the stability of the structure.

Investigation of 21 mutations in the two RBD-antibody complexes, namely RBD-COVA2-04 or RBD-COVA2-39, indicated an equal variability of DDG Affinity values ( Figure 6(A,B) ) as that found in the mutations of the SARS-CoV-2 RBD-ACE2 complex. These different values indicate a difference in the contribution of individual amino acids to protein-protein interaction. Equally, they suggest that not all mutations are equally important. Overall, as expected based on previous arguments, positive DDG Affinity values caused a decrease in the protein-protein binding affinity (Table S1 and S2).

3.2.1.2. Solvation free energy. Since most proteins operate under the aqueous solvation, information about the interaction of the protein and its domains with water can be useful for predicting the stability of the protein and its specific interaction with ligands or the nonspecific binding to surfaces such as the cell membrane (K€ onig et al., 2013; Martini et al., 2013) . In general, the information obtained from the investigation of the effect of side chains on the SFE can be important for understanding the protein folding and its stability (K€ onig et al., 2013) . Moreover, defining the most appropriate mode of interaction of antibodies with antigens is a prerequisite for an effective antibody design. For this purpose, physicochemical properties of wild and mutant type antibodies, including the nature of the interaction energy where water molecules by default play a key role, must be considered (Webster et al., 1994) . SFE (G sol ) values calculated for the wild-type and the variants of SARS-CoV-2 RBD and SARS-CoV RBD are shown in Figure S2(A,B) and Table 3 . By comparing these values, it can be noted that SARS-CoV RBD has a higher SFE and, thus, hydrophobicity than SARS-CoV-2 RBD. The higher charge of the SARS-CoV-2 RBD compared to that of SARS-CoV RBD can be one contributing factor to the lower SFE of the former protein domain (Table 3 ).

In addition to this crude comparison, three groups of mutations in the SARS-CoV-2 RBD were considered for a more selective study of the SFE. They included mutations converting (i) a charged amino acid to a neutral one (E484A, E484Q, K417N), (ii) a neutral amino acid to a charged one (N493K, G485R, S477R, T478K, G504D), and (iii) a non-hydrophobic amino acid to a hydrophobic one (Q493L, Y453F, S477I, S494L) (Table 4 ). Here, two mutations particularly stand out because they had a most dramatic effect on the SFE. They include K417N, which increased the SFE by nearly 50 kcal mol À1 , and S477R, which reduced it by over 30 kcal mol À1 (Table 3) . With K417N substituting a neutral residue for a charged one and S477R substituting a charged residue for a neutral one, the key effect of charge on the SFE and the hydrophobicity, the latter of which is inversely proportional to charge, becomes obvious. Overall, the comparison of the wild-type SARS-CoV-2 RBD with the variants from the groups (i) and (ii), including S477R and K417N (Figure 7) , shows the effects different mutations can have on charge, SFE, hydrophobicity, and electrostatic potential of the protein. What is common to all these single-point mutations, including S477R and K417N, is that they are most likely to affect the SFE by altering the folding and the stability of the protein structure. Needless to add, considering these changes made to the RBD due to mutations, their effect on the function and the properties of this domain can be harnessed for various therapeutic aims.

3.2.1.3. Free energy surface. The plots in Figure 8 provide a comparative view of the free surface energy landscape as a function of RMSD and R g for the RBDs of SARS-CoV-2 and SARS-CoV. The general trend applying to both RBDs is that the higher the value of the two parameters (RMSD and R g ), the lower the free energy surface. Also, given that the more stable structure possesses a lower free energy, the SARS-CoV RBD structure appears to be more stable than that of SARS-CoV-2 RBD (Figure 8(A,B) ).

Electrostatic interactions between protein residues and the solvent are one of the key determinants of protein folding and stability (Strickler et al., 2006; Zhou & Pang, 2018) . The electrostatics also regulates protein interactions with other proteins as well as with molecules other than proteins, such as small-molecule drugs (Voet et al., 2013) . The electrostatic features of the protein are determined by the distribution of whole and partial charges across the 3D protein structure (Vascon et al., 2020) . Here, Coulomb's law is not appropriate for describing electrostatic effects in proteins because it applies to a system with uniform dielectric properties, whilst proteins dispersed in a medium display a gradient dielectric constant, with a hydrophobic core surrounded by the hydrophilic surface and covered by the solvent. Hence, electrostatic calculations for proteins were carried out using the Poisson-Boltzmann equation (Vascon et al., 2020) .

As in agreement with the earlier demonstrated effect of charged residues on the SFE, electrostatic interactions between the solvent (water or physiological media) and the protein were evidenced as determinants of the degree of hydrophobicity. The distribution of the electrostatic potential on the surfaces of the RBDs of SARS-CoV-2 and SARS-CoV was evaluated too. Figure 9 (A,B) shows this distribution, with blue and red surfaces indicating electropositive and electronegative surfaces, respectively. Investigation of the amino acid sequence of RBDs showed that the net charge for SARS-CoV-2 and SARS-CoV is 7 and 4, respectively. This difference in the net charge is explained by the larger number of charged amino acid residues in the RBD of SARS-CoV-2 than in the RBD of SARS-CoV. Charged residues that particularly stand out are K445 in the SARS-CoV-2 RBM and R426 and D463 in the SARS-CoV RBM. Since these residues are involved in binding to ACE2, they are important to examine in more detail.

With the electrostatic potential having an essential role in various biochemical reactions, such as the binding of a ligand to a specific protein target, electrostatic effects of mutations are important to consider (Figure 9(C-G) ). Here we consider five different mutations, namely N439K, G485R, S477R, G504D and T478K, and demonstrate that the distribution of charges on the RBD profoundly changes with each of these single-point mutations. This can be evidenced by noticing intense changes in the two colors representing positive and negative surface charges entailing each of the mutations and affecting most segments in the domain. Still, the distribution of charges is most drastically affected at the site of substitution and its immediate vicinity on the surface of the folded protein. In view of this, the effects of mutations must be manifold and may critically change the specific interactions with drugs, antibodies or ACE2. These structural changes allow the virus to potentially evade the therapy through timely mutation. These findings may be of benefit not only for delineating the most critical amino acids for the RBD function, but also for finding out how specific mutations may lead to disorder or improvement in the function of the RBD and its binding affinity through alterations in the electrostatic potential.

R426 is known to be one of the key residues for the binding of antibodies to the SARS-CoV RBD (Amin et al., 2020) . This residue is positively charged under the physiological conditions, forming a salt bridge in interaction with receptors and antibodies. As a result, mutations at this site can have a large effect on protein-protein interactions. The R426G change, for example, causes the removal of this positively charged amino acid in the SARS-CoV-2 RBD, meaning that no salt bridge can be created by the RBD at this position. The significant decrease in the binding affinity of SARS-CoV-2 RBD to antibodies may be due to the lack of this salt bridge formation. The N439K mutation swapping an uncharged group for a positively charged one may elicit the opposite effect, possibly inducing the formation of a salt bridge at a site that previously did not form them. Here, it is important to add that N439K is the mutation with the second highest frequency of occurrence (Figure 3(A) ). The effects of this mutation may be such that it induces the creation of a salt bridge as in the SARS-CoV RBD, thus increasing the strength of the binding of the RBD to the receptors or antibodies. This increased affinity, of course, is a double-edged sword because it may improve the interaction with drugs and antibodies, thus increasing their efficacy, but it also strengthens the interaction of the vital particle with the host cell.

MD simulations were performed for the SARS-COV-2 RBD, the SARS-CoV RBD, and the SARS-CoV-2 RBD-ACE2 complex in order to understand the difference in the structure of the two RBDs, alongside the implications for the stability of the SARS-CoV-2 RBD induced by these changes. To determine the stability and convergence of structure, RMSD and R g were calculated (Sargsyan et al., 2017) . RMSDs for the SARS-CoV-2 RBD and the SARS-CoV RBD were found to be $ 0.5 and $0.36 nm, respectively. The RMSD plot indicates that both RBDs reached equilibrium after 40 ns and that the SARS-CoV RBD is more stable than the SARS-CoV-2 RBD (Figure 10(A) ). The RMSD of the RBD-ACE2 complex during 10 ns of the MD simulation time displayed the lowest RMSD ($0.25 nm) (Figure 10(C) ), which can be due to the presence of the stabilizing interactions in the protein-protein complex.

R g was calculated from the 50 ns trajectory to gain insight on the compactness and rigidity of the RBDs. The largest deviation was observed in the SARS-CoV-2 RBD, which can indicate functional movements and structural 'breathing' due to the lesser compactness of the structure (Figure 10(B) ). R g was also calculated for the SARS-CoV-2 RBD-ACE2 complex (Figure 10(D) ), showing deviation within $4.4 ns ($3.17 to $3.12 Å). In addition, the residue-based root mean square deviation (RMSF) plot was constructed from the 50 ns data in order to understand the deviation of each RBD. It was shown that the RMSF of the backbone for the SARS-CoV-2 RBD displays more flexible residues at residue numbers 477 to 487 as compared to the SARS-CoV RBD (Figure 11 ). This conformational flexibility of the binding domain may be one of the key factors in ensuring the high infectivity of SARS-CoV-2, as in analogy with other active sites in proteins (Obadan et al., 2019; Penfold et al., 2004) .

In order to analyze H-bond interactions during the 50 ns and 10 ns simulation times, intramolecular H-bond plots were constructed for the RBDs of SARS-CoV-2 and SARS-CoV (Figure 12(A) ). From the comparative plots, it could be deduced that the SARS-CoV-2 RBD has a similar number of intramolecular H-bonds as the SARS-CoV RBD (Figure 12(A) ). In addition, Figure 12 (C) indicates the number of H-bond in the SARS-CoV-2 RBD-ACE2 protein complex. The high number of H-bonds here indicates the relatively high stability of the interaction between the RBD and its receptor. The solvent-accessible surface area (SASA) as the surface region of the protein that is accessible to the solvent was calculated for the RBDs of SARS-CoV-2 and SARS-CoV and for the SARS-CoV-2 RBD-ACE2 complex as a function of time and the resulting plots are shown in Figure 12(B,D) . The results obviously indicate that the SARS-CoV-2 RBD has a higher SASA ($134 nm 2 ) than the SARS-CoV RBD ($126 nm 2 ). In addition, based on these plots, the SASA for the SARS-CoV-2 was $370 nm 2 . With SASA being a measure of the degree to which an amino acid is exposed to the environment, the higher SASA for the RBD of SARS-CoV-2 indicates its more open and diffused structure than that of the more compact RBD of SARS-CoV. Any increase or decrease in the SASA in response to structural changes indicates alterations in the conformation of the protein.

PCA is a method for the analysis of the MD trajectory that extracts dominant modes in the overall molecular motion. Here, it was used to determine the dominant modes of molecular flexibility and functional motions in the RBDs. To perform the PCA, eigenvectors were obtained by diagonalization of the covariance matrix of the coordinate fluctuations of the RBDs of SARS-CoV-2 and SARS-CoV. To understand the total motion of these structures in the phase space, the first two eigenvectors (1 and 2) were projected onto it. The comparison has shown that the motion properties described by these first two eigenvectors are different in the SARS-CoV-2 RBD as compared to the SARS-CoV RBD (Figure 13 ).

The RIN corresponding to the receptor-bound form of the SARS-CoV-2 RBD is presented in Figure 14 . The greater number of clusters and the correspondingly higher MCODE score was obtained for the SARS-CoV-2 RBD complex with ACE2 than for the SARS-CoV-2 RBD without the receptor. The greatest overlap in cluster members in both the free and the complex form of the RBD was noted for the residues I358, C361, C336, F456, and F515. Among these residues, F456 appears particularly significant because the RIN in Figure 14 indicates its direct role in mediating the interaction between the RBD and ACE2. In addition, this residue was identified as one of the mutations (F456L) in the RBD of SARS-CoV-2 ( Figure 3) . Critical residues at the interface between the SARS-CoV-2 RBD and ACE2 in their complex are indicated in the RIN presented in Figure 14 . Also, Table 5 shows the seemingly most important residues of the SARS-CoV-2 involved in the interaction between its RBD and ACE2 via H-bonds, van der Waals bonds, ionic and p-p stacking interactions. Among these residues, K417 (!417 N), Y453 (!453 F), F486 (!486 L), A475 (!475 V), and F456 (!456 L) are of utmost importance because their mutations could alter the ability of the RBD to bind to ACE2, which could lead to the loss of efficacy of antibodies/drugs targeting the RBD or competing with it for ACE2.

The betweenness analysis identified S349, A352, W436, and L452 as residues of interest in the SARS-CoV-2 RBD without ACE2 (Figure 15(A) ). The stress metric, in contrast, identified N422, W436, Y486, W353, F347, F400, and L461 as residues of interest in the SARS-CoV-2 RBD without ACE2 ( Figure  15 (C)). As far as the SARS-CoV-2 RBD with ACE2 is concerned, the betweenness analysis identified L452, S349, and A352 as residues of interest ( Figure 15(B) ), while the stress metric identified residues Y453, W436, Y505, R403, Y508, F342, and F338 as those of interest ( Figure 15(D) ). Different residues emerging as the results of the stress and betweenness analyses applied to the free RBD and the RBD-ACE2 complex may be indicative of significant structural alterations occurring upon the binding of the RBD to the receptor. Interestingly, by appearing in 3 out of 4 of these distinct analyses, W436 appears to be a key residue in both the free and the receptor-bound structures of the SARS-CoV-2 RBD and may play an essential role in ensuring both the stability and the binding of the RBD to its receptor.

Dipole moment can play an important role in the mechanism of interaction of proteins with receptor and drugs (Lien et al., 1982) . The computed dipole moments of the SARS-CoV-2 RBD, the SARS-CoV RBD, and the SARS-CoV-2 RBD-ACE2 were 720 Debye (D), 204 D, and 760 D, respectively. As shown in Figure 16 , it appears that the RBDs of SARS-CoV-2 and SARS-CoV, albeit having similar functions, have quite distinct dipole moments. This difference in the dipole moment reiterates the different charge distribution on the surface of the two RBDs and is expected in view of the fact that the distribution of charges in the protein defined the dipole moment. Further, the dipole moment of the SARS-CoV-2 RBD-ACE2 complex is more inconstant than the relatively stable, albeit naturally fluctuant dipole moments of the RBDs of SARS-CoV-2 and SARS-CoV alone. Overall, the dipole moment values are significant because large dipole moments may promote binding and also control other aspects of the protein function (Felder et al., 2007) . Large dipole moments lead to stabilization of the protein-protein complexes and are also of critical importance for achieving the correct orientation prior to the formation of the given complexes (Pichierri, 2013) . For this reason, many studies consider the dipole moment as a key parameter governing the drug-receptor interaction (Lien et al., 1982) . Here, it should be kept in mind that this large of the dipole moment value was calculated for the RBD only and may differ from the dipole moment of the spike protein as a whole.

This computational study consisted of the analysis of the major and the minor structural factors of the SARS-CoV-2 RBD playing a role in determining the virulent function of the SARS-CoV-2 spike protein. The protein dynamics (Hollingsworth & Dror, 2018) , electrostatics (Zhou & Pang, 2018) and residue networking (Gordon et al., 2020) were examined to elucidate this relationship between the protein structure and function. All-atom molecular dynamics (MD) simulations were carried out in the course of 50 ns at 310 K and a number of physiochemical parameters were subjected to analysis. A special emphasis was put on the assessment of the impact of amino acid mutations on the structure and thermodynamic properties of 42 variants of the SARS-CoV-2 RBD discovered so far. It goes beyond the scope of this paper to comment on the continuously evolving interrelation of different structural parameters investigated here, but the major findings pertaining to individual parameters could be compiled. First of all, the binding affinity of SARS-CoV-2 RBD for ACE2 is shown to be higher than that of SARS-CoV RBD. Next, the binding of COVA2-04 antibody to SARS-CoV-2 RBD is more energetically favorable than the binding of COVA2- Asp38 VDW Tyr505

Lys353 VDW 39, but also less favorable than the formation of the SARS-CoV-2 RBD-ACE2 complex, hinting at a possible thermodynamic advantage that the virus has over the immune system. The net charge, the dipole moment and hydrophilicity of SARS-CoV-2 RBD are higher than those of SARS-CoV, producing lower solvation and surface free energies and thus lower stability. This lower stability of SARS-CoV-2 was confirmed in the MD simulations of RMSD, RMSF, radius of gyration and H-bond distribution. It is paralleled by the higher flexibility of SARS-CoV-2 RBD compared to more rigid SARS-CoV RBD, owing in part to two less b sheets in the RBD of SARS-CoV-2 than in that of SARS-CoV. The structure of SARS-CoV-2 RBD is correspondingly more open, with the solvent-accessible surface area larger by 8 nm 2 than that of SARS-CoV RBD. It also undergoes intense structural alterations upon binding to ACE2. Mutations have a dramatic effect on distribution of charges, most prominently at the site of substitution and its immediate vicinity. Charge alterations upon single-point mutation have a prominent effect on the free energy landscape, increasing the solvation free energy by $50 kcal mol À1 for a charge-promoting substitution such as K417N and reducing it by over 30 kcal mol À1 for a charge-neutralizing substitution such as S477R. Meanwhile, all X!F mutations exhibit a stabilizing effect on the RBD structure through p stacking. Considering the effect of different mutations on the thermodynamics of SARS-CoV-2 RBD, some of them have imposed considerable changes to DDG Destabilizing and DDG Affinity . These changes proceed in such a way that with an increasing variety of mutations, the value of DDG Affinity has tended to increase and the value of DDG Destabilizing to decrease. Finally, F456 and W436 emerged from the RIN analysis of betweenness and stress as two key residues governing the stability and affinity of the spike protein for its ACE2 receptor.

Notwithstanding that the effect that varying physiological conditions and the presence of different drugs can have on the mechanism of interaction of the SARS-CoV-2 RBD with ACE2, the network analysis can provide a valuable tool for predicting changes to the interaction potential and the rate of interaction between the RBD and its receptor in response to structural changes, such as mutations. It can also help identify key mutations that may occur in the future and thus facilitate the antiviral drug discovery from the fundamental standpoint of protein-protein interaction between the virus and the host (Ackerman et al., 2018; Brito & Pinney, 2017) . This research route is worth examining in the future when the number of mutants reaches a statistically acceptable level, which would allow the determination of the direction in which the RBD mutations proceed with satisfactory levels of confidence.

SARS-CoV-2 has gripped the world in an ongoing pandemic and timely investigations of various aspects of its nature are of paramount importance. MD simulations were performed to structurally compare the RBDs on the spike proteins of SARS-CoV-2 and SARS-CoV, with and without binding to ACE2 and also in response to different mutations. Given the significant recent work reported on SARS-CoV-2, an attempt was made to examine parameters that were considered to a lesser extent in previous studies. The physicochemical properties of the RBD analyzed included the free energies of binding, destabilizing and solvation, the electrostatic potential and hydrophobicity; the effects of mutations on each of these properties was assessed. Some of the most critical mutations were highlighted, such as K417 and N439, and their effect on the physicochemical properties of the RBD discussed, along with the corollaries for the design of antibodies and drugs, most of which utilize the given RBD as the target. The study suggests that changes induced by these mutations may play a role in reducing or increasing the number of protein interactions and protein surface charges, consequently affecting the stability and function of the RBD and thus the infectivity of SARS-CoV-2. Therefore, the structural details and physicochemical parameters of the RBD reported here may help in the development of antibodies or drugs at the molecular level. Also, previous studies indicated that many disease-associated amino acid mutations are located at the protein-protein interface, for which reason we have also focused on understanding the structural changes occurring at the interface between SARS-CoV-2 and ACE2. Mutations are also known to affect protein-protein interaction by changing the binding affinity, which was analyzed in this study in detail. Information on the impact of mutations on protein-protein interaction may aid in classifying mutations into pathogenic and harmless. From there on, this information could prove essential for the accelerated development of appropriate therapeutic and drug design strategies. This study suggests that next to investigating the SARS-CoV-2 RBD structure, researchers dealing with drugs and antibodies field could benefit from the fast and precise computational tools to predict the binding affinity and its changes in response to mutation.

Network-guided discovery of influenza virus replication host factors

Software news and updates gabedit-A graphical user interface for computational chemistry softwares

Comparing the binding interactions in the receptor binding domains of SARS-CoV-2 and SARS-CoV

Hþþ 3.0: Automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations

An automated method for finding molecular complexes in large protein interaction networks

Protein-protein interactions in virushost systems

Potent neutralizing antibodies from COVID-19 patients define multiple targets of vulnerability

Interaction with the surrounding water plays a key role in determining the aggregation propensity of proteins

Structural and thermodynamic investigations on the aggregation and folding of acylphosphatase by molecular dynamics simulations and solvation free energy analysis

Effect of mutation on structure, function and dynamics of receptor binding domain of human SARS-CoV-2 with host cell receptor ACE2: A molecular dynamics simulations study

A server and database for dipole moments of proteins

CoRINs: A tool to compare residue interaction networks from homologous proteins and conformers

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

Cryo-electron microscopy structures of the SARS-CoV spike glycoprotein reveal a prerequisite conformational state for receptor binding

Molecular dynamics simulation for all

Residue interaction network analysis of Dronpa and a DNA clamp

Structural and functional properties of SARS-CoV-2 spike protein: Potential antivirus drug development for COVID-19

Neutralizing antibodies against SARS-CoV-2 and other human coronaviruses

The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin

Mutations at protein-protein interfaces: Small changes over big surfaces have large impacts on human health

Improvements to the APBS biomolecular solvation software suite

In silico vaccine strain prediction for human influenza viruses

Absolute hydration free energies of blocked amino acids: Implications for protein solvation and stability

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

PROCHECK: A program to check the stereochemical quality of protein structures

AQUA and PROCHECK-NMR: Programs for checking the quality of protein structures solved by NMR

Structural biology: Structure of SARS coronavirus spike receptor-binding domain complexed with receptor

Use of dipole moment as a parameter in drug-receptor interaction and quantitative structureactivity relationship studies

Composition and divergence of coronavirus spike proteins and host ACE2 receptors predict potential intermediate hosts of SARS-CoV-2

Molecular architecture of early dissemination and massive second wave of the SARS-CoV-2 virus in a major metropolitan

SARS-CoV-2 receptor ACE 2 and TMPRSS 2 are primarily expressed in bronchial transient secretory cells

RING: Networking interacting residues, evolutionary information and energetics in protein structures

Water-protein interactions: The secret of protein dynamics

A general method applicable to the search for similarities in the amino acid sequence of two proteins

Flexibility in vitro of amino acid 226 in the receptor-binding site of an H9 subtype influenza A virus and its effect in vivo on virus replication, tropism, and transmission

SAAMBE-3D: Predicting effect of mutations on protein-protein interactions

Flexibility in the receptor-binding domain of the enzymatic colicin E9 is required for toxicity against Escherichia coli cells

UCSF Chimera-A visualization system for exploratory research and analysis

Dipole-dipole interactions in protein-protein complexes: A quantum mechanical study of the ubiquitin-Dsk2 complex

The RING 2.0 web server for high quality residue interaction networks

EMBOSS: The European molecular biology open software suite

DynaMut: Predicting the impact of mutations on protein conformation, flexibility and stability

Dynamic asymmetry exposes 2019-nCoV prefusion spike

How molecular size impacts RMSD applications in molecular dynamics simulations

Protein stability and surface electrostatics: A charged relationship

Why have nanotechnologies been underutilized in the global uprising against the coronavirus pandemic?

Contacts-based prediction of binding affinity in protein-protein complexes. eLife, 4, 1-15

Protein electrostatics: From computational and structural analysis to discovery of functional fingerprints and biotechnological design

Electrostatic similarities between protein and small molecule ligands facilitate the design of protein-protein interaction inhibitors

Quality control of protein models: Directional atomic contact analysis

Structural and functional basis of SARS-CoV-2 entry by using human ACE2

Antibody-antigen interactions

ProSA-web: Interactive web service for the recognition of errors in three-dimensional structures of proteins

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

An alternative binding mode of IGHV3-53 antibodies to the SARS-CoV-2 receptor binding domain

I-TASSER server: New development for protein structure and function predictions

Key residues of the receptor binding motif in the spike protein of SARS-CoV-2 that interact with ACE2 and neutralizing antibodies

MERS, SARS and other coronaviruses as causes of pneumonia

Mutation effect estimation on protein-protein interactions using deep contextualized representation learning

Electrostatic interactions in protein structure, folding, binding, and condensation

The support and resources from the Center for High Performance Computing at Shahid Beheshti University of Iran are gratefully acknowledged.

No potential conflict of interest was reported by the authors.

http://orcid.org/0000-0003-0538-0252 Vuk Uskokovi c http://orcid.org/0000-0003-3256-1606