Introduction

COVID-19 pandemic which is the result of infection by SARS-Coronavirus-2 (CoV-2), had claimed over 6 million lives as of June 2023, and continues impacting worldwide health [1]. CoV-2 expresses S (Spike) protein [2, 3], which is responsible for binding to the host cell receptor followed by fusion of the viral and cellular membranes, respectively [4]. To engage a host cell receptor, the receptor-binding domain (RBD) of the S protein undergoes hinge-like conformational movements [5], thus mediating interaction with the receptor angiotensin-converting enzyme 2 (ACE2) [6–9] (the RBD and ACE2 binding structure is shown in Fig 1(A). The binding of RBD with ACE2 is believed to be the critical initial event in the infection cascade. The interface residues on RBD thus play key roles in the spike protein’s binding with ACE2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

The binding structure of ACE2 (in blue) with RBD (in green) (a), with the mutations (shown in red sticks) on RBD labeled; and the comparison of the initial and final structures of ACE2 bound with RBD-trimutant variant (b). In b), the initial structures of ACE2 and RBD are shown in wheat, and the final structure of ACE2 is shown in raspberry and magenta for RBD. The mutations of E484K, K417N, and N501Y are shown in blue sticks in both the initial and final structures.

https://doi.org/10.1371/journal.pone.0289432.g001

Coronaviruses are enveloped viruses that possess a positive-sense single-stranded RNA genome 26–32 kb in length [10], and have a moderate mutation rate that allows them to evolve in a way that perhaps by increasing their transmissibility [11–14] and/or their resistance to protective immunity induced by previous infection or vaccines [15–22]. Up to now, different variants of concern (VOCs) have shown up, including Alpha, Beta, Gamma, Delta, and Omicron [23–30] (Some mutations involved are shown in S1 Table in S1 File). The Alpha variant was found in UK in September of 2020, while Beta variant was found in South Africa in May of 2020. Gamma variant was found in Brazil [23–25, 27]. The Omicron variant which contains 11 mutants on RBD, is still the dominant cause of the resurgence of infection in many countries currently [29]. Most of the mutations (around 58% of mutations) on RBD are located at the receptor binding motif which is the binding interface with ACE2 [31].

N501Y on RBD, which appeared in Alpha, Beta, Gamma and Omicron variants, increased the binding affinity of RBD with ACE2 [32, 33] and enhanced the viral resistance to neutralizing antibodies [34]. S477N which appeared in Omicron, and E484K which appeared in both Beta and Gamma variants, were reported to enhance the binding affinity of S protein with ACE2 [35, 36]. K417N found in both Beta and Omicron variants was reported to decrease the binding affinity of RBD and ACE2 [37], and was thought to cause immune escape [38]. The triple mutants of N501Y-E484K-K417N appeared in Beta variant, has been shown to increase the infectivity of virus considering the effect of each single mutants [37]. N439K, which was found in multiple regions in the world and a quite common mutation, can enhance the binding of RBD with human ACE2 receptor while evade antibody mediated immunity [39]. P479S, which also happened with high frequency at most countries based on SARS-CoV-2 genome samples deposited at GISAID in early 2020 [40], was predicted to make the virus less infective based on algebraic topology-based machine learning modelling [31]. Besides that, based on the report of most frequently occurring mutations in RBD in January 2021 (https://www.gisaid.org/hcov19-mutation-dashboard), T478I which appeared in Alpha variant, however did not influence the binding of RBD and ACE2 [41]. Those mutants all appeared in RBD binding interface.

Besides those experimental work, in order to find out the mechanism of CoV-2 variants binding with ACE2 and to design effective therapeutics to deal with CoV-2 variants, molecular dynamics simulation methods have been applied to study RBD binding with ACE2 [35, 36, 42–49]. Different simulation methods were used to predict the binding free energy change (ΔΔG) [50, 51], including free energy perturbation (FEP) method [43, 45], molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method [44], molecular mechanics-generalized Born surface area (MM/GBSA) method [47] and so on. Li et al. [51] and other researchers found that FEP method can predict ΔΔG reasonably [43, 50, 52] and is better than other methods [51]. Thus FEP method was applied to predict the point mutation on RBD to the binding of RBD variants with ACE2 in this work. Besides that, dynamic correlation was calculated to predict the correlation matrix and dynamic network in RBD and ACE2 in order to understand how the S spike mutations on RBD binding interface can affect the binding of RBD with ACE2. Since point mutations happened quite often on CoV-2 RBM [31], a comparative study on multiple single mutations on RBD binding with ACE2 is necessary. In total 7 RBD single mutants which happened naturally and can potentially affect the binding and dynamics of RBD with ACE2 were studied including N501Y, E484K, K417N, N439K, P479S, T478I, S477N, and one triple mutant N501Y-E484K-K417N. Their effects on the binding of RBD with ACE2 were analyzed. Besides all-atom NAMD simulations, FEP method was applied to predict ΔΔG of RBD with ACE2 because of the point mutation. The simulation predictions helped to clarify the contradicting available experimental results. Calculating the dihedral angle covariance matrices and doing dynamic network analysis, an intrinsic network community was predicted. Our results revealed the molecular mechanism that underlies the dynamics and energy change of some RBD variants binding with ACE2. Such kind of information will be valuable for the development of future vaccines and neutralizing antibodies against mutant forms of the SARS-CoV-2 virus.

Materials and methods

1). Based on the survey up to now, 7 naturally occurring single RBD mutants were studied: E484K, N439K, P479S, T478I, N501Y, K417N, S477N, besides one triple mutants of RBD which has E484K, N501Y and K417N. The selection of those mutations is based on different SARS-CoV-2 variants shown as listed in S1 Table in S1 File. The mutation sites on ACE2 and RBD binding interface are shown in Fig 1(A). The crystal structure of the RBD domain of the Spike protein is available in complex with ACE2 at 2.45 Å resolution in the PDB with ID 6M0J [53], which was used as the starting structure to build the RBD mutants bound with ACE2. VMD program [54] was applied to mutate residues on RBD to generate RBD mutants. The corresponding RBD and ACE2 structures in [53] were used in reference simulations of RBD and ACE2 in free forms respectively. The bound and free forms of simulations conducted in this project are listed in S2 Table.

