A computational evaluation of structural stability of omicron and delta mutations of SARS-CoV-2 spike proteins and human ACE-2 interactions.

Several more infectious SARS-CoV-2 variants have emerged globally since SARS-CoV-2 pandemic and the discovery of the first D614G variant of SARS-CoV-2 spike proteins in 2020. Delta (B.1.617.2) and Omicron (B.1.1.529) variants have proven to be of major concern out of all the reported variants, considering their influence on the virus' transmissibility and severity. This study aimed at evaluating the impact of mutations on these two variants on stability and molecular interactions between the viral Spike protein and human angiotensin converting enzyme-2 (hACE-2). The spike proteins receptor binding domain (RBD) was docked with the hACE-2 using HADDOCK servers. To understand and establish the effects of the mutations on the structural stability and flexibility of the RBD-hACE-2 complex, molecular dynamic (MD) simulation of the docked complex was performed and evaluated. The findings from both molecular docking analysis and binding free energy showed that the Omicron (OM) variant has high receptiveness towards hACE-2 versus Delta variant (DT), thereby, responsible for its increase in transmission. The structural stability and flexibility evaluation of variants' systems showed that mutations on DT and OM variants disturbed the stability of either the spike protein or the RBD-hACE-2 complex, with DT variant having greater instability impact. This study, therefore, assumed this obvious instability observed in DT variant might be associated or responsible for the reported severity in DT variant disease over the OM variant disease. This study provides molecular insight into the effects of OM and DT variants on stability and interactions between SARS-CoV-2 protein and hACE-2.

Introduction

