In silico analysis of SARS-CoV-2 spike protein N501Y and N501T mutation effects on human ACE2 binding.

The SARS-CoV-2 is an RNA-based virus and the most vital step of its survival is the attachment to hACE2 through its spike protein. Although SARS-CoV-2 has the ability to maintain high accurate replication and it can be accepted as a low mutation risked virus, it already showed more than nine thousand mutations in spike protein, of which 44 mutations are located within a 3.2 Å interacting distance from the hACE2 receptor. Mutations on spike protein, N501Y and N501T raised serious concerns for higher transmissibility and resistance towards current vaccines. In the current study, the mutational outcomes of N501Y and N501T on the hACE2-SARS CoV-2 spike protein complexation were analyzed by employing all-atom classic molecular dynamics (MD) simulations. These simulations revealed that both N501Y and N501T mutations increased the binding strength of spike protein to the host hACE2, predicted by binding free energy analysis via MM/GBSA rescoring scheme. This study highlights the importance of energy-based analysis for identifying mutational outcomes and will shed light on handling long-term and effective treatment strategies including repurposing anti-viral drugs, anti-SARS-CoV-2 antibodies, vaccines, and antisense based-therapies.

Introduction

The SARS-CoV-2 (Severe Acute Respiratory Syndrome Coronavirus 2) is an RNA-based virus, which uses RNA as a genetic material instead of DNA [1,2]. The life cycle of these viruses begins with the attachment of viruses to the human host through receptor-ligand interactions. The viral surface receptor, receptor binding domain (RBD) of spike protein, comes in contact with the angiotensin-converting enzyme 2 (hACE2) on the surface of type 2 alveolar epithelium cells of the human upper respiratory tract (Fig. 1 ). After the binding of spike protein to the hACE2, SARS-CoV-2 virions begin to enter host cells either, via endosomes, and or membrane fusion. Following the viral entry, viral RNA is allowed to spread inside the host. Host machinery, ribosomes on the rough ER, translate viral genetic material directly as it is a positive-stranded RNA. Virion precursors are then assembled in the host, which is followed by the transport to the cell surface via small vesicles. SARS-CoV-2 virions released from the host are following the same life cycle, leading to high viral loads in the upper respiratory tract [[3], [4], [5], [6], [7], [8]].

Fig. 1

X-ray structure of hACE2-SARS CoV-2 spike protein RBD complex (PDB ID: 6M0J). The active site Zn2+ ion of hACE2 was shown with magenta sphere. Residues comprising the hot spot focused on the study were labeled.

Generally, RNA viruses are more open to new mutations as they lack exonuclease-proofreading activity of the virus-encoded RNA polymerases. The only exceptions are Nidoviruses, which.

Comprise the coronaviruses as well [9,10]. SARS-CoV-2, a member of the coronavirus family, has the ability to maintain high accurate replication. Thus, SARS-CoV-2 can be accepted as a low mutation risked virus [1]. However, COVID-19 caused by SARS-CoV-2 infection has affected millions of people around the world. Thereby, the COVID-19 pandemic provided a good environment for the creation of novel mutations. Recent studies have shown the presence of more than nine thousand mutations have raised in spike protein, of which 44 mutations located within a 3.2 Å interacting distance from the hACE2 receptor [9,11,12]. However, there is undoubtedly an increase in the number of single mutations that occurred during pandemic.

Previous studies on MERS-CoV and SARS-CoV and another member of the coronaviruses have demonstrated that single amino acid mutations make them resistant to neutralizing antibodies [[13], [14], [15]]. Furthermore, several mutations on the SARS-CoV-2 spike protein changed the SARS-CoV-2-hACE2 binding affinities and also transmissivity as well [14,16,17]. The bioinformatic analysis also showed that the spike protein-encoding region, RBD in the spike protein in particular, which generates proteins to be utilized in the host attachment, are more prone to new mutations compared to other genomic regions [9,18].

In this study, mutational outcomes of N501Y and N501T, from structural and energetics point of view, on the human hACE2-SARS-CoV-2 spike protein complex will be analyzed using molecular dynamics (MD) simulations. MD simulations will provide valuable data on the dynamics of each complex and binding free energy calculations will shed light on the change on the binding strength of spike protein upon mutations.