Different simulation systems were set up using VMD program and CHARMM36m forcefield [55]. The system was solvated with enough amount of TIP3P water molecules with the minimum distance between protein and water box edge being at least 12 Å, and was neutralized with 0.15 M NaCl solvent. After a brief energy minimization using the conjugate gradient and line search algorithm, 4 ps of dynamics was run at 50 K, and then the system was brought up to 310 K over an equilibration period of 1 ns using NAMD program [56]. All-atom molecular dynamics simulation was performed on each system using NAMD program and CHARMM36 forcefield [57]. The temperature was 310 K and a desired pressure of 1 atm using standard thermo- and barostats. The sampling run was 100 ns with a time step of 1 fs. The trajectory output frequency was 1 ps. To do comparisons on the binding of RBD mutant with ACE2, the RBD bound with ACE2 wildtype systems were also set up and simulations were performed for 100 ns for RBD mutant & ACE2. The simulation systems, set up, the number of atoms and box size information are shown in S2 Table.

To analyze the trajectories, the Root Mean Square Deviation (RMSD) and Fluctuations (RMSF) of the proteins were calculated using the VMD program and an in-house analysis script based on the coordinates of the backbone Ca atoms after aligning the trajectories respectively, to the original crystal structure of the RBD, ACE2, and to the initial complex structure of the RBD and ACE2 from Lan et al. [53]. The buried surface area (BSA) for the complex was calculated in two steps using the VMD program and a script using the Richards and Lee method with the water probe size of 1.4 Å [58]. First, the total solvent accessible surface area of the complex (ASA_complex) was calculated based on the complex’s trajectory. Second, the accessible surface area of each protein in the complex (ASA_RBD, ASA_ACE2) was calculated for each protein individually. Then, the buried surface area (BSA) is calculated using Eq (1): (1)

The number of hydrogen bonds between the RBD and ACE2 were calculated using the VMD program with the heavy atom distance cutoff of 3.0 Å and the angle cutoff of 20 degrees deviation from H-bond linearity. The time a particular H-bond is formed over the course of the simulation is monitored and is expressed as % occupancy. In order to find out the residues on the binding interface, the closest distance between every residue atom (including hydrogen) between the RBD and ACE2 was calculated and averaged over the trajectory run. The average distances between each residue on RBD and on ACE2 are shaded by proximity on a red to white color-scale and were used to build the distance maps.

2). Dihedral Angle Correlation, network community calculation. To quantitatively describe the motional correlations of the backbone in RBD, we calculated the ϕ and ψ cross-correlation matrix for all the residues except the first or last one on the sequence. ϕ describes the rotation of the N−Cα bond and involves the CO−N−Cα−CO bonds. ψ describes the rotation of Cα and the CO bond and involves the N−Cα−CO−N bonds. Both ϕ and ψ were calculated for residues at different times using the wordom program [59, 60]. Because both ϕ and ψ are angular variables, to avoid the periodicity problem, the circular correlation coefficient, which is a T-linear dependence, was calculated in this project following the method of Fisher [61, 62] and Mardia and Jupp [63], and in the same way as in our previous work [64, 65]. The ϕ−φ cross-correlation coefficients (matrices) were calculated of all combinations of dihedral angles. In order to simplify the data representation, the ϕ−ϕ, ϕ−φ, and φ−φ correlation matrices were combined and the element with the largest magnitude at residue pair position was always chosen to build the final matrix.

In the dynamic network community calculation, the network model was built by the NetworkView plug-in of VMD [66] and the program Carma [67]. In the network, the nodes could be considered as a single atom or a cluster of atoms. In this project, each Cα atom in one amino acid was treated as one node. The edge between two nodes was defined with the cutoff distance of 4.5 Å for at least 75% of molecular dynamics simulated trajectory.

3). In order to find out the contribution of mutant residues on RBD to its binding with ACE2, FEP method has been applied in the same way as in [52]. The change of free energy of each mutation to the binding was predicted, including E484K, N439K, N501Y, P479S, S477N, T478I, K417N on RBD binding with ACE2.

