key: cord-1042131-xucpbsbm authors: Lai, Hien T. T.; Nguyen, Ly H.; Kranjc, Agata; Nguyen, Toan T.; Nguyen-Manh, Duc title: Elucidating the differences in the molecular mechanism of receptor binding between 2019-nCoV and the SARS-CoV viruses using computational tools date: 2020-04-21 journal: bioRxiv DOI: 10.1101/2020.04.21.053009 sha: 8635dca8f394aae9df6294b571f948dcdacec568 doc_id: 1042131 cord_uid: xucpbsbm The outbreak of the 2019-nCoV coronavirus causing severe acute respiratory syndrome which can be fatal, especially in elderly population, has been declared a pandemic by the World Health Organization. Many biotechnology laboratories are rushing to develop therapeutic antibodies and antiviral drugs for treatment of this viral disease. The viral CoV spike (S) glycoprotein is one of the main targets for pharmacological intervention. Its receptor-binding domain (RBD) interacts with the human ACE2 receptor ensuring the entry of the viral genomes into the host cell. In this work, we report on the differences in the binding of the RBD of the previous coronavirus SARS-CoV and of the newer 2019-nCoV coronavirus to the human ACE2 receptor using atomistic molecular dynamics techniques. Our results show major mutations in the 2019-nCoV RBD with respect to the SARS-CoV RBD occurring at the interface of RBD-ACE2 complex. These mutations make the 2019-nCoV RBD protein backbone much more flexible, hydrophobic interactions are reduced and additional polar/charged residues appear at the interface. We observe that higher flexibility of the 2019-nCoV RBD with respect to the SARS-CoV RBD leads to a bigger binding interface between the 2019-nCoV RBD and ACE2 and to about 20% more contacts between them in comparison with SARS-CoV. Taken together, the 2019-nCoV RBD shows more stable binding interface and higher binding affinity for the ACE2 receptor. The mutations not only stabilize the binding interface, they also lead to overall more stable 2019-nCoV RBD protein structure, even far from the binding interface. Our results on the molecular differences in the binding between the two viruses can provide important inputs for development of appropriate antiviral treatments of the new viruses, addressing the necessity of ongoing pandemics. The COVID-19 disease that has rapidly spread from the Chinese Wuhan city to the rest of the world is caused by the 2019 novel coronavirus (2019-nCoV) named also the severe acute respiratory syndrome corona virus 2 (SARS-CoV-2). 1 The sequencing of its genome revealed that it is closely related to other coronaviruses, such as severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV), which have killed hundreds of people in the last two decades. 2, 3 Coronaviruses are enveloped, positive sense, single-stranded RNA viruses that carry on their surface spike-like projections giving it a crown-like appearance under the electron microscope; hence the name coronavirus. 4 These spikes (S) enable the fusion between viral and host membranes and are essential for the beginning of the enveloped virus infection. 5, 6 They are composed of a large ectodomain, a trans-membrane anchor and a short intracellular tail ( Fig. 1) . The ectodomain consists of a receptor-binding subunit S1 and a membrane-fusion subunit S2 which are crucial for binding of the virus to the host cell surface and for entry of the viral genomes into the target cell, respectively. [7] [8] [9] [10] [11] The S1 subunit is composed of N-terminal domain (S1-NTD) and C-terminal domain (S1-CTD), which can both function as the Receptor Binding Domain (RBD), recognizing and binding to the host protein receptor. The cellular entry receptor for the SARS-CoV is Angiotensin Converting Enzyme 2 (ACE2). 12 The RBD of SARS-CoV and SARS-CoV-2 share significant sequence similarity, therefore several research groups investigated if SARS-CoV-2 uses the same cellular entry receptor. 2, [13] [14] [15] [16] And they all confirmed that this is indeed the case. The ACE2 enzyme is a part of the renin-angiotensin-aldosterone system (RAAS) that controls blood pressure, fluid and electrolyte balance, systemic vascular resistance, tissue damage etc. 17 The primary physiological role of the ACE2 is to convert a peptide hormone angiotensin II into angiotensin (1) (2) (3) (4) (5) (6) (7) , that have protective role as they diminish the inflammation and fibrosis. 17, 18 The SARS-CoV and SARS-CoV-2 viruses bind to ACE2 preventing the metabolism of the angiotensin II, which accumulates and can cause deleterious effects. There are many controversial debates about the beneficial or harmful role of the RAAS inhibitors in patients with COVID-19. [18] [19] [20] [21] [22] Number of experts defend their beneficial role, but the qualified clinical studies will have to be done to answer this important question. The SARS-CoV and SARS-CoV-2 are closely related, nonetheless, significant differences were Figure 1 : A schematic view of the coronavirus spike protein. The spike protein bears a large ectodomain that is composed of a subunit S1 which has at the top the Receptor Binding Domain (RBD) and the subunit S2. It further features the transmembrane anchor (TM anchor) and the intracellular tail (IC tail). ACE2 is the Angiotensin Converting Enzyme 2, the main host entry receptor for SARS-CoV and SARS-CoV-2. reported for their protein structures, importantly the mutations in the SARS-CoV-2 RBD domain with respect to the SARS-CoV that can impact the binding affinity for the host receptor. 14, 23 In our study we aimed to compare the structural and energetic differences in the binding of the RBDs of the SARS-CoV-2 and SARS-CoV to the ACE2 receptor by Molecular Dynamics (MD) simulations. The viral RBD represents one of the main possible targets for the development of antiviral drugs, therefore, understanding the binding between the viral RBD and ACE2 can be very useful for rational drug design. We exploited the fact that different 3D protein structures of the ACE2 in complex with the RBD of SARS-CoV and SARS-CoV-2 are available in the Protein Data Bank: (i) SARS-CoV RBD in complex with ACE2, PDB code: 2AJF 24 (ii) SARS-CoV-2 RBD in complex with ACE2, PDB code: 6VW1 25 (iii) SARS-CoV-2 RBD in complex with ACE2, PDB code: 6M0J. 26 By using two SARS-CoV-2 samples and one SARS-CoV sample, we can effectively compare the variations among the viruses, as well as identify the important and unimportant mutations and interactions. Molecular dynamics is performed to investigate the characteristics of the binding mechanism between viral RBDs and its receptor. Sequence, structural and dynamical analyses are carried out to identify similarities and differences between the new SARS-CoV-2 viruses and the previous SARS-CoV virus. Our results provide a molecular understanding of the binding complex and show that the new viruses bind more strongly to the receptor. Despite variations between the sequences of the two SARS-CoV-2 viruses, the two samples of the new viruses behave very similarly, and have similar binding properties for the human receptor. By eliminating unimportant differences in the sequences of the new viruses, and comparing them to the sequence of the old SARS-CoV virus, four important mutations between SARS-CoV-2 and SARS-CoV are identified: Y to L455, VP to EI471-472, -PP to GVE482-484, and Y to Q498 (the indices are that of SARS-CoV-2 virus, which is shifted higher by 13 compared to SARS-CoV virus numbering). They are shown structurally to be at or close to the binding interface with the receptor protein. The mutations lead to higher flexibility of the backbones by substitution of Prolines and insertion of Glycine. The mutations also reduce number of hydrophobic interactions at the interface since hydrophobic residues are substituted with charged or polar residues. With higher backbone flexibility, new residues can move closer to the receptor interface, increasing the number of non-bonded interactions between the proteins. Specifically, the mutations eliminate one hydrogen bonding with the receptor, but increase the coordination number (number of contacts) by 20%. This leads to stronger non-bonded interactions (electrostatics and van der Waals) with the receptor compared to the SARS-CoV virus. The overall results suggest that the binding energy in the new SARS-CoV-2 viruses is slightly higher relative to the SARS-CoV virus. This is in reasonable agreement with experimental results reporting that the new corona virus has 10-20 times higher binding affinity for the ACE2 receptor than the old SARS-CoV. 27 Structurally, stronger binding is shown to not only stabilize the binding interface, but also to stabilize the overall structural fold of the viral RBD of the new SARS-CoV-2 viruses. This paper is organized as follows. After the introduction in Section 1, the detail of the computational procedure is presented in Section 2. In Section 3, the results are presented and discussed. We conclude in Section 4. The RBD -ACE2 complexes were obtained from the RCSB PDB database with ID: 6VW1, 25 6M0J 26 for the new SARS-CoV-2 viruses, and 2AJF 24 for the SARS-CoV virus. For easy identification, these systems are named 6VW1, 6M0J and 2AJF respectively. The 6VW1 and 2AJF systems contain two complexes, chains A and E, chains B and F. The 6M0J contains only one complex, chains A and E. We keep only the complex of chain A (ACE2 receptor) with chain E (viral RBD) in these samples to better compare them to one another. This is also understandable since the complex of chain B and F is a near 180 • rotation symmetry of the complex of chain A and E. It is an artifact of protein crystallization, not actual arrangement of the complexes in real biological system. There are some notable structural elements in these complexes that need special care when setting up the simulation systems such as: the disulfide bonds between various pairs of Cysteine residues; the N−linked glycosylations of various Asparagine residues; and the ZnGlu 2 His 2 zinc−finger structure that needs proper protonation states of the associated amino acids. 28 Details of these structural elements are listed in Table 1 for better identification in later discussion of the simulation results. . Table 1 : Detail information of various special structural elements of the SARS-CoVs RBD -ACE2 complexes that need proper care when setting up the system. For alignment of the primary sequences of the proteins, ClustalW web-server 29 with BLO-SUM matrix, 30 BioEdit software 31 are used. Missing residues in the central protein structure were added using homology modeling method. [32] [33] [34] The initial systems for simulation were created using CHARMM-GUI web-server 35 and then modified manually to properly describe some of above structural elements. Molecular Dynamics (MD) simulations 36 are performed on the systems using the GROMACS/2018.6 software package. 37 Charmm-36 force-field 38 was chosen for parametrization of the proteins and ions. GLYCAM06 force-field 39 was used to parametrize the glycans and TIP3P 40 model is used for water in the explicit solvent simulation model. After solvation, sodium and chlorine ions are added to the system to neutralize the total charges and to set physiological electrolyte concentration of solution at 150 mM NaCl. The simulation box size was chosen so that the proteins in neighbor periodic box are at least 3 nm apart from each other. Since the electrostatic screening length at 150mM NaCl concentration is about 7Å, this 3 nm distance is more than enough to eliminate the finite size effect due to long range electrostatic interactions among proteins in neighboring simulation boxes, yet small enough to keep the size of the system manageable with our current computational resources. The total numbers of molecules, residues and atoms for the three systems simulated are listed in Table 2 : The temperature of 310 K and the pressure of 1 atm Standard analyses such as mean deviation, mean fluctuation, hydrogen bonding are performed using GROMACS software. Visual Molecular Dynamics (VMD) program 42 is used for visualization, spatial inspection and some surface analyses. Some collective variable analyses are performed using PLUMED software. 43 In-house Python scripts are used to integrate various tools for Principal Component Analysis (PCA), setting up the systems, process simulation data, graphical plottings. SciPy python package is used for calculating probability density, Gaussian quadrature from which the difference in the binding free energy of the proteins is estimated. Before going into detailed discussion of simulation results, let us analyse the sequences and structures of the three RBD−ACE2 complexes of the experimental samples, 2AJF for SARS-CoV virus, 6VW1 and 6M0J for SARS-CoV-2 viruses. This helps to clarify common features as well as differences between SARS-CoV and SARS-CoV-2 viruses. This information will be valuable for later interpretation and discussion of results of our MD simulations. In the experimental X−ray crystal structures of 2AJF, 24 6VW1 25 and 6M0J, 26 sequences of the human ACE2 receptor are identical but glycosylations of the human ACE2 protein in these complexes show some small variations (see Table 1 ). For the 2AJF sample, the N-link glycosylations occur at the four Asparagine residues, N53, N90, N322, and N546. They link to chains of two or three N-Acetyl-D-Glucosamine (NAG) and may end with β-D-Mannose (BMA) molecules. For the 6VW1 sample, the same residues are also N-link glycosylated, albeit with slightly different length of the oligosacharides. Additionally, residue N103 is also glycosylated with a single NAG residue. For the 6M0J sample, all three Asparagine residues, N90, N322, and N546 are N−link glycosylated with a single NAG sugar molecule. The N−link glycosylation plays an important role in the protein structure stability. Despite different number of glycosylated sites in studied ACE2 proteins, we will see later in our MD simulation results that the three ACE2 proteins have nearly identical structural and dynamical properties in low energy modes. Only the 6M0J sample with less glycosylation degrees shows some differences at localized high energy dynamical modes, but these differences do not influence the binding interface of the viral RBD and the ACE2 receptor. The viral receptor-binding domain Figure 2 : Sequence alignments of the viral RBD of two SARS-CoV-2 stuctures (6M0J and 6VW1) and the SARS-CoV structure (2AJF). The *, :, or . notations below the sequence denote a fully, strongly, or weaker conserved mutation respectively. A blank space below the sequences means a non-conserved mutation. Different mutations emerge between the compared RBDs, but four groups of mutations (highlighted in yellow) show aminoacids fully conserved between SARS-CoV-2 samples and non-conserved in the SARS-CoV sequence (see text for more discussion). The grayed−out residues at the end of the sequences are flexible N-and C-termini loops that are missing in the PDB crystal structure. Unlike the human receptor, the sequences of viruses' receptor-binding domain (RBD) of the coronaviruses in all three structures, 2AJF, 6VW1 and 6M0J are different to varying degrees of conservations. In Fig. 2 , these three sequences of the viral RBD are aligned using ClustalW web-server. 29 The sequence identity between the two SARS-CoV-2 coronaviruses is 85%. The sequence identity between SARS-CoV (2AJF) with two SARS-CoV-2 samples (6M0J and 6VW1) are 72% and 85%, respectively. Thus, these RBD of coronaviruses are very similar in sequence and structure. This is expected because they are closely related in the coronavirus family. Since we are looking for major differences between SARS-CoV and SARS-CoV-2, we only look for aminoacids that are the same in the two SARS-CoV-2 viruses but are not-conserved with respect to SARS-CoV virus. Using these criteria, we identify four locations of important mutations. These are highlighted in yellow in Fig. 2 . Remarkably, upon spatial inspection, all these four locations are at or close to the binding interface with the human ACE2 receptor. This is a strong hint that these mutations can be critical for the stronger receptor binding of the new viruses. We will focus more on these mutations in later discussions of results of MD analyses. Before discussing the MD results, a quick scan on the physical changes at some of these mutations can already be made: • Y to L455: this mutation in SARS-CoV-2 leads to the loss of the aromatic ring and especially the hydroxyl group of Tyrosine amino acid. So it removes the possibility for hydrogen bond with the receptor counterpart residue Glu35. In MD analyses, we will see that the loss of hydrogen bond binding at this location is compensated with other non-bonded interactions. • VP mutated to EI at index 471-472 of SARS-CoV-2: this mutation substitutes the rigid Proline residue and adds the negatively charged Glutamic Acid E471. This increases the flexibility of the protein backbone. • -PP to GVE at index 482-484 of SARS-CoV-2: The insertion of Glycine and the substitutions of two Proline amino acids clearly make this segment much more flexible. This flexibility allows the segment to move closer to its receptor and enables the negatively charged E484 residue to form electrostatic bond with the positively charged residue K31 of the human ACE2 receptor. • Y to Q498 mutation in SARS-CoV-2 turns the T-stacking hydrophobic interaction between the two tyrosines at the RBD-ACE2 interface into electrostatic interactions with the opposing Q42 residue of the ACE2 receptor. In summary, these mutations increase the flexibility of the backbone atoms at the viral binding interface and they add to this interface polar or charged residues. Higher flexibility allows them to adapt better to the receptor binding interface, increasing the number of contacts between the two proteins. The increase in the non-bonded electrostatic and van der Waals interactions then leads to higher binding affinity as found later by MD simulation. Residues at N-and C-termini of SARS-CoV-2 viral RBD Fig. 2 shows, the viral RBD of SARS-CoV-2 samples have some additional residues at the Nand C-termini compared to the SARS-CoV sample. However, Fig. 3 shows that these RBD residues are located far away from the binding interface with the human ACE2 receptor. Because our main goal is to discern the differences between SARS-CoV and SARS-CoV-2 receptor binding, whenever appropriate when aligning the structures for analyses, we use only the common "core" part of the three viral RBD proteins (from amino acid sequence CPFGE to VLSFE in Fig. 2 ). Other differences among the viral RBD structures are the missing residues at the N-and C-termini. In Fig. 2 , the grayed−out residues at the two ends of SARS-CoV-2 sequences are missing in the protein experimental crystal structures since the precise 3D structures of the termini loops cannot be determined by crystallography because of their flexibility. To find out if these residues may play an important role for the stability of the structure, we simulate two different 6VW1 systems, the first with original PDB structure and the second where the missing termini residues are added to the PDB structure by homology modelling. [32] [33] [34] Both systems are simulated for 200 ns and statistics are collected from 50 ns onward. By comparing to the systems with and without added missing residues, we see no noticeable differences in the structures of the binding interface. In Fig. 4 , the mean fluctuations of C α atoms in the core region of the viral RBD are plotted for the systems with or without terminal residues added (blue and red curves, respectively). As one can see, the difference Figure 3 : Compared to the experimentally-resolved viral RBD structure of SARS-CoV, the SARS-CoV-2 viral RBD structure has additional residues at the N-and C-terminals. However, these residues (shown in sticks) are far from the binding interface and will not be used for alignment of structures. In this figure, the 6VW1 sample is shown. The 6M0J structure is very similar. Figure 4 : Root-mean-square fluctuations of the "core" region of the viral RBD structure of the 6VW1 complex. The blue curve is for the complex where missing residues at the termini are added by homology modeling method. The red curve is for the original complex. among the systems is minimal, mostly less than 0.2Å. Quantities such as radius of gyrations, moment of inertia along principle axes of systems with and without missing residues also show identical values within standard uncertainty as shown in Table 3 . Thus, one can conclude that adding missing terminal residues does not alter significantly the structure of the core binding region of the viral RBD. On the other hand, since these flexible and disordered terminal loops extend far into the aqueous solution, they significantly increase the size of the simulation box. As a result from these considerations, for the rest of this work, we simulate only the residues present in the crystal structure without the missing residues at the termini to save computational resources. The ACE2 protein of the 6VW1 system misses a D615 residue at the C−terminal compared to the ACE2 protein in the 2AJF and 6M0J systems. Since this is just one residue, costing minimal computing resource, we add this missing residue in the simulated systems. The three systems simulated thus have identical ACE2 protein sequences. Most of the differences are due to changes in the viral RBD sequences. Let us now present and discuss the results of our molecular dynamics (MD) simulations of the three systems and compare the differences between SARS-CoV and SARS-CoV-2 samples. As a standard procedure, one first looks at the deviations of the structural proteins from their native experimental structures. In Fig. 5 , the root-mean-square deviations (RMSD) of backbone atoms in the three systems are plotted as function of time. Fig. 5(a) is plotted for the human ACE2 receptor, while Fig. 5(b) is plotted for the viral RBD protein. The viral RBD, being smaller (∼ 200 residues) shows faster saturation time (at around 50ns) than the ACE2 receptor (∼ 600 residues) whose saturation time is at around 100 ns or longer. Among SARS-CoV versus SARS-CoV-2 viruses, the viral RBD is remarkably more stable in the new viruses as they both deviate only 1.5Å from the native structure while that of the old SARS-CoV virus shows about 3Å deviation. Similarly to the RBD, the human ACE2 in the SARS-CoV-2 systems are more stable and show lower RMSD values than in the SARS-CoV complex. The ACE2 receptor in 2AJF system also indicates longer relaxation time. Nevertheless, both proteins in all systems are stable with RSMD less than 4Å deviation from native structure. As a result of this RMSD analysis, in all subsequent equilibirium statistical analyses, the first 100 ns of the MD trajectories will be dropped. Next, let us look at the root mean square fluctuations (RMSF) of the atoms. This fluctuation is due to thermal effect of finite temperature (at 310 K) compared to the low temperature of experimental crystallized structure where proteins are frozen. It is a good measurement of the backbone stability of the proteins. In Fig. 6(a) , RMSF of the backbone C α atoms of the human ACE2 receptor in the three systems are plotted for individual residues. One can see from this figure that this receptor behaves similarly across the three systems, with variations less than 0.5Å from one system to another. Our MD simulations show as well that the different ACE2 glycosylation states (see Table 1 for more details) have minimal influence on the rigidity of the protein backbones. One noteworthy difference is the slightly higher fluctuations of ACE2 residues 134-177, 262-278, 489-500 in the 6M0J system compared to the other two. As it turns out, these residues surround a Cl − ion in the 2AJF and 6VW1 systems (see Fig. 6(b) ). For unknown reasons, this chlorine ion is absent in the experimental crystal structure of 6M0J system. These residues are at the other extreme of (a) (b) Figure 5 : The root-mean-square deviation of backbone atoms (backbone RMSD) of the human ACE2 receptor (a) and the viral RBD protein (b) in MD simulations from the initial crystal structure. Both systems of SARS-CoV-2 (green and cyan curves) show smaller deviations with faster time to saturation than that of SARS-CoV system (red curve). Both proteins are stable with deviations less than 4Å from the native structure. the ACE2 receptor, far from the binding interface, so we believe that the slightly higher fluctuations in this region have no important consequences. (a) (b) Figure 6 : Root-mean-square fluctuations of the backbone C α atoms of the human ACE2 receptor as function of the residues along the protein backbone (a) for different viruses. Although very similar to each other, the 6M0J system shows slightly higher fluctuations at certain ranges of residues. Figure (b) show these residue ranges in blue color. They are all located near a Cl − ion (shown in green). See text for more details. The RMSF for the viral RBD core region is plotted in Fig. 7 . Remarkably, Figure 7 also shows that, at the N-terminal part, around the indices 340, Figure 7 : Root-mean-square fluctuations of the backbone C α atoms of the viral RBD as function of the residues along the protein backbone for different viruses. Only the core binding region is shown. The residues are aligned and shifted according to Fig. 2 for better comparison between the two viruses. 365, and 390, the SARS-CoV sample is much less stable than the two SARS-CoV-2 samples. Close inspection of the secondary structure in this region (not shown) show that the various β−sheets here are longer in SARS-CoV-2 viruses than those of SARS-CoV virus, indicating a more ordered and compact structure. Stronger receptor binding seems to stabilize not only the ACE2-RBD binding interface (see the next Section), but also the overall structure of viral RBD of the new viruses. After performing several sequence and structural investigations, let us move to the interaction picture among amino acids at the binding interface of the viral RBD and its human ACE2 receptor. In Fig. 8(a) , the distribution of the number of hydrogen bonds between the viral Fig. 2 , it is observed that the first mutation from Tyrosine to Leucine turns off hydrogen bonds at this location due to the lack of side chain hydroxyl (−OH) group in Leucine as expected. It should be noted that the number of hydrogen bonds can only be accurately determined from a full quantum mechanical calculation. However, such calculations are beyond our computing capacity for such a large system. Here, simple geometric criteria are used to determine the hydrogen bonding, namely the distance between the donor and acceptor atoms is less than 3.5 Å, and the angle of the three atoms making up the hydrogen bond is less than 30 • . This may overestimate the total number of hydrogen bonds since one atom can participate in more hydrogen bonds at once. Nevertheless, even with this classical definition, the loss of one hydrogen bond due to this mutation is quite clear in our figure. In addition, in Fig. 8(b) , the coordination number between the C α atoms of the viral RBD and its receptor is plotted for the three systems in order to estimate the number of interactions between the studied RBDs and ACE2 receptor. Here, the PLUMED driver version 2.5.1 43 is used for the calculation, for which the default contact function for a pairs of atoms i and j is given by: Here, r ij is the distance between the atoms, r 0 is the neighbor cutoff distance (set to 8Å in our calculation), n = 6, and m = 12. The meaning of the contact function is similar to inverted step function: when r ij < r 0 , s ij 1, when r ij > r 0 , s ij 0. Therefore, the contact function denotes whether the two atoms i and j are within distance r 0 from each other. The Principal component analysis and the free energy landscape in col- Principal component analysis is a useful method to analyze dynamics of proteins. 44 By calculating the covariance matrix among the atoms, and diagonalize the resultant matrix we obtain eigenvectors of dynamical motions of the protein (similar to the concept of normal modes). This analysis is very helpful for biomolecules since the long−time dynamics will be dominated by only a few lowest energy modes. All the strong oscillatory high energy modes are screened out efficiently using this analysis, leading to a huge reduction in dimensionality of the systems. By our own inspection, three lowest modes are enough to show distinct clusters of configurations for the viral RBD. In Table 4 , various information about the dynamics of the proteins using PCA method are listed. Here, the P max (a i , a j ) is the maximum probability density in the plane of the projections (a i , a j ) on the principal eigenvectors i-th and j-th respectively. The trace of the covariance matrix is a measure of the overall flexibility of the proteins. As one can see from this table, the trace for ACE2 for SARS-CoV and SARS-CoV-2 6VW1 are large and similar to each other. This is expected, since the human ACE2 receptor is a big protein, with same sequence and almost identical native structures in these systems. The ACE2 protein of SARS-CoV-2 6M0J sample shows even a higher value. This may be related to the missing Cl − ion in the experimental structure, and less degrees of glycosylations of ACE2 protein in this sample (see the preliminary Section). The difference between SARS-CoV and SARS-CoV-2 viruses are apparent from the PCA information for the viral RBD domain. Here, the SARS-CoV 2AJF shows much larger trace of 1.8 nm 2 while both SARS-CoV-2 samples shows a trace of only 1.2 nm 2 . This is again in line with previous analyses that SARS-CoV-2 viral RBD are more stable. In Fig. 9 , the normalized histogram of the projections of the protein structure on their first two lowest energy modes are shown. The projection on the third mode follows a simple Gaussian distribution (not shown). This map is the probability density of the configuration Figure 9 : The probability density in the plane of the projections of the structure of the backbones of proteins on its two lowest energy modes. Top row is for the ACE2 receptor and bottom row is for the viral RBD. Note that, the colorbar scale is different for ACE2 versus viral RBD. See text for more discussion. of the protein in the phase space spanned by these collective variables. Differences in the probability density for the viral RBD are once again significant. The SARS-CoV-2 virus systems show two smooth, Gaussian-like clusters, while that of SARS-CoV shows weak and scattered profiles. The maximum probability densities in Fig. 9 together with projection on the third eigenvector are listed in Table 4 for both ACE2 and viral RBD backbone structures. Using the relation, ∆G ∼ −k B T ln p, these result suggest that the depth of the free energy landscape of RBD for the SARS-CoV-2 viruses is lower by k B T ln(p SARS2 max /p SARS max ) = 0.5 ± 0.1 kcal/mol. Clearly, the lower fluctuations of the backbone of the viral RBD in the SARS-CoV-and are beyond the scope of this work. Here, by many different sequence, structural and dynamical analyses, we have shown that the new SARS-CoV-2 viruses have stronger binding affinity for the human ACE2 receptor than SARS-CoV virus. In this paper, we present a MD study of the viral receptor-binding domains from SARS-CoV-2 and SARS-CoV viruses in complex with the human ACE2 receptor. By using various sequence, structural and energetic analyses, it is shown that the SARS-CoV-2 viruses bind stronger to its receptor than the SARS-CoV virus. The stronger binding not only creates more stable binding interface, it can also causea more ordered structure of the viral RBD in SARS-CoV-2 virus. Despite very similar sequence identities, four important mutations in the new viruses play important role in this stronger binding process, namely Y to L455, VP to EI471-472, -PP to GVE482-484, and Y to Q498 (the indices are that of SARS-CoV-2 virus). These mutations all lead to a reduced internal rigidity of the viral protein backbone near its binding interface and add polar and charged residues to the interface. The higher flexibility in the backbone of the new viruses then allows it to move closer to the human receptor protein surface, and to bind stronger to the receptor using non-bonded interactions. Overall, the coordination number between the proteins increases by 20% in the cases of the new SARS-CoV-2 viruses. This increase in different non-bonded interactions (electrostatic, van der Waals, hydrophobic) overcome the loss of one hydrogen bond between the proteins, leading to a stronger binding free energy. There are some variations in the primary sequences of the two samples of SARS-CoV-2 viruses, with about 85% identity. However, our results show that these variations do not lead to significant differences in the physical properties of the binding complex. Both samples of the new viruses show very similar behaviours. These mutations are therefore non critical for the binding process. Only the non-conserved mutation differences with SARS-CoV are important. Preliminary estimate of the binding free energy agrees reasonably with the experimental data that SARS-CoV-2 viruses show 10 to 20 folds higher binding affinity for ACE2 receptor as compared to SARS-CoV virus. A more thorough sampling method such as metadynamics, umbrella sampling or steered molecular dynamics is requested for more quantitative comparison. Similarly, the presence of several other mutations in the sequence of new viruses, its variations, stability/conservation are other interesting aspects that needs more thorough investigations. However, all these studies go beyond the current scope of this paper and will be presented in a near future work. Within this work, our sequence, structural and dynamical analyses have provided strong support that SARS-CoV-2 viruses bind stronger to the human ACE2 receptor than SARS-CoV virus. Our molecular description of the virus-receptor binding interface, particularly the insight we obtained in the role of mutations, can help in the development of vaccines and antiviral drugs by rational drug design addressing the important needs in ongoing pandemics. 2020) A novel coronavirus outbreak of global health concern A pneumonia outbreak associated with a new coronavirus of probable bat origin A new coronavirus associated with human respiratory disease in China Clinical virology manual Mechanisms of coronavirus cell entry mediated by the viral spike protein Dissecting virus entry via endocytosis Structure, function, and evolution of coronavirus spike proteins Pre-fusion structure of a human coronavirus spike protein Cryo-electron microscopy structure of a coronavirus spike glycoprotein trimer Architecture of the SARS coronavirus prefusion spike Conformational states of the severe acute respiratory syndrome coronavirus spike protein ectodomain Angiotensinconverting enzyme 2 is a functional receptor for the SARS coronavirus Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor Receptor recognition by the novel coronavirus from Wuhan: an analysis based on decade-long structural studies of SARS coronavirus Characterization of the receptor-binding domain (RBD) of 2019 novel coronavirus: implication for development of RBD protein as a viral attachment inhibitor and vaccine SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor Sound Science before Quick Judgement Regarding RAS Blockade in COVID-19 Interaction between RAAS inhibitors and ACE2 in the context of COVID-19 2020) COVID-19 and the cardiovascular system Are patients with hypertension and diabetes mellitus at increased risk for COVID-19 infection? SARS-CoV2: should inhibitors of the reninangiotensin system be withdrawn in patients with COVID-19 The proximal origin of SARS-CoV-2 Structure of SARS coronavirus spike receptor-binding domain complexed with receptor Structural basis for receptor recognition by the novel coronavirus from Wuhan Crystal structure of the 2019-nCoV spike receptor-binding domain bound with the ACE2 receptor Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation Overcharging of the Zinc Ion in the Structure of the Zinc-Finger Protein Is Needed for DNA Binding Stability Multiple sequence alignment using ClustalW and ClustalX Using BLOSUM in sequence alignments BioEdit: an important software for molecular biology Comparative Protein Modelling by Satisfaction of Spatial Restraints Homology modeling Homology modeling of mouse NLRP3 NACHT protein domain and molecular dynamics simulation of its ATP binding properties CHARMM-GUI: a web-based graphical user interface for CHARMM Molecular dynamics simulations GROMACS: a messagepassing parallel molecular dynamics implementation CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data GLYCAM06: a generalizable biomolecular force field Hydrophobic solvation of methane and nonbond parameters of the TIP3P water model LINCS: a linear constraint solver for molecular simulations VMD: visual molecular dynamics PLUMED 2: New feathers for an old bird Three-dimensional cluster analysis identifies interfaces and functional residue clusters in proteins11Edited by Energy Landscape and Transition State of Protein-Protein Association Escaping free-energy minima Chasing the Full Free Energy Landscape of Neuroreceptor/Ligand Unbinding by Metadynamics Simulations Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling Nonequilibrium Equality for Free Energy Differences