{'fraction_non_alphanumeric': 0.034904688304997426, 'fraction_numerical': 0.022926326635754765, 'mean_word_length': 5.062939247227862, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Most frequently harboured missense variants of hACE2 across different populations exhibit varying patterns of binding interaction with spike glycoproteins of emerging SARS-CoV-2 of different lineages

Since the emergence of SARS-CoV-2 at Wuhan in the Hubei province of China in 2019, the virus has accumulated various mutations, giving rise to many variants. According to the combinations of mutations acquired, these variants are classified into lineages and greatly differ in infectivity and transmissibility. In 2021 alone, a variant of interest (VoI) Mu (B.1.621), as well as, variants of concern (VoC) Delta (B.1.617.2) and Omicron (BA.1, BA.2) and later in 2022, BA.4, BA.5, and BA.2.12.1 have emerged. Since then, the world has seen prominent surges in the rate of infection during short periods of time. However, not all populations have suffered equally, which suggests a possible role of host genetic factors. Here, we investigated the strength of binding of the spike glycoprotein receptor-binding domain (RBD) of the SARS-CoV-2 variants: Mu, Delta, Delta Plus (AY.1), Omicron sub-variants BA.1, BA.2, BA.4, BA.5, and BA.2.12.1 with the human angiotensin-converting enzyme 2 (hACE2) missense variants prevalent in major populations. In this purpose, molecular docking analysis, as well as, molecular dynamics simulation was performed of the above-mentioned SARS-CoV-2 RBD variants with the hACE2 containing the single amino acid substitutions prevalent in African (E37K), Latin American (F40L), non-Finnish European (D355 N), and South Asian (P84T) populations, in order to predict the effects of the lineage-defining mutations of the viral variants on receptor binding. The effects of these mutations on protein stability were also explored. The protein-protein docking and molecular dynamics simulation analyses have revealed variable strength of attachment and exhibited altered interactions in the case of different hACE2-RBD complexes. In vitro studies are warranted to confirm these findings which may enable early prediction regarding the risk of transmissibility of newly emerging variants across different populations in the future.


Introduction
The novel coronavirus, named Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), also referred to as nCoV-19, was first identified in the Wuhan, Hubei province of China in late December 2019 as the causative agent of the Coronavirus Disease 2019 (COVID-19) [1]. Since then, the virus has spread all over the world, affecting more than 440 million people and causing more than 5.97 million deaths worldwide, as of March 2022. These numbers only represent the clinically diagnosed cases. The real number of infections is much higher because the virus can cause asymptomatic infections [2]. SARS-CoV-2 belongs to the Betacoronavirus genus under the Coronaviridae family. The single-stranded positive-sense RNA genome of the SARS-CoV-2 is approximately 29.9 kb in length. It is replicated by an RNA-dependent RNA polymerase (RdRP), encoded in the viral genome [3].
The virus spreads via respiratory droplets, targets the respiratory tract of humans, and infects the human angiotensin-converting enzyme 2 (hACE2) positive lung cells, which also express transmembrane protease, serine 2 (TMPRSS2) e.g: nasal and bronchial epithelial cells, also alveolar type II pneumocytes [4]. Upon entering the body it interacts with hACE2 of cells via its spike glycoprotein (S). The spike protein has two distinct subunits, S1 and S2, which are responsible for receptor binding and membrane fusion, respectively. Receptor-binding domain (RBD) of the spike protein interacts directly with the hACE2 and is a part of the S1 subunit, spanning residues 319-541. Cathepsin L and TMPRSS2 mediated priming of the spike protein also facilitates membrane fusion of the virus, through cleavage at the S1/S2 junction and S2' site [5].
Many generic therapeutic regimens implemented in the case of previous outbreaks of coronaviruses, i.e. SARS-CoV and MERS-CoV have been repurposed for SARS-CoV-2 infection. Moreover, with advanced computer-aided drug design, many specific therapeutic components in the forms of small molecules, monoclonal antibodies, siRNA and miRNA have been proposed since the emergence of the novel coronavirus [6][7][8].
In this short period of time, several vaccines have been produced with considerable efficacy and many are currently in clinical trial [9][10][11][12]. In addition, multi-epitope vaccine constructs have been designed through in silico approaches as potential prophylactics [13].
Like other RNA viruses, the SARS-CoV-2 has a high mutation rate, and this has caused the emergence of many variants of this virus since its identification [14]. The mutations that are necessarily present in all isolates of a particular variant of the virus are called lineage-defining mutations for that variant. An isolate is assigned to a variant clade based on the presence of certain lineage-defining mutations [15]. Based on the nature of the mutations, the variants may exhibit alteration to many different characteristics of the virus such as transmissibility, disease severity, immune escape, increase in virulence, change in chemical manifestations of disease, and diagnostic or therapeutic escape. These mutations may also contribute to community transmission and bring about detrimental changes in COVID-19 epidemiology.
According to the CoVariant portal of GISAID (https://covariants. org/), 27 variants have been identified so far. The World Health Organization (WHO) has declared 5 of them (Alpha, Beta, Gamma, Delta, Omicron) as Variant of Concern (VoC) and 2 of them (Mu, Lambda) as Variant of Interest (VoI). The VoCs are named according to their time of emergence with respect to the identification of the Wuhan strain, which is considered the reference for SARS-CoV-2. The Alpha variant emerged in the UK (reported in December 2020), Beta in South Africa (December 2020), Gamma in Brazil (January 2021), Delta and Delta plus in India (May 2021), and Omicron in South Africa (November 2021). The VoI Lambda emerged in Peru (December 2020) and Mu in Columbia (January 2021) [16]. The mutations cause changes in the respective protein's structures [17]. Mutations occurring in the spike protein have conferred the virus higher transmissibility and allow them to escape from the immune response generated by the monoclonal neutralizing antibody or vaccination [18]. There are nine mutations in the spike protein of Delta and Mu variant, 36  The host receptor hACE2 harbours many polymorphisms, of which, some are missense variants that result in changes in the amino acid residue of the protein, and subsequent conformational change. These may alter the structure of the hACE2 which may affect the attachment with the spike protein of SARS-CoV-2. The polymorphisms are inherited within and can be attributed to a population in a region, where one may prevail in a given population, and others may prevail in other populations inhabiting different regions. The population frequency of a particular polymorphism might appear as a small number in percent value, but technically it includes a very large number of people. For example, the P84T missense variant in the hACE2 occurs at a frequency of 3.58% in the South Asian population, according to gnomAD exomes r2.1.1 dataset, and is estimated to be harboured by almost 56 million people in India and Bangladesh alone [19]. Therefore, it is important to take into consideration the effect of hACE2 polymorphism while analysing the virulence of SARS-CoV-2. The mutations occurring within the viral genome as well as natural variants within its receptor, hACE2 may contribute to the binding affinity followed by the entry and transmission of the virus in a population or region-specific manner [20]. Thus, the unique combinations of population-specific hACE2 polymorphisms and SARS-CoV-2 spike protein mutations may elucidate the region-dependent emergence and preference of each SARS-CoV-2 variant for one or more particular populations.
Molecular docking analysis provides a prediction of the binding affinity and stability of a complex. Therefore, it is possible to detect alterations in binding affinity between the reference and the variant complexes, resulting from the minute deviations of the interacting protein structures. An increase or decrease in binding affinity depends on the molecular interactions between the interacting molecules. Amino acid substitutions, depending on their physicochemical properties, often disrupt existing interactions and establish new interactions. Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) is end-point free energy calculation methodology which is vigorously used in the prediction of binding free energy (BFE) and identification of appropriate protein-protein complexes. The technique in conjunction with perresidue energy decomposition of the binding free energy of a proteinprotein complex is particularly useful in highlighting the key residues present in the binding interface. Furthermore, Molecular dynamics (MD) simulation provides insights into the changes in a protein's structure and behaviour in physiological circumstances with respect to time. Initially, proteins were conceptualized to be rigid entities as the X-ray diffraction method produces structures that are fixed in a crystal lattice [21]. However, just like all other macromolecules in physiological systems, proteins exist in a dynamic equilibrium. Advances in computational technologies have enabled us to simulate a physiological system in a computer system and observe the dynamics of any biological macromolecule in that system [22]. Through MD simulation, the stability and dynamics of protein-protein interactions can be observed, thus providing insights into the interactions between proteins, and the behaviour of each of the protein's constituent amino acids during the course of simulation runtime.
Several in silico and experimental studies have shown promising associations of hACE2 polymorphisms with population-wise disease susceptibility and prevalence [23][24][25][26]. Additionally, in silico and experimental studies have also been conducted to evaluate the binding affinity of the novel SARS-CoV-2 variants to hACE2 [27][28][29]. These approaches generate observations regarding variations in receptor-binding interactions based on either receptor hACE2 or the viral spike protein RBD. We aim to bridge the gap in the comparative assessment of binding interactions between polymorphic hACE2 and emerging SARS-CoV-2 variants.
This study, firstly, aims to evaluate the structural effects of the mutations occurring in the spike glycoproteins of the recently emerged SARS-CoV-2 variants. Secondly, considering the structural deviations of the emerging SARS-CoV-2 variants induced by these mutations, we sought to investigate the role of polymorphisms of hACE2 on the variable susceptibility of these variants among different populations.