In total, two kinds of systems were set up. One is the RBD in wildtype form in solvent (unbound state). The other is the RBD mutant bound with ACE2 (bound state). For both kinds of systems, mutation on one residue (E484K, N439K, K417N, T478I, P479S, S477N, N501Y) on RBD (mutant systems) was conducted using VMD program [54]. After solvating the unbound and bound systems in both the wildtype and mutant forms with enough amount of water and neutralizing them with 0.15 M NaCl, NAMD all-atom molecular dynamics simulations were performed to equilibrate both kinds of systems using an NVT ensemble after a brief energy minimization. The dual-topology methodology [68–70], and the soft-core potentials [71, 72] were applied to overcome endpoint singularities [73] in the FEP simulation, following the same method as in our previous work [52]. The free energy change in mutation for residue i in thermodynamic cycle can be calculated using Eq (2): (2)

Here, ΔΔG_i is the relative change of alchemical free energy of residue i during the binding of RBD with ACE2 as the thermodynamic cycle shown in S1 Fig. In each FEP simulation, 3 ns simulations were performed in each window for both forward and backward directions and for both the bound and free states after 0.5 ns simulations were conducted to equilibrate the system. In this project, the interaction calculation includes 10 windows, ranging from a thermodynamic coupling parameter λ values of zero to 1 (with a gap of 0.1) for a total of 20 simulations per mutation. The time step was 1.0 fs.

Results

1). RMSD, RMSF, ΔRMSF, DSSP results for ACE2 and RBD variants

RMSD result.

Based on all-atom NAMD simulations on different RBD variants bound with ACE2, the RMSDs of the complex, RBD, ACE2 were calculated after aligning the trajectories to the RBD wildtype bound with ACE2, RBD wildtype, and ACE2 individually. RBD mutants, ACE2, and the complexes are very stable during 100 ns simulations and are very consistent with those from ACE2 bound with RBD wildtype (with the RMSD vs simulation time result shown in S2 Fig in S1 File), as the comparison of the initial and final structures for RBD-trimutant&ACE2 shown in Fig 1(B). Calculating the average RMSD of the complex, the RBD, and the ACE2, the results are shown in Table 1. Only ACE2&RBD-N501Y complex has a RMSD (2.09±0.2 Å) lower than that of the ACE2&RBD wildtype (2.21±0.4 Å). Although RBD-N501Y has the highest RMSD (1.54±0.1 Å), ACE2 bound with this RBD variant has the lowest RMSD (1.84±0.2 Å). That suggests that RBD-N501Y can fit the binding pocket of ACE2 with the minimum structure change on ACE2 and the maximum structure adaptation of RBD-N501Y. On the other hand, the complexes of ACE2 with RBD tri-mutant, RBD-T478I, and RBD-N439K showed the largest structural deviation from the wildtype(3.05±0.4 Å, 3.06±0.4 Å, 3.34±0.5 Å). When bound with RBD-N439K, ACE2 has the largest structural deviation (2.60±0.3 Å) comparing to the ACE2 wildtype (1.88±0.2 Å). Overall, the RMSD result suggests that point mutation did not affect the binding structure of ACE2 with RBD significantly.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. The average and standard deviations of RMSDs of ACE2&RBD complex, RBD wildtype and variants, and ACE2, the number of hydrogen bonds (Hbonds) formed, and the buried surface area (BSA) between ACE2 and RBD.

The standard deviations were calculated based on the second-half of NAMD simulations.

https://doi.org/10.1371/journal.pone.0289432.t001

RMSF and ΔRMSF result.

Calculating the Root Mean Squared Fluctuation (RMSF) of RBD and ACE2 from different simulations, the results are shown in Fig 2(A) and 2(B). Residues from 455 to 505 on RBD bound with ACE2 at the region of RES19 to RES81, which are highlighted in light green for RBD on Fig 2(A) and for ACE2 on Fig 2(B). Within the binding region, all the mutants have a smaller RMSF than the wildtype, especially S477N, N439K, T478I and the tri-mutant. However, at other regions such as RES360 to RES390, E484K showed a larger structure fluctuation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Comparison of RMSF from RBD variants and wildtype (a) and ACE2 (b). The variant results are shown in black solid line while the wildtype shown in red dashed line. The binding regions on RBD and ACE2 are shaded in light green, and the secondary structures of RBD and ACE2 are shown on the top of RMSF figures.

https://doi.org/10.1371/journal.pone.0289432.g002

