Cooperative and structural relationships of the trimeric Spike with infectivity and antibody escape of the strains Delta (B.1.617.2) and Omicron (BA.2, BA.5, and BQ.1)

Abstract

Herein, we conducted simulations of trimeric Spike from several SARS-CoV-2 variants of concern (Delta and Omicron sub-variants BA.2, BA.5, and BQ.1) and investigated the mechanisms by which specific mutations confer resistance to neutralizing antibodies. We observed that the mutations primarily affect the cooperation between protein domains within and between protomers. The substitutions K417N and L452R expand hydrogen bonding interactions, reducing their interaction with neutralizing antibodies. By interacting with nearby residues, the K444T and N460K mutations in the Spike^BQ.1 variant potentially reduces solvent exposure, thereby promoting resistance to antibodies. We also examined the impact of D614G, P681R, and P681H substitutions on Spike protein structure that may be related to infectivity. The D614G substitution influences communication between a glycine residue and neighboring domains, affecting the transition between up- and -down RBD states. The P681R mutation, found in the Delta variant, enhances correlations between protein subunits, while the P681H mutation in Omicron sub-variants weakens long-range interactions that may be associated with reduced fusogenicity. Using a multiple linear regression model, we established a connection between inter-protomer communication and loss of sensitivity to neutralizing antibodies. Our findings underscore the importance of structural communication between protein domains and provide insights into potential mechanisms of immune evasion by SARS-CoV-2. Overall, this study deepens our understanding of how specific mutations impact SARS-CoV-2 infectivity and shed light on how the virus evades the immune system.

Graphical abstract

Introduction

The rapid evolution of SARS-CoV-2 is dominated by the accumulation of mutations in the viral genome, impacting many aspects of the biology of the virus, such as its speed of propagation, pathogenicity and infectivity [1]. SARS-CoV-2 variants whose emergence is associated with important waves of infection are often designated as variants of concern (VOCs). In order to allow monitoring of the spread of their genetic variants, VOCs are classified into Nextstrain clades and Pango lineages [2]. To facilitate discussions for both technical and non-scientific audiences, the World Health Organization (WHO) recommends the use of letters from the Greek alphabet to match the order of appearance of the SARS-CoV-2 VOCs, such as Alpha, Beta, Gamma, Delta, and Omicron [3]. Thus, VOCs have been classified as Alpha (pango lineage, B.1.1.7; nextstrain clade, 20I), Beta (pango lineage, B.1.351; nextstrain clade, 20H), Gamma (pango lineage, P.1; nextstrain clade, 20 J), Delta (pango lineage, B.1.617.2; nextstrain clade, 21A, 21I, 21 J), and the globally predominant strain, the Omicron variant (pango lineage, B.1.1.529; nextstrain clade 21 M). The Omicron VOCs include several sub-lineages such as BA.1 (nextstrain clade 21 M), BA.2 (nextstrain clade 21L), BA.4 (nextstrain clade 22A), BA.5 (nextstrain clade 22B), BA.2.12.1 (nextstrain clade 22C), BA.2.75 (nextstrain clade 22D), BQ.1 (nextstrain clade 22E), and XBB (recombination between BA.2.10.1 and BA. 2.75, nextstrain clade 22F) [2]. Although the Delta variant has a high degree of transmissibility and infectivity, the Omicron variant quickly became predominant worldwide, thus proving to be even more transmissible than Delta [1]. BA.1 and BA.2 are less lethal in vaccinated individuals [4], but are able to escape the immune system [5, 6]. BA.4, BA.5, and BQ.1 are currently the most predominant strains worldwide and they are related to the most recent surges in cases of COVID-19 [7, 8].

When SARS-CoV-2 variants are compared, mutations can be found in any gene, but a higher number of mutations concentrate in the gene region that encodes the Spike protein [1]. Spike is a homotrimer anchored in the viral membrane composed by S1 and S2 subunits, each performing a different role in the viral infection [1, 9]. The S1 subunit contains the N-terminal domain (NTD), the receptor binding domain (RBD) that play key roles in recognition of the hACE2, SD1 and SD2 domains [1, 9,10,11]. The RBD may be present in two conformational states, up or down. In its up-conformational state, the RBD is capable of binding to hACE2. Conversely, the S2 subunit mediates the fusion between virus and host membranes [1, 12]. Spike may be recognized and cleaved by the furin protease in the junction connecting S1 with S2 (S1/S2 site), at position R682-R685 (RRAR₆₈₅↓) (Figures S1 and S2) [1, 13, 14]. In addition, transmembrane Serine Protease 2 (TMPRSS2) and Cathepsin isoforms also recognize and cleave the S2’ site (KPSKR₈₁₅↓) [9, 15,16,17,18,19]. After hACE2 interaction and S2’ cleavage, the S2 subunit undergoes a major structural change to expose the fusion loop, which mediates the fusion of virus envelope and host cell membrane, allowing the formation of a pore for the release of viral RNA into the cell cytoplasm [20].

Previous experimental studies have shown that, in addition to high binding affinity, other biochemical variables are related to the enhanced viral infection and increased transmissibility of VOCs [1, 9], including: (1) increased structural flexibility that would make it difficult to be recognized by neutralizing antibodies in Spike^BA.1 [6, 21]; (2) independence of TMPRSS2 cleavage, leading to changes in cell tropism [21,22,23]; (3) increased spike density on the viral surface particle [1, 24, 25]; (4) improved exposure of RBD to interact with hACE2 [21, 26,27,28]. Although the Omicron variant is more transmissible than the Delta variant, the evolution of the mortality rates for these variants suggests that the Omicron variant is less lethal than the Delta variant [29, 30].

Our previous studies showed that the Spike^VOCs have similar binding affinities to hACE2, but these are higher than that of Spike^WT [1, 9]. Previously results have shown that substitutions K417N, S477N, T478K, E484K and N501Y cause an increase in binding affinity between RBD and hACE2 [1, 9, 31, 32]. The mutations present in the Omicron variant of SARS-CoV-2 have elicited concerns due to their potential implications for the spike protein structure and its recognition by neutralizing antibodies. Structural investigations have revealed that mutations occurring in the receptor binding domain (RBD) and other regions of the Spike protein impact inter-domain and inter-subunit interactions, resulting in a stable open conformation of the trimeric spike protein. Notably, specific amino acid substitutions have been identified in the spike protein of the Omicron variant, such as from Asn856 to Lys853 and from Asn764 and Thr547 to Lys761 and Lys544, respectively. These mutations have the capability to modify polar interactions and the conformation of the trimeric Spike protein, potentially affecting antibody recognition [32, 33]. In a broader context, mutations in the NTD may significantly alter the antigenic epitope region, conferring resistance to neutralization by antibodies targeting the NTD in the Omicron variant [32]. Indeed, cryo-EM studies conducted on the Spike protein of the Omicron subvariant BA.2 have demonstrated a higher proportion of Spike proteins with the up-RBD position when compared with the BA.1 [6, 33]. These mutations in the BA.2 Spike protein induce a more compact packing of the RBDs, within and between each RBD, thereby contributing to enhanced thermal stability of the BA.2 Spike protein [6]. Furthermore, these structural changes influence the inter-protomer and intra-protomer communication, thereby impacting the transition between down- and up-RBD [5, 6, 32, 33]. Both the BA.1 and BA.2 variants have the potential to disrupt the interaction with neutralizing antibodies, enabling evasion of neutralization by various classes of antibodies directed against the RBD and NTD [34,35,36,37,38]. In general, mutations in the RBD and NTD may induce conformational changes that attenuate interactions with specific antibodies, consequently leading to immune escape and reduced efficacy of neutralizing antibodies [34,35,36,37,38].

The mutations present in Omicron’s Spike protein are diverse and play important roles in the evolution of the virus. For example, the del143-145 results in the removal of a crucial epitope for the 4A8 antibody, reducing its effectiveness in neutralizing the virus [37]. Conversely, the S477N, Q493R, Q496S, Q498R and N501Y substitutions establish additional interactions for the formation of the RBD/hACE2 complex, increasing the binding affinity between the virus and the host cells [37]. Other mutations, such as S371L, S373P, G339D and S375F, directly impact RBD conformation and stabilize its up form, favoring binding to hACE2 [37]. These structural changes give the virus an advantage in entering cells and spreading. The L452R mutation, also present in the Delta variant, interferes with the recognition of the virus by the JMB2002 antibody, making it less effective in neutralizing the virus [37]. In addition, mutations K417N, S371L, S373P and S375F can lead to structural alterations in the Spike protein, decreasing recognition of the S304 antibody [38]. These mutations can affect the ability of the S304 antibody to bind and neutralize the Omicron variant [38]. The Omicron variant is resistant to most clinically approved neutralizing antibodies or induced by vaccines or previous infection [39]. Studies indicate that mutations identified in the Omicron Spike protein result in changes in epitopes recognized by neutralizing antibodies such as AZD1061, REGN10987, REGN10933 and LYCoV555 [37]. Thus, these structural alterations in the Spike protein, together with other mutations, may impact the recognition of neutralizing antibodies. Thus, this immunological evasion may decrease the effectiveness of vaccines directed at neutralizing epitopes present in the S1 region of the Spike protein. Omicron variants BA.4 and BA.5 have additional mutations in the Spike protein, such as L452R and F486V, which compete with the earlier subvariant BA.2 in terms of infectivity and antibody evasion [40, 41]. In general, all these mutations may confer greater capacity for immune escape than BA.2. In general, mutations in the Spike protein, especially in BA.1, BA.2, BA.4 and BA.5, have the potential to alter binding epitopes in the Spike protein, making it difficult to recognize and effectively neutralize the virus [42]. This ability to immune escape raises concerns about the effectiveness of existing antibody therapies and vaccines against these variants. In this regard, it is critical to continue closely monitoring mutations and adapting therapeutic and vaccination strategies to combat Omicron variants and ensure an effective response against SARS-CoV-2.

