Investigation of the binding and dynamic features of A.30 variant revealed higher binding of RBD for hACE2 and escapes the neutralizing antibody: A molecular simulation approach.

With the emergence of Delta and Omicron variants, many other important variants of SARS-CoV-2, which cause Coronavirus disease-2019, including A.30, are reported to increase the concern created by the global pandemic. The A.30 variant, reported in Tanzania and other countries, harbors spike gene mutations that help this strain to bind more robustly and to escape neutralizing antibodies. The present study uses molecular modelling and simulation-based approaches to investigate the key features of this strain that result in greater infectivity. The protein-protein docking results for the spike protein demonstrated that additional interactions, particularly two salt-bridges formed by the mutated residue Lys484, increase binding affinity, while the loss of key residues at the N terminal domain (NTD) result in a change to binding conformation with monoclonal antibodies, thus escaping their neutralizing effects. Moreover, we deeply studied the atomic features of these binding complexes through molecular simulation, which revealed differential dynamics when compared to wild type. Analysis of the binding free energy using MM/GBSA revealed that the total binding free energy (TBE) for the wild type receptor-binding domain (RBD) complex was -58.25 kcal/mol in contrast to the A.30 RBD complex, which reported -65.59 kcal/mol. The higher TBE for the A.30 RBD complex signifies a more robust interaction between A.30 variant RBD with ACE2 than the wild type, allowing the variant to bind and spread more promptly. The BFE for the wild type NTD complex was calculated to be -65.76 kcal/mol, while the A.30 NTD complex was estimated to be -49.35 kcal/mol. This shows the impact of the reported substitutions and deletions in the NTD of A.30 variant, which consequently reduce the binding of mAb, allowing it to evade the immune response of the host. The reported results will aid the development of cross-protective drugs against SARS-CoV-2 and its variants.

Introduction

COVID-19 is a disease instigated by a newly emerged, vastly infectious coronavirus, SARS-CoV-2, first reported in the city of Wuhan, China. The common symptoms of SARS-CoV-2 are fever, headache, shortness of breath, and cough. The spike protein acts as an entrance factor for the SARS-COV-2 infection in humans. There are six open reading frames in the coronavirus genome which encode for proteins, ordered from 5′ to 3′, M (membrane), E (envelope), S (spike), and N (nucleocapsid), as well as nonstructural proteins like protease and RNA-dependent RNA polymerase (RDRP) [1,2].

The coronavirus infection in cells begins when the spike protein binds to the host angiotensin-converting enzyme 2 (ACE2). Human ACE2 (hACE2) is mostly expressed in the lungs, kidneys, and small intestine, which might result in significant sickness [3]. After binding with ACE2, the host cell proteases split the SARS-COV-2 S-protein into the S1-ectodomain and S2 membrane-anchored domain, which are located at the N-terminal and C-terminal respectively. The S1 subunit aids in the recognition of cell surface receptors as well as facilitating virus entry into the cell [[4], [5], [6]]. The Mac-I domain of the SARS-COV-2 NSP3 protein has been identified as a promising therapeutic target due to its function in altering the innate immune response and increasing virulent qualities [7]. The Mac-1 domain of NSP3 was shown to enhance the response of interferons for viral neutralization and to play a critical role in viral infectivity [8]. Further studies revealed Mac-1 hijacking of the host immune response, which is accomplished by deactivation of the IFN pathway and suppression of STAT1 [9,10].

The SARS-CoV-2 are evolutionary gamblers that lacks a proofreading capabilities thus prone to enormous mutation frequency (million times > DNA-containing cells). Many variants of the SARS-CoV-2 virus have now been detected, as part of the ongoing worldwide pandemic, and have been associated with rapid transmissibility, evasion of immune responses, high occurrence of indisposition, and recurrence [11]. There are several variants of concern (VOC) including the alpha (B.1.1.7) variant, with 40–80% enhanced transmission, as well as beta (B.1.1.351), gamma, delta, and B.1.1.529 variants [12,13]. K417T, E484K, and N501Y mutations were detected in the receptor binding domain (RBD) of the P.1 variant reported in early 2021 in Brazil, which showed enhanced lethality, by 50%, and transmission, by 38% [14]. Another strain with significantly increased lethality was discovered in India in late 2020, officially known as B.1.617.2, which exhibited transmissibility enhanced by 87% [15].