Mapping the change of RMSF (ΔRMSF = RMSF_mutant-RMSF_wt) on the secondary structure of RBD and ACE2, the result for N501Y is shown in Fig 3(A) and the tri-mutant bound with ACE2 in Fig 3(B). Comparing to the wildtype, RBD N501Y is slightly more rigid at RES501 site and RES484 site, although other binding region has a similar rigidity to the free form. However, ACE2 is more rigid on the binding interface except the loop region which becomes more flexible than the ACE2 when bound with RBD wildtype. Fig 3(B) shows that part of the binding interface on RBD becomes more rigid, part unchanged while the left region becomes more flexible comparing to the RBD wildtype, while ACE2 also has part of the interface more rigid, part unchanged, and other part more flexible comparing to the ACE2 bound with RBD wildtype. Those are consistent with RMSF result in Fig 2. Similarly, the ΔRMSF of other RBD mutants and ACE2 were calculated and mapped to the secondary structure of RBD and ACE2 (shown in S3–S5 Figs). In summary, the binding interface of ACE2 molecule becomes more rigid when bound with RBD variants except with S477N and N439K which have more flexible region than rigid region on the binding interface. Comparing to the RBD wildtype, RBD variants only become slightly more rigid except RBD N439K variant which becomes mostly more flexible. Faraway regions of ACE2 bound with RBD variants usually become more flexible (such as RBD-S477N, RBD-N439K, and RBD-K417N, RBD-T478I), but not for the ACE2 bound with E484K, which overall become more rigid. Those can agree with results in RMSF result shown in Fig 2, and the result suggests that point mutations on RBD can change the overall structure flexibility of RBD and ACE2 comparing to the wildtype. Such kind of change of flexibility should contribute to the adaptation of RBD variants to fit the binding pocket of ACE2 which also can adjust its flexibility accordingly.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

The ΔRMSF of RBD-N501Y and ACE2 (a) from wildtype to variant forms of RBD mapped to the initial bound structure of ACE2&RBD-N501Y; and the ΔRMSF of RBD-trimutant and ACE2 (b) mapped to the initial bound structure of the complex. Colors range from blue to white to red (with white representing no change in RMSF, blue representing RMSF decrease, and red representing RMSF increase). The ΔRMSF is in the range from -0.5 (totally blue) to 0.5 (totally red).

https://doi.org/10.1371/journal.pone.0289432.g003

DSSP result.

In order to track the structure change of ACE2 and RBD contributed by the point mutation, the evolution of the secondary structure was predicted based on DSSP [74] calculation over the simulation time using the wordom program. The DSSP program assigns the protein secondary structure to one of the nine states defined, including three helix types (3₁₀, G; α, H; π, I), two extended sheet types (antiparallel and parallel, E and X), two types of turns (including the helix turn (T) and β bridges (B)), bend (S), and others (L). The secondary structure was calculated for each residue at each time frame during the simulation time. Based on that, the average proportion of helix (including G, H and I), sheet (including E), turn (including B and T) and other (including X, S and L) structures on RBD and ACE2 over the simulation time were calculated, and the results in comparison with available literature data [41] are shown in Fig 4 for RBD and in S6 Fig in S1 File for ACE2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparison of RBD DSSP with available literature data.

https://doi.org/10.1371/journal.pone.0289432.g004

Overall, the proportions of RBD secondary structure elements are very consistent with available experimental data [41], with around 13%-15% of helix, 30%-32% of sheet (including both parallel and anti-parallel), and 13%-14% of turns, and 40%-43% of others. That implies the RBD structures during NAMD simulations are reasonable. In Fig 4, comparing the RBD mutants with RBD wildtype from the experimental data (the last 6 columns which only have RBD in free state available), it shows that the average percentage of other structure is higher in RBD mutants than the RBD wildtype. That is consistent with the simulation prediction for RBD wildtype and variants (the first 9 columns which are all in bound form). That suggests that point mutation can increase undefined structure (other structure) in RBD. Although DSSP of RBD in wildtype free state (wt, free) has similar helix and turn ratios to the experimental data of RBD in wildtype free state (wt, free, ex), its sheet ratio is higher than the experimental data by around 6% while its other structure ratio is lower than the experimental data by around 6%. Such kind of discrepancy should come from the different definitions of sheet from Anderson [73] and the experimental method. Because of that, the average percentage of sheet structure from simulations are overall higher than experimental data while the average percentage of other structures are overall lower than experimental data.

Comparing results from simulation only or from experiments only, overall point mutation only changes the secondary structure distribution of RBD slightly in both bound and free states. It almost does not change the secondary structure distribution of ACE2 as shown in S6 Fig. So, overall, the mutants and wildtype have similar secondary structure distribution and point mutation should only affect the secondary structure locally.

2). Hydrogen bonds formed on the binding interface

Forming intermolecular hydrogen bonds is one of the major driving forces for the binding of RBD with ACE2. Calculating the number of hydrogen bonds formed between ACE2 and RBD mutants over the simulation time (shown in S7 Fig), overall RBD mutants formed similar amount of hydrogen bonds with ACE2 to the RBD wildtype during the simulation time. Calculating the average and standard deviation of the number of hydrogen bonds formed, the results are shown in Table 1. S477N and P479S formed more hydrogen bonds with ACE2 than the wildtype, while T478I formed the least, and the tri-mutant formed the second to the least amount of hydrogen bonds with ACE2. Tracking the residue pairs on the binding interface, the results are shown in S3 Table in S1 File. There are 8 pairs of residues forming hydrogen bonds on the RBD wildtype and ACE2 interface, including Gly502-Lys353, Tyr505-Glu37, Asn487-Tyr83, Thr500-Tyr41, Lys417-Asp30, Gln493-Glu35, Tyr449-Asp38, Gln498-Gln42. With the triple-mutation, three pairs (Lys417-Asp30, Tyr449-Asp38, Gln498-Gln42) were broken and only the other five still formed hydrogen bonds in the tri-mutant&ACE2 interface. However, a new hydrogen bond formed between Thr500-Asp355, and the interface structure generated from ligplot program [75] in comparison with RBD wildtype is shown in Fig 5. The details of hydrogen bonds formed between ACE2 and RBD wildtype and variants are shown in S3 Table. Overall, similar binding interface formed in RBD wt&ACE2 and RBD variant&ACE2, although the exact residues forming hydrogen bonds changed and the number of hydrogen bonds formed also changed by certain amount. Such kind of small modification on the binding interface contributes to the different binding free energy of ACE2 and RBD which mainly comes from hydrogen bonds formation and electrostatic and hydrophobic interactions [47], and affects the infectivity of CoV-2 variants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Hydrogen bonds formed on the ACE2&RBD interface for the wildtype (a) and tri-mutant (b) generated from ligplot program. Residues formed hydrogen bonds are connected using dashed green lines. Residues on RBD are labelled in light green with bonds shown in magenta, while residues on ACE2 are labelled in blue with bonds shown in orange.