Mutations may give rise to new variants of the coronavirus, potentially triggering new pandemics. It is crucial to investigate the impact of these mutations on the Spike protein to understand the virus’s evolution and devise strategies to mitigate their exacerbating effects. An effective approach in this regard would involve developing vaccines that specifically target each variant. Generally, novel coronavirus variants exhibit higher mutation rates compared to the wild-type (WT), particularly in the NTD and RBD regions of the Spike protein [34,35,36,37,38]. Additionally, mutations in the S2 subunit of Spike are also observed, which can confer both evolutionary advantages and disadvantages to SARS-CoV-2. Notably, the D614G mutation (present in all VOCs), P681R (found in the Delta variant), and P681H (present in the Omicron variant) are significant. The D614G mutation in the Spike protein provides several evolutionary advantages to the virus, including enhanced entry into host cells, improved structural flexibility, increased hACE2-binding capacity, and a higher density of viral spikes [25, 43,44,45,46,47]. Conversely, the P681R mutation facilitates the cleavage of the Spike protein and enhances the fusogenicity of SARS-CoV-2, leading to greater pathogenicity compared to the wild-type [48,49,50]. The P681H mutation, found in both the Alpha and Omicron variants, results in reduced fusogenicity compared to the P681R mutation [48] and exhibits similarities to the wild-type in certain cell lines [49]. In the Omicron variant and a pseudovirus carrying the P681H mutation, there is an observed independence from TMPRSS2 in viral infection, suggesting a potential change in cell tropism and a decrease in virus infectivity [48]. Although these mutations significantly impact the virus’s infectivity, the precise relationship between these mutations, the structure of the Spike protein, and their effects on the entry mechanism of SARS-CoV-2 VOCs and immune system escape remains unclear.

Despite the abundance of information available on Variants of Concern (VOCs) and their Spike protein, our understanding of how the Omicron variant evades the immune system remains limited. Additionally, the relationship between these mutations and the reduced infectivity of the virus in host cells is still unclear. Our previous study revealed a rapid succession of variant dominance, with BA.1 swiftly replacing the Delta variant [1], followed by BA.2. Subsequently, BA.5 emerged as the dominant variant over BA.2, and eventually, BQ.1 predominated over BA.5 [2]. However, the precise connection between the Delta variant’s higher infectivity, the lower infectivity observed in the Omicron variant and its sub-lineages, increased immune evasion, the prevalence of these current variants, and the evolution of the Spike protein remains uncertain. To gain further insight into these questions, it is crucial to investigate the direct association between these factors and the structure of the Spike protein. Herein, we focused on examining how mutations contribute to conformational changes among Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1. We conducted molecular dynamics simulations of the trimeric ectodomains of these Spike variants for 300 ns. By calculating hydrogen bonding occupancies, we analyzed the impact of specific substitutions that confer resistance to neutralizing antibodies. To evaluate the effect of these mutations on the structural conformations of the Spike proteins, we calculated root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) values and performed principal component analysis (PCA). These analyses enabled us to identify significant variations in conformational flexibilities, particularly in the N-terminal domain (NTD) and the receptor-binding domain (RBD). Furthermore, we examined the dynamical cross-correlation matrix (DCCM) of the simulated trajectories, which revealed alterations in inter- and intra-protomer communication that may contribute to a reduced efficiency of the fusogenic process. We also established a connection between loss of sensitivity to neutralizing antibodies and inter-protomer communication using a multiple linear regression model. Given the potential increased pathogenicity, transmissibility, and ability to evade the immune system exhibited by SARS-CoV-2 VOCs, it is crucial to gain a detailed understanding of how Spike protein mutations impact its structure and function. This knowledge will guide predictions of fitness gains and facilitate the rational development of updated vaccines for effective COVID-19 therapy.

Materials and methods

Homology modeling

We make available all input and/or output files containing parameters and codes for homology modeling, molecular dynamics, principal component analysis, statistical analysis script (in .R), and dynamical cross-correlation matrix in the Github (https://github.com/anacletosouza/molecular-dynamics-B.1.617.2-BA.2-BA.5-BQ.1) and in supplementary files S1-S6. In order to investigate how molecular dynamics may influence transmissibility and infectivity, we simulated the trimeric Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1 ectodomains in a water-filled box. The available structures have mutations and/or unresolved regions. In order to correct these problems, we used cryogenic electron microscopy (Cryo-EM) of the Spike^BA.1 structure (PDB ID 7WK2) as a template to obtain the corrected structures in the Swiss-Model web server [50]. Initially, the genomic sequences of the B.1.617.2 (GISAID, accession code EPI_ISL_15728122), BA.2 (GISAID, accession code EPI_ISL_12905607), BA.5 (GISAID, accession code EPI_ISL_14439279) and BQ.1 (GISAID, accession code EPI_ISL_15731012) were translated in Translate tools (available on the Expasy web server) to obtain the primary sequence of the corresponding Spikes proteins. All primary sequences were submitted in the Swiss Model web service [50,51,52], ignoring the transmembrane regions. Furthermore, we considered the unglycosylated state of Spike because the difference between glycosylated and unglycosylated Spike is smaller than seen between different Spike variants [53]. In order to compare all molecular dynamics simulations, we used only the structures in the down conformation.

Setup for molecular dynamics simulations

In order to study the structural features of the Spikes, initially we determined the protonation states of ionizable residues in an implicitly aqueous environment at pH 7.0 using the PROPKA module implemented in the Maestro program (academic version v. 2021-4, by Schrödinger) [54]. Thus, (1) all glutamic and aspartic residues were represented as non-protonated; (2) arginine and lysine residues were assumed to have a positive charge; (3) in the N- and C-terminal regions, the amino and carboxyl groups have been converted to charged groups; (4) the histidines were protonated according to prediction of PROPKA and manual inspection, such as: Spike^B.1.617.2, δ-tautomer (H75, H664, H1067, H1073, H1092, H1097, and H1110) and ε-tatomer (H58, H78, H216, H254, H528, H634, H1057) and charged histidine in H155 residue; Spike^BA.2, δ-tautomer (H72, H251, H511, H960, H1054, H1064, H1070, 1089, 1094, 1107) and ε-tatomer (H55, H75, H152, H213, H525, H631, H687); Spike^BA.5, δ-tautomer (H72, H509, H629, H958, H1052, H1062, H1068, H1087, H1092, H1105) and ε-tatomer (H55, H150, H211, H249, H523, H685); Spike^BQ.1, δ-tautomer (H72, H149, H508, H628, H957, H1051, H1061, H1067, H1086, H1091, H1104) and ε-tatomer (H55, H210, H248, H522, H684). In Spike^B.1.617.2, the disulfide bridges were considered between residue pairs: C140/C175, C300/C310, C345/C370, C388/C441, C400/C534, C489/C497, C547/C599, C626/C658, C671/C680, C747/C769, C752/C758, C1041/C1052, C1091/C1135. In Spike^BA.2, the disulfide bridges were considered between residue pairs: C137/C172, C297/C307, C342/C367, C385/C438, C397/C531, C486/C494, C544/C596, C623/C655, C668/C677, C744/C766, C749/C755, C1038/C1049, C1088/C1132. In Spike^BA.5, the disulfide bridges were considered between residue pairs: C135/C170, C295/C305, C340/C365, C383/C436, C395/C529, C484/C492, C542/C594, C621/C653, C666/C675, C742/C764, C747/C753, C1036/C1047, C1086/C1130. In Spike^BQ.1, the disulfide bridges were considered between residue pairs: C135/C169, C294/C304, C339/C364, C382/C435, C394/C528, C483/C491, C541/C593, C620/C652, C665/C674, C741/C763, C746/C752, C1035/C1046, C1085/C1129. All structures were previously minimized in Chimera tools, using steepest descent steps of 500, steepest descent step size of 0.02 Å, conjugate gradient steps of 10 and update interval of 10.

GROMACS is a widely used molecular simulation package that employs physical equations to describe the behavior of proteins of interest. The OPLS-AA force field formalism for the total energies of ionic liquid systems is determined by summing up various energy components. These components include harmonic bond stretching and angle bending energies, torsional energetics represented as a cosine series, and nonbonded interactions described by Coulombic and Lennard–Jones potentials. The parameters used in these equations include force constants (k), equilibrium bond and angle values (r₀ and θ_o), Fourier coefficients (V) for torsion, partial atomic charges (q), and Lennard–Jones parameters (σ and ε, radii and well-depths, respectively). The subscripts b and θ represent the binding term and angle term, respectively [55].

$$\begin{aligned} E_{Total} = & \sum\limits_{{\text{i}}} {k_{b,i} \left( {r_{i} - r_{o,i} } \right)^{2} } + \sum\limits_{{\text{i}}} {k_{\theta ,i} \left( {\theta_{i} - \theta_{o,i} } \right)^{2} } \frac{1}{2}\sum\limits_{{\text{i}}} {\left[ {V_{1,i} \left( {1 - cos\varphi } \right)} \right]} \\ & + V_{2,i} \left( {1 - cos2\varphi } \right) + V_{3,i} \left( {1 - cos3\varphi } \right) + V_{4,i} \left( {1 - cos4\varphi } \right)] \\ \end{aligned}$$

$$\sum\limits_{i} {\sum\limits_{i > j} {\left\{ {\frac{{q_{i} q_{j} e^{2} }}{{r_{ij} }} + 4\varepsilon_{ij} \left[ {\left( {\frac{{\sigma_{ij} }}{{r_{ij} }}} \right)^{12} - \left( {\frac{{\sigma_{ij} }}{{r_{ij} }}} \right)^{6} } \right]} \right\}} }$$