As of August 24, 2022, about 598, 180, 048 cases of SARS-CoV-2 infections have been recorded (https://coronavirus.jhu.edu/map.html). This huge spread of the virus could be associated with emerging changes/variation in the virus genetic content because of mutation, resulting to emergence of several variants of the virus. Study has shown that the virus RNA, rapidly mutate in host cells [1], thereby produces the diverse variants of the virus. Of the different classifications of SARS-CoV-2 variants by WHO, variants of concern (VOC) have been the most infectious and virulent variants [2]. WHO defined VOC as variants with a) modification in their genetic code that alters SARS-CoV-2 characteristics, b) variants that cause notable viral transmission, and severe COVID-19 epidemiology, and c) reduced drugs/vaccines efficiency [2].

Presently, five SARS-CoV-2 variants fall under the VOCs, namely, Alpha (lineage B.1.1.7) (Domingo & Benito, 2021), Beta (lineage B.1.351) (Reincke et al., 2022), Gamma (lineage P.1) (Nonaka et al., 2021), Delta (lineage B.1.617.2) [3], and Omicron (lineage B.1.1.529) [4] VOCs. Out of these, Delta (B.1.617.2) and Omicron (B.1.1.529) variants have proven to be of major concern, considering their influence on the virus’ transmissibility and severity [5]. Delta variant (lineage B.1.617.2) was discovered in countries within few months and WHO labeled the B.1.617.2 sublineage on May 11, 2021, as a VOC delta variant [6]. By the end of third quarter of 2021, delta variant has become the global dominant VOC [2,4,7]. The predominant mutation in the spike protein of the B.1.617 sublineage includes D614G, L452R, and P681R (Fig. 1 and Table 1 ). A study by Mlcochova et al. showed that the variant is eight-fold less susceptible to antibodies produced by vaccines and six times less sensitive to serum antibodies from recuperating persons when compared to its wild type [3].

Fig. 1

A comparison of Delta and Omicron variant spike mutation [8]. OM spike protein's variation is determined by 30 mutations (3 small deletions and 1 insertion).

Table 1

Comparism of Spike protein mutations present in Delta and Omicron variants [8].

	Variants
	Wild type (WT)	Delta Variant (B.1.617.2)	Omicron (B.1.1.529)
Sequence	NCBI ID:P0DTC2	NCBI: QWK65230.1	GSAID ID: R40B60_BHP_3321001247/202
Amino acid/position of Mutation	–	T19R, G142D, Δ156–157, R158G, Δ213–214, L452R, T478K, D614G, P681R, D950 N	A67V, Δ69–70, T95I, G142D, Δ143–145, N211I, L212V, ins213-214RE, V215P, R216E, G339D, S371L, S373P, S375F, K417 N, N440K, G446S, S477 N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681H, N764K, D796Y, N856K, Q954H, N969K, L981F

Omicron variant (lineage B.1.1.529) is another VOC. It was detected by the Network for Genomic Surveillance in South Africa on November 24, 2021 [9] and was named by the WHO (WHO, 2021d). The three sublineages of Omicron include B.1.1.529.1 or (BA.1), B.1.1.529.2 or (BA.2), and B.1.1.529.3 or (BA.3) (Lineage B.1.1.529.2021). Within few months, there has been rapid spread of the viral variant in more than 35 countries around the world [9]. When compared to other VOCs, Omicron's genome was discovered to be a distinctive and predominant variant, thereby causing approximately 55%–60% cases of new SARS-CoV-2 infection cases. It was uncovered to have six distinct mutations which differ from the wild-type (2019 Wuhan variant) resulting to its high transmissibility [[10], [11], [12]] and possible evasion of antibodies.

More than 55 mutations were observed in the genome of the omicron variant. Thirty-one (31) are situated in the S1 domain of the spike glycoprotein [2,13]. 15 of the mutations in the variants situated in the spike protein hACE2 RBD as shown in Fig. 1 and Table 1 [[14], [15], [16]].

This study is focused on the investigation and evaluation of the impact of mutations on these two variants of concerns (Delta and Omicron variant) on the structural stability and flexibility of the RBD-hACE-2, binding affinity and molecular interactions between the viral Spike protein and hACE-2, employing computational techniques. To understand SARS-CoV-2 viral infectivity and pathogenesis, it vital to have better knowledge of the virus' mechanism of recognizing its host receptor (hACE-2). Molecular docking analysis is a widely held and less costly approach to gain understanding of this interaction between SARS-COV-2 and it host receptor, hACE-2. This approach has been employed in several studies that involves computational investigation of either protein-protein or protein-ligands integrations [[8], [17], [18], [19], [20]].

Methods

Proteins acquisition, preparation and molecular docking

The X-ray crystal structures of the complex of SARS-CoV-2 spike glycoprotein with hACE-2 (PDB ID: 7DF4) was obtained from the RSCB Protein Data Bank [21]. The SARS-CoV-2 spike glycoprotein and hACE-2 structures were separated, and the respective mutations for each variant were inserted on the spike protein using rotamer library [22] plugin on the UCSF Chimera software package [23]. The structures of the proteins were prepared by removing water molecules, nonstandard naming, protein residue connectivity. Molecular docking of the structures of Spike RBD for each variant with hACE-2 was done using Haddock molecular docking server [24], with default docking parameters. The docked complexes were then subjected to molecular dynamics simulations.

Molecular dynamic (MD) simulations

The MD simulation was performed as described by Idowu et al. with little modification [19]. The simulations were performed using the GPU version provided with the AMBER package (AMBER 18) [25], in which the FF18SB variant of the AMBER force field [26] was used to describe the systems.

The Leap module of AMBER 18 allowed for the addition of hydrogen atoms and Cl− and Na + counter ions to complexes, to neutralize all systems. The systems were then suspended implicitly within an orthorhombic box of TIP3P water molecules such that all atoms were within 10 Å of any box edge [27]. An initial minimization of 2000 steps were carried out with an applied restraint potential of 500 kcal/mol for both solutes. They were performed for 1000 steps using the steepest descent method followed by 1000 steps of conjugate gradients. An additional full minimization of 1000 steps were further carried out using the conjugate gradient algorithm without restraint. A gradual heating MD simulation from 0 K to 300 K was executed for 50 ps, such that the systems maintained a fixed number of atoms and fixed volume. The systems' solutes were imposed with a potential harmonic restraint of 10 kcal/mol and collision frequency of 1.0 ps. Following heating, an equilibration estimating 500 ps of each system was conducted; the operating temperature was kept constant at 300 K. Additional features such as several atoms and pressure were also held constant, mimicking an isobaric-isothermal ensemble. The system's pressure was maintained at 1 bar using the Berendsen barostat [28,29].

The total time for the MD simulations conducted was 50 ns. In each simulation, the SHAKE algorithm was employed to constrict hydrogen atoms' bonds [30]. The step size of each simulation was 2fs, and an SPFP precision model was used. The simulations coincided with the isobaric-isothermal ensemble (NPT), with randomized seeding, the constant pressure of 1 bar maintained by the Berendsen barostat [29], a pressure-coupling constant of 2 ps, a temperature of 300 K and Langevin thermostat [31] with a collision frequency of 1.0 ps.

Post-dynamic analysis

Analysis of Root Means Square Deviation (RMSD), Radius of Gyration (RoG), and Solvent Accessible Surface Area (SASA) was done using the CPPTRAJ module employed in the AMBER 18 suit. All raw data plots were generated using python 3.9 on Anaconda3 software.

Binding free energy calculations

To estimate and compare the systems' binding affinity, the free binding energy was calculated using the Molecular Mechanics/GB Surface Area method (MM/GBSA) [32]. Binding free energy was averaged over 50000 snapshots extracted from the 50ns trajectory. The free binding energy (ΔG) computed by this method for each molecular species (complex, ligand, and receptor) can be represented as: Egas denotes the gas-phase energy, which consists of the internal energy Eint, Coulomb energy Eele and the van der Waals energies Evdw. The Egas was directly estimated from the FF14SB force field terms. Solvation free energy, Gsol, was estimated from the energy contribution from the polar states, GGB, and non-polar states, G. The non-polar solvation energy, SA. GSA was determined from the solvent-accessible surface area (SASA), using a water probe radius of 1.4 Å. In contrast, the polar solvation, GGB, the contribution was estimated by solving the GB equation. S and T denote the total entropy of the solute and temperature, respectively.

Results and discussion

In this study, the result of molecular docking analysis of the two variants and wild type of SARS-CoV-2 with hACE-2 showed that the OM variant has highest docking score of −134.4 kcal/mol, while DT variant showed a lower score of −59.2 kcal/mol (Fig. 2 ). A molecular docking score is therefore a measurement of the fitness of a molecule/ligand into the catalytic active/binding site pocket of an enzyme or protein, and the more negative the value, the better the fitness of the molecule [33]; Shode et al., 2021). This scoring functions allowed the estimation and prediction of the binding affinities of individual molecules [34]. This result corroborates the finding of Kumar et al. that earlier reported OM exhibited highest docking sore of −539.81 kcal/mol, higher than the docking score reported for the DT variant (−529.62 kcal/mol) [8]. Furthermore, the trend in our result is similar to the trend reported in the finding of Hwang et al. that reported docking energy of 340.0 kcal/mol for the WT, 25.2 kcal/mol for the DT, and −382.1 kcal/mol for the OM variant [35]. This indicates that the OM variant is more receptive to the host receptor, hACE2 than the DT. This result might be the first clue to explain the higher infectivity/transmission of SARS-CoV-2 observed in OM variant over the DT variant.