Materials and methods

Systems construction

X-ray structure of SARS CoV-2 spike protein RBD complexed with human hACE2 receptor (PDB ID: 6M0J, resolution: 2.45 Å) [19] was used as the initial model. Crystal water molecules were deleted and the most likely protonation state of each titratable residue at physiological pH (∼7.0) was predicted by the PROPKA server [20]. Asn501Tyr and Asn501Thr mutations were performed in YASARA Structure software [21], where the energy minimized conformations of Tyr and Thr residue side chains were obtained.

Molecular dynamics (MD) simulations and analysis

GROMACS 5.1.4 software [22,23] was utilized for MD simulations and GROMOS 53A6 force field was used [24]. The force field also includes the parameters for Zn2+. Each system, hACE2-WT spike protein, hACE2-N501Y mutant spike protein, and hACE2-N501T mutant spike protein, was placed into the center of the box with dimensions 10.0 x 10.0 × 15.0 nm. These dimensions were chosen to ensure that proteins stay in the box (unit cell) throughout the simulations and does not interact with its copy image dictated by periodic boundary conditions (PBC). The protein amino acid distances at the edges, calculated by the last atom to the solvent atoms at the edge of the unit cell were 1.7 nm, 1.3 nm, 1.3 nm for x, y, and z dimensions, respectively. These distances were much larger than chosen cut-off of 1.2 nm. The boxes were filled with extended simple point charge (SPC/E) model water molecules [25] and 58 Na+, 34 Cl- ions were added to neutralize systems and mimic the physiological environment. WT system comprised of ∼144,577 atoms, N501Y system of 144,590 atoms and N501T of 144,566 atoms. Each starting system was subsequently energy-minimized using the steepest descent method for 50,000 steps with energy minimization tolerance of 2.39 kcal mol−1 nm−1 and step size of 0.01 nm. Then, energy-minimized structures were taken for the production phase. All-atom MD production simulations were run for 100 nano-seconds (ns) for each system with the constant number of particles (N), pressure (P), and temperature (T), i.e. NPT ensemble, for each system. The SETTLE algorithm [26] was employed to constrain the bond length and bond angle of the solvent molecules. LINCS algorithm [27] was utilized to constraint bond lengths of the amino acid residues were constrained using. Long-range electrostatic interactions were treated with Particle-mesh Ewald (PME) method [28]. A constant pressure of 1 bar was applied with a coupling constant of 1.0 pico-seconds (ps) and water molecules/ions were coupled separately to a bath at 300 K with a coupling constant of 0.1 ps. Leap-frog algorithm [29] was used for integrating the equation of motion, which was integrated at 2 femto-seconds (fs) time steps. Tools available in GROMACS and VMD 1.9.1 [30] programs were utilized to analyze trajectories. The most representative structure was obtained by cluster command implemented in GROMACS 5.1.4. It adds a structure to a cluster when its distance to any element of the cluster is less than a cut-off value of 0.3 nm. The center of a cluster is the structure with the smallest average RMSD from all other structures of the cluster.

All results and discussion were based on the data from the entire 100 ns trajectory for each system. For some cases, the most representative structure was mentioned throughout the manuscript, which was obtained from GROMACS software based on equilibrated part (60–100 ns) of RMSD value graphs.

Binding free energy calculations and solvation energy estimation

Many free energy calculation methods have been developed to calculate the protein-protein binding free energies and determine the hot spots at interface surfaces thus far. Chen et al. studied the rescoring capability of Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) [31]. It was revealed that the use of the MM/GBSA scheme is a more efficient way to predict binding affinities for protein-protein interactions compared to alternative rescoring functions such as MM/PBSA. Therefore, we have decided to implement MM/GBSA scheme as a free energy calculation method. Binding free energy calculations were performed on the HAWKDOCK server (using default parameters), which utilizes MM/GBSA method [32]. As HAWKDOCK server calculates binding free energies both for the total complex and per residue, which helped us analyze contributions of each residue on the binding interface in detail. We provided snapshots (a total of 21 snapshots), obtained from every one ns of the equilibrated time frame of the simulations (80–100 ns) for hACE2-WT and mutant spike proteins to the HAWKDOCK server (MM/GBSA). For WT system, calculations were also repeated with snapshots for each 500 ps of the simulations, for the total of 41 snapshots and the results were similar (Supporting information 2).