https://doi.org/10.1371/journal.pone.0289432.g005

3). BSA and distance map result

BSA comparison.

Calculating the buried surface area (BSA) between ACE2 and RBD variants over the simulation time, the results in comparison with the RBD wildtype are shown in Fig 6. RBD wildtype and variants can bind with ACE2 consistently during the simulation time. Calculating the average and standard deviations of BSA based on the second half of the simulation time, the results are shown in Table 1. Overall, the RBD variants and ACE2 form a similar BSA as the RBD wildtype and ACE2. The RBD-T478I formed a slightly smaller BSA than other mutants, while RBD-E484K, RBD-N439K and RBD-S477N form a larger BSA with ACE2 than the RBD wildtype; the RBD-trimutant forms the smallest BSA among all the simulations especially during 20 to 100 ns period in the simulations as shown in Fig 6. Based on the findings from Chan et al. [45], a larger BSA suggests a stronger binding of RBD with ACE2, thus it is expected that E484K, N439K and S477N variants should form a stronger binding with ACE2 than other variants, while the tri-mutant should form a less stable binding with ACE2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Comparison of the BSA between RBD variants and ACE2 with that of the ACE2&RBD wildtype started from the crystal structure.

https://doi.org/10.1371/journal.pone.0289432.g006

Distance map result comparison.

Calculating the distance map between residues on ACE2 and RBD, the results of RBD-E484K bound with ACE2 in comparison with RBD wildtype bound with ACE2 are shown in Fig 7 as an example. The other RBD variants bound with ACE2 distance maps are shown in S8 and S9 Figs in S1 File. Those results show that different RBD variants bound with ACE2 at the interface consistent with that of RBD wildtype and ACE2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7.

Distance map comparison of a). RBD-E484K mutant bound with ACE2 and b). RBD wildtype bound with ACE2. The binding region (RES455 to 505) on RBD are highlighted in light green. The color bars are shown on the right side of the distance maps, with a residue pair distance being no larger than 4 Å shown in red, while a residue pair distance larger or equal to 14 Å shown in white.

https://doi.org/10.1371/journal.pone.0289432.g007

Even with the small modification on the residue pairs forming interactions on the binding interface, a similar BSA and distance maps were formed for RBD variants with ACE2 comparing to RBD wildtype. Thus it is not surprising to see that overall the binding structures of RBD mutants&ACE2 complexes are very consistent with the wildtype which is the crystal structure of RBD&ACE2 from Lan et al. [53] (the comparison of the initial and final structures for RBD-trimutant&ACE2 shown in Fig 1(B) as an example). The structures of RBD mutants bound with ACE2 showed limited deviation from the original structure, which can agree with RMSD, RMSF, BSA and hydrogen bonds results.

In summary, RBD-mutants bind with ACE2 consistently and at a consistent interface. Because the point mutation did not change the binding interface of RBD and ACE2 significantly, it is reasonable to start from the crystal structure of ACE2&RBD to do FEP simulation and predict the point mutation effect on ΔΔG.

4). FEP result

Analyzing FEP simulation results, ΔΔG was predicted for RBD mutants. The results in comparison with Upadhyay et al. [41] are shown in Fig 8. The comparison of simulation prediction with four more available literatures [36, 37] are shown in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. ΔΔG for RBD variants binding with ACE2 predicted from FEP (shown in blue bars) in comparison with literature data from Upadhyay et al. (shown in orange bars).

Reference ΔΔG was calculated based on available data in the table as ΔΔG = ΔG_mutant-ΔG_wt.

https://doi.org/10.1371/journal.pone.0289432.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. ΔΔG predicted from FEP in comparison with available literature data.

https://doi.org/10.1371/journal.pone.0289432.t002