The recently emerged B.1.1.529 strain of SARS-CoV-2 arose in November in South Africa and has spread to other countries, posing a serious public health threat. The spike protein of the B.1.1.529 strain has 30 different alterations including two deletions and seven substitutions found in the N-terminal domain (NTD), and fifteen substitutions found in the RBD. Several other strains, including new variants, have been identified to have the B.1.1.529 spike protein, which might affect SARS-CoV-2 transmission and infection. Recently, using structural modelling and simulation approaches, it has been reported that the binding affinity of B.1.1.529 spike protein RBD is enhanced for ACE2 and reduced for mAb [16].

Strain A.VOI.V2, also known as the A.30 variant, was first reported in several patients in Sweden and Angola in 2021, and reportedly originates from Tanzania. According to a recently published report on the cellular and molecular immunology of the A.30 variant, it evades vaccine-induced antibodies and has enhanced extra pulmonary spread, which is a great concern for public health response. The A.30 variant has 5 deletion and 10 substitution mutations in the N terminal domain of the S1 protein, which acts as an antigenic epitope and can be bound by neutralizing antibodies [17]. Furthermore, three mutations were discovered in the RBD, which interacts with the human ACE2 receptor and is a primary target of antibodies for neutralization. T478R and E484K are two of these mutations, occurring at the ACE2 binding site; however, the letter is responsible to develop resistance to antibody-mediated neutralization [18].

The application of molecular modelling and simulations in deciphering the molecular mechanism of pathogenesis for different variants of SARS-CoV-2, and therapeutics development has been invaluable. Here, we have also employed bio-molecular and simulation tools to comprehend the impact of these mutations on the binding of RBD to the host receptor ACE2, and of monoclonal antibodies to the NTD of A.30 variant. Protein coupling, simulation, and post-simulation investigations demonstrated the binding and dynamic features of this variant in contrast to the wild type. We explored the binding variations between the wild type and A.30 variant RBD and NTD interaction with the host receptor ACE2 and mAb. The current study will aid the development of prophylactic therapeutics using these key features for structure-based and rational drug design.

Methodology

Structure retrieval and mutant modelling

The spike protein, which is required for viral interaction with the host cell, was reported to have mutations in emerging strains. The recently submitted SARS-CoV-2 amino acid sequence was retrieved from the UniProt database (accession number P0DTC2) for identification of mutation loci in newly evolved variants [19]. The wild type protein structure of SARS-CoV-2 spike protein (6M0J) was retrieved from the PDB database and AlphFold2 was used to model the predicted mutations in the wild type protein structure [20,21]. The NTD and mAb structures were also collected from RCSB using 7C2L accession number [20].

Binding network and dissociation constant (KD) determination