We also used PISA server to evaluate the solvation-free energy gain (ΔiG) upon binding and interface surface [33]. These calculations were performed on cluster structures for each system, obtained with clustering method explained above.

Results and discussion

Even though the SARS-CoV-2, which is a member of the coronavirus family, has the ability to maintain replication with high accuracy, the current pandemic affecting all the continents, and regions on earth provided a huge chance for the natural selection of favorable mutations. In this study, we analyzed two mutations on the critical virus-binding hotspot (N501Y, and N501T) following computational methodologies and revealed structural details on how the mutations affect the binding energy of spike protein.

Validation of MD simulation protocol

We validated our MD simulation protocols by considering root mean square deviation (RMSD) values for the hACE2-spike protein complex during the simulation. RMSD values can also validate the MD simulation protocols: the smaller the RMSD, the higher accuracy for the reproduction of the crystal structure [34,35]. For this purpose, we obtained RMSD values throughout MD simulations (0–100 ns), which revealed that the protein complex reached equilibrium after ∼55 ns of the simulation (Fig. 2 ).

Fig. 2

Root mean square deviation (RMSD) values for WT, N501Y and N501T mutant ACE2-spike protein complexes throughout the MD simulation.

Moreover, to validate the methodology of the MD simulations, the most representative structure of the hACE2-WT spike protein complex from MD simulations and the X-ray structure (PDB ID: 6M0J) was superimposed. The RMSD value was 0.19 nm, which revealed good alignment on the overall complex and the hACE2 active site Zn2+ binding moiety (Fig. 3 ).

Fig. 3

Superposition of X-ray ACE2 (red) - SARS CoV-2 spike protein (blue) complex with ACE2 (orange) - SARS CoV-2 spike protein (cyan) the most representative structures from MD simulations. Active site Zn2+ ion coordination shells for both structures were provided in the inlet figures. The overall RMSD of the alignment was calculated as 0.19 nm. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

As a consequence, we concluded here that our MD simulations produced data in line with the experimental structures. After the validation of the MD simulation protocol, we approached the results of this study from two different perspectives, binding free energy and structural features leading to these energy values.

Analysis of N501Y, and N501T mutations on the receptor-binding domain (RBD)

First, we analyzed the binding free energy for hACE2 and WT, N501T, and N501Y mutant spike proteins by the MM/GBSA method to predict the experimental strength of hACE2-spike protein interactions [32]. From 80 to 100 ns time frame, snapshots were taken for every 1 ns, providing 21 structures for MM/GBSA analysis. The mean calculated binding energies from 21 snapshots were −61.94 ± 8.69 kcal/mol for WT, −94.16 ± 9.45 kcal/mol for N501Y mutant, and −82.91 ± 10.18 kcal/mol for N501T mutant, where more negative free energies indicated stronger binding. We applied t-test statistics to reveal if the means of binding free energies for WT, N501Y, and N501T were significantly different from each other. The time evolution of binding free energies and t-test statistics were provided in the Supporting Information Fig. 1. N501Y mutant provided much lower free energy than the WT, whereas N501T mutant also yielded much lower binding energies than the WT, but slightly higher free energy than N501Y mutant throughout last 20 ns. This result clearly showed that mutations produced remarkably lower binding free energies.

We have analyzed the individual amino acid contributions to the total binding free energy from both binding interface of hACE2 and RBD of spike protein. 12 hACE2 and RBD amino acids that yielded higher energies were listed in Table 1 and Table 2 , respectively. Upon mutations, amino acids’ contribution profile was significantly altered on the RBD residues, while the same profile stayed similar on the hACE2 residues. However, the sum of energy contributions from both 12 RBD and hACE2 residues increased significantly for both N501Y and N501T mutants compared to those of WT.