Based on Fig 8, N501Y, N439K, K417N contribute to the binding of RBD with ACE2, T478I and S477N have negligible effect, while neither E484K nor P479S contributes to the binding. Upadhyay et al. [41] did ITC experiments on the binding of RBD mutants with ACE2. Calculating ΔΔG using equation: ΔΔG = ΔG_mutant-ΔG_wt based on available data, ΔΔG data from Upadhyay are shown in Fig 8 and Table 2 above. As can be seen from Upadhyay et al., N501Y can increase the binding of RBD with ACE2, T478I has negligible effect while E484K can make the binding weaker. Our prediction of those three mutations can agree with Upadhyay et al., and the effect of N501Y also agrees with Luan et al. [43] who conducted all-atom molecular dynamics simulations and FEP on RBD N501Y mutant binding with ACE2, and predicted the ΔΔG = -0.81kcal/mol as shown in Table 2. Comparing to the FEP result from Pavlova et al. [76], our prediction on N501Y can agree with theirs well, while the effects of E484K and K417N reversed although neither of them contributed significantly to the free energy change comparing to that of N501Y. However, the effects of E484K and P479S predicted can agree with findings from both earlier experimental data by Linsky et al. [77] and computational prediction by Chen et al. [78].

Our FEP calculation predicted that S477N should have a negligible effect on RBD binding with ACE2, while a stronger binding of RBD-S477N variant with ACE2 by Upadhyay et al. Although our prediction deviated marginally from Upadhyay et al., the effect of S477N from Upadhyay et al. was opposite to the binding of RBD and ACE2 by Barton et al. The FEP prediction on K417N is opposite to Upadhyay et al.’s work, although the FEP prediction can agree with the experimental work by Laffeber et al. [79] and Barton et al. [80] as shown in Table 2. The ΔΔG result from Laffeber et al. [79] as shown were calculated after converting experimental data K_D for wildtype and mutants to ΔG using the relationship as: ΔG_exp = -RTlnK_D at T = 298K.

N439K, which was found in multiple regions in the world and was a quite common mutation, can enhance the binding of RBD with human ACE2 receptor while evade antibody mediated immunity as found by Thomson et al. [39]. FEP predicted that N439K can decrease the binding free energy of RBD with ACE2 thus should increase the binding affinity between them. Our prediction can agree with Thomson et al.’s result [39].

Based on our FEP prediction, both N501Y and K417N can increase the binding affinity of RBD and ACE2, but E484K has a negative effect on the binding. Because of that, the triple mutants of K417N-E484K-N501Y should contribute to the binding of RBD with ACE2 if assuming a linear relationship between the effect of each mutation and the total effect. Because of that, it is expected that ACE2 should bind with the triple mutant stronger than the wildtype. That agrees with the infectivity of Beta variant. The FEP prediction on N501Y suggests that Alpha variant should be more infectious than the original SARS-CoV-2. That also agrees with experimental and clinic findings.

5). Dihedral angle correlation covariance result

In order to find out how point mutations affect the dynamic correlation network of RBD, dihedral angle ϕ, φ combined covariance matrix for RBD was calculated in both the wildtype and variant forms, and in both bound and free states. A set of highly correlated residues in the dihedral angle covariance matrices are revealed as in Fig 9(A) for RBD-N501Y and in Fig 9(B) for RBD-tri-mutant. The RBD wildtype bound with ACE2 covariance matrix is shown in the lower-right triangle below the diagonal line while the covariance matrices of RBD-N501Y variant and tri-mutant variant are shown in the upper-left triangle above the diagonal line in Fig 9(A) and 9(B). RBD-N501Y variant lost the original correlation in RBD wildtype between RES501 and other residues as pointed out by the orange arrows, but other correlations in the range of RES470 to RES483 became stronger in RBD-N501Y variant (pointed out by red arrows) which were weak in the original RBD wildtype. Another region around RES370 (pointed out by red arrows) which is also faraway from RES501 site, shows a strong correlation in RBD-N501Y variant, although such kind of strong correlation was not present in the RBD wildtype. Overall, a stronger correlation shows up in RBD-N501Y variant comparing to the RBD wildtype. Similarly, RBD-trimutant (three mutations as N501Y, E484K, K417N) (Fig 9(B)) has more discrete correlations (pointed out by red arrows) which were not originally in the RBD wildtype, while the original correlations between N501, E484 and K417 and other residues (pointed out by orange arrows) were lost. RES501 region in RBD-tri-mutant is strongly correlated with other regions as pointed out by red arrows.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9.

Comparison of the averaged mixed dihedral angle covariance matrix for RBD-N501Y (upper-left triangle above the diagonal line) and RBD in wildtype (lower-right triangle) in bound state (a), and the averaged mixed dihedral angle covariance matrix for RBD-tri-mutation (upper-left triangle above the diagonal line) in comparison of RBD in wildtype (b). The color bars of covariance map are shown on the right of the figures from blue to yellow corresponding to the covariance coefficient in the range of -0.1 to 0.1. The secondary structure of RBD is shown on top of the covariance matrices, and the residue numbers of RBD are labelled on both the x and y-axes.

https://doi.org/10.1371/journal.pone.0289432.g009