The HADDOCK (high ambiguity-driven protein-protein docking) online server was used to investigate the binding differences between the wild type and novel strain A.30 RBD and NTD [22,23]. A guru interface was used for the docking where all the available features are executed for the best docking. A restraint docking by defining the interface residues i.e. 21, 24, 27, 28, 30, 35, 38, 79, 80, 82, 83, and 353 for ACE2, 449, 453, 455, 456, 486, 487, 489, 493, 496, 498, 500, 501, 502, and 505 for RBD, 25–32, 25–32, 51–58, and 100–116 for mAbs, while a region between 145 and 150 was defined for NTD. This server generate different docking conformations, which are then clustered according to the related orientations, and the final conformation is based on higher number of cluster size and lowest Z score. The PDBsum online server (https://www.ebi.ac.uk/) was used to determine the different binding networks (hydrogen bonds, salt bridge and non-bonded contacts) between the wild type and mutant proteins [24]. Subsequently, to provide more accurate information about the strength of the binding complex between the spike protein and ACE2 receptor of both mutant and wild type strains, we used the online server PROtein binDIng enerGY prediction (PRODIGY) to calculate the dissociation constant (KD) and the binding affinity [25].

Molecular dynamics simulation

The molecular dynamics behavior of both wild type and mutant complexes was analyzed with the AMBER20 package, which uses the FF19SB force field [[26], [27], [28]]. The system was neutralized by 23 sodium ions added to each system and the system solvation was performed by using the OPC (optimal point charge) water box model (9180 water molecules). Afterward, the system energy was minimized for the removal of bad clashes. 3000 cycles of conjugate gradient algorithm and 6000 cycles of steepest descent algorithm were used and, after heating to 300 K, 1 atm of constant pressure with weak restraint was used for system equilibration [29,30]. Finally, the 300 ns production was run [31]. A time step 0.02 was used for the simulation. The CPPTRAJ package of AMBER20 was used to analyze trajectories and CUDA was used to run the simulation [32]. For structural stability, root mean square deviation (RMSD) analysis, as a function of time, was performed using the following equation:Where:

di = position difference between atoms,

= reference and superimposed structure.

The residues flexibility was indexed by estimating the root mean square fluctuation (RMSF) employing B-factor [33], which is the most imperative constraint for calculating the flexibility of all residues in a protein. Numerically, the flexibility can be calculated with the given equation:

The radius of gyration (Rg) measures an overall size of the protein during simulations. For the calculation of radius of gyration, the following equation was used.Where.

mi = mass of the atom i,

ri = distance of atom i.

Analysis of total binding free energies

MM/GBSA method, which is reported to be the most accepted protocol for the estimation of the complex free energy in molecular interactions. The protocol calculated the BFE for the WT and mutant RBD in association with hACE2 [[33], [34], [35], [36], [37]].

Different contributing components of total binding energy were calculated by the following equation:

It has a wide range of applications i.e. used to estimate the binding energy for proteins in different studies including SARS-CoV-2 and neurological disorders [[38], [39], [40], [41], [42], [43]].

Results and discussion

Structural modelling and analysis

The emergence of more lethal variants of SARS-CoV-2, that can also have increased transmissibility, chance of re-infection, and evasion of the immune response, has been ongoing since the global pandemic was declared in March 2020 and unceasing reports of emerging variants has aggravated global public health. Many variants, including B.1.1.7, B.1.1.529, P.1, B.1.351, have been reported in different geographic locations around the world with a diverse set of mutations that help the virus to improve its survivability. Recently, a heavily mutated strain of SARS-CoV-2, formally known as A.30, has been reported to have significantly increased transmissibility and immune evasion. This variant has R346K, T478K and E484K mutations in the RBD of the spike protein that help the virus to increase binding affinity and, therefore, infectivity. Additionally, D80Y, D215G, L249 M, and W258L mutations, and deletions at loci 144, 246, 247 and 248, have been reported in the NTD. The distribution of these mutations in the spike protein is shown in Fig. 1 A.

Fig. 1

(A) Mutational landscape of A.30 variant RBD and NTD. (B) Superimposed structure of the wild type RBD and A.30 RBD with R346K, T478K and E484K mutations. (C) Demonstrates the superimposed structures of the wild type NTD and A.30 NTD.

Although previous studies have classified A.30 as an antibody escaping strain, however, atomic-level insights are needed to understand the key differences between the wild type and A.30 variants to guide structure-based drug design. In this study, we investigated the effect of these unique mutations on the structure, function, and binding of the RBD to ACE2, and of mAbs to the spike protein NTD, in great detail. To illustrate the interaction variations and associate them with the infection rates of the A.30 variant, the present work combines integrative protein-protein modelling and molecular simulation techniques. Prior to protein-protein docking, AlphFold2 was used to model the 3D structures of A.30 RBD and NTD for comparison with wild type. A RMSD difference of 0.225 Å was seen between the wild type and A.30 RBDs, while a 0.662 Å RMSD difference was recorded between NTDs. This indicates that the variants have undergone secondary structural reconfiguration, resulting in a new strategy for binding and infection. The reference and modelled structures were superimposed (wild type-RBD, A.30-RBD, wild type-NTD, and A.30-NTD) as shown in Fig. 1B and C.

Protein-protein docking of the wild type and A.30 variant RBDs with human ACE2

Prediction of protein complexes for bound ACE2 and RBD was carried out using HADDOCK web docking server. The binding interfaces for RBD-ACE2 and NTD-mAb were defined for constraint docking and accuracy. Nine residues of hACE2 and 13 residues of the RBD directly interacted with each other. The regions 140–150 were defined for the NTD and the three essential regions, CDR1-3 (25–32), (51–58), and (100–116), from mAb were selected as the binding interface residues for mAb and NTD. Previously these residues have been used to study the interaction of RBD from different variants such as Alpha, delta, delta+, and Omicron variants [20,44,45].

For the wild type RBD-ACE2 complex, the docking score was predicted to be −111.8 ± 1.5 kcal/mol, which has been reported previously by other studies [16,35,46,47]. In contrast, the docking prediction for the A.30 RBD-ACE2 complex reported a value of −129.6 ± 1.8 kcal/mol ( Table 1 ). The variations in the docking scores demonstrate the impact these substitutions have on the function, with a higher binding affinity seen for the A.30 RBD towards ACE2. This report strongly aligns with the already published data on other SARS-CoV-2 variants, including B.1.1.7, B.1.351, P.1, B.1.617, and B.1.618, where the spike mutations, particularly in the RBD, decreased the docking score compared to wild type [35,44,[47], [48], [49]].

Table 1

Estimated docking affinities and KD (dissociation constant) for wild type RBD, wild type NTD, A.30 RBD and A.30 NTD.

Parameters	Wild Type-RBD	A.30-RBD	Wild Type-NTD	A.30-NTD
HADDOCK score	−111.8 ± 1.5	−129.6 ± 1.8	−63.6 ± 5.2	−51.4 ± 1.3
Cluster size	51	160	9	8
RMSD (Å)	14.6 ± 0.2	14.5 ± 0.4	13.0 ± 0.3	21.6 ± 0.3
Van der Waals energy	−48.1 ± 1.3	−54.0 ± 4.1	−55.2 ± 4.6	−34.8 ± 5.6
Electrostatic energy	−169.7 ± 13.2	−248.1 ± 31.6	−131.8 ± 17.0	−165.5 ± 23.1
Desolvation energy	−30.0 ± 3.4	−26.3 ± 5.2	−1.6 ± 1.7	−4.8 ± 2.8
Restraint's violation energy	1.9 ± 1.1	3.1 ± 4.0	195.5 ± 40.9	26.6 ± 24.7
Buried surface area (A2)	1661.1 ± 57.7	1829 ± 72.9	1341.8 ± 21.7	1105.2 ± 80.7
Z-score	−1.6	−1.7	−1.5	0.7
KD (dissociation constant)	3.2E−09	1.3E−10	1.9E−09	5.4E−08
Ref	[16]		[16]

We further compared the van der Waals (vdW) and electrostatic energies for both complexes, where wild type reported a vdW of −48.1 ± −1.3 kcal/mol, while the A.30 complex demonstrated a vdW of −54.0 ± −4.1 kcal/mol. When determining the electrostatic energy, the wild type complex reported −169.7 ± −13.2 kcal/mol and for the A.30 complex, −248.1 ± −31.6 kcal/mol. This suggests that the higher binding affinity of the A.30 variant is due to the enhanced contribution by vdW and electrostatic energy. Increased electrostatic energy has been seen in other variants and is regarded as the main cause of stronger interactions in these complexes [16,35,44,47,49].

Further study of these complexes showed both had 11 hydrogen bonds, however, wild type RBD-ACE2 had only one salt bridge with hACE2, while the A.30 variant complex was reported to have three. Both the complexes reported some conserved interactions, while additional interactions were also seen in the A.30 RBD-ACE2 complex only. As given in Fig. 2 A, the interaction Tyr83-As487 was seen in the wild type but altered in the A.30 complex where Tyr83 was shown to interact with Gln498 and Glu496. The interaction of Lys417 in the RBD for both wild type and A.30 variant was strongly conserved for its interaction with Glu30. A single salt bridge between Glu30-Lys417 was seen in both complexes, while two additional salt bridges between Glu38-Lys484 were reported in the A.30 complex. Additionally, the interaction Glu35-Arg493 was only reported in the A.30 variant and has been shown to be accountable for the anchor locking and correct orientation of the RBD when binding to ACE2 [35,50]. Likewise, the hydrogen bond established by Lys353 in each complex determines the role of these residues to recognize and bind to ACE2 [35,50]. The conserved salt bridges indicate the functional significance of Glu38 and the substitution Lys484 that essentially increase the binding of wild type and various variants [35,47,49]. The binding variations displayed by both the complexes, and particularly for the A.30 variant, shows the importance of key substitutions and consequent binding alteration, forming an altered approach towards binding and host cell entrance. The 3D-interactions of wild type RBD and A.30 RBD complexes are shown in Fig. 2A and B. The docking scores for each complex are given in Table 1.

Fig. 2

Interaction profiling of wild type and A.30 RBD-ACE2 complexes. (A) Shows the binding mode of wild type-RBD with ACE2. (B) Shows the binding mode of A.30-RBD with ACE2.

Protein-protein docking of the wild type and A.30 NTD variant with mAb

Evaluation of the binding variations between the wild type and A.30 NTDs with mAb was also performed using the HADDOCK docking approach. The wild type NTD-mAb complex demonstrated a docking score of −63.6 ± 5.2 kcal/mol, while the A.30 complex reported −51.4 ± 1.3 kcal/mol. This variation in the HADDOCK docking scores demonstrated conformational changes persuaded by the deletions and replacement of key amino acids, which consequently, altered binding. Investigation into the binding differences discovered the impact of the mutations and deletions on the changes in binding conformations. Comparative investigation of the wild type and A.30 NTD-ACE2 complexes show that the key interactions required for recognition and processing by mAbs are lost in the A.30 complex. The variations in the A.30 NTD induced different binding dynamics in the CDR regions of mAb compared to wild type. In the key binding regions of NTD (amino acids 140–150) and CDR 1–3 of the mAbs, A.30 was shown to have lost interactions that are sustained in the wild type. Earlier studies have shown that the deletion of two key residues, Tyr145 and, to a greater extent, His146, reduces the identification and targeting of the B.1.618 variant, corroborate the hypothesis that the deletion of key residues helps the virus to evade neutralization by escaping the antibody response [47]. The interaction patterns of the wild type and the A.30 NTD in complex with mAb are shown in Fig. 3 A and B. The docking scores for each complex are given in Table 1.

Fig. 3

The comparative binding analysis of wild type and A.30 NTD to the mAb. (A) Shows the binding mode of wild type NTD with mAb while (B) demonstrates the binding of A.30 NTD with mAb.

Determining the dissociation constant (KD) for RBD and NTD

The binding strength of the two A.30 complexes, RBD-ACE2 and NTD-mAb, was estimated through KD prediction which has been previously used for other variants such as B.1.1.7, B.1.351, P.1, B.1.617, and B.1.618 [35,47,49]. This study, which employed the same approach as these previous reports, revealed a KD of 1.3 E−10 for the A.30 RBD-ACE2 complex and 1.3 E−10 and 5.4 E−08 for the NTD-mAb complex. When compared, the RBD of A.30 binds more strongly than the wild type, while the NTD of the wild type binds more strongly than that of A.30. The results strongly corroborate the docking score analysis as well as earlier studies based on other variants [35,47,49]. The docking scores and KD values for the wild type RBD and NTD, and the A.30 RBD and NTD complexes are given in Table 1.

Conformational dynamics of the wild type and A.30 RBD

Dynamics stability analysis of RBD

The conformational dynamics of wild type and A.30 variant were investigated through root mean square deviation (RMSD) based on carbon alpha atoms. This analysis is fundamental in understanding the binding affinity of wild and variant spike proteins to the ACE2 receptor and is vital for explaining the effects of mutations on the structure, function, and overall binding characteristics of the spike protein. This is also important in deciphering the virus infection ability in terms of attachment, replication, and transmission. The RMSD was estimated for a 300 ns (ns) period and is presented in Fig. 4 A. As shown in Fig. 4A, the wild type displayed constant equilibrium and revealed stable dynamics with a net RMSD of ∼1 Å (Å) throughout the simulation time. The maximum RMSD reported for the wild type spike system was ∼1.5 nm at 150 ns? This suggests that the intermolecular interactions were strong, keeping the binding conformation stable. For the A.30 RBD, similar dynamics were seen, however, a short period of instability occurred between 150 and 200 ns? These variations corresponded to structural adjustments to achieve better binding conformation between the molecules and involved the binding/unbinding of chemical interactions. However, good binding stability between the variant RBD and ACE2 was seen towards the end of the period, showing stable dynamics for the system.

Fig. 4

Structural and dynamic stability analysis of wild type/A.30 RBD-ACE2 predicted by RMSD analysis. (A) shows the RMSD of wild type/A.30 complexes, (B) Rg plot for wild type/A.30 RBD variants, (C) hydrogen bonds analysis of the wild type and variant.

Our results are consistent with previous research that shows that RBD global stability leads to a greater ACE2 binding capacity [51]. Furthermore, earlier research has shown a close link between the RBD structural stability and binding affinity, with mutations that improve structural stability and stiffness accompanied by increases in binding affinity [52,53]. Other investigations have shown that the mutation C432D in the RBD reduced ACE2 mediated entrance through destabilization of the RBD structure [51]. Recent research on the variants B.1.1.7, B.1.351, P.1, B.1.617, and B.1.618 demonstrated an improvement in structural stability that is significantly linked to variant stability and stronger binding affinity in these variants [35,[47], [48], [49]]. These data substantially support our findings, which demonstrated that A.30 structural stability has increased along with its binding affinity, following the global pattern of greater stability and binding affinity shown in other variants.

The change in the dynamics of A.30 compared to that of wild type might be an adaptation for evading neutralizing antibodies. This is supported by a very recent study demonstrating that humoral antibody response is impaired due to A.30 mutations [54]. These dynamics may also be helpful for the A.30 virus in retaining high transmissibility due to the presence of P681H and D614G mutations. Previously, said mutations were shown to enhance the strain's ability to infect and transmit more aggressively than the parent strain [54]. The A.30 mutations are further reported to allow close interactions between the HR1 and HR2 domains of the virus spike, resulting in the virus fusion with the host cells and initiation of a strong infection cycle [55].

Structural compactness analysis of RBD

To validate the RMSD findings, radius of gyration (Rg) analysis was performed to examine the system's compactness over time (Fig. 4B). As with RMSD, Rg illustrated that the wild type system is more compact and did not experience any major relaxation event in the structure during the simulation time. The mean and maximum Rg of the wild type spike-ACE2 complex was 30.8 Å, while for the A.30 variant; the complex underwent deviation at the same points seen in RMSD. Again, these changes were due to unstable binding in the start and middle period of the simulation time but binding stability was gained later. The mean and maximum Rg of the variant spike-ACE2 complex was 32.5 Å.

Hydrogen bonding analysis of RBD

The protein-protein complexes are holded by hydrogen and hydrophobic bonds. of the interface is always surrounded by water molecules that also compete for hydrogen bonds [56,57]. Considering the higher importance of hydrogen bonds in protein coupling, it is worthy to desipher hydrogen bonding landscape of protein interfaces. For instance, the H-bonds understanding has been previously implemented for SARS-CoV-2 variants to reveal the binding differences. Hydrogen bonding is a very important stabilizing factor in biological systems. These forces are formed when a hydrogen atom is shared between strong electronegative molecules. Fig. 4C shows the number of hydrogen bonds over time, with the average the number of hydrogen bonds throughout reported to be 390 for the A.30 variant. Throughout the simulation time, both studied systems showed the maximum number of hydrogen bonding. This reveals that the interactions between wild type and A.30 variant spike proteins and ACE2 are enriched with strong and close distance hydrogen bonds, keeping both system dynamics stable. Post-simulation hydrogen bonding occupancy for each complex is given in Table 2 .

Table 2

Post-simulation hydrogen bonding occupancy analysis for each complex.

Index	Interaction		Wild type RBD %	A.30 RBD%
1.	LEU19-ALA475		34%	39%
2.	ASP355-THR500		8%	3%
3.	TYR83-ASN487		66%	69%
4.	LYS353-GLY502		71%	52%
5.	GLU38-GLY496		26%	41%
6.	Tyr41-THR500		37%	2%
7.	LYS353-TYR493		64%	12%
8.	LYS353-TYR446		47%	39%
9.	GLU38-LYS417		42%	78%
10.	TYR83-GLY496		4%	64%
11.	TYR83-GLN498		12%	51%
Index	Interactions	Wild type NTD %	Interactions	A.30 NTD%
1.	GLY56-LEU249	28%	TYR27-ASN147	29%
2.	PRO53-LYS147	36%	ASP55-ASN98	32%
3.	ALA101-LYS147	54%	GLU31-ASN147	46%
4.	ALA103-ASN148	69%	GLY104-LYS148	20%
5.	THR105-ASN148	33%	ASP55-ASP176	33%
6.	GLU31-LYS150	57%	GLY56-ASP176	15%
7.	GLU31-LYS150	42%	TYR111-LYS180	31%
8.			GLU54-LYS185	12%
9.			ASP55-ASN186	17%

Residue flexibility analysis of RBD

For understanding flexibility on the residue level, root mean square fluctuation (RMSF) was applied to the simulation trajectories of the systems. It is crucial to examine the deviations of a system's residues to conclude which residues are vital for holding the interaction between spike and ACE2 and providing stabilization to the overall complex. As shown in Fig. 5 A, the mean RMSF for the wild type system and A.30 system were 2.20 and 2.26 Å, respectively. Similarly, for the wild type spike-ACE2 complex, the maximum RMSF observed was at residues 100–200, 300–340, and 450–550, while for the A.30 complex, highly flexible residues also included the region of 650–680, alongside those seen for wild type. Additionally, we calculated the residue flexibility index for the three important loops in the RBD of wild type and A.30 variants, which come in direct contact with ACE2. As given in Fig. 5C-D, the three loops demonstrated differential values for residue flexibility index in each complex. These loops (residues 484–505) were more flexible in the wild type, with minimal fluctuation in the A.30 variant. This shows the A.30 loops flexibility is stabilized by binding with ACE2 and by the mutation induced variation in conformational dynamics. Hence, A.30 shows greater stability in binding ACE2 than the wild type.

Fig. 5

(A) Residue flexibility of wild type and A.30 variants calculated as RMSF. (B–D) represent the flexibility index for the three important loops from residues 484–505.

Conformational dynamic of the wild type and A.30 NTD

Dynamics stability analysis of NTD

We further assessed the stability variations between the wild type and A.30 NTD in complex with mAb. As given in Fig. 6 A, the RMSD of both complexes remained comparable until 150 ns, where an increase was seen for the A.30 NTD-mAb complex. For wild type, the RMSD remained ∼0.85 Å during the first 150 ns, while during this period the RMSD for A.30 was also reported to be over 0.80 Å. After 150 ns, the wild type NTD gained more stability and demonstrated a uniform RMSD for the remaining period. The average RMSD for the wild type NTD-mAb complex was reported to be 0.9 Å. The RMSD for the A.30 NTD-mAb complex continued to increase gradually but reported more structural perturbation than the wild type. The average RMSD increased during the last 150 ns and was calculated to be 1.10 Å. Previously, similar findings were reported for other variants, such as B.1.1.7, B.1.1.617, B.1.1.618, and B.1.1.529, suggesting that the mutations and deletions which change the protein dynamics help the A.30 variant to escape the neutralizing antibodies [16,35,45,46].

Fig. 6

Structural and dynamic stability analysis of wild type/A.30 NTD variant complexes with mAb, predicted by RMSD analysis. (A) shows the RMSD of wild type/A.30 NTD complexes, (B) Rg plot for wild type/A.30 NTD variants, (C) hydrogen bonds analysis of the wild type and A.30 variant NTD.

Structural compactness analysis of NTD

The radius of gyration for both the complexes reported a strong agreement with the RMSD results. As can be seen in Fig. 6B, the Rg graphs for the wild type and A.30 variant reported a similar pattern to RMSD. Wild type and A.30 complexes demonstrated comparable Rg during the first 150ns which then increased during the last 150ns for the A.30 complex only. This shows that the A.30 variant reported significant binding and unbinding events, thus destabilizing the neutralizing antibody, and consequently helping the virus to escape the immune response. These findings strongly align with the previous reports of the other variants [16,35,45,46]. The average Rg for the wild type complex was calculated to be 32.80 Å while for A.30 it was calculated to be 34.2 Å.

Hydrogen bonding analysis of NTD-mAb

Estimation of hydrogen bonding during the simulation provided information about the binding stability incurred by the hydrogen bonds. To calculate the total number of hydrogen bonds in the simulation trajectory, hydrogen bonding analysis was performed. As given in Fig. 6C, the wild type reported more hydrogen bonds than the A.30 NTD-mAb complex. In the wild type complex, the average number of hydrogen bonds was reported to be 248 while in the A.30 complex the average number of hydrogen bonds was reported to be 234. This shows that the mutations and deletions in the NTD of A.30 induce conformational changes that enable the virus to lose important contact with the mAb. Post-simulation hydrogen bonding occupancy for each complex is given in Table 2.

Residue flexibility analysis of NTD

We also predicted the local level residue flexibility for the NTD of wild type and A.30 variants in complex with mAb. As shown in Fig. 7 , the wild type NTD is comparatively more stable than the A.30 NTD. In particular, the A.30 NTD demonstrated higher fluctuation between residues 1–240, while a less substantial fluctuation was displayed by the wild type NTD between residues 420–450. The difference in the flexibility index shows how conformational optimization by these mutations affects the change in binding strength. The RMSF (flexibility) graph for the wild type and A.30 NTD in complex with mAb are given in Fig. 7.

Fig. 7

Residual flexibility of wild type and A.30 NTDs in complex with mAb, calculated as RMSF.

Binding free energy calculation for the wild type and A.30 complexes RBD-ACE2 & NTD-mAb

The MM/GBSA technique for computing the BFE of biological partners is a widely used method for examining the putative docking configuration. This method, which is less costly than the alchemical free energy methods, displays the binding stability of interacting key regions and the BFE. It is also regarded to be more precise than any rational scoring function. We employed the MM/GBSA technique because it allows us to see how the mutations in the spike RBD influence the binding with hACE2 and NTD with mAb. The BFE results are given in Table 3 .

Table 3

MM/GBSA results, which show the binding free energy for each complex. All the energies are represented in kcal/mol.

Complex	vdW	Electrostatic	GB	SA	Total Binding Energy	Ref
Wild Type RBD	−87.75	−616.79	658.08	−11.79	−58.25	[16]
A.30 RBD	−99.23	−1168.78	1215.92	−13.5	−65.59
Wild Type NTD	−78.78	−212.72	236.2	−10.46	−65.76	[16]
A.30 NTD	−84.81	−828.43	875.93	−12.03	−49.35

Binding free energy for RBD-ACE2 complexes

As given in Table 3, the vdW for the wild type RBD-ACE2 complex was estimated to be −87.75 kcal/mol while for the A.30 complex the vdW was reported to be −99.23 kcal/mol. Moreover, for the wild type complex, the electrostatic energy was estimated to be −616.79 kcal/mol, whereas the A.30 complex was estimated to have electrostatic energy of −1168.78 kcal/mol. Previously, a higher electrostatic energy for the variant has also been reported and was considered the prime factor in its enhanced binding to ACE2 [35,47,49]. Finally, BFE for the wild type complex was −58.25 kcal/mol in contrast to the A.30 complex, which reported −65.59 kcal/mol. The higher BFE for the A.30 complex implies vigorous binding with the ACE2 than for wild type, allowing the variant to bind and spread more promptly.

Binding free energy for NTD-mAb complexes

The BFE for NTD-mAb complexes revealed contrasting results to the above RBD-ACE2 complexes. The BFE for the wild type NTD-mAb complex was computed to be −65.76 kcal/mol while for the A.30 complex it was −49.35 kcal/mol. This shows the impact of the reported substitutions and deletions in the NTD of A.30 variant which reduce the binding affinity of mAb towards NTD and helps in evading the immune response instigated by the host. Our findings corroborate with an earlier study that showed that mutations and deletions in the NTD of B.1.618 variants led to immune evasion [47].

Conclusions

In conclusion, the current study employed protein-protein coupling and molecular simulation approaches to decipher the key features required for stronger interaction with the ACE2 and escaping the neutralizing antibodies for the A.30 variant. We reported that the key substitutions altered the binding modes of RBD and NTD towards ACE2 and mAb. The only limitation of the current study is long run replicated simulations are missing which could further provide deep understanding on the reproducible results and basis for accurate interactions. The explored key features can be considered while designing novel therapeutics against this variant. These data will aid in the development of cross-protective drugs against SARS-CoV-2 and its variants.

Funding

Dong-Qing Wei is supported by grants from the Key Research Area Grant 2016YFA0501703 of the of China, the (Grant No. 32070662, 61832019, 32030063), the (Grant No.: 19430750600), as well as JiRLMDS Joint Research Fund and Joint Research Funds for Medical and Engineering and Scientific Research at (YG2021ZD02).

Availability of data and material

All the data is available on RCSB, UniProt and any simulation data would be provided on reasonable demand. The accession numbers to access this data are given in the manuscript.

Ethics approval and consent to participate

N/A.

Consent for publication

N/A.

Conflicts of interest

Declared none.

Declaration of competing interest

Authors declare there is no declaration of interest.