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.

Recent studies have pointed the role of angiotensin-converting enzyme-II (ACE2) in mediating the entry of SARS-CoV-2 to the host cell by binding to the receptor-binding domain (RBD) of viral spike protein, and successive priming by cellular proteases initiates the infection. SARS-CoV replication rate and disease severity is controlled by the binding affinity of RBD with ACE2. To understand, how mutations in the conserved residues of RBD affect the molecular interaction with ACE2, we generated five alanine mutants i.e. Y449A, N487A, Y489A, N501A and Y505A in the receptor binding motif (RBM) of the ACE2-RBD SARS-CoV-2 complex (PDB: 6M0J). Computational site directed mutagenesis induced dynamics in wild-type and mutant complexes were extensively studied through all-atoms molecular dynamics (MD) simulations of 150 ns. In silico mutational analysis revealed loss of important intermolecular hydrogen bonds and other non-bonded contacts, critical for molecular recognition of SARS-CoV-2 RBD to ACE2, which is well supported by saturation mutagenesis study of binding interface residues. MD simulations results showed that RBM motif is flexible, where mutant residues are relatively more mobile than corresponding wild-type residues. Global motion analysis through principal component studies revealed that RBD exhibits protuberant in-ward motion towards the human ACE2 binding interface which may be crucial for molecular interaction. Conclusively, the present finding are in congruence with previous experimental reports and provides detailed information on the structural basis of receptor binding by human SARS-CoV-2, which will crucial for the development of novel inhibitors or drugs to combat against SARS-CoV-2.Communicated by Ramaswamy H. Sarma.

Introduction

During the end of 2019, multiple patients with symptoms of novel coronavirus designated as 2019-nCoV (now termed as SARS-CoV-2) infection were reported in relation to the Huanan Seafood Wholesale Market in Wuhan, South China (Wu et al., 2020). Within a short span of time, the virus SARS-CoV-2 has spread across the globe affecting more than 250 countries and according to the latest report, it has infected more than 10.5 million people and caused deaths of more than 511,851 individuals worldwide. The high fatality rate of SARS-CoV-2 has raised global concerns and now declared as a global pandemic by the WHO. The SARS-CoV-2 impose a continuous threat to human health as well as the economy as they emerge periodically, unexpectedly, have a broad host range, spread easily, and lead to catastrophic consequences (Li et al., 2020).

Coronaviruses are positive-sense single stranded RNA viruses with the largest genome among all RNA viruses, ranging from 26 to 32 Kilobases in length (Malik et al., 2020). The viral genome is packed inside a helical capsid made up of nucleocapsid protein and further surrounded by an envelope which is formed of three major structural proteins, namely spike protein, small envelope protein, and membrane protein (Li et al., 2020). These four important glycoproteins majorly contribute to the structure of all coronaviruses and play vital roles in pathogenesis. Particularly, the spike protein mediates coronavirus entry into the host cells by binding to host cell receptor which triggers a cascade of events leading to fusion of both viral and host membranes (Li, 2016). Additionally, the spike protein is a critical determinant of wide host range and viral tissue tropism.

Electron microscopy studies have revealed the three domain (i.e. N-terminal large ectodomain, middle transmembrane domain and C-terminal short intracellular tail) architecture of the spike protein of coronavirus. The ectodomain is comprised of two functional subunits i.e. S1-subunit having a receptor binding domain (RBD) and a membrane-fusion S2-subunit. The spike protein mediates viral genome to enter into host cell by first binding to a host receptor through the RBD domain in the S1 subunit and subsequently fusing the host and viral membranes with the help of the S2 subunit (Li et al., 2005; Liu et al., 2004; Tai et al., 2020).

Sequence analysis of RBD region of SARS-CoV-2 and SARS-CoV spike proteins showed high conservation and homology with no deletion or insertion except for one amino acid insertion on a loop away from the functional receptor binding domain, speculating that SARS-CoV2 shares the same receptor, ACE2 (angiotensin-converting enzyme 2) as SARS-CoV (Wan et al., 2020). Consistent with high sequence and structural similarity, protein-protein binding assays corroborated the important role of ACE2 in mediating SARS-CoV-2 viral attachment and entry (Tian et al., 2020; Walls et al., 2020; Zhou et al., 2020). A recent report on the structure of SARS-CoV-2 spike showed that the RBD domain of spike protein binds to the peptidase domain of ACE2 (Wrapp et al., 2020). It has been also observed that host susceptibility to SARS-CoV infection is primarily determined by the affinity between the viral RBD domain of spike protein and the host receptor ACE2, and there are some specific amino acid residues that are potentially involved in the interaction and viral binding (Huentelman et al., 2004; Wan et al., 2020; Zhang et al., 2005). Recently, various computational studies using sequence alignment and phylogeny methods have reiteratively applied this prior structural knowledge on SARS CoV-ACE2 binding to map the similar residues at equivalent positions that are putatively involved in SARS CoV-2 and ACE2 interaction and the likely host range of SARSCoV-2 based on the properties of the conserved amino acids located in the biding interface (Othman et al., 2020; Wan et al., 2020; Yan et al., 2020).