The correlation matrices for RBD-E484K variant and RBD wildtype in free form were analyzed with results in comparison with RBD wildtype in bound state shown in Fig 10(A) and 10(B), and results for other RBD variants in S10–S12 Figs. It was interesting to see that RBD-E484K has strong correlation in the head (RES334-374) and tail (RES520-526) regions which are faraway from each other, comparing to the RBD wildtype. The RBD wildtype in free state has a stronger correlation at its tail region but weaker correlation at its head region than the RBD wildtype in bound state. RES372 region has strong correlation with other residues. Since the RBD free state simulation starts from the crystal structure of RBD (by removing the bound ACE2 in PDB ID of 6M0J) which is in open state, while binding of RBD with ligand can make the binding region more rigid [49], that suggests that binding with ACE2 weakened the correlation of RBD in its tail region which is the binding interface with ACE2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10.

Comparison of the averaged mixed dihedral angle covariance matrices for RBD-E484K (a) and RBD wildtype in free state (b) (upper left triangle above the diagonal line) with RBD wildtype in bound state (lower right). The color bars of covariance map are shown on the right of the figures from blue to yellow corresponding to the covariance coefficient in the range of -0.1 to 0.1. The secondary structure of RBD is shown above the covariance matrices, and the residue numbers of RBD are labelled on both the x and y-axes.

https://doi.org/10.1371/journal.pone.0289432.g010

Similarly, other RBD variants have their correlation matrix changed comparing to the RBD wildtype as results shown in S10–S12 Figs.

Since RES373 on RBD wildtype in bound form, RES372 on RBD wildtype in free form, RES480 on RBD-N501Y, and RES501 on RBD-trimutant have strong correlation with other residues as shown in Figs 9 and 10, those residues were picked as base for RBD in different simulations correspondingly. Mapping the covariance coefficients of the base residue with other residues on the b-factor of the secondary structure of RBD, the results for RBD wildtype in bound form, RBD wildtype in free form, RBD-N501Y, RBD-trimutant are shown in Fig 11(A)–11(D) correspondingly (The secondary structures mapped with covariance coefficients for other RBD variants not shown). In Fig 11, residues having a covariance coefficient higher than 0.05 were labelled, and in red for a positive correlation coefficient while in blue for a negative correlation coefficient. Overall, RBD wildtype in free form, RBD-N501Y, RBD-trimutant all have more correlations than the RBD wildtype in bound form. In the secondary structure of RBD-N501Y with the covariance coefficients mapped as shown in Fig 11(C), RES480 has strong correlation with P479, G476, F456, S469, N440, I441, T500, Y501, P507, G485, H519, C525, G526, L335, N334. Besides that, other residues labelled in red and blue also have correlation with E484. Checking all the secondary structures mapped with covariance coefficients (picking different base residues) for RBD in different forms, some common correlations exist among residues, such as S366/A363, S373/A372, I402/R403, N440/L441, S443/K444, E465/F464, A475/G476, G485/RES484, E516/L517, H519, P521/A522, G526/C525 in different RBD covariance matrices. Those residues are located at the head, binding interface and tail regions, and are on either α-helices or loops on RBD. Thus it suggests that key residues on the helices and loops at the head, binding interface and tail regions have stronger intrinsic dynamic correlations than the β sheet region.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11.

Covariance coefficients between residue 373 (a)/residue 372 (b)/residue 480 (c) and residue 501(d) and other residues mapped on the structure of RBD wildtype in bound state (a), RBD wildtype in free state (b), and the N501Y variant (c), and tri-mutant variant (d). Positive covariance is shown in red while negative covariance in blue and the covariance equals to zero in white. The range of covariance is in the range of -0.1 to 0.1. Residues have positive correlation with the picked residue (or base residue) are labelled in red while in blue for negative correlation.

https://doi.org/10.1371/journal.pone.0289432.g011

Because there are interactions between atoms inside a protein, the protein can be viewed as a network. Treating the Cα atom in each residue as a node, the interaction between nodes inside the network can be analyzed. Multiple nodes with strong interactions can form one community, which is a subnetwork that partitions the original network in the protein. The communication communities can be calculated using VMD program. The dynamical network and communication communities in RBD wildtype and variants were analyzed based on residue−residue interactions over time. It was found that usually 7 to 9 network communities formed in RBD in different forms and at different states. The results for RBD-N501Y, RBD wildtype bound form, RBD wildtype free form are shown in Fig 12 and in S13 and S14 Figs for other RBD variants. Interestingly, the largest community (shown in blue) always extends from the tail region of RBD to the binding interface in the community maps. Such kind of extensive community showed up in all other RBDs (as their network community maps shown in Figs 12 and S13 and S14). That suggests that there should be an intrinsic dynamic network community inside RBD distributed from the tail region to the central, and to the binding interface, which is consistent with the dihedral angle covariance result shown in Fig 11. Thus, the dynamics in RBD binding interface region should be transmitted to the tail of RBD thus to the other sections of S protein. It is expected that disrupting of this intrinsic network can affect the binding of RBD with ACE2 and with antibodies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12.

Community networks formed in the RBD-N501Y(a), RBD wildtype in bound form (b) and RBD wildtype in free form (c) based on MD simulation and dynamical network analysis. The RBD network structures are oriented in the same direction as the RBD shown in Fig 1.

https://doi.org/10.1371/journal.pone.0289432.g012