Fig. 2

Molecular Docking scores SARS-CoV-2 Spike glycoprotein Variants a) Wild-type (WT), b) Delta-type (DT) and c) Omicron type (OM) toward hACE-2.

Prediction of the binding energy between variant spike RBD and hACE-2

Molecular dynamics is a popular computational method for calculating the binding affinity and probing energy interaction. In this study, molecular mechanics/generalized born surface area (MMGBSA) computational technique was employed to estimate the binding free energies (ΔGbind) between the variants and wild type spike proteins toward hACE-2. The molecular dynamic simulation of the RBD−hACE-2 docking complexes was used for the predictions of binding affinity between the hACE2 receptor and two variant and WT. Table 2 showed the thermodynamic binding energy profiles for the variants toward hACE-2. The estimated ΔGbind result showed the binding free energy for RBD and hACE-2 is stronger in OM variant (−43.591 kcal/mol) than the DT variant (−19.743 kcal/mol) and WT- (−29.700 kcal/mol). Our finding is similar to the reports of previous studies that showed the binding affinity between hACE-2 and the RBD of the OM variant to be stronger than that of the WT and DT variant [[35], [36], [37]]. For instance, Hwang et al. reported the binding energy was−894.4 kcal/mol for the WT, −980.8 kcal/mol for the DT variant, and −1444.5 kcal/mol for the OM variant [35]. The result of our binding energy together with the molecular docking result, further suggests that the OM variant is more receptive to hACE2 than the DT, and might be associate or responsible for the greater transmission of the virus observed in OM variant than DT variant and WT.