Over the last six months, a number of cryo-EM and crystal structures of important drug-targets of SARS-CoV2 have been reported in RCSB PDB (Protein Data Bank). Recently, a crystal structure of the SARS-CoV-2 spike RBD bound to the ACE2 cell receptor (PDB: 6M0J) at 2.45 Å resolution has been resolved (Lan et al., 2020). Furthermore, Yan et al. (2020) have reported the full length cryo-EM structure SARS-CoV-2 RBD-ACE2 complex (PDB ID; 6M17). These structural data supported the earlier findings stating that the SARS CoV-2 RBD and SARS CoV RBD shares nearly identical binding interface and the underpinning residues are highly conserved with similar side chain properties (Lan et al., 2020).

The highly conserved receptor binding motif of RBD of SARS CoV-2 is crucial for mediating interactions with host cellular receptor ACE2. However, how mutation in these conserved residues of RBD influence the molecular recognition and the underlying effect of these mutations on the structure, dynamics and interactions with ACE2 is not well understood. Molecular dynamics simulations have been a pivotal approach to investigate mutation-induced changes in protein structure, function and molecular interaction at atomic scale (Dehury et al., 2015, 2017; Pan et al., 2019, Dehury et al., 2014). Using the experimental ACE2-RBD complex (PDB: 6M0J), we combined site-directed computational mutagenesis along with all-atoms molecular dynamics simulation to provide more insights into the structure of the SARS-CoV-2 spike RBD and ACE2 complex. The key non-bonded interactions mediated by RBD of SARS-CoV-2 with ACE2 interactions have been explored comprehensively, and enlightened how mutation in the receptor binding motif affects the dynamics of recognition. The results from our study provides depth understanding of the stability and flexibility of RBD spike protein binding to host receptor at atomistic level which could be useful for structure-based design of potent drugs or inhibitors to combat novel coronavirus SARS-CoV-2 infection.

Materials and methods

System preparation

In SARS-CoV-2, viral entry mediated by the trimeric spike protein is most important step in its lifecycle and represents an attractive therapeutic intervention point by blocking the virus-cell membrane fusion event or the co-receptor interactions. SARS-CoV-2 and other human coronaviruses have similar infection mechanisms which share the same human ACE2 receptor crucial for viral entry. Therefore, the experimental crystal structure of ACE-RBD complex (PDB ID: 6M0J) was used as the reference structure for our study. The missing side chains in the complex were modelled using WHAT-IF server (Hekkelman et al., 2010). Subsequently, various mutants were generated using site directed mutagenesis protocol in PyMOL version 2.3 (The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC.). In this study five mutations were introduced in three conserved and two variable sites of RBD of SARS CoV-2 randomly. The resultant mutant models were cross-checked using BIOVIA Discovery Studio Visualizer version 4.5 (Dassault Systèmes BIOVIA, BIOVIA DSV, 4.5, San Diego: Dassault Systèmes, 2019) to remove steric clashes. Finally, a total of six RBD-ACE2 complexes including the wild-type, were then subjected to all-atoms MD simulations using GROMACS version 2018.4 (Abraham et al., 2015).

MD simulation

The six systems were simulated using CHARMM36 force field (Huang et al., 2017) using TIP3P water model. Each complex was solvated by ∼76,666 (approximately) waters and neutralized a strength of 0.15 M counter Na+ and Cl− ions in a cubic box of the initial side length of 90 Å (for details see Table S1). After electro-neutralization, each system was subjected to energy minimization using the steepest descent algorithm until the maximum force of 1000 kJ mol−1 nm−1 has been achieved. NVT ensemble via the Nose-Hoover method at 303.15 K was employed to maintain the temperature of each system, while the pressure was equilibrated using an NPT ensemble at 1 bar with the Parinello–Rahman algorithm. Finally, a production runs of 150 ns were performed for each system starting from different random initial velocities. The non-bonded interactions were treated using Verlet cut-off scheme, where, the particle mesh Ewald (PME) method was employed to treat long-range electrostatic interactions, while the short-range electrostatic and van der Waals interactions were calculated with a cut-off of 12 Å. Periodic boundary conditions were applied to all simulations, and bonds involving hydrogen atoms were constrained using the linear-constraint-solving (LINCS) algorithm. Trajectory data were saved at time intervals of 200 ps.

MD trajectory analysis