This force field employs geometric combining rules to determine the Lennard–Jones coefficients, where σ_ij is calculated as the square root of the product of σ_ii and σ_jj, and ε_ij is determined as the square root of the product of ε_ii and ε_jj.

All molecular dynamics simulations were performed using GROMACS, v. 2022.1 [56,57,58,59], using the OPLS-AA force field [60]. All systems were then explicitly solvated with TIP3P water models in a triclinic box and neutralized, maintaining the concentration of 150 mM NaCl and minimized until reaching a maximum force of 10.0 kJ/mol or a maximum number of steps in 50,000. The systems were consecutively equilibrated in isothermal-isochoric (NVT) ensembles by 2000 ps (number of steps and intervals of 1.000,000 and 2 fs, respectively) and isothermal-isobaric (1 bar, 310 K, NpT) by 2000 ps (number of steps and intervals of 1,000,000 and 2 fs, respectively). All simulations were then performed in a periodic triclinic box considering the minimum distance of 1.2 nm between any protein atom and the walls of the box. Molecular dynamics runs were performed in isothermal-isobaric (1 bar, 310 K, NpT) by 300 ns (number of steps and intervals of 150,000,000 and 2 fs, respectively) in the NpT ensemble. We use the leap-frog algorithm to integrate Newton equations. We selected LINCS (LINEar Constraint Solver) that satisfies the holonomic constraints, whose number of iterations and order were 1 and 4, respectively. We use neighbor search grid cells (Verlet clipping scheme, frequency to update 20-step neighbor list, and clipping distance for 12 Ų short-range neighbor list). In the van der Waals parameters, we smoothly shift the forces to zero between 10 and 12 Å. In electrostatic Coulomb, we use Particle-Mesh Ewald (PME) fast and smooth electrostatics for long-range electrostatics. In addition, we set the distance for the Coulomb cut to 12 Å, order of interpolation for PME to value 4 (cubic interpolation) and grid spacing for Fast Fourier Transform (FFT) to 1.6 Å. In temperature coupling, we use velocity rescheduling with a stochastic term (V-resale, modified Berendsen thermostat). After obtaining two coupling groups (protein and water/ions), we consider the time constant (0.1 ps) and 310 K as the reference temperature. In the pressure coupling (NpT ensembles), we use the Parrinello-Rahman barostat (isotropic type, time constant of 2 ps, reference pressure of 1 bar and isothermal compressibility of 4.5 × 10^–5 bar⁻¹). We saved the compressed coordinates every 10.0 ps. In molecular dynamics calculations, we use periodic boundary conditions in xyz coordinates (3D space). We then calculated the percent hydrogen bond occupancy (all frames, considering cut-off of distance and angle of 4 Å and 20°, respectively), root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) using the GROMACS modules and Visual Molecular Dynamics (VMD) [61]. Analysis of C_α dynamical cross-correlation matrix (DCCM) was performed using the R-based package Bio3d [62]. The correlations were projected in the first frame (t = 0 ns from MD trajectory corresponding to each SARS-CoV-2 VOC). The scripts used to calculate DCCM, statistical test t, and average number of correlations and anti-correlations of intra-protomers and inter-protomers of the Spike proteins of all strains, in R files, are available in the supplementary material S5.

Principal component analysis

PCA (Principal Component Analysis) is a statistical technique used to reduce the dimensionality of a dataset by identifying the main linear combinations of variables that capture the majority of the data’s variation. This enables simplified visualization and analysis of the data while preserving the most relevant information. Principal component analysis was calculated using the GROMACS modules, covar and anaeig. Using simulated molecular dynamics trajectories, we determined the average position of the Spike backbone atoms and calculated covariance matrix fitting non-mass weighted. Covariance matrices were constructed using the backbone of the Spike^B.1.617.2 (30,591 × 30,591), Spike^BA.2 (30,537 × 30,537), Spike^BA.5 (30,483 × 30,483), and Spike^BQ.1 (30,456 × 30,456). From covariance matrices, we diagonalized them to obtain the eigenvalues (whose sum of all eigenvalues were of 510.3, 572.0, 561.8, and 385.8 nm²) and, afterwards, orthogonal eigenvectors. Using the first eigenvalue, we determined the eigenvector 1 (named principal component 1, PC1) that represents the first main movement detected along molecular dynamics trajectories. We projected molecular dynamics trajectories of the Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1 in PC1, skipping 300 frames, to obtain 100 snapshots and to visualize the structural changes of these proteins. In order to quantify these conformational changes of the main movements of the Spike in 310 K, we calculated the entropy by Schlitter method [63] (ΔS^Schlitter), using anaeig of the GROMACS. All input parameters are available in supplementary file S3.

Multiple linear regression model

One approach to establish the association between structural data and the functional activity of SARS-CoV-2 variants of concern (VOCs) is through the utilization of multiple linear regression. In order to increase the amount of data, we utilized the correlation matrices from Spike^WT and Spike^BA.1, which were published in our previous study [21]. The scripts used to calculate the correlations and anti-correlations within and between the Spike proteins of all strains, in R files, are available in the supplementary material S5. Subsequently, we searched for publications containing IC₅₀ values of Sotrovimab tested against all strains studied in this research. IC₅₀ is defined as the antibody concentration in ng/mL required to inhibit 50% of SARS-CoV-2 replication [64,65,66,67,68]. The independent variables considered in our analysis included the average frequency of intra-protomer correlations (corr_intra), the average frequency of intra-protomer anti-correlations (anti_corr_intra), the coefficient of variation for intra-protomer correlations (corr-CV_intra), the coefficient of variation for intra-protomer anti-correlations (anti-corr-CV_intra), the average count of interprotomer correlations (corr_inter), the average frequency of interprotomer anti-correlations (anti-corr_inter), the coefficient of variation for interprotomer correlations (corr-CV_inter), and the coefficient of variation for interprotomer anti-correlations (anti-corr-CV_inter). Based on the structural information and aiming to relate the structure and function throughout the evolutionary process of SARS-CoV-2 VOCs, we developed various models of multiple linear regression (the scripts used to build the models, in R files, are available in the supplementary material S5). Specifically, we focused on models with F-statistics greater than 20, indicating strong overall significance. These models were used to explore the relationship between the dependent variable, IC₅₀ values of Sotrovimab [64,65,66,67,68], and the independent variables, representing the structural data. This approach enabled us to investigate the interplay between structural characteristics and the functional activity of SARS-CoV-2 VOCs. Thus, we are interested in obtaining the model:

$$\begin{aligned} {\text{IC}}_{{{\text{5}}0}} \, = & \,\beta _{0} \, + \,\beta _{{\text{1}}} *{\text{corr}}\_{\text{intra}}\, + \,\beta _{{\text{2}}} *{\text{anti}}\_{\text{corr}}\_{\text{intra}}\, + \,\beta _{{\text{3}}} *{\text{corr}} - {\text{CV}}\_{\text{intra}}\, + \,\beta _{{\text{4}}} *{\text{anti}} - {\text{corr}} - {\text{CV}}\_{\text{intra}} \\ \, + \,\beta _{{\text{5}}} *{\text{corr}}\_{\text{inter}}\, + \,\beta _{{\text{6}}} *{\text{anti}} - {\text{corr}}\_{\text{inter}}\, + \,\beta _{{\text{7}}} *{\text{corr}} - {\text{CV}}\_{\text{inter}}\, + \,\beta _{{\text{8}}} *{\text{anti}} - {\text{corr}} - {\text{CV}}\_{\text{inter}} \\ \end{aligned}$$

where, β₀, β₁, β₂,…, β₈ are the regression coefficients. Assess the goodness of fit of the model by examining statistical measures such as p-values, F-statistics, adjusted r², and r². These measures provide insights into the significance and explanatory power of the independent variables.

Results

Sequence and structural variations among Spike^VOCs

Several mutations in the Spike^VOCs are related to increased binding affinity, increased infectivity, immune system escape and increased transmissibility of SARS-CoV-2 [1]. From amino acid sequences (Figure S1), we performed the multiple sequence alignment of Spike^WT, Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1 proteins (Figure S2), which shows that Spike^B.1.617.2 presents mutations in NTD (T20R, G75R, D80Y, T95I, G142D), RBD (L452R, T478K), in the SD1 subdomain (T547K, D614G, P681R), the SD2 subdomain (D614G, H655Y), and in S2 subunit (D950N and V1264L). Conversely, Spikes of BA.2, BA.5 and BQ.1 share the following mutations when compared to Spike^WT: NTD (T20I, del25/27, A28S, G339D), RBD (S371F, S373P, S375F, T376A, D405N, R408S, K417N, S477N, T478K, E484A, Q498R, N501Y, Y505H), SD1 subdomain (D614G, H655Y, N679K, P681H, N764K, D796Y). In Spike^BA.2, residues H69, V70, N440, L452, F486 are the same as the residues in the wild type, but there is the substitution Q493R only found in this Omicron sub-variant. From the multiple alignment of amino acid sequences, we built a phylogenetic tree to investigate which proteins are more closely related and share similarities. The phylogenetic tree shows that the Delta variant is more similar to the WT than to the Omicron sub-variants (Figure S3). On the other hand, among the Omicron sub-lineages, BA.5 and BQ.1 share more similarities with each other than with BA.2. Based on the analysis of the aforementioned mutations, we can observe that there are few differences between BA.5 and BQ.1. Spike^BA.5 and Spike^BQ.1 have the L452R substitution and deletions of the residues H69 and V70. Comparing the sequences of Spike^BQ.1 and Spike^BA.5, we observed few differences such as a deletion at Y145 in BQ.1 and the substitutions K444T and N460K. These three mutations are only found in Spike^BQ.1. Nevertheless, these two mutations resulted in the dominance of BQ.1 over BA.5 [2].