Table 1

Binding energy contributions of top 12 residues of hACE2 domain of spike proteins. (Energies provided in kcal/mol with their mean values ± standard deviation).

hACE2
WT		N501Y		N501T
Residue	Energy	Residue	Energy	Residue	Energy
Y41	−5.54 ± 0.76	Y41	−4.38 ± 0.59	Y41	−5.25 ± 0.93
Q24	−4.92 ± 0.92	Q24	−3.78 ± 0.70	Y83	−4.29 ± 1.01
Y83	−2.10 ± 0.56	Q42	−3.51 ± 0.95	D38	−3.63 ± 2.01
T27	−2.00 ± 0.80	L45	−3.35 ± 0.52	K31	−3.2 ± 1.34
Q42	−1.68 ± 0.479	Y83	−2.46 ± 0.18	Q24	−3.15 ± 0.67
H34	−1.49 ± 0.66	L79	−2.45 ± 0.30	L79	−2.77 ± 0.20
L45	−1.39 ± 0.21	K31	−2.39 ± 1.55	L45	−2.27 ± 0.42
F28	−1.32 ± 0.34	F28	−2.21 ± 0.23	Q42	−2.18 ± 0.33
D38	−1.15 ± 1.02	T27	−2.20 ± 0.27	M82	−2.11 ± 0.35
L79	−1.06 ± 0.65	M82	−1.95 ± 0.28	F28	−1.98 ± 0.15
M82	−0.84 ± 0.26	D38	−1.14 ± 1.99	T27	−1.91 ± 0.28
K31	−0.75 ± 0.62	N61	−1.09 ± 0.96	H34	−1.01 ± 0.40
Sum	−24.2		−30.9		−33.8

Table 2

Binding energy contributions of top 12 residues of RBD domain of spike proteins. (Energies provided in kcal/mol with their mean values ± standard deviation).

Spike protein (RBD)
WT		N501Y		N501T
Residue	Energy	Residue	Energy	Residue	Energy
Y505	−5.97 ± 0.32	Q493	−9.07 ± 1.71	F486	−5.46 ± 0.75
F486	−3.30 ± 0.66	Y501	−5.82 ± 1.55	Y489	−4.8 ± 0.32
Y449	−3.00 ± 1.70	Y449	−5.23 ± 0.75	Q498	−4.71 ± 0.54
G476	−2.98 ± 1.72	F486	−5.21 ± 0.50	Y449	−3.8 ± 0.83
N501	−2.79 ± 0.64	Y489	−4.99 ± 0.44	Q493	−3.51 ± 0.38
Y489	−2.36 ± 0.47	S494	−4.09 ± 1.95	V445	−3.18 ± 0.59
N487	−2.09 ± 1.07	V445	−3.59 ± 1.50	Y505	−2.82 ± 0.97
F456	−1.99 ± 0.30	L455	−2.99 ± 0.37	F456	−2.76 ± 0.25
A475	−1.73 ± 0.94	F456	−2.46 ± 0.24	L455	−2.59 ± 0.59
S477	−1.72 ± 1.13	Y505	−2.12 ± 0.52	S494	−1.97 ± 0.29
V445	−1.58 ± 1.05	Y495	−1.99 ± 0.40	N487	−1.95 ± 0.91
Q498	−1.48 ± 0.39	G496	−1.66 ± 0.64	T501	−1.73 ± 1.52
Sum	−31.0		−49.2		−39.3

When Asn 501 mutated to Tyr, the energy contribution of this residue increased by 3.03 kcal/mol (lower binding free energy), indicating a stronger interaction with hACE2, whereas Thr mutation significantly decreased the contribution by 1.06 kcal/mol (higher binding free energy). The reason for these alterations was the new environments around these residues. Tyr being an aromatic amino acid, brings the aromatic amino acid residues Tyr 41 (hACE2), Tyr 449 and Tyr 505 (RBM) together bound into highly hydrophobic binding pocket (Fig. 4 ). hACE2 residues Leu 45 and Leu 351 lies in this hydrophobic pocket as well.