The molecular dynamics stability parameters were analysed using toolkits of GROMACS and VMD version1.9.1 (Humphrey et al., 1996). The 2D plots were generated using Grace Version5.1.23 (http://plasma-gate.weizmann.ac.il/Grace/) program. PyMOL, BIOVIA DSV and LigPlot+ were used for ACE2-RBD interaction analysis and image rendering. Principal component analysis (PCA) was used to determine the correlated motions of the residues to a set of linearly uncorrelated variables called principal components (David & Jacobs, 2014; Dehury et al., 2014; 2015, 2017). The positional covariance matrix main chain coordinates and its top two eigenvectors (EVs) were used. This method is based on the construction of the covariance matrix of the main chain-atoms of coordinate fluctuations of the simulated complexes. To explore the conformational heterogeneity in the ensemble of ACE2-RBD complex structures i.e. conformations that visited most frequently along the trajectory, we employed clustering approach employed in GROMOS clustering algorithm with a cut-off of 0.15 nm to extract representative structure. PCA and clustering were performed by using the last 100 ns MD trajectories. To determine the evolution of the secondary structural elements in the complex during MD simulations, the gmx do_dssp program was used. Number of distinct hydrogen bonds formed between specific amino acids residues between ACE2 and RBD was determined utilizing the gmx hbond utility with the donor-acceptor set at a maximum of 0.35 nm.

Electrostatic potential surface

The electrostatic potential surfaces were computed using the APBS-Electrostatics plug-in tool of PyMOL. The linearized Poisson–Boltzmann (PB) equation with a bulk solvent radius of 1.4 Å and a dielectric constant of 78 was used to compute the solvent-accessible surface area of the Wt and Mt type representative ACE2-RBD cluster structures. The electrostatic positive and negative surfaces in each structure were viewed using a contour (kT/e) value of 1.

Results and discussion

Sequence-structure similarities of RBDs

The crystal structures of ACE2-RBD complex from SARS-CoV-2 (PDB ID: 6M0J) and ACE2-RBD from SARS-CoV (2AJF) were used as reference framework in the present study for further in depth structural analysis underpinning molecular interactions. The sequence-structure alignment of the RBD regions in SARS-CoV and SARS-CoV-2 spike protein was performed to identify the major similarities and differences between them (Figure 1(A)). Among the variable region, residues falling in beta sheets (β5-β9) exhibited maximum degree of divergence. The electrostatic surface potential of both the complexes displays striking resemblance in the charge distribution of amino acids at the receptor binding motif (Figure 1(B)) in both the RBDs from SARS CoV-2 and SARS CoV. The overall Cα-of RBDs of SARS-CoV-2 with SARS-CoV revealed high conservation with RMSD value of 0.455 Å, whereas the overall Cα-RMSD of whole ACE2-RBD protein complex is 0.458 Å which signifies that the RBDs of SARS-CoV-2 and SARS-CoV binding fashion is quite similar (Figure 1(C)). Sequence-structure analysis of the residues at the ACE2-RBD interface for SARS-CoV-2 and SARS-CoV were found to be conserved with least variation (Figure 1(A)), and particularly, the hydrophobic residues known to be indispensable for interactions were highly conserved (Lan et al., 2020; Li et al., 2005). The binding site residues responsible for interactions in both the RBDs are labelled in stars (Blue: SARS-CoV and Green: SARS-CoV-2) where majority of the amino acids are either are highly conserved or share similar side chain. The receptor-binding motif (RBM) of RBD prefers to binds to the external surface of the claw-like structure of ACE2 (Figure 1(C)) and share similar side chain properties with those in the SARS-CoV RBD. A number of natural mutations occur near RBM region of RBD spike protein widely termed as ‘hot spots’ as it is mainly responsible for viral-host binding, are presumed to determine the host range of SARS-CoV. In SARS-CoV several amino acids i.e. Asp442, Leu472, Asn479, Asp480, and Thr487 underwent natural selections, which are found to be critical for host receptor recognition, cell entry, and host range of SARS-CoV (Li, 2008; Wu et al., 2012). Therefore, it is crucial to understand how mutation affects the strength of affinity of these complexes which will provide structural basis mediating SARS-Cov-2 entry into human cells.

Figure 1.

Sequence-structure comparison of receptor binding domain (RBD) from SARS CoV-2 and SARS CoV and structural features of ACE2-RBD complex (6M0J: ACE2-RBD complex of SARS CoV-2 and 2AJF:ACE2-RBD complex of SARS CoV. (A) Pair-wise sequence structure alignment of RBD from of SARS Cov-2 and SARS CoV (B). Electrostatic properties of ACE2-RBD complexes using the experimental structure 6M0J and 2AJF (Blue, red, and white colors represent positively charged, negatively charged, and neutral surfaces, respectively). (C) Structural superimposed view of ACE2-RBD complexes displaying the receptor binding motif marked in pink (green: 6M0J; SARS-CoV-2 and cyan: 2AJF; SARS-CoV).

Here, based on the experimental structure of ACE2-RBD complex from SARS-CoV-2, the conserved residues of the RBM was mutated to alanine followed by molecular dynamics simulations experiments (Figures 2 and S1) to gauge the stability of ACE2-RBD residual interactions.

Figure 2.

Schematic representation of the present work where point mutations are induced at Tyr449, Asn487, Tyr489, Asn501 and Tyr505 of the RBM in ACE2-RBD (SARS Cov-2) complex and subjected to all-atoms MD simulation.

Molecular dynamics

MD approaches have been proven useful to predict how biomolecules will respond to changes such as mutation and deciphering functional mechanisms of proteins and their interactions with other molecules in uncovering the structural basis underlying molecular mechanisms (Hollingsworth & Dror, 2018; Karplus & McCammon, 2002; Maximova et al., 2016). To understand the dynamic changes in the conformations of the wild type and mutant ACE2-RBD complexes, molecular dynamic simulation was carried out for 150 ns. The conformational changes observed during the 150 ns simulation for the ACE2-RBD complexes are discussed below. To determine the stability and mechanistic aspects of the wild type and mutant complexes, backbone RMSD and radius of gyration (Rg) of all systems were determined. RMSD is frequently used to depict the dynamic stability of systems as it measures the global fluctuations of proteins or complexes. The root mean square deviation (RMSD) of backbone Cα atoms from the starting structures of production dynamics are computed and portrayed in Figure 3(A). The RMSD of all starting configurations slightly increased during the initial equilibration phase but quickly converged after 10 ns (Figure 3(A)). Figure 3(A) shows the backbone RMSD of the wild type and mutant systems ACE2-RBD complex over the entire 150 ns of MD simulations. The Wild type and Y505A system of the ACE2-RBD complex displayed lowest RMSD of 0.25 ± 0.03 nm and 0.25 ± 0.03 respectively, followed by Y449A, N487A, Y489 and N501A mutant systems. Among the mutants systems, N501A displayed higher RMSD of 0.33 ± 0.09 nm as compared to the other complex systems. Overall most of the systems deviate to a quite similar small extent from their starting structures resulting in a backbone RMSD of ∼0.27 ± 0.05 nm during 150 ns time scale of molecular dynamics simulations. The above data indicates that all of the systems reached equilibrium after few ns of equilibration. Comparing to the wild type 6M0J complex (solved at 2.45 Å resolution), the Cα-RMSD of representative structure of each systems were found to be 1.87, 1.68, 1.66, 2.66, 1.77 and 2.28 Å respectively for Wild type, Y449A, N487A, Y489A, N501A and Y505A systems indicates that complex systems maintained their structural integrity throughout MD simulation.

Figure 3.

Stability parameters of the wild type and mutant ACE2-RBD complexes over the course of the simulation. Backbone root mean square deviation (RMSD), root mean square fluctuation (RMSF) of Cα-atoms RBD and ACE2, and radius of gyration (Rg) of each system over the time scale of 150 ns.

The root means square fluctuation (RMSF) which imitates the mobility of Cα atoms of each residue around its mean position was employed to study the dynamics flexibility of each system. The results are shown in Figure 3 for both ACE2 and RBD in wild type and mutant systems. Close observation of the RMSF plot analysis revealed that the conserved RBM (437-508 in SARS-Cov-2) is more flexible as compared to other regions of RBD. Among the most flexible regions we observed higher flexibility in the loop regions positioned at 472–490, and 495-506 near the binding site. In addition, we also observed that mutations have influence the RMSF of residues in the conserved regions of RBM. As compared to wild type (0.16 nm), the average RMSF values for the RBM site was computed to be 0.28, 0.18, 0.37, 0.40 and 0.13 nm respectively for the mutant Y449A, N487A, Y489A, N501A and Y505A systems respectively. The largest flexibility was seen at the characteristic inherent high fluctuating β-strand loop feature of RBM. It is well known fact the flexibility determines the binding as it may not only alter the binding interface between two interacting partners but also a vital contributor to the entropy penalty upon binding (Tuffery & Derreumaux, 2012). Large-scale mobility in the binding site residues implies a greater magnitude of flexibility, therefore, it can be argued that mutation in the RBM motif distress the flexibility of this binding site which in turn affect the strength of binding of RBD to host receptor ACE2.

The radius of gyration (Rg) depicts the overall compactness and dimension of protein-protein complex systems during MD simulations. It also elucidates how consistent secondary structure elements are closely packed into 3D structure of complex. The average Rg was within range of 3.17 to 3.37 nm which signifies that all systems were remained compact during MD simulation. Among the mutants, N501A displayed somewhat little higher Rg of ∼3.37 nm (Figure 3). The solvent accessible surface area (SASA) of each system was also computed which has been illustrated in Figure S2 with little variation in their solvent accessible surface area profile. The intrinsic dynamics stability statistics of ACE2-RBD complex systems during 150 ns MD has been summarized in Table 1.

Table 1.

Intrinsic dynamics stability statistics of wild type and mutant ACE2-RBD complexes during 150 ns MD (values within brackets represents the standard deviation).

System	RMSD (nm)	Rg (nm)	SASA (nm2)	RMSF
				ACE2 (nm)	RBD (nm)	RBM (nm)
WT	0.25 (0.03)	3.18 (0.02)	363.27 (4.16)	0.14 (0.04)	0.17 (0.06)	0.16
Y449A	0.26 (0.03)	3.21 (0.02)	371.37 (3.15)	0.21 (0.06)	0.25 (0.05)	0.28
N487A	0.26 (0.03)	3.17 (0.03)	365.27 (6.26)	0.14 (0.04)	0.21 (0.09)	0.18
Y489	0.30 (0.06)	3.20 (0.03)	368.92 (3.66)	0.29 (0.09)	0.32 (0.07)	0.37
N501A	0.33 (0.09)	3.37 (0.31)	369.20 (4.08)	0.34 (0.06)	0.34 (0.08)	0.40
Y505A	0.24 (0.04)	3.20 (0.02)	368.11 (4.37)	0.15 (0.04)	0.16 (0.06)	0.13

Evolution of secondary structure elements during MD

To further confirm the stability of the complex systems, we also monitored the changes in secondary structure during MD simulations using DSSP algorithm (Figure S3). The results showed that there were no noteworthy changes in structural elements i.e. helical and β-sheet content observed during the entire simulation time, which further confirmed the stability of our studied systems. To obtain the structural representative from each simulation systems, we performed ensemble RMSD based cluster analysis through a stringent cut-off of 0.15 nm. Based on the RMSD from structural ensembles, the representative snapshot was superimposed with the experimental ACE2-RBD complex (6M0J) (as shown in Figure 4) using PyMOL. As evidenced from Figure 4, it can be clearly observed that the MD snapshots superpose well with the experiential complex with a Cα-RMSD < 2.5 Å which signifies that all systems retained the structural integrity and maintain the same fold with minor changes in the RBM of SARS-Cov-2.

Figure 4.

Structural superimposition of the top ranked cluster representative of wild type and mutant ACE2-RBD complexes obtained from MD (shown in cyan) as compared to the experimental complex 6M0J (green).

Essential dynamics of ACE2-RBD complex systems

PCA or essential dynamics is one of the dimensionality reduction techniques, has been widely used to analyse molecular simulation data (Amadei et al., 1993; Kitao & Go, 1999). It also allows identifying the dominant modes of molecular flexibility in a rigorous manner, and displays the information in the form of variations in the values of a small number of collective coordinates. In order to investigate the significant motions in wild type, Y449A, N487A, Y489A, N501A and Y505A ACE2-RBD complexes, PCA was carried out. The eigenvalues obtained from the diagonalization of the covariance matrix of the atomic fluctuations of each system has been displayed in Figure 5(A). From the eigenvalue plot, it can be observed that five eigenvalues are relative to concerted motions rapidly declined in amplitude to reach a number of constrained, more localized fluctuations. PCA also suggested that first 20 principal components (PC) or EVs account for more than 87% of the total motions observed in the last 100 ns of the trajectories for ACE2-RBD complex systems. To know the global motion of these complexes (in phase space), we projected the first two EVs, EV1 and EV1 into the phase space. It can be seen from Figure 5(B) that the properties of the motions described by the first two EVs are different in mutant systems as compared to the wild type ACE2-RBD complex. In all systems, we observed least differential scattering of main-chain atoms, which postulate the occurrence of least conformational changes in the mutant complexes as compared to wild type which is in agreement with RMSD analysis. Among the mutant systems, Y489A and N501A occupied more conformational space with high trace value (trace value of 108.07 and 104.57 nm2 respectively) as compared wild type (trace value of 60.89 nm2). These higher scattering of atoms with high trace values in the mutant systems indicates overall increased flexibility than the wild type which is in agreement with thermal-fluctuation. In order to find understand the movement of ACE2 and RBD in the wild type and mutant complexes motions described by EV1 and EV2, we plotted the porcupine plots (as show in in Figures 5(C–H) and S4). Interestingly, we consistently observed prominent in-ward motion of RDB towards the ACE2 binding interface in both wild type and mutant complexes. We speculate these consistent in-ward motion displayed by the RBD towards host receptor initiatives the molecular interaction, where mutations in the conserved region of RBM affects the dynamics of these systems. While, in ACE2 we also observed significant inward motions in Y449A, Y489A and Y505A systems towards the RBD (SARS-Cov-2).

Figure 5.

Principal component analysis of the wild type and mutant RBD2-ACE2 complex systems using the last 100 ns trajectories. (A) Eigenvalues of each system (B) 2-Dimensional projection of top two eigenvectors (EV1 vs EV2) of each system in phase space. (C–H) Porcupine plot of the EV1 of each system displaying the movement of main-chain atoms where the cone represents the movements and length of cone displays the strength of movement.

Dynamics of inter-molecular hydrogen bonds

The number of inter-molecular hydrogen bonds and their residence period was monitored during the last 100 ns of MD simulations using the gmx hbond tool in GROMACS. The dynamics of hydrogen-bonds for each system are displayed in Figure 6. In wild type system, the hydrogen bonds were found to be maintained during the MD simulations, while, in mutant displayed decline in the H-bonds with least number in N487A, and Y489A systems. A number of crucial bonds are found to be broken during MD simulation of mutant systems, which signifies that replacement of the conserved amino acids with alanine results in unfavourable interaction at the RBD-human ACE2 interface thereby affects the stability of inter-molecular H-bonds. The crucial information on hydrogen bonds and other non-bonded contacts has been discussed elaborately in the subsequent section.

Figure 6.

The inter-molecular hydrogen bond dynamics of wild type and mutant RBD-ACE2 complexes during the last 100 ns MD.

Inter-molecular contact analysis

To elucidate how mutation affects the SARS-CoV-2 RBD and ACE2 interaction at molecular level, we chose to compare the structural representative of each simulation and compared the experimental complex structure of SARS-CoV-2 RBD bound with ACE2, 6M0J (Table 2). In addition, we also computed the electrostatic surface of representative structure of wild type and mutant type and compared with experimental complex (as shown in Figure S5). The SARS-CoV-2 RBM comprised of β5-β6 strands, α4-α5 helices and loops containing most of the contacting residues of SARS-CoV-2 for ACE2 binding in the wild type ACE2-RBD complex. The four pairs of i.e. Cys336-Cys361, Cys379-Cys432, Cys480-Cys488 and Cys391-Cys525 of disulfide bridges which stabilize the β-sheet structure were intact throughout MD. Comparative analysis of the MD simulated binding mode of the SARS-CoV-2 RBD to the ACE2 was found to be more or less very similar to that of experimental SARS-CoV RBD-ACE2 complex structure (2AJF, Li et al., 2005).

Table 2.

Comparative analysis of hydrogen-bond and electrostatic contacts of wild type SARS-Cov-2 RBD with ACE2 (Experimental vs. MD simulation).

6M0J (RBD-ACE2 experimental complex)			MD simulated complex
Interacting Pairs	Distance	Category	Interacting Pairs	Distance	Category
B:LYS417:NZ - A:ASP30:OD2	2.90	Hydrogen Bond; Electrostatic	A:LYS31:NZ - B:GLU484:OE1	4.33	Electrostatic
A:TYR41:OH - B:THR500:OG1	2.71	Hydrogen Bond	A:LYS31:HZ2 - B:GLN493:OE1	2.03	Hydrogen Bond
A:GLN42:NE2 - B:GLY446:O	3.24	Hydrogen Bond	A:GLN42:HE21 - B:GLN498:OE1	1.95	Hydrogen Bond
A:GLN42:NE2 - B:TYR449:OH	2.79	Hydrogen Bond	A:TYR83:HH - B:ASN487:OD1	1.78	Hydrogen Bond
A:TYR83:OH - B:ASN487:OD1	2.79	Hydrogen Bond	B:GLN493:HE22 - A:GLU35:OE1	2.01	Hydrogen Bond
A:LYS353:NZ - B:GLY496:O	3.08	Hydrogen Bond	B:GLN498:HE21 - A:GLN42:OE1	2.26	Hydrogen Bond
B:TYR449:OH - A:ASP38:OD2	2.70	Hydrogen Bond	B:THR500:HG1 - A:ASP355:OD2	1.91	Hydrogen Bond
B:ASN487:ND2 - A:GLN24:OE1	2.69	Hydrogen Bond	B:GLY502:HN - A:LYS353:O	3.02	Hydrogen Bond
B:GLN493:NE2:B - A:GLU35:OE1	3.13	Hydrogen Bond	B:TYR505:HH - A:GLU37:OE1	1.72	Hydrogen Bond
B:GLY502:N - A:LYS353:O	2.78	Hydrogen Bond	–	–	–

A: ACE2 B: RBD.

Previous studies have also shown that a positive-charged patch on the SARS-CoV-2 RBD contributed by Lys417 which is absent on the SARS-CoV RBD may be contributing differentially to the binding affinity of the SARS-CoV-2 (4.7 nM) and SARS-CoV (31 nM) to ACE2 receptor (Lan et al., 2020). Along this line, we also monitored changes in electrostatic surface for mutant complexes which clearly displayed mutations affect the accessible surface area and charge distribution at RBM (Figures S5 and S6). Most importantly, our study identified mutation in the conserved RBM affects the distribution of charged residue at the binding interface which in turn contributes to the differential interaction of RBD with ACE2 receptor.

The important hydrogen bonds and salt bridges presented at SARS CoV-2 RBD-ACE2 interfaces involve a network of hydrophilic interactions which are summarized in Table 2. Several Tyrosine residues like Y449, Y489 and Y505 and few other residues including two Asparagine residues viz., N487 and N501 in SARS-CoV-2 spike RBD region is seen to be important for favourable binding to receptor. Therefore, point mutation and subsequent MD simulation was focused on these residues. We first compared the experimental reference structure with the MD simulated systems and the retainment of strong hydrogen bonds observed for the aforementioned residues throughout the simulation analysis are in good agreement with the experimental data. However, such a stable interaction is not observed for the mutants with loss in strong Hydrogen bonds and binding affinity to the receptor (Table 3). In addition, we also measured the inter-atomic distance profile of the important interacting atom pairs (either hydrogen bond or electrostatic contact: for details see Table 3) of the RBD and ACE2 during the last 100 ns MD (Figure 7). Our results show that mutation caused a significant decrease in binding interactions for five residues of the interface corresponding to RBM segment. This was corroborated by the binding affinity of the residues with receptor as calculated by PRODIGY, where mutation decreased the binding affinity of the vital residues. Indeed, the RMSF variability of each of these mutated amino acid increases compared to the reference structure.

Table 3.

Important hydrogen-bonds interaction observed in the mutant systems (the cluster representative of each system was used to compute the hydrogen bonds using BIOVIA DSV).

Y449A			N487A
Interacting Pairs	Distance	Category	Interacting Pairs	Distance	Category
B:LYS417:HZ2 - A:ASP30:OD2	1.65	Hydrogen Bond; Electrostatic	B:LYS417:HZ3 - A:ASP30:OD2	1.64	Hydrogen Bond; Electrostatic
B:ARG403:NH1 - A:GLU37:OE1	4.10	Electrostatic	A:GLN24:HE21 - B:ALA475:O	2.05	Hydrogen Bond
B:LYS458:NZ - A:GLU23:OE1	4.82	Electrostatic	B:THR500:HG1 - A:ASP355:OD2	1.61	Hydrogen Bond
A:LYS31:HZ3 - B:GLN493:OE1	2.16	Hydrogen Bond	B:GLY502:HN - A:LYS353:O	2.10	Hydrogen Bond
A:HIS34:HD1 - B:TYR453:OH	1.77	Hydrogen Bond
A:TYR83:HH - B:ASN487:OD1	1.96	Hydrogen Bond
A:LYS353:HZ1 - B:GLY496:O	2.07	Hydrogen Bond
A:LYS353:HZ2 - B:GLN498:OE1	1.64	Hydrogen Bond
B:ASN487:HD22 - A:GLN24:OE1	2.57	Hydrogen Bond
B:GLN493:HE22 - A:GLU35:OE2	1.96	Hydrogen Bond
B:THR500:HG1 - A:ASP355:OD2	1.63	Hydrogen Bond
B:GLY502:HN - A:LYS353:O	1.93	Hydrogen Bond
B:TYR505:HH - A:GLU37:OE2	2.96	Hydrogen Bond
N489A			N501A
A:GLN24:HE22 - B:GLN474:OE1	2.19	Hydrogen Bond	B:LYS417:HZ2 - A:ASP30:OD2	1.61	Hydrogen Bond; Electrostatic
A:LYS31:HZ1 - B:GLN493:OE1	2.65	Hydrogen Bond	A:TYR83:HH - B:ASN487:OD1	1.82	Hydrogen Bond
B:SER477:HG1 - A:THR20:OG1	2.88	Hydrogen Bond	B:ASN487:HD21 - A:GLN24:OE1	2.14	Hydrogen Bond
B:ASN487:HD22 - A:TYR83:OH	2.20	Hydrogen Bond	B:GLN493:HE22 - A:GLU35:OE1	2.01	Hydrogen Bond
B:GLN493:HE22 - A:GLU35:OE1	2.29	Hydrogen Bond	B:THR500:HG1 - A:ASP355:OD2	1.81	Hydrogen Bond
B:THR500:HG1 - A:ASP355:OD2	2.46	Hydrogen Bond	B:GLY502:HN - A:LYS353:O	2.26	Hydrogen Bond
B:GLY502:HN - A:LYS353:O	2.46	Hydrogen Bond	–	–	–
Y505A
A:LYS31:HZ3 - B:GLU484:OE1	3.19	Hydrogen Bond; Electrostatic
B:LYS417:HZ2 - A:ASP30:OD2	1.80	Hydrogen Bond; Electrostatic
A:LYS31:HZ1 - B:GLN493:OE1	1.74	Hydrogen Bond
A:TYR83:HH - B:ASN487:OD1	1.76	Hydrogen Bond
B:ASN487:HD21 - A:GLN24:OE1	2.70	Hydrogen Bond
B:THR500:HG1 - A:ASP355:OD2	1.85	Hydrogen Bond
B:GLY502:HN - A:LYS353:O	1.82	Hydrogen Bond

A: ACE2 B: SARS CoV-2 RBD.

Figure 7.

The distance between the key interacting residue pairs i.e. ACE2 with RBD (SARS Cov-2) in wild type and mutant complexes during last 100 ns MD as a function of time. We only computed the distance between the interacting atom pairs of ACE2 and RBD forming electrostatic and hydrogen bonds.

Though Y449A system displayed a greater number of non-bonded contacts which may be attributed due to mutations, on the other hand it has lost few crucial H-bonds crucial for stability of the ACE2-RBD complex. Further a previous study by Wan et al. (2020), showed residue N501 in SARS-CoV-2 (corresponding to residue 487 in SARS-CoV) enhances viral binding to receptor ACE2 and plays a vital role in human to human transmission. We observed the same interaction in the simulated complex also but mutating the residue501 led to loss of favourable interaction which signifies that it has important implications for the pathogenesis, tropism and transmission of SARS-CoV-2.

The complete interaction between wild type and mutant SARS-CoV-2 RBD-ACE2 systems have been illustrated in Figure 8. Though significant difference between the binding affinity for the SARS-CoV2 and SARS-CoV RBD for ACE2 has been observed experimentally, it is unlikely that binding affinity alone explain the unusual transmissibility of SARS-CoV-2. It is noteworthy to mention here that other critical factors including furin cleavage site located that at the S1/S2 periphery of the spike protein also might be playing a crucial role in expediting the swift human-to-human transmission. Moreover, recent studies have also shown that viral entry also depends on TMPRSS2 protease activity and cathepsin B/L activity may be able to substitute for TMPRSS (Sungnak et al., 2020).

Figure 8.

Molecular interaction of wild type and mutant type complexes ACE2-RBD systems from MD trajectory. The top ranked snapshots extracted from the clustering analysis of the molecular dynamics trajectories were used for interaction analysis. (A) Interaction analysis of wild type and Y449 complex. (B) Intermolecular contacts of mutant N487, Y489A, N501 and Y505A complex systems. The hydrogen bond forming amino acids are colored in blue (ACE2) and green (RBD), respectively. The H-bonds are marked in green dotted lines along with their distance. Other non-bonded hydrophobic contact forming residues are labelled in semi-circles.

Saturation mutagenesis study of RBD-ACE2 interface

The entry of SARS-CoV to the invading cell depends on the binding of the spike protein to the specific host cell receptor and successive S protein priming by cellular proteases. The binding affinity of the spike protein and ACE2 was reported to be a chief factor for rate of replication and severity of the disease (Zhou et al., 2020 and Hoffmann et al., 2020). Therefore, to understand the binding affinity for SARS-CoV and SARS-CoV-2 RBD towards ACE2, we employed mCSM-PPI2 (Rodrigues et al., 2019), an integrated bioinformatics tool to assess the effects of non-synonymous mutations in protein-protein affinity (ΔΔGaffinity). The mCSM-PPI2 tool employs machine learning graph-based signature to predict the molecular effect of mutation in protein-protein binding affinity. Using both the experimental complexes of 2AJF (SARS-Cov-RBD-ACE2) and 6M0J (SARS-Cov-2-RBD-ACE2), saturation mutagenesis studies were performed to assess the effect of mutation on the interface residues of both the interacting partners (as summarized in Figures 9 and S7). All total 817 and 855 mutations were assessed for 2AJF and 6M0J complexes respectively. The average ΔΔGaffinity for all possible 19 mutations at the interface residues of 2AJF complex was computed to be −0.65 kcal mol−1 and for the 6M0J complex it was −0.55 kcal mol−1, which indicates mutation tends to decrease the ΔΔGaffinity at the interface for both the complexes. In addition, we also computed the binding affinity of two complexes using PRODIGY (Xue et al., 2016), a tool which uses contact-based prediction of binding affinity in protein-protein complexes (as summarized in Table S3). We observed slightly higher binding affinity for SARS-CoV-2-RBD towards ACE2 protein as compared SARS-CoV-RBD which can be correlated with rate of infection and disease severity by SARS-CoV-2 strains (Zhou et al., 2020 and Hoffmann et al., 2020). Consensus, our modeling study using both the tools predicted, mutation in the RBM of RBD affects the binding affinity and provides some crucial information on how inter-residue non-covalent interaction is crucial for protein-protein interaction.

Figure 9.

Saturation mutagenesis study depicting the effect of mutations at the interface of SARS-CoV-RBD-ACE2 (A) and SARS-CoV2-RBD-ACE2 (B). The interface residues are colored according to the average ΔΔGaffinity predictions values.

Altogether our MD simulation and saturation mutagenesis study of the interfacing residues has provided new insights to complex virus-receptor interactions which provide better avenues to understand SARS-CoV-2 disease biology. A series of structural studies on SARS-CoV have revealed that that receptor recognition by SARS-CoV is one of the utmost significant elements of its cross-species and human-to-human transmissions (Li et al., 2006 and Li, 2013). In this illustration, using experimental solved SARS CoV-2 RBD-ACE2 complex, multi-scale molecular dynamics and saturation mutagenesis were conducted to understand the host-virus interactions at atomic scale. The result from our study strongly supports how binding affinity between host cell receptor and virus spike protein is essential for infection and also provides some useful insights to the disease severity. The present study will certainly be useful as a framework to undertake further structure-based study on other hosts ACE2 receptors to provide novel insights into the host receptor-SARS CoV-2 RBD mediated interactions that may help to battle this novel COVID-19 outbreak.

Conclusion

Recent studies after the outbreak of COVID-19 have revealed that SARS-CoV-2 is highly homologous to human SARS-CoV (2005) and attributes to the human host cells through the binding of the spike protein to the angiotensin-converting enzyme II (ACE2) (Zhou et al., 2020). However, the molecular mechanisms of recognition and role of conserved amino acids at the RBM of RBD from SARS-CoV-2 binding to human ACE2 are under explored. In this study, we have comprehensively explored the wild type ACE2-RBD complex as well as five randomly selected mutants complexes (located in the conserved and variable sites of RBM of SARS-CoV-2) through all-atoms MD simulations.

Pair-wise sequence-structure alignment and structural superposition of ACE2-RBD complex from SARS-CoV-2 and SARS-CoV portrayed that ACE2 binds in the same fashion to the conserved receptor binding motif and share significant similarity in the side chain conformation, buried surface area of interacting residues, and other number of non-bonded contact networks with minor changes both in and outside the RBM (Figure 9).

Most of the systems reached equilibrium after 10 to 20 ns of equilibration and maintained their structural integrity like the experimental ACE2-RBD (SARS-CoV-2) complex till 150 ns. The RBM of spike protein of SARS-CoV-2 exhibits high degree of flexibility and most importantly induced mutations in these regions affects the mobility of RBM towards the ACE2 binding interface.

Principal component analysis suggested that RBD induces an inward motion towards the binding interface of host receptor ACE2, which is speculated to play a crucial role in the molecular recognition process besides other mechanisms behind the molecular interaction. The strength of binding depends on the perfect super-positioning side chains of both the interacting protein pairs at the binding interface mostly driven by the conserved binding residues located at RBM, where mutation tends to affect the binding affinity.

Site-directed mutagenesis and computational binding affinity studies showed that that mutation in the conserved RBM affects the structural-dynamics of the complex, affects the charge distribution and disturbs the inter-molecular non-bonded contacts thereby perturbs the strength of binding to host cell receptor ACE2. Recently a study has pointed out that SARS-CoV-2 may be more stable and can endure at higher temperature than SARS-CoV which points to the bat origin of SARS-CoV-2, as bats possess higher body-temperature as compared to humans (Chin et al., 2020; Qi et al., 2020). Hence, further studies are needed to understand the ACE2-RBD complex interactions at different temperature, which can provide crucial link of its infection ability at different temperature zones across the globe to formulate intervention therapeutic strategies for preventing and controlling SARS-CoV-2 in near future.

Overall, our results have shown encouraging findings in viral-host mediated receptor-interactions which can be very useful of development of novel inhibitors to combat against SARS-CoV-2. We strongly believe that our conclusions provide a suitable starting point for further development in understanding the highly contagious aetiology of SARS-CoV-2 infection in human.

Click here for additional data file.