In order to investigate how mutations are related to structural changes of different Spike of VOCs, we compared the backbone root-mean-square deviation (RMSD) of the protomers of Spike^B.1.617.2 (from Q14 to S1146), Spike^BA.2 (from Q14 to S1144), Spike^BA.5 (from Q14 to S1142) and Spike^BQ.1 (from Q14 to S1141) (Figure S4). The average RMSD values of Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1 converged with ~ 5.0, ~ 5.1, ~ 5.0 and ~ 4.9 Å, respectively. Our previous study resulted in RMSD values of ~ 5 Å for Spike^BA.1 and ~ 4 Å for Spike^WT [21]. Therefore, Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1 visit a greater number of different structural conformations, similar to Spike^BA.1, when compared with the wild type.

Next, we calculated the backbone root-mean-square fluctuation (RMSF) of the amino acid residues of Spike^VOCs (Fig. 1). We attributed a RMSF cut-off value of 3 Å and classified the residues with high conformational flexibility as those with an average RMSF greater than 3 Å. To facilitate the location of amino acid residues in the Spike (Fig. 1a–c), we have mapped the residue numbering based on the domain architecture of the Spike^WT protein, as shown in Fig. 1b. Spike^B.1.617.2 mutations are located in NTD (T20R, G75R, D80Y, T95I, G142D), RBD (L452R, T478K), SD1 (T547K, D614G, P681R), SD2 (D614G, H655Y), on HR1 (D950N) and HR2 (V1264L) (Figure S1 and S2). These mutations influence the structural flexibility of Spike^B.1.617.2, as shown by the high RMSF value in NTD (14–23, 71–74, 113, 145–155, 164–165, 167, 212–218, 248–257), in RBD (474–488), in SD2 subdomain (622–626), in proximity of the S1/S2 cleavage site (679–690), in proximity of the S2’ cleavage site (810–812), between the FP and HR1 domains (836–846), in HR1 domain (941–942), and in HR2 domain (1140–1146) (Fig. 1d). The residues with high RMSF values in Omicron sub-variants are: (1) BA.2, in NTD (14–20, 66–76, 142–150, 161–163, 210, 243–255), in RBD (468–486), in SD2 subdomain (620–623), in proximity of the S1/S2 site (675–685), in proximity of the S2’ site (806–810), between the FP and HR1 domains (842), in HR1 domain (938, 940–941), and in HR2 domain (1138–1144) (Fig. 1e); (2) BA.5, in NTD (14–21, 67–73, 139–146, 152, 175–180, 207–210, 240–253), in RBD (469–482), in proximity of the S1/S2 site (673–683), in proximity of the S2’ site (803–808), in FP (827), between FP and HR1 domains (830–841), and in HR2 domain (1134–1142) (Fig. 1f); (3) BQ.1, in NTD (14–21, 68–71, 140–145, 177, 206–208, 241–251), in RBD (470–479), in SD2 subdomain (616–620), in proximity of the S1/S2 site (673–680), in proximity of the S2’ site (803–807), and in HR2 domain (1135–1141) (Fig. 1g).

Fig. 1

Structural flexibility of the trimeric Spike^{B.1.617.2 (Delta)}, Spike^BA.2, Spike^BA.5, and Spike^BQ.1. a S1 and S2 subunits of the Spike trimer (PDB ID 7W94 [69]). b Spike^WT amino acid sequence showing the S1 and S2 subunits and the cleavage sites S1/S2 and S2’ [1]. c RBD and NTD of the Spike. Backbone root-mean-square fluctuation (RMSF) calculated from the molecular dynamics (MD) simulations of the d Spike^B.1.617.2, e Spike^BA.2, f Spike^BA.5, and g Spike^BQ.1

In our previous study, Spike^BA.1 showed high RMSF values for residues in NTD (14–21, 70–75, 142–149, 175–179, 208–215, 243–255), in RBD (442–444 and 468–485), in SD2 subunit (599–603, 635 and 654), in proximity of the S1/S2 site (672–690), in FP (828–838), in HR1 domain (938–944), and in HR2 domain (1139–1144) [21, 70]. When compared with Spike^BA.1 [21], our results for Spike^BA.2, Spike^BA.5, and Spike^BQ.1 show that these proteins are less flexible in the S2 subunit, but more flexible in the S1 subunit (mainly in the NTD and RBD). In our previous study, we defined two regions of the Spike RBD interaction interface with human receptor Angiotensin-converting enzyme 2 (hACE2), E1 (residues K417, L455, F456, and 470–490) and E2 (residues 446–453 and 493–505), which have different physicochemical features that impact the interaction with hACE2 [11]. E1 makes mostly hydrophobic interactions, while E2 interacts mostly with hydrophilic interactions with hACE2 [11]. The Spikes of BA.2, BA.5 and BQ.1 share the same mutations in NTD (T20I, del25/27, A28S, G339D), in RBD (S371F, S373P, S375F, T376A, D405N, R408S, K417N, S477N, T478K, E484A, Q498R, N501Y, Y505H), and SD1 (D614G, H655Y, N679K, P681H, N764K, D796Y). It is known that substitutions in E1 region (K417N, S477N, T478K, E484K) and E2 region (N501Y) contribute to increase the binding affinity between RBD and hACE2 [1, 9, 31]. All results are in total agreement with several previous studies, which show high RMSF values in the residues 470–480 (in E1 region), interacting or not with hACE2 [9, 11, 21, 70]. Several studies have shown that there are high RMSF values in the E1 region of the RBD of the WT, Alpha, Beta, Gamma, Delta and Omicron variants between residues 460 and 490 [9, 11, 71,72,73]. We compare our results with the molecular dynamics study of Zhou et al. and we noted that the RMSF patterns of the WT glycosylated Spike protein and the Beta, Delta and Omicron variants [74, 75] were similar to those obtained in our study. Interestingly, even when interacting with hACE2, the glycosylated Spike protein showed the same pattern of RMSF values also when not interacting with this enzyme [70, 74,75,76,77]. Molecular dynamics studies of the Spike protein from other VOCs also showed the similar RMSF pattern, indicating that the structural data are, on average, convergent to the same result, regardless of whether or not it is glycosylated or even interacting or not with hACE2 [70, 74,75,76,77]. These results are also in line with our previous studies that show the similar pattern of RMSF values for subvariant BA.1 [21].

Overall, we did not observe significant RMSF differences in the S1 subunit of the Spike^BA.1 [21], Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1. However, we noted a significant decrease of the RMSF values in the proximity of the S1/S2 and S2’ cleavage sites of the Spike^BA.5 and Spike^BQ.1 when compared with other Spike VOCs studied in this paper. These results suggest that, for both BA.5 and BQ.1, after hACE2 interaction and S2’ cleavage, the S2 subunit of Spike tends to decrease its structural changes to expose the fusion loop, when compared with previous variants, decreasing the fusion of virus envelope and host cell membrane. Indeed, the Omicron variant has lower fusogenicity than WT and the other previous variants, indicating that the Spike protein may be losing its ability to make important structural changes in the membrane fusion process [22, 78, 79]. Park and collaborators investigated the reverse mutation in specific positions of the BA.1 variant, in order to evaluate how these mutations affect the fusogenicity [79]. In this regard, eleven reverse mutations close to the S1/S2, S2’ and HR1 region cleavage sites were analyzed [79]. Reverse mutations near the S2’ site and in the HR1 region had a stronger effect on fusogenicity, where BA.2, BA.2.12.1 and BA.5 exhibited comparable levels of fusogenicity [79]. Overall, these results combined with our data show how mutations are significantly affecting the Spike structure and, consequently, the fusogenicity of Omicron sub-variants.

In order to identify key movements within the Spike trimers, a principal component analysis (PCA) was conducted by diagonalizing the covariance matrix of Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1. The eigenvalues obtained from PCA allowed to calculate Schlitter entropy (ΔS^Schlitter) of the Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1 conformers, resulting in values of 74.8, 74.6, 75.3, and 74.3 kJ.mol⁻¹.K⁻¹, respectively. The previously reported entropy of trimeric Spike^BA.1 was 75.1 kJ.mol⁻¹.K⁻¹ [21]. Thus, the entropy values of the analyzed variants were similar but significantly higher than that of Spike^WT from a prior study (72.5 kJ.mol⁻¹.K⁻¹) [21]. This demonstrates that the Spike proteins of VOCs, in general, exhibit significant structural changes, which are quantified by ΔS^Schlitter. Next, utilizing the first eigenvector (PC1, first principal component), it was observed that the major conformational changes in Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1 primarily occurred in the NTD and RBD regions (Figs. 2 and 3). Cryo-EM studies have revealed that Omicron sub-variants exhibit reduced affinity for antibody interactions due to the mutations [6, 32,33,34, 38, 80]. Antibodies targeting the N-terminal domain (NTD) and receptor-binding domain (RBD) display limited binding to these sub-variants, particularly those specifically designed to recognize the RBD [5, 6, 21, 32, 81,82,83].

Fig. 2

MD trajectories projected in eigenvector 1. The MD trajectory projected in eigenvector 1 (also called principal component 1, PC1) represents the main movements filtered from MD trajectory. The main movements were detected in the NTD and RBD. The panels represent structure superpositions of 100 snapshots projected in the eigenvector 1 (PC1) for a Spike^B.1.617.2 and b Spike^BA.2. The backbone conformational changes were quantified by calculating entropy using Schlitter method (ΔS^Schlitter) at absolute temperature of 310 K. Spike^B.1.617.2, ΔS^Schlitter = 74.8 kJ.mol⁻¹.K⁻¹; Spike^BA.2 = 74.6 kJ.mol⁻¹.K⁻¹. The small rectangle is located at the S1 subunit while the big rectangle is located at the S2 subunit

Fig. 3