Table 2

Thermodynamic Binding Free Energy Profiles for Sgp Variants toward hACE-2.

Energies	WT	DT	OM
Δ EvdW	−98.725 ± 6.270	−57.436 ± 19.141	−88.452 ± 16.263
ΔEelec	−951.678 ± 88.028	−1047.226 ± 63.408	−1778.178 ± 67.000
ΔGgas	−1050.403 ± 86.036	−1104.662 ± 76.974	−1866.63 ± 74.596
ΔGsolv	1003.848 ± 87.233	1093.919 ± 65.055	1836.167 ± 70.992
ΔGbind	−29.700 ± 6.383	−19.743 ± 5.187	−43.591 ± 9.938

ΔEele: electrostatic energy, ΔEvdW: van der Waals energy, ΔGsol: solvation free energy, ΔEgas: gas-phase free energy, ΔGbind: total binding free energy.

Protein-protein molecular interactions

A typical receptor-ligand interaction examined the molecular interactions between the bound ligand and the amino acid residues at the binding sites of the protein [38,39]; Obakachi et al., 2022). However, in this study we evaluate the protein-protein molecular interactions between the interacting amino acid residues of the spike RBD and the hACE-2. Fig. 3 showed the interaction poses of the RBD-hACE-2 complexes for the variants and WT after 50 ns MD simulation. Table 3 summarized the number and type of interactions that exist within each complex. For the WT complex, a total of 21 interactions (17 hydrogen bonds and 5 hydrophobic bonds) was observed. However, the presence of mutations on the DT spike protein lowers the number of interactions in the DT spike protein-hACE-2 complex to a total of 16 bonds (9 hydrogen and 7 hydrophobic bonds). Furthermore, the DT mutations significantly lower the hydrogen bond in the complex from 17 in WT to 9 in DT. The reduction in interaction, undoubtedly, is responsible for the lowered docking score and binding affinity reported for DT in this study. For the OM variants, insignificant reduction in the number of interactions, 19 (10 hydrogen and 9 hydrophobic bonds) was observed.

Fig. 3

Interacting Amino residues between the Variants' Spike proteins and hACE-2.

Table 3

Interacting Amio residues between the Variants Spike protein and hACE-2.

WT	DT	OM
Phe1070-Met64 (Pi-alkyl), Gln6-Ser1061 (Hb), Gly1060-Gln6 (Hb), Tyr65-Ala1059 (Hb), Lys13-Tyr1073 (Pi-sigma), Lys13-Gln1077 (Hb), Lys13-Gln1077 (Hb), Gln1077-Glu17 (Hb), Lys1001-Asp12 (Hb, Salt Bridge), Hie16-Tyr1037 (Pi-pi T-shaped), Tyr1037-Hie16 (Hb), Hie16-Ser1078 (Hb), Arg987-Hie16 (Hb), Tyr1033-Asp20 (Hb), Lys335-Tyr1079 (Hb), Tyr1089-Glu19 (Hb), Tyr1089-Lys335 (Pi-akyl), Tyr23-Gln1082 (Hb), Asn1085-Lys335 (Hb), Gly1086-Lys335 (Hb), Hie16-Ser1078 (Hb, 2.79)	Thr306-Ser1061 (Hb), Asn1071-Ala366 (Hb), Tyr367-Asn1071 (Hb), Ala369-Asn1071 (Hb), Phe338-Ala1059 (Pi-alkyl), Ala368-Ala1059 (Akyl), Gly336-Ala1059 (Hb), Tyr1057-Lys335 (Hb), Phe1040-Lys335 (Pi-alkyl), Lys335-Leu1039 (Hb), Lys335-Leu1039 (Hb), Gln1077-Hie16 (Hb), Lys1001-Asp20 (Hb, Electrostatic), Arg987-Glu17 (Hb; Electrostatic), Lys50-Arg992 (Unfavorable Positive-positive), Arg987-GLU17 (Hb; electrostatic)	Lys1062-Lys335 (Hb), Tyr1073-Ala369 (Hb), Hie16-Leu1039 (pi-alkyl), Arg1077-Glu17 (Hb; electrostatic), Arg1077-Glu17 (Hb; electrostatic), Arg1077-Asp12 (Hb), Arg1077-Asp12 (Hb), Arg1077-Asp12 (Hb), Arg1077-Lys13 (Hb), Lys13-Glu990 (Hb; Electrostatic), Lys13-Glu990 (Hb; Electrostatic), Ser1-Asn1032 (Hb), Ser1-Asn1032 (Hb), Ser1030-Gln6 (Hb), Ser1030-Gln6 (Hb), Arg1082-Tyr65 (Hb; Electrostatic), Tyr1085-Leu61 (Pi-alkyl), Hie1089-Leu61 (Pi-alkyl), Met64-Tyr1085 (Pi-sulfur)
Total = 21 (17 hb, 5 Hydrophobic Bonds)	Total = 16 (9 hb, 7 Hydrophobic Bonds)	Total 19 (10 Hb, 9 Hydrophobic Bonds)