The largest community includes the following residues in common: L335, A348, S373, D428, N437, E471, F490, Q493, Y495, G496, E516, H519, P521-G526. Those residues are at the N-term, central, binding region, and tail of RBD. It is expected that those residues could be the mutation sites in the future pandemics. Up to now, among those residues, S373, Q493, G496 are the Omicron RBD mutation sites as S373P, Q493R, G496S.

Discussion

Based on the findings from Chan et al. [45], a larger BSA suggests a stronger binding of RBD with ACE2, thus it is expected that E484K, N439K and S477N variants should form a stronger binding with ACE2 than other variants, while the tri-mutant should form a less stable binding with ACE2 based on results shown in Fig 6. However, those cannot agree with the ΔΔG result predicted from FEP as shown in Fig 8, which shows that N439K should decrease the binding free energy of RBD and ACE2 (cause stronger binding), but not E484K. S477N almost does not change the binding free energy of RBD and ACE2. Thus, the size of BSA does not correlate to the binding affinity of RBD and ACE2 based on our test in this project.

Up to now, different results on point mutation effects on the binding affinity of RBD with ACE2 were reported. Tian et al. [81] found that N501Y mutation can increase the binding affinity of RBD with ACE2 receptor using both experimental and simulation methods, thus contribute to a higher rate of transmission of SARS-CoV-2 variant. However, it was found that E484K did not contribute to the binding of RBD with ACE2. That is different from Wang et al. [82], who found that E484K mutant on RBD can increase the binding affinity of RBD with ACE2 receptor, while decreasing the binding affinity of RBD with antibodies. Upadhyay et al. [41] found that mutants differ in their expression levels. Of the eight best-expressed mutants, two (N501Y and K417T/E484K/N501Y) showed stronger affinity to ACE2 compared with the wildtype, whereas four (Y453F, S477N, T478I, and S494P) had similar affinity and two (K417N and E484K) had weaker affinity than the wildtype. Compared with the wildtype, four mutants (K417N, Y453F, N501Y, and K417T/E484K/N501Y) had weaker affinity for the CC12.1 antibody, whereas two (S477N and S494P) had similar affinity, and two (T478I and E484K) had stronger affinity than the wildtype. Mutants also differ in their thermal stability, with the two least stable mutants showing reduced expression. Taken together, these results indicate that multiple factors contribute toward the natural selection of variants, and all these factors need to be considered to understand the evolution of the virus. In addition, since not all variants can escape a given neutralizing antibody, antibodies to treat new variants can be chosen based on the specific mutations in that variant. Wang et al. [83] indicated that T478I and N501Y substitutions were two S mutations important for receptor adaption of SARS-CoV-2, potentially contributing to the spillover of the virus to many other animal hosts. Therefore, more attention should be paid to SARS-CoV-2 variants with these two mutations. In this project, it has been found that point-mutations on RBD do not change the binding interface between ACE2 and RBD comparing to the wildtype. However, each point-mutation can make a specific contribution to the binding free energy of RBD with ACE2. N501Y can increase the binding of RBD with ACE2 significantly based on FEP prediction. That agrees with some other researchers as shown in Table 2. Such kind of discrepancy emphasized the importance of theoretical prediction and systematical comparative study on the SARS-CoV-2 point-mutation variants.

Based on the dihedral angle covariance matrix calculation, it was found that there are some common correlations in RBD, no matter in variant or wildtype form, in bound or free state. Those residues are mostly located in the largest network community predicted by dynamic network analysis with results shown in Figs 12 and S13 and S14. Since the community is intrinsic in RBD, and it extends from the tail region of RBD to the binding interface region, it is expected that mutation on those residues can affect the binding of RBD with ACE2 and with antibodies.

Since the spike protein of SARS-CoV-2 is a trimer and each monomer has more than 1000 residues [84], more work should be done in order to understand how the point mutations on S protein affect its structure, dynamics and function thus to supply insight on designing novel drugs against future pandemics.

Conclusions

In this project, all-atom molecular dynamics simulations using NAMD program were conducted on seven RBD single mutants and one RBD triple mutant binding with ACE2. It was found that point mutations on RBD do not affect the binding interface of RBD with ACE2 significantly. RBD and ACE2 still form similar number of hydrogen bonds on the binding interface after point mutation, although some of the exact residues forming hydrogen bonds changed. The distance map result shows that RBD mutant and ACE2 form the binding interface consistent with the RBD wildtype, and they also form a similar buried surface area on the interface. Thus, the structures of RBD mutants bound with ACE2 are very consistent with RBD wildtype with ACE2. Calculating ΔΔG of point mutation using FEP method, it was found that N501Y, N439K, K417N can enhance the binding of RBD with ACE2; E484K and P479S could weaken the binding of RBD with ACE2; S477N and T478I do not affect the binding of RBD with ACE2. The point mutation changed the dynamic correlation network of residues on RBD based on the dihedral correlation calculation, but an intrinsic network community exists in RBD which extends from the tail region to the central region, then to the binding interface region. The result in this project can supply the methodology and molecular insight on studying the molecular structure and dynamics of possible future pandemics and design novel drugs.

Supporting information

Acknowledgments

The simulations were performed using supercomputer time from XSEDE via an award to Zhang L. (MCB160041), and with some trajectory analysis partially using the high performance computers at UTC, and one simulation performed on Unity computers of URI.