MD trajectories projected in eigenvector 1. The MD trajectory projected in eigenvector 1 (also called principal component 1, PC1) represents the main movements filtered from MD trajectory. The main movements were detected in the NTD and RBD. The panels represent structure superpositions of 100 snapshots projected in the eigenvector 1 (PC1) for a Spike^BA.5 and b Spike^BQ.1. The backbone conformational changes were quantified by calculating entropy using Schlitter method (ΔS^Schlitter) at absolute temperature of 310 K. Spike^BA.5, ΔS^Schlitter = 75.3 kJ.mol⁻¹.K⁻¹; Spike^BQ.1 = 74.3 kJ.mol⁻¹.K⁻¹. The small rectangle is located at the S1 subunit while the big rectangle is located at the S2 subunit

In general, the significant conformational changes in the RBD and NTD result in substantial decreases in antibody neutralization. One contributing factor is the reconfiguration of the RBD loop (residues 371–376), which disrupts the binding of neutralizing antibodies, ultimately hindering antibody neutralization [6]. The wide range of conformations in these regions is crucial as it may diminish recognition by neutralizing antibodies against SARS-CoV-2 [84]. Notably, high fluctuations in nearby loops 144–155 and adjacent loops 246–260, observed in different VOCs, contribute to decreased neutralizing antibody activity [85]. Structural alterations in the NTD and RBD significantly decrease binding affinity to neutralizing antibodies, facilitating immune system evasion. SARS-CoV-2 VOCs can evade neutralizing antibodies from convalescent and/or immunized individuals [1, 81, 82, 86]. Additionally, deletions 242–244 in Spike induce structural changes in the NTD, further reducing the efficiency of neutralizing antibodies [85]. Importantly, NTD and RBD contain crucial epitopes for neutralizing antibody recognition, and mutations in these domains lead to reduced recognition by vaccine-induced neutralizing antibodies [1, 5, 6, 40, 41, 87]. Numerous studies have demonstrated significant decreases in neutralization by antibodies against the Delta variant and Omicron sub-variants [1, 5, 6, 40, 41, 87]. Intriguingly, previous research has shown that reinfection risk with Omicron sub-variants is six times higher than other SARS-CoV-2 VOCs [88]. The mechanism of action of antibodies that recognize RBD is based on blocking the Spike/hACE2 interaction, preventing viral infection [89, 90]. However, although some antibodies interact with RBD, they do not prevent interaction with hACE2 but destabilize the homotrimer [91, 92]. Considering that the first principal component reveals the main movements during molecular dynamics simulation and that such movements are related to the structural changes of NTD and RBD, our results suggest a potential association between the loss of sensitivity to neutralizing antibodies and these structural changes. Our hypothesis is that the structural changes in these domains, which are the primary targets of neutralizing antibodies, could diminish the efficacy of vaccines against SARS-CoV-2. Therefore, the PCA results suggest that the pronounced conformational changes in the RBD and NTD of these variants are associated with diminished recognition by neutralizing antibodies.

Dynamical cross-correlation matrix reveals different profiles of cooperativity in the evolution of the Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1

Cooperativity is a fundamental property observed in proteins, which facilitates coordinated conformational changes in response to various stimuli, including environmental factors, ligand interactions, and interactions with other macromolecules [93,94,95,96,97,98]. This property plays a critical role in the biological and functional functions of proteins. Specifically, the Spike protein exhibits communication between different domains [21, 96, 99,100,101], and mutations in this protein impact its ability to infect host cells or evade neutralizing antibodies [21, 23, 32, 40, 41, 99]. To investigate this cooperative phenomenon, we employed Pearson’s correlation analysis to examine the movement of amino acid pairs within the Spike protein of the Delta and Omicron variants (BA.2, BA.5, and BQ.1) over time. Our objective was to verify the presence of these cooperative effects. Through our correlation analysis, we successfully identified a cooperative effect within the Spike protein. This effect arises from interactions among amino acid residues that span both short and long distances.

We calculated the dynamical cross-correlation matrix (DCCM) between all residue pairs i and j using molecular dynamics trajectories. DCCM is defined as the matrix ${DCC}_{ij} = {\langle \Delta R}_{i}.{\Delta R}_{j}$⟩, where ${\Delta R}_{i}={R}_{i}-{\langle R}_{i}\rangle$ is the difference between the Euclidean vector R_i, for position of the C_α of residue i, and the average position of this residue, ⟨R_i⟩. Next, we calculated the value for the Pearson correlation matrix after normalization as:

$${nDCC}_{ij} =\frac{{DCC}_{ij}}{{[{DCC}_{ii} .{DCC}_{jj} ]}^\frac{1}{2}}$$

where ${nDCC}_{ij}$ can range between − 1 and 1. The ${nDCC}_{ij}=0$ means no correlations, ${nDCC}_{ij}=1$ indicates complete correlation and ${nDCC}_{ij}=-1$ shows complete anticorrelation. We used DCCM for investigating cooperativity between protomers of the trimeric structure of the Spike. Due to the high number of uncorrelated residue pairs, we applied a cut-off value of 0.6 to the nDCC_ij values.

In order to investigate the distribution of local interactions and long-distance interactions in each protomer, we conducted a statistical analysis to examine changes in the average number of intra-protomer and inter-protomer correlations (above 0.6) and anti-correlations (below -0.6) (Tables S5 and S6). The correlation matrices for the Delta variant and the sub-variants BA.2, BQ.1, and BA.5 of the Omicron are displayed in panels a) of Figures S5-S8. These correlations were visualized by mapping them onto the Spike protein structure of each variant (Figs. 4, 5 and Figures S5-S8). Additionally, we performed a Student’s t-test to determine if the mutations had a significant impact on the average correlations at a significance level of 0.10. Regarding the average intra-protomer correlations above 0.6, the Delta variant (B.1.617.2) exhibited 51,473 ± 3687, which was comparable to BA.2 (62,620 ± 9462) (p-value = 0.17) and BA.5 (55,745 ± 12,324) (p-value = 0.62), but higher than BQ.1 (34,991 ± 2519) (p-value = 0.048) at a significance level of 0.10 (Table S5). BA.2 showed similarity to BA.5 (p-value = 0.47) but had a higher average than BQ.1 (p-value = 0.005), with a significance level of 0.10. Furthermore, the mean number of intra-protomer correlations in BA.5 was greater than BQ.1 (p-value = 0.095). In general, intra-protomer correlations are associated with local intermolecular interactions responsible for maintaining the structural flexibility of the trimeric Spike protein (Fig. 4). We also calculated the coefficient of variation (CV) for these correlations, resulting in values of 7.2, 15.1, 22.1, and 7.2 for Delta, BA.2, BA.5, and BQ.1, respectively (Table S5). The CV is a statistical measure that describes the relative variability of a dataset compared to its mean, providing insight into the dispersion of the data as a percentage. It is commonly used to compare dispersion among different datasets, particularly when the means differ, allowing for an assessment of data stability or consistency regardless of the measurement scale. In this context, the CV enables us to evaluate how correlations and anti-correlations are distributed within and between protomers. A low CV close to zero indicates correlations (or anti-correlations) with less variation relative to the mean, indicating an even distribution within (or between) protomers. Conversely, a high CV suggests a relatively large variation compared to the mean, indicating greater dispersion of these correlations and highlighting their uneven distribution within (or between) protomers. The Delta variant and BQ.1 exhibited the lowest CV values, while BA.2 and BA.5 had the highest, with BA.5 showing the greatest coefficient of variation. This demonstrates that Delta and BQ.1 have similar correlation distributions within each protomer, while in BA.2 and BA.5, the internal correlations are unevenly distributed. Therefore, these results indicate that, for Delta and BQ.1, the interactions responsible for maintaining the structural integrity of the Spike protein are similarly distributed in each protomer, while in BA.2 and BA.5, they are unevenly distributed.

Fig. 4

DCCM analysis for Spike^{B.1.617.2 (Delta)}, Spike^BA.2, Spike^BA.5, Spike^BQ.1. Correlations above 0.8 are depicted in the initial frame (t = 0 ns) of trimeric Spike of a B.1.617.2 (Delta), b BA.2, c BA.5, d BQ.1. These data are available in supplementary files S6, in which they can be visualized in the Pymol software

Fig. 5

DCCM analysis for Spike^{B.1.617.2 (Delta)}, Spike^BA.2, Spike^BA.5, Spike^BQ.1. Anti-correlations below -0.6 are depicted in the initial frame (t = 0 ns) of trimeric Spike of a B.1.617.2 (Delta), b BA.2, c BA.5, d BQ.1. These data are available in supplementary files S6, in which they can be visualized in the Pymol software

Table S5 also presents the average number of intra-protomer anti-correlations below − 0.6: 591 ± 586, 413 ± 405, 437 ± 602, and 58 ± 54 for B.1.617.2 (Delta), BA.2, BA.5, and BQ.1, respectively. The coefficients of variation (CV) for Delta, BA.2, BA.5, and BQ.1 were 99.2, 98.1, 137.8, and 93.1%, respectively. These high CV values, along with the standard deviation (SD), indicate a highly uneven distribution of intra-protomer anti-correlations compared to intra-protomer correlations. Figure 5 shows, in these variants, that the anti-correlations are directly linked to long-range interactions, such as the communication between domains within each S1 and S2 subunit, as well as between the S1 and S2 subunits. The results from Table S5 reveal both the unequal distribution and different average values of these anti-correlations. For instance, the Delta variant exhibits a higher average value of intra-protomer anti-correlations compared to all sub-variants of Omicron. Among the Omicron sub-variants, BA.2 has the highest value, while BQ.1 has the lowest. Considering the distribution of intra-protomer anti-correlations, as indicated by the CV values, BA.5 displays the greatest discrepancy among all analyzed VOCs. This indicates that the distributions of intra-protomer anti-correlations are uneven within each protomer. Furthermore, our results suggest a decreasing intra-protomer communication between domains in the S1 and S2 subunits from Delta to Omicron. Moreover, the communication between domains remains unequal, with BA.5 showing the highest level of inequality among the sub-variants. Although not yet fully elucidated, we hypothesize that this loss of communication between domains may affect cellular tropism during the virus’s evolutionary process [22, 78, 79] and reduce dependency on protease cleavages, such as Cathepsin L and TMPRSS2 [22, 102,103,104].