However, a significant increase in strong hydrophobic interaction/bond was recorded compared to the WT (5 hydrophobic bond). This study suggests that the increase in the number of strong hydrophobic bonds in OM variant might be the reason for high binding energy and docking score reported in this study and other related studies.

Structural stability and flexibility evaluation of variants’ systems

The structural stability of the protein complexes was measured following experimental simulation of the spike protein RBD (for both variants and WT) together with the hACE-2. To confirm the dynamic of the complexes, 50 ns MD simulations was performed. Binding of a molecule (either drugs or protein) to a specific biological target is usually associated with structural and conformational changes, which in most cases influence the biological activity of that target [[40], [41], [42]]. To establish the stability and accurate equilibration of the investigated complexes, RMSD, RoG, RMSF and SASA of alpha carbon (Ca) atoms were monitored and analyzed along with the entire duration of 50 ns of the MD simulation for the WT, DT and OM. These parameters were calculated for the RBD site only (RBD), the whole spike glycoprotein only (SPIKE) and the spike RBD-hACE-2 (RBD-ACE-2) complex.

The measurement of the systems convergence and stability is referred to as RMSD [43]. Fig. 4, Fig. 5, Fig. 6 showed the result of the RMSD plot that measures the complexes’ convergence and stability for the RBD, SPIKE and RBD-hACE-2 complex, respectively. Table 4 showed the average values of each parameter used to interpret structural stability of each system. For the RBD (Fig. 4a), the RMSD plot, after maintaining convergence at approximately 8 ns, all the complexes exhibited favorable stability throughout the MD simulations, except for the RBDOM with raised RMSD plot at 28 ns. The average RMSD value for RBDOM (2.95 Å) is higher that the estimated values for RBDDT and RBDWT with average values of 2.54 Å and 2.47 Å, respectively (Table 4). The mutations on the OM variant increase the RMSD plot, thereby altering its stability. However, for the whole spike systems, after all the systems maintained convergence at approximately 5 ns, SPIKEOM and SPIKEDT systems exhibited unfavorable stability throughout the MD simulations (Fig. 5a). With the DT variant showing the highest average RMSD value of 8.95 Å, and SPIKEOM exhibited the lowest value of 7.28 Å. A similar trend was observed in the RBD-

Fig. 4

Comparative a). RMSD b). RoG c). RMSF, and d). SASA profile plots of C-a atoms of RBD calculated throughout 50 ns molecular dynamics simulation.

Fig. 5

Comparative a). RMSD b). RoG c). RMSF, and d). SASA profile plots of C-a atoms of Whole Spike proteins calculated throughout 50 ns molecular dynamics simulation.

Fig. 6

Comparative a). RMSD b). RoG c). RMSF, and d). SASA profile plots of C-a atoms of Spike-hACE-2 complex calculated throughout 50 ns molecular dynamics simulation.

Table 4

Calculated Average Values of Parameter used to Interpret Structural stability and Flexibility of Variants’ Systems.