Retrieval of amino acid sequences and sequence variants of SARS-CoV-2 spike glycoprotein and human ACE2
The complete reference sequence of SARS-CoV-2 surface glycoprotein (accession: YP_009724390) was retrieved from the NCBI virus database  Table 1).
The reference sequence of human angiotensin-converting enzyme 2 (hACE2) (accession: Q9BYF1) was obtained from the UniProt database (https://www.uniprot.org/uniprot/Q9BYF1 Population frequencies of each variant were obtained from the gnomAD v2.1.1 dataset (https://gnomad.broadinstitute.org/). The most highly prevalent missense variant within a population was chosen as the representative missense variant of that population. For more than one variant classified into the same population, co-occurrence was tested for the variant pair using the Variant Co-Occurrence tool (https://gnomad. broadinstitute.org/variant-cooccurrence). Between the variants occurring in different haplotypes, the one with the higher population frequency was selected, as the representative variant of the respective population. Four SNPs were thus listed for four populations: African (AFR), Latin American (AMR), non-Finnish European (NFE), and South Asian (SAS); as reported in the gnomAD v2.1.1 dataset.

Construction and validation of mutant models of the viral spike glycoprotein receptor-binding domain (RBD) and human ACE2 protein using homology modelling
The three-dimensional structure of SARS-CoV-2 spike glycoprotein RBD bound with hACE2 (PDB ID: 6M0J) was retrieved from the Protein Data Bank (https://www.rcsb.org/structure/6m0j). The lineagedefining mutations were incorporated into the amino acid sequences of the respective spike glycoprotein variants (Delta, Mu, and Omicron). Chain E (RBD) of the hACE2-RBD complex (PDB ID: 6M0J) served as the template for the construction of the RBD models of the SARS-CoV-2 variants on the Swiss-Model server (https://swissmodel.expasy.org/) [37]. Chain A (hACE2) of the hACE2-RBD complex (PDB ID: 6M0J) was entered along with the single amino acid substitutions specified for the populations into the Missense 3D server (http://missense3d.bc.ic.ac. uk/missense3d/) [38] to generate the population-specific mutant models of hACE2 protein.

Molecular docking assay and analysis of protein-protein interactions
Template-based protein-protein docking between the hACE2 and RBD was performed using the HDOCK server (http://hdock.phys.hust. edu.cn/) [42]. Protein-protein docking utilises the measure of intermolecular interactions to predict the strength of a molecular complex. The protein chains were prepared for docking experiments by removing water, heteroatoms, and ligand groups. The interacting residues for hACE2 were set as reported in the UniProt database as well as those detected by the PDBSum server (http://www.ebi.ac.uk/thornton-srv/ databases/pdbsum) [43], using the hACE2-RBD complex (PDB ID: 6M0J) as reference. Using the same complex as the reference, the interacting residues of the RBD were set in a similar manner. The interacting residues of both hACE2 and RBD are presented in Table 1. The docking energy scores were used for further analysis.
Protein-protein interactions across hACE2 and RBD interfaces were visualised by BIOVIA Discovery Studio Visualizer v21.1.0.20298. Details regarding the two-dimensional interactions between the RBD variants and hACE2 (wildtype and selected missense variants) molecules were generated using the PDBsum server.

Prediction of binding energy and decomposition of free energy contributions
The binding energy of the reference and variant hACE2 -RBD complexes were calculated using an end-point free energy calculation methodology, Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) on the HawkDock server (http://cadd.zju.edu.cn/h awkdock/) [44]. The server analyses protein-protein interactions by integrating the ATTRACT docking algorithm, the HawkRank scoring function and the MM/GBSA free energy decomposition analysis. MM/GBSA is calculated based on the ff02 force field, the implicit solvent model and the GB OBC1 model (interior dielectric constant = 1). The whole system was minimized for 5000 steps with a cutoff distance of 12 Å for van der Waals interactions (2000 cycles of steepest descent and 3000 cycles of conjugate gradient minimizations).
In MM/GBSA analysis, the binding free energy (ΔG bind ) of the system can be expressed as a change in free energy, which can be calculated by Equation (1) ΔG bind = ΔE internal + ΔE electrostatic + ΔE vdw + ΔG GB + ΔG SA -TΔS (1) Here, ΔE internal denotes change in internal MM energy, ΔE electrostatic indicates change in electrostatic energy, ΔE vdw depicts change in van der Waals energy, ΔG GB denotes the electrostatic solvation energy, ΔG SA depicts the deviation of solvation free energy combined with nonelectrostatic solvation, and TΔS is the total system binding entropy.

Molecular dynamics simulation
To evaluate the changes in stability and binding patterns of the docked protein-protein complexes with time, we performed a 125 ns molecular dynamics (MD) simulation of the Omicron BA.1 and BA.2 complexes, as well as, the highest and lowest scoring complexes of Delta and Mu variants. The simulation was run using the Desmond v3.6 program in the Schrödinger suite [45]. A single-point charge (SPC) water model was used to solvate the system and an orthorhombic periodic boundary box with the dimension of 10 Å * 10 Å * 10 Å was employed to maintain a stable volume. Na + and Cl − ions were used at a concentration of 0.15 M to electronically neutralize the system. The OPLS3e force field was employed to optimize the system. Ensembles of the NPT (constant number of particles, pressure, and temperature) were maintained at 300.0 K temperature and 1.01325 bar pressure. After an initial relaxation, the system was recorded at 125 ps intervals with an energy of 1.2.
The results of the MD simulation are presented as the evaluation of the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) data.

RMSD analysis
RMSD is a measurement of the average displacement of selected atoms in a particular timeframe with respect to a reference frame [46]. During the initialization of the simulation, the protein structures were aligned with the backbone frame, and the RMSD was calculated using Equation (2).
Here, N is the number of selected atoms; t ref is the reference time (usually the first frame is considered the reference frame, and time is taken as t = 0); Frame x is the location of the chosen atoms when superimposed with the reference frame. Every simulation frame underwent this calculation to provide the RMSD data.

RMSF analysis
RMSF is useful to observe the behavior of each amino acid residue during the simulation period. It is an indicator of local conformational change around the amino acid considered. The calculation is done by using Equation (3).
where T = trajectory time interval used to calculate the RMSF values, t ref = the reference time interval, r i = the location of residue i, r' I = the location of the atoms in residue i after their overlap upon that reference, and the angle brackets (<>) indicates that the square distance is averaged across the residue's atoms.

Retrieved high-frequency population-specific hACE2 missense variants within the hACE2-RBD binding interface
Six unique SNPs, specifying the hACE2-RBD binding interface residues, were initially retrieved. The African and the Non-Finnish European population each had two missense variants within the specified regions. When tested for co-occurrence, both of these variant pairs appeared to occur on different haplotypes. For the rest of the populations, each had only one missense variant. The four hACE2 SNPs thus obtained according to the highest population frequency are provided in Table 2.
All the modelled RBD structures produced QMEAN "Z-score" for absolute quality estimates of protein models well above − 4.0 on Swiss-Model structure assessment indicating good quality models [47]. All RBD models scored an overall quality factor for non-bonded atomic interactions within a protein above 90 on ERRAT which satisfies the criteria for good quality model. The Ramachandran plots of the RBD models obtained from PROCHECK are presented in Supplementary  Fig. 1. 90% or more residues of all the modelled RBD structures were in energetically most favoured region, hence, qualify as acceptable for further analyses [48]. The respective Z-scores of all RBD models, produced by ProSA-web server, were consistent with those determined experimentally as depicted in the Z-score plots ( Supplementary Fig. 2).

Effects of spike protein mutations on protein stability and hACE2 binding affinity
Lineage-defining amino acid substitutions of the SARS-CoV-2 variants are presented in Fig. 1.
Predictions on the stability of 48 amino acid substitutions found in the spike protein of SARS-CoV-2 variants we investigated, provided by different tools, are presented in Fig. 2.
For the spike glycoprotein of the Delta variant, three mutations were destabilising out of the six single amino acid substitutions analysed. Two of these mutations (L452R and T478K) occupy the receptor-binding domain. In addition to these three destabilising mutations found in the Delta variant, a stabilising mutation (A222V) was identified outside the RBD of the Delta plus variant. In the case of the Mu variant, eight mutations were analysed and four of them were destabilising. The RBD of      mutation G496S. 13 out of 27 analysed mutations in Omicron BA.4, were shown to be destabilising, whereas, 12 were shown to be stabilising. For the most recent Omicron subvariant BA.2.12.1, 13 out of 27 mutations analysed were shown to be destabilising and 13 to be stabilising. Considering the emerging sequence of the variants, mutations are accumulating rapidly within the spike protein (Fig. 3), particularly, among RBD residues.

Docking energy score analysis of the hACE2-RBD complexes
The docking energy scores of thirty-five complexes of spike RBD variants with the reference and missense hACE2 variants prevalent in four major populations are depicted in Fig. 4

Comparative analysis of the intermolecular interactions for different hACE2-RBD complexes
The summary of the interactions derived from the docked complexes of the hACE2 and RBD has been presented in Fig. 5.
The depiction compares the RBD residues which undergo mutations in at least one variant investigated and/or forms a salt bridge. The residues are Arg403, Lys417, Gly446, Glu484, Gln493, Gly496, Gln498, Asn501, Tyr505 of the spike protein RBD. The complex formed by the reference hACE2 and the spike RBD of the Wuhan-Hu-1 variant has been referred to as the 'reference' complex. The rest of the complexes are 'variant' complexes, in which either hACE2 or RBD is a variant.
The notable H-bond and salt bridge interactions found within different variant hACE2-RBD complexes compared with the reference hACE2-RBD complex are depicted in Fig. 7.

Prediction of free energy of protein-protein complexes and free energy decomposition of protein-protein interactions
The binding free energies predicted for the hACE2 -RBD variant complexes are summarized in Supplementary Table 4. The binding free The free energy contributions at the level of amino acid residues were calculated using free energy decomposition analyses. The differences in residual free energy contributions of wild type hACE2 residues and the mutated residues (E37K, F40L, D355 N, and P84T) caused by SNPs, highly prevalent in African (AFR), Latin American (AMR), Non-Finnish European (NFE), and South Asian (SAS) populations are depicted in Fig. 9. The numerical values are provided in Supplementary  Table 5.
All mutated residues except F40L, were shown to shift the free energy towards a more positive value. Although it seems to decrease the binding affinity at a single residue level, the overall the binding free energy does appear to follow a similar trend. There is less deviation on the MD trajectory of the AFR hACE2 -Omicron BA.1 RBD complex at the initial stage, indicating more stable binding that was maintained for the runtime. The SAS hACE2 -Omicron BA.1 RBD complex also showed stronger initial interaction as it had lower RMSD at the beginning compared to the Ref hACE2      RMSF analysis of the protein backbone of the SARS-CoV-2 RBD variants revealed significant variations in RMSF values in the Omicron BA.1 RBD variants interacting with different hACE2 variants. The Asn360 was stable while interacting with AFR and SAS hACE2, however, exhibited fluctuation while interacting with AMR hACE2. The wild type Asn370 residue and the mutated Ser371Leu residue in the RBD were stabilised while interacting with AFR hACE2 but fluctuated higher in the rest of the cases. Leu387 and Asp389 showed high fluctuations while interacting with AMR hACE2 but were most stable in case of AFR and SAS hACE2. The mutated Gly446Ser residue exhibited greater stability while interacting with AFR and SAS hACE2 compared with the rest. The RBD region spanning residues 496-508 of Omicron BA.1 harbours four mutations and they have stabilised while interacting with all hACE2 variants compared with the reference hACE2.

RMSD analysis
In the Omicron BA.2 RBD variant complexes with different hACE2, greater fluctuations were observed at the RBD residues Asn360, Thr393, Gly413, Asp428, Asn460, Ile468, and Thr500 in interaction with the SAS hACE2 compared with the rest. The RBD residues of the AFR hACE2 -Omicron BA.2 RBD complex were most stabilised in terms of RMSF. Interestingly, we observed a high fluctuation in the Gly447 residue of the Omicron BA.2 RBD while interacting with reference hACE2. However, greatly reduced fluctuations were observed for interactions with SAS, NFE, and AFR hACE2. While interacting with AMR hACE2, the residue stabilised only moderately. On the other hand, the RBD residues of the Delta variant were observed to undergo major fluctuations in complex with reference hACE2, whereas, were substantially stabilised in complex with SAS hACE2.
The N-terminal and C-terminal residues of all complexes show higher RMSF value due to their positions at the ends of the protein molecule which fail to take part in any secondary structure formation, resulting in greater flexibility.

Discussion
Among the SARS-CoV-2 variants of concern, we especially investi-  As hACE2 serves as the primary receptor for the entry of SARS-CoV-2, it is the first probable 'genetic gateway' that might determine disease progression [49]. The genetic polymorphisms of hACE2 are not evenly distributed across populations. Rather, the frequencies of these polymorphisms considerably differ among populations. Hence, a population may preferentially harbour one polymorphism at a much higher frequency, compared with the rest. Moreover, both experimental and in silico studies demonstrated that the binding interactions of SARS-CoV-2 vary in presence of hACE2 polymorphisms [50]. Furthermore, previous experiments have shown that the SARS-CoV-2 variants also interact differently with hACE2 [51]. The spread of the viral variants, thus, is likely to differ extensively among the host populations harbouring distinct hACE2 polymorphisms. However, none of the currently available studies have considered the RBD mutations and the hACE2 polymorphisms simultaneously. We combined both the recently emerged SARS-CoV-2 variants and the hACE2 single missense variants prevalent in major populations in a combinational approach to investigate the differences in the receptor-binding interactions between different RBD -hACE2 complexes. Our results indicate notable differences in receptor-binding interactions of the spike RBD of different SARS-CoV-2 variants within a population represented by a hACE2 single missense variant. The receptor-binding affinity of a particular SARS-CoV-2 variant also varies across populations carrying distinct single missense variants of hACE2. The findings may provide useful insights into the differential transmission rate of the emerging SARS-CoV-2 variants over individual populations.
This study focuses particularly on the missense mutations of hACE2 protein, present within its binding interface with the RBD of the spike glycoprotein of SARS-CoV-2. Substitutions of these residues are more likely to impact the binding affinity of the hACE2-RBD complex. Accordingly, instead of the whole spike glycoprotein, our focus remained on its RBD, which directly interacts with hACE2. Furthermore, the hACE2 missense variants with the highest frequencies within a population were preferred as they represent a larger target population, broadening the scope of this study. The hACE2 missense variants assigned to the African (AFR), Latin American (AMR), non-Finnish European (NFE), and the South Asian (SAS) populations represent 5.12%, 0.001%, 1.75%, and 3.58% of the respective populations.
Protein stability and protein-protein interactions are greatly affected by the accumulation of both stabilising and destabilising mutations. It is hence worthwhile to study the mutational pattern of the SARS-CoV-2 variants. This is particularly necessary for the evaluation of the efficacy of protective treatments against the disease, as well as, for predicting the fate of the ongoing pandemic. Stabilising mutations in the spike protein might increase its longevity under unfavourable conditions, whereas, destabilising mutations could significantly reduce the fitness of the protein. However, destabilising mutations tend to enhance protein flexibility, which often considerably affects the infectivity of a viral strain through an increase in receptor-binding affinity [52]. On the other hand, stabilising mutations render structural rigidity which may confer increased stability of the protein-protein complex [53]. In this study, the lineage-defining mutations of SARS-CoV-2 variants: Delta, Delta plus, Mu, Omicron sub-variants BA.1, BA.2, BA.4 and BA.2.12.1 were analysed. The destabilising mutation K417 N in the RBD region of the Omicron variants was reported to increase receptor-binding affinity [54]. This mutation along with five other notable mutations (N440K, G446S, E484A, Q493R and N501Y) have shown to markedly reduce the activity of monoclonal antibodies: bamlanivimab, etesevimab, casirivimab, and imdevimab [55][56][57] and moderate reduction in the activity in tixagevimab and sotrovimab [57]. Another destabilising mutation, L452R, present in the Delta variant was found to resist the activities of several monoclonal antibodies in therapeutic use [58]. The destabilising mutation E484K, present in the RBD of the spike glycoprotein of the Mu variant, exhibited an alteration in the complementarity of binding plasma polyclonal neutralizing antibodies [59]. The stabilising mutation N501Y, shared by the Mu, Omicron sub-variants BA.1, BA.2, BA.4 and BA.2.12.1 was found to increase the hACE2-binding affinity of the RBD [60]. This may also present a possible route of immune escape [61].
We sought to explore the effects of the SARS-CoV-2 variant spike glycoprotein RBD mutations on the binding with hACE2 containing the most prevalent single amino acid substitutions in the RBD binding regions for the African, Latin American, non-Finnish European and South Asian populations using molecular docking. A more negative docking energy score should indicate higher stability between receptor and ligand, hence, lower binding energy, increased affinity, and stronger interaction [62]. Furthermore, MM/GBSA calculations are theoretically more rigorous than scoring functions as commonly used in many of the docking algorithms, and can be used to validate the docking energy scores generated in this manner [63].
Previous studies have confirmed that Lys417, Tyr449, Phe456, Phe486, Gln493, Gly496, Gln498, Asn501 and Tyr505 are some key residues essential in the formation of a stronger complex with hACE2 [64,65]. Among which, a stable salt bridge between Lys417 of RBD and Asp30 of hACE2, along with three sustained H-bonds between RBD residues: Tyr449, Gln493 and Gln498 and hACE2 residues: Asp38, Glu35 and Lys353, distinguishes SARS-CoV-2 from SARS-CoV [66]. Hydrogen bonding plays a key role in hACE2-RBD binding among other intermolecular interactions. Mutations in the primary interface of hACE2-RBD induce stronger binding by increasing H-bonds as well as expanding buried solvent accessible surface [67]. Hence, we compared two types of interactions in terms of their number: salt-bridge, and hydrogen bonds between the reference and variant complexes. Increase and decrease in the quantity of these bonds fairly correlate with the strength, and to some extent, stability of the reference and variant complexes under investigation.
Molecular Dynamics simulation was performed to further evaluate the stability of interactions between the RBD and hACE2 variants over time using Root Mean Square Deviation (RMSD) values of the hACE2-RBD complexes. The stability of a protein-protein complex decreases as the RMSD value increases [68]. MD simulations also revealed differences in fluctuation among the RBD variant residues while in complex with different hACE2 variants through Root Mean Square Fluctuation (RMSF) values. The RMSF values represent the degree of fluctuation for each of the amino acid residues of the molecule of interest. High RMSF values may be correlated with the flexibility of the protein which may alter a protein-protein binding [69].
The N501Y substitution, first emerged in the Mu variant, had been reported to increase the affinity of the hACE2-RBD complex by forming new protein-protein interactions that modify internal structural dynamics [60,70]. Preferential selection of the mutation N501Y in more transmissible variants were also observed in an in vitro evolution study of the spike glycoprotein RBD [71] Our result is consistent with this finding as the binding free energy drastically moves towards the negative in case of the Mu variant in complex with the reference hACE2 (− 80.43 kcal/mol) compared with that of the Delta variant (− 68.77 kcal/mol) which lacks the given mutation, indicating higher binding affinity of the Mu variant for the reference hACE2. In comparison with the NFE hACE2, the Mu RBD manifests lower average RMSD (3.5072 Å) in complex with AMR-hACE2. Furthermore, key RBD residues of this variant, except Phe486, show lower RMSF while in complex with AMR hACE2 and much higher RMSF while in complex with NFE hACE2, when compared with reference hACE2. This finding aligns with the failure of the variant of South American origin to propagate into other populations, especially the Non-Finnish Europeans, as it revealed a less stable interaction with the NFE hACE2, compared with the AMR hACE2 [72].
The Lys353-Tyr501 [N-H⋯O] bond was also present in the NFE hACE2 -Omicron BA.1 RBD complex which exhibited a much higher affinity (− 77.53 kcal/mol) compared with the variant complex with reference hACE2 (− 74.51 kcal/mol). However, a substantially high RMSD (3.3771 Å) was observed in the NFE hACE2 -Omicron BA.1 RBD complex compared with the reference hACE2 (2.7918 Å), which indicates the binding was not maintained properly during the simulation runtime. The Gly496 residue is mutated in the RBD of the Omicron BA.1,   The higher RMSF at residues 496-508 in both Omicron BA.1 and BA.2 complexes may relate to previous reports of reduced antibody response [74]. Although a lower RMSF indicates a more rigid behaviour of the residue in concern, it does not guarantee a more stable protein-protein interactions as protein-protein interactions largely rely on the conformational flexibility [75].
The exceedingly low affinity of the reference hACE2 -Delta RBD variant complex, compared with a much higher affinity of the variant complex with the SAS hACE2 might be explained by comparing their patterns of interactions. The Delta RBD forms much fewer H-bonds with the reference hACE2, involving only four (Gln493, Gln498, Asn501, Phe486) of the key interacting residues mentioned above. Each of the residues exhibit high RMSF values in complex with the reference hACE2. Only Gln493 forms H-bonds with Lys31 and Gln35 of hACE2. Whereas, the Delta RBD residues manifest much lower RMSF values while interacting with SAS hACE2, of which Gln493 and Gln498 form H-bonds and Lys417 forms a salt bridge with hACE2 residues Lys31 and His34, Asp38, and Asp30, respectively. During the course of 125 ns MD simulation, the Delta variant maintained a much stable binding in complex with the SAS hACE2 compared with that of reference hACE2. This outcome is also notably consistent with the devastating spread of the variant in South-East Asia, as a large population of around 49.7 million in India alone would represent the SAS hACE2 harbouring the P84T substitution [76].
Although stronger binding to the receptor yields a higher probability of receptor endocytosis of the virus and subsequently successful replication, it is not the only determinant of the rate of transmissibility of the virus and the severity of the disease. A furin-like protease, transmembrane serine protease 2 (TMPRSS2), plays an equal role in the virus's ability to infect a permissive cell. The protease cleaves the viral spike protein at the S1/S2 junction, enabling its fusogenic activity and allowing subsequent fusion with the host cell [77]. An earlier report has suggested the Delta variant outcompetes the Omicron BA.1 variant in infecting TMPRSS2 overexpressing VeroE6 cells in culture [78]. In our result, the Omicron BA.1 and BA.2 RBD variant exhibited much higher receptor affinity, compared with the Delta RBD variant. However, the Delta variant may overcome the lower receptor-binding affinity with the heightened ability of infection using TMPRSS2 activity. Nevertheless, Omicron variants have been shown to employ a TMPRSS2-independent pathway for infection [79]. Hence, the higher receptor-binding capacity of the Omicron variants may still imply higher transmissibility of these emerging variants.

Conclusions
Our study suggests variation in the transmissibility of the emerging SARS-CoV-2 variants across populations harbouring specific hACE2 missense variants. Although viral susceptibility at the population level depends on many factors related to individual hosts, this study may provide insights into transmissibility of emerging SARS-CoV-2 variants among different populations. Furthermore, our results also reinforce the available evidence supporting that the Omicron subvariant BA.2 shows the highest transmissibility as well as the increasing incidences of Omicron BA.2.12.1, among the recently emerged SARS-CoV-2 variants. The severity and outcome of disease, however, can only be predicted upon considering host immune responses and related factors. Nevertheless, predictions on variant transmissibility with accurate interpretations may encourage effective implementation of preventive measures to contain community transmission and the rate of infection caused by the newly emerging variants.

Author contributions
AHMNN conceived the idea and supervised the research work. AT, RA and PB retrieved the sequences. AT and RA performed modelling of the proteins. MA, AAS and SC performed the mutational analyses. AT, RA and MA performed the docking analyses. AT performed the MMGBSA analyses. RA performed the MD simulations and subsequent data analysis. AHMNN, TY, AT and RA prepared the first draft. AT, RA, PB, MA and TY prepared the revised manuscript. AHMNN, AKMMH and TY then curated the revised draft critically. All authors read and approved the revised manuscript.

Declaration of competing interest
The authors declare that there is no conflict of interest.