We also investigated the inter-protomer correlations and anti-correlations (Table S6) to evaluate the potential impact of SARS-CoV-2 mutations on inter-protomer communication. The analysis revealed average inter-protomer correlation values exceeding 0.6 for the different variants: Delta (1656 ± 1358), BA.2 (1922 ± 711), BA.5 (3389 ± 1611), and BQ.1 (478 ± 123). The corresponding coefficient of variation (CV) values were 82, 37, 47.5, and 25.7%, respectively. These results indicate that inter-protomer communication is highest in BA.5, which exhibits the second highest CV value among the variants. This finding suggests that the communication among protomers in this variant is significantly higher but unevenly distributed in the Spike protein. On the other hand, the Delta variant demonstrates the second highest average number of inter-protomer correlations and also has the highest CV value, indicating high but completely uneven communication between protomers. Although BQ.1 shows the lowest average number of inter-protomer correlations, the CV values are more evenly distributed compared to the other variants. These findings imply that Spike protein mutations play a role in influencing the communication network among protomers, thereby affecting the stability and overall functionality of the homotrimer structure. Conversely, BA.2 and BQ.1 exhibited lower CV values, suggesting a relatively consistent distribution of inter-protomer correlations and indicating a more stable inter-protomer communication network in these variants.

Additionally, during our analysis, we observed significant impacts of mutations on inter-protomer anti-correlations for the Delta, BA.2, BA.5, and BQ.1 variants. The data revealed the following average numbers of anti-correlations: (406 ± 300) for Delta, (1257 ± 109) for BA.2, (1950 ± 2469) for BA.5, and (33 ± 31) for BQ.1. These results indicate notable differences in the CV values of inter-protomer anti-correlations among the different variants, with values of 79.9% for Delta, 8.7% for BA.2, 126% for BA.5, and 93% for BQ.1. Interestingly, the average number of inter-protomer anti-correlations was higher in BA.5, followed by Delta, and lower in BQ.1 compared to the other variants. Both BA.5 and Delta exhibited high CV values, indicating a significant decrease in inter-protomer anti-correlations from Delta and BA.5 to BQ.1. However, in all cases, the communication remains highly dispersed. BA.2 displayed the second-highest average number of inter-protomer anti-correlations and the lowest CV value, suggesting a more similar communication pattern between protomers compared to the other variants. Overall, the mutations in these variants affect the average number of inter-protomer anti-correlations, showing different magnitudes. With the exception of the BA.2 sub-variant of Omicron, all inter-protomer communications are unequal.

Overall, these findings have important implications for understanding the structure and function of the trimeric Spike protein in these SARS-CoV-2 variants. Mutations significantly influence intra- and inter-protomer correlations and anti-correlations, thereby impacting the structural stability of the protein and its ability to interact with hACE2 and neutralizing antibodies. Moreover, these results highlight the complexity of the virus in intra- and inter-protomer domain communication, shedding light on its fundamental role in Spike protein functionality and the virus’s ability to infect host cells. Understanding these cooperative effects is of utmost importance for the development of effective therapeutic strategies against SARS-CoV-2 variants, as they provide valuable insights into the interactions and dynamics of the Spike protein. These findings can contribute to the search for more targeted treatments and the creation of more effective vaccines against different virus variants.

Cooperative behavior of D614G substitution in Spike^B.1.617.2, Spike^BA.2, Spike^BA.5 and Spike^BQ.1

According to cryo-EM studies, the D614G substitution does not cause high conformational changes in the Spike structure, except for the loss of the hydrogen bonding interaction between D614^WT and K854^WT, favoring the exposure of RBD to interact with hACE2 [45]. Our previous study has shown that D614G is predominant in the Spikes of SARS-CoV-2 VOCs [1]. This mutation also is essential for the virus because (1) it causes the disruption of the interprotomer contact, increasing up-RBD state [26]; (2) modulates structural conformations that promote membrane fusion [45]; (3) increases virion spike density and infectivity [24, 25]; and (4) increases the cleavage rate by furin [105]. Our previous molecular dynamics study of the Spike^BA.1 has shown that D614G mutation contributes to exposure RBD due to loss of hydrogen bonding interactions between _SCD614^BA.1 the residue pairs _SCK854^BA.1, _SCY837^BA.1, _SCK835^BA.1, _SCT859^BA.1 [21]. In order to understand the impact of the D614G mutation in the Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1, we investigated the relation between D614G with gain or loss of anti-correlations and correlations. In this regard, we then focused only on the correlations between residues equivalent to D614^WT (G614^B.1.617.2, G611^BA.2, G609^BA.5 and G608^BQ.1) and all other amino acid residues (Figures S9-S12). For all Spike sub-variants, intra-protomer correlations between residues equivalent to D614^WT and P681^WT and their adjacent residues correspond to local hydrogen bonding interactions. In Spike^B.1.617.2, we observed, for the protomers B and C, an intra-protomer correlation between G614 and the residues 316–319 (region between NTD and RBD) (Figure S9). In contrast, for the Spike^BA.2, both intra- and inter-protomer correlations were observed. Spike^BA.2 has correlations between G611_C and residues G32-V33_A, N162_A, and D284-V286_A, all of which are located in the NTD, the subscripts correspond to the A, B or C chains of trimeric Spike structure (Figure S10). However, Spike^BA.2 also shows intra-protomer correlations for the G611 residue of all chains and residues located between the NTD and RBD (S313-V317) and anticorrelations for residues A369 and P488 (Figure S10). In Spike^BA.5, there is an intra-protomer correlation of the G609_A with residues between NTD and RBD (F313-S320_A) and between G609_C and residues I307-Q309_C and S313_C (Figure S11). Spike^BQ.1 also showed only two intra-protomer correlations with the region between NTD and RBD, involving residue G608_A and residues N311-V314_A and residue G608_B and F312_B (Figure S12). In our previous study, Spike^WT showed few cooperative relationships, including a correlation of D614_B and the residues S71_A (in NTD) and I312_B, which is located between NTD and RBD [21]. On the other hand, the D614G mutation in Spike^BA.1 gave origin to a direct communication between G611^BA.1 and residues belonging to NTD and between NTD and RBD [21]. Overall, the D614G substitution seems to play a key role in the emergency of communication between D614G with residues in the β-sheet located between the NTD and the RBD (residues 306–334) [105]. Taken together, our results are in agreement with previous studies [105], where small structural changes in this region could be critical to generate significant cooperative effects that will interfere with the availability of RBD for the interaction with hACE2 or affect the efficiency of neutralizing antibodies.

Cooperativity supports a role of P681R mutation in high lethality the SARS-CoV-2 Delta variant

In Figure S2, we show that all Omicron sub-variants carry a P681H substitution while the Delta variant has a P681R substitution. In order to understand the effect of these substitutions in the evolution of the Delta and the Omicron sub-variants, we investigated the impact of these substitutions in the network of correlated residues. The detailed results of the DCCM analysis for the substitutions P681R in Spike^B.1.617.2 and P681H in Spike^BA.2, Spike^BA.5, and Spike^BQ.1 are shown in Figures S13-S16. In Spike^B.1.617.2, the P681R substitution causes a significant gain of correlations between distant residues, thus implying higher intra- and inter-protomer communication between S1 and S2 subunits than in the wild type (Figure S13).

In contrast, we clearly noted that the P681H substitutions of the Omicron sub-variants tend to decrease the inter-protomer correlations between the S1 and S2 subunits. A previous study showed that Spike^B.1.1.539 is significantly less fusogenic than the Spike^B.1.617.2 [78]. In addition, Suzuki and co-authors also showed that hamsters infected with SARS-CoV-2^Omicron exhibited fewer respiratory disturbances than those infected with SARS-CoV-2^B.1.617.2, suggesting that Omicron is less virulent than Delta [78]. Our results suggest that the P681R mutation may favor intra- and inter-protomer communication between S1 and S2 subunits. Therefore, given that Omicron is less fusogenic than Delta, we hypothesize that the abundant long-distance communication may be a key factor to the higher fusogenicity of the Delta variant when compared to Omicron sub-variants.

In contrast, it is already known that the Omicron variant spreads faster than Delta [106] and is more resistant to neutralizing antibodies than other SARS-CoV-2 VOCs, including Delta [22, 107,108,109,110]. Conversely, all Omicron sub-variants are highly transmissible and carry the P681H substitution (in Spike^B.1.1.539) [78]. Although still unclear, the low lethality of the Omicron sub-variants may be explained in part due to the decreased virulence of Omicron, which might be in part a consequence of the wide coverage of the global vaccination rate in humans [29]. In our computational analysis, a decrease in the cooperative behavior is strongly associated with the P681H substitution and probably affects Spike function via impairment of protomer-protomer interactions (Figures S13-S16). We propose that the low lethality of Omicron sub-variants and the associated decrease of virus fusogenicity could be explained by the loss of cooperativity caused by mutations such as the P681H substitution. In contrast, the P681R substitution present in the Delta variant is correlated with an elevated number of intra- and inter-domain correlations and, critically, this mutation contributes to enhanced communication among subdomains responsible for initiation of the hACE2 interactions and those associated with membrane fusion in a mechanism not well described. P681R, therefore, probably contributes greatly to enhancing fusogenicity, thus leading to higher infectivity and lethality.

Mutations in K417, L452, N444, and N460 make them less exposed to interact with neutralizing antibodies