Parameters	RBD			Spike Protein			Spike-hACE-2
	WT	DT	OM	WT	DT	OM	WT	DT	OM
RMSD	2.47 ± 0.30	2.54 ± 0.20	2.95 ± 0.35	8.75 ± 2.30	8.95 ± 3.46	7.28 ± 1.26	9.23 ± 2.60	10.58 ± 4.81	7.81 ± 1.88
RMSF	3.85 ± 0.77	4.82 ± 0.63	3.56 ± 0.89	4.62 ± 1.95	7.49 ± 3.10	4.12 ± 1.37	6.11 ± 2.90	9.66 ± 3.74	7.05 ± 3.32
ROG	18.43 ± 0.10	18.75 ± 0.13	18.57 ± 0.14	44.93 ± 1.32	47.69 ± 0.63	45.16 ± 1.00	58.41 ± 1.62	60.94 ± 1.71	60.58 ± 1.41
SASA	8680.43 ± 213.87	9548.53 ± 376.15	9372.75 ± 285.18	56596 ± 906.08	59078.34 ± 604.38	58152.34 ± 731.50	80911.81 ± 773.16	82721.99 ± 1209	82051.91 ± 921.76

ACE-2 complexes, where the RBD-ACE-2DT showed the highest instability as evidenced by its high average RMSD value of 10.58 Å compared to the RBD-ACE-2WT and SPIKE-ACE-2OM complexes with average RMSD of values of 9.23 Å and 7.81 Å, respectively (Fig. 6a). The RMSD plots of the RBD, SPIKE and RBD-ACE-2 showed that both OM and DT mutations affect the stability of the protein but, with the DT variant exhibiting pronounced instability in its protein complexes.

To understand the compactness of the alpha carbon backbones of the protein complexes, the RoG values for each complex were examined. The RoG value is a measure of the extent of compactness of the alpha carbon backbones of the proteins. An increase in RoG values implies a decrease in protein structure compactness, thereby suggesting decreased stability [44]. Fig. 4, Fig. 5, Fig. 6b showed the RoG plots for all the RBD, SPIKE and RBD-ACE-2 systems, respectively and Table 4, showed the average RoG values for all the systems. The result revealed that RoG plots of RBDDT and RBDOM exhibited decrease in structural compactness suggesting high mobility and less stability as evidenced by the higher average RoG values of 18.75 Å and 18.57 Å for RBDDT and RBDOM, respectively when compared to the RBDWT (18.43 Å). For the whole spike (SPIKE) systems, a similar trend was observed. The RoG plots of.

SPIKEDT and SPIKEOM demonstrated decline in structural compactness suggestive of high mobility and less stability as supported by the higher average RoG values of 47.69 Å and 45.16 Å for SPIKEDT and SPIKEOM, respectively when compared to the SPIKEWT (44.93 Å). Likewise, just as observed in both RBD and SPIKE systems, the estimated RoG plots of RBD-ACE-2DT and RBD-ACE-2OM also demonstrated decline in structural compactness indicative of high mobility and less stability as evidenced by the higher average RoG values of 60.94 Å and 60.58 Å for RBD-ACE-2DT and RBD-ACE-2OM, respectively when compared to the RBD-ACE-2WT (58.41 Å). The results of the RoG plots of the RBD, SPIKE and RBD-ACE-2 showed that both OM and DT mutations disturb the stability of the interaction between hACE-2 and RBD, however, the DT mutations exhibit more obvious instability in its protein complexes than the OM mutations.

Fig. 4, Fig. 5, Fig. 6c showed the RMSF plots for all the RBD, SPIKE and RBD-ACE-2 systems, respectively and Table 4, showed the average RMSF values for all the systems. RMSF is a measure of the impacts of the binding of molecule on the behavior of the active residue [45]. And high RMSF values indicated increase flexible movements, and in contrast, lower values meant restricted fluctuations. For the RBD systems (Fig. 4c), RBDDT and RBDOM showed higher flexible movements as evidenced by higher overall average values of 18.75 Å and 18.57 Å, compared to the RBDWT. Nevertheless, at amino acid residues 1060–1070 (the residues on RBD involved in interaction with hACE-2), all the RBD systems showed a noticeable increase in fluctuation, with the RBDOM exhibiting more flexibility than others. This finding is suggestive of how much activity and flexibility required of the interacting amino acid residues in interacting with hACE-2. A look at the whole spike structure, we observed the SPIKEDT exhibited higher fluctuation (7.49 Å) than both the SPIKEWT and SPIKEOM. This finding further corroborates our previous findings that showed that the DT mutations mostly disturbed the protein stability than the SPIKEOM mutations. Similarly, the both the RBD-hACE-2DT (9.66 Å) and RBD-hACE-2OM (7.05 Å) induces more flexible movements on the amino acid residues, with both variants showing higher flexibility than the RBD-hACE-2WT (6.11 Å).