Fig. 4

2-D amino acid interaction maps around N/Y/T 501 positions of WT, N501Y, and N501T mutants. A chain represents hACE2 and B chain represents spike protein (ex. GLN B:506 represents spike protein residue Gln 506).

It was clear that Tyr 501 was very well accommodated in highly hydrophobic pocket, creating face to edge (π-π T-shaped) interaction with Tyr 41 (Fig. 5 ). This aromatic interaction was not observed in WT and in N501T mutant. Moreover, there were π-π stacking (face to face) interaction between Tyr 41 (hACE2) and Tyr 505 (Spike), and π-alkyl interaction between Leu 351 (hACE2) and Tyr 501 (Spike) for N501Y mutant (Fig. 6 ).

Fig. 5

Aromatic environments of residue 501 interaction site for WT, N501Y, and 501T mutants.

Fig. 6

Zoomed in amino acid interactions around N/Y/T 501 positions of WT, N501Y, and N501T mutants. Residue 501 was labeled in purple. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Tyr 501 created a strong H-bond with Asp 355 (hACE2). On contrary, Thr 501 in N501T mutant did not interact directly with any of hACE2 residues.

Furthermore, we created binding free energy contribution map through the hACE2-RBD interface, which also reflects the altered contributions shown in Table 1. Both mutations changed contributions mainly on four regions.

Fig. 7 highlights that N501Y mutation resulted in stronger hACE2 interaction in regions 1 (local area of mutation), 2, and 3. On the other hand, N501T mutation did not alter the interaction significantly in region 1, but it led to stronger interaction in regions 3 and 4.

Fig. 7

Binding free contribution maps for hACE2-Spike protein complexes (WT, N501Y, and N501T).

Although, N501T did not change the energy contribution in the local mutation area (region 1), it significantly increased energy contributions in other regions, which explained the overall increased binding free energy for this mutant.

We also calculated the solvation-free energy gain (ΔiG) upon binding by utilizing the PISA server [36]. The most representative (cluster) structures were used for these calculations. hACE2-N501Y spike protein complex provided 1.4 kcal/mol decrease in ΔiG, while N501T mutant provided 0.5 kcal/mol lower energy compared to hACE2-WT complex (Table 3 ).

Table 3

Solvation free energy gains and interface areas in Å2 calculated from cluster structures by PISA server. More negative ΔiG indicates stronger binding of proteins.

Complex	ΔiG (kcal/mol)	Interface surface (Å2)
WT	−9.1	941.2
N501Y	−10.5	1114.8
N501T	−9.6	956.4

This was expected because tyrosine comprised of highly hydrophobic side chain gained more energy upon complexation compared to threonine and asparagine, comprised of the polar side chain and Tyr 501 was accommodated better than Thr 501 and Asn 501 in region 1. Increased total binding interface surface area for the mutants also supports the observed higher binding free energies.

Conclusion

In this study, the binding energy-based consequences of two spike protein RBD mutations, N501Y and N501T, were analyzed by utilizing all-atom classical MD simulations. Binding free energy calculations revealed that N501Y mutant spike protein bound to hACE2 with higher affinity, 32.22 kcal/mol lower binding free energy, than the WT protein. Golubchick et al. [37] studied on the N501Y reported increased viral load for the lineage including N501Y mutation in-vitro. Additionally, in-silico studies via FEP simulations by Luan et al. provided increased binding affinity towards hACE2 upon N501Y mutation [38], which is parallel to our results. The calculated binding free energy for N501T mutant spike protein was 20.97 kcal/mol lower than the WT spike protein, but 11.25 kcal/mol higher than the N501Y mutant. Shang et al. also shared an in-vitro data [17] proposing reduced spike binding affinity towards hACE2 upon N501T mutations but, other studies conducted by Wang et al. [39] showed opposite results, leading increase in the RBD and ACE2 binding affinity. Our results for increased binding strength of N501T mutant spike protein was in line with data presented in the study by Wang et al. This study combines the structural and energy-based approaches together for the detailed analysis of N501Y and N501T mutations on SARS-CoV-2 spike protein.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.