One of the consequences of Spike mutations in the Omicron variant is the loss of correlations during the virus’s evolutionary process. This loss impacts the structure and function of the Spike protein, leading to a significant disruption in protomer-protomer and/or domain-domain communication. Changes in the pattern of structural communication result in structural changes in the Spike protein and, consequently, affect local intermolecular interactions. This negatively influences the viral infection process by reducing fusogenicity and altering cellular tropism. On the other hand, these changes in intermolecular and structural interactions may confer the virus with the ability to conceal regions of great importance for neutralizing antibody recognition. Indeed, SARS-CoV-2 VOCs have several mutations that confer resistance to vaccine-induced antibodies and antibodies from convalescent patients [1, 111, 112]. According to several studies, mutations at the sites K417, D420, and L455, all located in the RBD, may promote resistance to neutralizing antibodies [108, 113]. In addition, in BQ.1 and BQ.1.1, the substitutions K444T and N460K favor resistance to neutralizing antibodies [108]. Since these mutations involve substitutions of amino acid residues that may make novel hydrogen bonding interactions, we investigated how these substitutions could affect the recognition of neutralizing antibodies due to increased interactions with local amino acid residues. To verify this hypothesis, we calculated the hydrogen bonding occupations of Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1. Hereafter, occupancies are expressed as the proportion of simulation frames where the residues are close enough and in a suitable geometry for the formation of hydrogen bonding interactions.

In order to investigate the effects of substitutions K417N, L452R, K444T, and N460K in the hydrogen bonding networks, we calculated the hydrogen bonding occupancies (Tables S1-S4). In Spike^B.1.617.2, residue K417^B.1.617.2 side chain (SC) interacts with _SCN422^B.1.617.2 (occupancy = 35.96%), _SCE406^B.1.617.2 (76.32%), _SCY421^B.1.617.2 (27.40%), and main chain (MC) of _MCF374^B.1.617.2 (30.24%). In Spike^BA.2, the mutation K417N changed the interactions of the corresponding residue (_SCN414^BA.2) by adding hydrogen bonding interactions with _SCN367^BA.2 (33.41%) and _MCR451^BA.2 (62.75%). In Spike^BA.5, we observed that _SCN412^BA.5 interacts with _SCQ404^BA.5 (17.37%), _MCR449^BA.5 (66.12%), and _SCY416^BA.5 (10.72%). For Spike^BQ.1, the equivalent residue, _SCN411^BQ.1 presents increased hydrogen bonding interaction networks with residues _SCE409^BQ.1 (28.68%), _MCA375^BQ.1 (15.03%), _MCY372^BQ.1 (23.51%), _MCR457^BQ.1 (49.13%), _MCL458^BQ.1 (10.31%), and _SCY424^BQ.1 (12.36%). Overall, these results show a global increase of hydrogen bonding interactions and a more complex network of interactions between local residues. The K417N mutation has been demonstrated to be critical to evade the immune system [114]. Therefore, the K417N substitution expands the network of interactions to neighboring residues, possibly making this asparagine residue less exposed to interact with neutralizing antibodies.

We also investigated how the L452R substitution influences neutralizing antibody resistance. In Spike^B.1.617.2, we observed that _SCR452^B.1.617.2 makes hydrogen bonding interactions with residues _SCY351^B.1.617.2 (18.41%), _MCN450^B.1.617.2 (23.53%), _SCE484^B.1.617.2 (87.43%), and _SCS494^B.1.617.2 (9.43%). This network of interactions may be directly related to the loss of binding affinity with neutralizing antibodies, where _SCR452^B.1.617.2 may preferentially interact with these amino acids, weakening interactions with neutralizing antibodies. On the other hand, in Spike^BA.2, _SCL452^BA.2 residue does not mutate and thus preserved hydrophobic interactions with residues _SCY348^BA.2, _SCY448^BA.2, _SCL449^BA.2, _SCY450^BA.2, _SCC477^BA.2, _SCG479^BA.2, _SCV480^BA.2, _SCA481^BA.2, _SCY486^BA.2, _SCF487^BA.2, _SCL489^BA.2. In Spike^BA.5, _SCR447^BA.5 interacts with residues _SCY346^BA.5 (35.37%), _MCN445^BA.5 (29.54%), _SCS489^BA.5 (18.39%), and _MCS489^BA.5 (14.23%). In Spike^BQ.1, a wide network of hydrogen bonding interactions forms between _SCR445^BQ.1 and residues _SCD436^BQ.1 (86.92%), _MCS437^BQ.1 (18.56%), _SCF341^BQ.1 (35.88%), _SCN444^BQ.1 (43.39%), _SCY489^BQ.1 (20.50%), and _MCY489^BQ.1 (16.10%). L452R substitution has been reported to escape HLA-A24-restricted cellular immunity [115]. Therefore, the Spike^BQ.1 contains an increase in hydrogen bonding interactions due to the L452R substitution. Thus, our results show that the mutation L452R expands the hydrogen bonding interactions with neighboring residues, a change that may likely reduce the availability of this residue for interactions with neutralizing antibodies against SARS-CoV-2.

BQ.1 and BQ.1.1 also have two other substitutions, K444T and N460K, that are known to favor resistance to neutralizing antibodies [108]. We measured the hydrogen bonding interaction occupancies for residues _SCK439^BA.5 and _SCT438^BQ.1 with neighboring residues. Residue _SCK439^BA.5 makes hydrogen bonding interactions with _MCN443^BA.5 (11.29%), _MCG442^BA.5 (22.57%), and _SCN443^BA.5 (10.56%). However, _SCT438^BQ.1 interacts with _MCS432^BQ.1 (14.12%), _MCD436^BQ.1 (12.49%), _SCN442^BQ.1 (43.39%), and _SCR503^BQ.1 (53.53%). Thus, the N444T mutation (_SCT438^BQ.1) probably also leads to an escape of the immune system recognition, as suggested from a significant increase in hydrogen bonding occupancy for this threonine residue. Most notably, a large change in occupancy for the hydrogen bonding interaction between _SCT438^BQ.1 and _SCN442^BQ.1 and the emergency of interactions with _SCR503^BQ.1 highlights the differences between BA.5 and BQ.1. In addition, in Spike^B.1.617.2, hydrogen bonding interactions are observed between _SCN460^B.1.617.2 and the residues _SCD420^B.1.617.2 (80.47%), _SCK424^B.1.617.2 (62.39%), and _SCY421^B.1.617.2 (39.35%). For Spike^BA.2, there are interactions between _SCN457^BA.2 and _SCD417^BA.2 (74.96%), and _SCK421^BA.2 (63.93%). Spike^BA.5 presents hydrogen bonding interactions between _SCN455^BA.5 and residues _SCK419^BA.5 (38.71%), _SCD415^BA.5 (57.35%) and _SCY416^BA.5 (10.11%). For Spike^BQ.1, interactions are observed between _SCK454^BQ.1 and _SCD979^BQ.1 (35.07%), _SCD414^BQ.1 (94.89%) and _SCY415^BQ.1 (38.01%).

Our results suggest SARS-CoV-2 Omicron sub-variants are evolving through the acquisition of evolutionary advantageous mutations that enhance interactions of T438^BQ.1 and residues N442^BQ.1 and R503^BQ.1, thus decreasing the availability of these amino acids to neutralizing antibodies. Another fundamental contribution is the N460K substitution observed in Spike^BQ.1, which mediates a strong saline interaction between _SCK454^BQ.1 and _SCD414^BQ.1 and thus may also contribute to a low availability of the _SCK454^BQ.1 for interactions with neutralizing antibodies. Indeed, Qu and co-authors showed that K444T and N460K mutations in Spike^BQ.1 and Spike^BQ.1.1 promote resistance to neutralizing antibodies [83]. In addition, previous studies have reported that people previously infected with Omicron can be reinfected with different Omicron sub-variants in a shorter time than SARS-CoV-2 VOCs. In this case, the reinfection can happen after two months [88]. Thus, our results are important because they show how the Spike protein evolves into variants that hide the recognition region of neutralizing antibodies, increasing immune system evasion and consequently transmission of SARS-CoV-2.

Unveiling structural relationships: evasion of neutralizing antibodies revealed through multiple linear regression

One of the key considerations is how to relate simulation data from different variants to the evasion of neutralizing antibodies. One approach is to utilize the data on average number of correlations and anti-correlations, as well as coefficient of variation (CV), within-protomer and inter-protomer, as independent variables. We used data from our previous studies of WT and BA.1 (Tables S7 and S8) and complemented them in Table S9. In the multiple linear regression model, we calculated the coefficients of multiple linear regression. Among all the combinations of independent variables used to generate the models, we selected five models that exhibited F-statistic values greater than 20 (Table S10). As these five models showed satisfactory values of r² and adjusted r², we used a significance level of p-value < 0.05 to evaluate each of these models. This is because the p-value indicates the strength of evidence against the null hypothesis. If the p-value is smaller than the significance level of 0.05, we consider that there is sufficient statistical evidence to reject the null hypothesis and conclude that there is a significant relationship between the independent variables and the dependent variable. Based on this criterion, we selected Model 2, in Table S10, as the most appropriate one.

$${\text{IC}}_{{{5}0}} = { 1}00{49}.{96 } + { 3}.{4}0*{\text{corr}}\_{\text{inter }} - { 7}.0{1}*{\text{anti}}{-} {\text{corr}}\_{\text{inter }} - { 147}.{58}*{\text{corr}}{-}{\text{CV}}\_{\text{inter}}$$

$$N = { 6};{\text{ F}}{-}{\text{Statistic }} = { 56}.0{5};r^{2} \, = \, 0.{99};{\text{ adjusted }}r^{2} \, = \, 0.{97};{\text{ p}}{-}{\text{value }} = \, 0.0{12 } < \, 0.0{5}.$$