To measures the proteins exposure to solvent molecules, the SASA plots for all the systems were examined (Fig. 4, Fig. 5, Fig. 6). A high SASA value has been reported to be an indication of decrease in the exposure of buried hydrophobic residues which suggest decrease in systems stability [46,47]. The result of average SASA values for the systems (Table 4), showed both RBDOM and RBDDT increase the SASA plots and average values far higher than the RBDWT. This finding indicates that both RBDOM and RBDDT decreased the exposure of buried hydrophobic residues which means decrease in the RBD stability. Likewise, in the whole spike complexes, the two variants decreased the SPIKE stability as proven by high SASA values of 59078 Å2 (SPIKEDT) and 58152 Å2 (SPIKEOM) compared to the SPIKEWT (56596 Å2). Equally, as reported in both the RBD and SPIKE systems, in the RBD-hACE-2 systems the RBD-hACE-2OM (82051 Å2) and RBD-hACE-2DT (82721 Å2) increase the SASA average values far higher than the RBDWT, which indicates that both RBDOM and RBDDT decreased the exposure of buried hydrophobic residues which mean decrease in the RBD-hACE-2 complex stability.

The findings of this study from all the structural stability parameters uniformly showed that mutations from both DT and OM variants disturbed the stability of either the spike protein or the RBD-hACE-2 complex. From the examination of the mutations present on the RDBs of the two variants; OM (G339D, S371L, S373P, S375F, K417 N, N440K, G446S, S477 N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H) and DT (L452R, T478K), it could be inferred that the presence of 15 mutations on OM RBD increases its affinity to hACE-2. Only mutation T478K is common to both variants RBD. Furthermore, the study reveal instability impact of the mutations on DT variant is more pronounced than the OM variant as evidenced by higher average values of structural stability plots and average values. This study, therefore, assumed this obvious instability observed in DT variant might be associated or responsible for the reported severity in DT variant disease than the OM variant disease [48,49]. Furthermore, this study believes, as shown in our result (binding energy together with the molecular docking results), that high transmission of the OM variants is associated with higher receptiveness reported to exist between OM spike protein and hACE2.

Conclusion

After the emergence of Severe Acute Respiratory Syndrome Coronavirus type 2 (SARS-CoV-2), variants Delta (B.1.617.2) and Omicron (B.1.1.529) have proven to be of major concern out of all the reported variants, considering their influence on the virus' transmissibility and severity. Mutations on the variants’ spike protein have either increase the disease transmissibility or severity. Therefore, this study evaluates the impact of mutations on these two variants on stability and molecular interactions between the viral Spike protein and hACE-2. The findings of this study showed the OM variant has high receptiveness towards hACE-2, thereby, increase its transmission. The study further showed that DT and OM variants disturbed the stability of either the spike protein or the RBD-hACE-2 complex, with DT variant having more instability impact. This study provides molecular dynamic insight into the effects of OM and DT variants on the stability and molecular interactions between SARS-CoV-2 Spike protein and hACE-2.

Evaluation of a monomer of Spike protein instead of trimers over a molecular dynamic simulation period of 50 ns is the only limitation to this study. This is because molecular dynamic simulation of complex structures of protein-protein interactions is relatively difficult. However, this is not expected to change or affect the findings of this study. Furthermore, similar studies have also worked on spike protein monomers.

Funding source

This project is sponsored by the NIH-RCMI (U54MD007605) and CPRIT-CTTP (RP210043) Grants.

Declaration of competing interest

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