The estimated regression coefficients in the selected model provide valuable insights into the relationship between the independent variables and the dependent variable. In this study, we investigated the impact of inter-protomer communication on the resistance of neutralizing antibodies against SARS-CoV-2, as represented by the IC₅₀ value of Sotrovimab, a measure of the inhibitory concentration required to reduce viral replication by 50%. We found that the average number of inter-protomer correlations, the average number of inter-protomer anti-correlations, and the CV of inter-protomer correlations play crucial roles in determining the virus’s resistance to neutralizing antibodies. The regression coefficient for the ‘corr_inter’ variable was estimated at 3.40, indicating that a one-unit increase in the average number of inter-protomer correlations is associated with a 3.40 ng/mL increase in the IC₅₀ value of Sotrovimab, while holding other variables constant. Conversely, the regression coefficient for the ‘anti-corr_inter’ variable was estimated at − 7.01, suggesting that a one-unit increase in the average number of inter-protomer anti-correlations is associated with a 7.01 ng/mL decrease in the IC₅₀ value, assuming other variables remain constant. Furthermore, the regression coefficient for ‘corr-CV_inter’ was estimated at − 147.58, indicating that a one-unit increase in the coefficient of variation of inter-protomer correlations results in a 147.58 ng/mL decrease in the IC₅₀ value, holding other variables constant. These findings highlight the importance of inter-protomer communication in immune evasion and resistance to neutralizing antibodies. In particular, we observed a significant decrease in inter-protomer anti-correlations in the BQ.1 strain, which exhibits higher antibody resistance. This discovery suggests that inter-protomer communication plays a critical role in the virus’s ability to evade the immune response.

In summary, the estimated coefficients in the multiple regression model provide compelling evidence that the average number of inter-protomer correlations and anti-correlations, along with the CV of inter-protomer correlations, are key factors in neutralizing antibody resistance. These findings may contribute to the development of effective strategies to combat SARS-CoV-2 and enhance our understanding of the intricate relationship between virus structure and function throughout the evolutionary process.

Discussion

The virulence of SARS-CoV-2 is primarily driven by genomic mutations that disrupt crucial processes such as hACE2 recognition, protease cleavage, and fusogenic activity. In our study, we performed phylogenetic analysis of Spike^WT, Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1, revealing distinct groups with Spike^B.1.617.2 being more closely related to Spike^WT than the Omicron sub-variants (Figure S3). To explore the relationship between structural cooperativity of these mutants and their biological outcomes, we conducted 300 ns simulations of their trimeric ectodomains, employing RMSD, RMSF, PCA, DCCM, and a multiple linear regression model. Notably, these variants exhibited greater conformational changes compared to the wild-type Spike protein (Spike^WT). The mutations affected the structural flexibility of the protein, particularly in the N-terminal domain (NTD) and receptor-binding domain (RBD). These structural alterations have the potential to impact the Spike protein’s ability to interact with the cellular receptor (hACE2) and neutralizing antibodies. Principal component analysis (PCA) further revealed that the most significant conformational changes in the variants occurred within the NTD and RBD regions, which are known targets of neutralizing antibodies. Such structural changes may contribute to reduced vaccine efficacy and decreased sensitivity of neutralizing antibodies to these variants.

In addition, correlation analysis allowed us to investigate the dynamic behavior of amino acid pairs within the Spike protein of Delta variant and Omicron sub-variants over time, revealing a cooperative effect within the protein. We observed the distribution of local and long-range interactions within each protomer, highlighting differences in the distribution of correlations and anti-correlations. The mutations in these variants significantly affected intra- and inter-protomer correlations and anti-correlations, thereby influencing the structural stability of the protein and its ability to interact with hACE2 and neutralizing antibodies. Collectively, our findings underscore the complexity of the virus and its ability to infect host cells and evade the immune system through intricate intra- and inter-domain communication within the Spike protein.

We hypothesize that SARS-CoV-2 initially evolved to enhance interaction with hACE2, thereby increasing the infection rate (as seen in Spike^B.1.617.2). However, the Omicron sub-variants appear to have emerged to evade the immune system. A previous study demonstrated that the sub-variants BQ.1, BQ.1.1, BA.4.6, BF.7, and BA.2.75.2 are resistant to neutralizing antibodies from individuals who received triple vaccination [83]. However, all Omicron sub-variants, particularly BA.2.75.2, retained their reduced infectivity in the human lung cancer cell line Calu-03 [83]. Additionally, transmissibility and immune system evasion play significant roles in determining the severity of infections, hospitalizations, and fatalities. Among the variants of concern (VOCs), Omicron has been identified as having reduced virulence compared to the wild-type (WT) strain and previous variants [22, 116,117,118]. The prevalence of BA.5 and BQ.1 globally, including B.1.617.2 and BA.2, motivated our investigation into how mutations impact in the structure of the Spike and, consequently, in transmissibility and antibody escape of these SARS-CoV-2 VOCs.

Mutations may give rise to new variants of the coronavirus, potentially triggering new pandemics. It is crucial to investigate the impact of these mutations on the Spike protein to understand the virus’s evolution and devise strategies to mitigate their exacerbating effects. An effective approach in this regard would involve developing vaccines that specifically target each variant. Generally, novel coronavirus variants exhibit higher mutation rates compared to the wild-type (WT) form, particularly in the NTD and RBD regions of the Spike protein [34,35,36,37,38]. Additionally, mutations in the S2 subunit of Spike are also observed, which can confer both evolutionary advantages and disadvantages to SARS-CoV-2. Notably, the D614G mutation (present in all VOCs), P681R (found in the Delta variant), and P681H (present in the Omicron variant) are significant. We observed that the cooperativity between the three protomers within the trimeric Spike structure is gained in BA.2 but predominantly lost in the Omicron sub-variant BQ.1. We hypothesized that P681H also may exhibit a strong correlation with low lethality but P681R may be associated with increased lethality and enhanced fusogenicity. Therefore, our results support the hypothesis that cooperativity, including mutations D614G, P681R, and P681H, impact both the infectivity and lethality of SARS-CoV-2.

Indeed, D614G mutation in the Spike protein provides several evolutionary advantages to the virus, including enhanced entry into host cells, improved structural flexibility, increased hACE2-binding capacity, and a higher density of viral Spikes [25, 43,44,45,46,47]. P681R mutation facilitates the cleavage of the Spike protein and enhances the fusogenicity of SARS-CoV-2, leading to greater pathogenicity compared to the wild-type [48,49,50]. The P681H mutation, found in both the Alpha and Omicron variants, results in reduced fusogenicity compared to the P681R mutation [48] and exhibits similarities to the wild-type in certain cell lines [49]. In the Omicron variant and a pseudovirus carrying the P681H mutation, there is an observed independence from TMPRSS2 in viral infection, suggesting a potential change in cell tropism and a decrease in virus infectivity [48].

Our computational molecular dynamics simulations revealed significant structural changes in the S1 subunit of Spike^B.1.617.2, Spike^BA.2, Spike^BA.5, and Spike^BQ.1, particularly in the NTD and RBD regions. We also propose that substitutions K417N, L452R, N444T, and N460K enhance resistance to neutralizing antibodies by strengthening their hydrogen bonding interactions with adjacent residues, thereby obstructing interactions with neutralizing antibodies. These structural alterations provide a basis for understanding the loss of sensitivity to neutralizing antibodies in VOCs. Furthermore, our results suggest potential mechanisms by which certain mutations contribute to immune system evasion in the Delta and Omicron variants, thereby informing the development of updated vaccines against novel Omicron sub-lineages. In this regard, our study using DCCM, along with the multiple linear regression model and hydrogen bond analysis, suggests that the mechanism by which the virus escapes the immune system is associated with the loss of communication between protomers and reduced exposure of potential residues that could interact with neutralizing antibodies. Indeed, when compared to other VOCs [80], the Omicron sub-variants exhibit an increased ability to evade the immune system. The mutations observed in Spike^BQ.1, as analyzed in this study, seem to favor loss of communication between protomers and domains as a means of immune evasion. Remarkably, these mutations do not significantly affect the interaction with hACE2, as reported in published studies [33].

In summary, our study revealed that correlations in the Spike protein of each variant are primarily influenced by local hydrogen interactions and coordinated movements of atoms resulting from these interactions. Conversely, anti-correlations are more associated with long-range interactions. Further analysis of the anti-correlations showed that the furin cleavage site region in Spike^{B.1.617.2 (Delta)} is more closely associated with the S2 subunit than the S1 subunit, indicating that cleavage in this region may facilitate membrane fusion during viral infection. In Spike^BA.2 and Spike^BA.5 variants, we observed stronger anti-correlations between structural changes in the RBD and NTD of one protomer with the RBD and NTD of the neighboring protomer, particularly within the S1 subunit. Notably, mutations in Spike^BQ.1 result in a significant decrease in anti-correlations, impairing communication between the S1 and S2 subunits, as well as between the corresponding domains of these subunits. This loss of correlations may hinder the viral infection process, aligning with previous studies that reported lower fusogenic capacity of Omicron compared to the original WT strain and the Delta variant [22, 78, 79]. Our hypothesis is that this loss of cellular invasion capacity may be attributed to decreased communication between the S1 and S2 subunits. However, it is plausible that the virus is developing immune system escape mechanisms, which could compensate for this evolutionary disadvantage. These findings underscore the importance of inter-domain communication within the Spike protein and provide valuable insights into potential mechanisms driving virus evolution. The observed disparities in inter-protomer communication underscore the complexity and adaptability of SARS-CoV-2. The virus’s remarkable ability to evolve rapidly and adapt to changing environments is crucial for its survival and transmission. Understanding the intricate interplay between mutations, intra- and inter-protomer interactions, and their effects on viral fitness is paramount for comprehending viral dynamics and devising effective therapeutic strategies.