Computation of the Binding Energies between Human ACE2 and Spike RBDs of the Original Strain, Delta and Omicron Variants of the SARS-CoV-2: A DFT Simulation Approach.

The receptor binding domain (RBD) of SARS-CoV-2 binds to human ACE2 leading to infection. In this study, the complexes that are formed by the attachment of the SARS-CoV-2 spike RBDs of the original strain, delta and omicron variants to the human ACE2 are investigated via density functional theory (DFT) simulations to obtain binding energies. The DFT computations are performed without fragmenting the interfaces to involve longer-range interactions for improved accuracy, which is one of the primary features of the approach used in this study. Basis set superposition error corrections and van der Waals dispersions are also included in the DFT simulations. The binding energies of the SARS-CoV-2 spike RBDs of the original strain, delta and omicron variants to the human ACE2 are computed as -4.76, -6.68, and -11.77 eV, respectively. These binding energy values indicate that the binding of the omicron variant to the ACE2 is much more favorable than the binding of the original strain and the delta variant, which constitute a molecular reason for the takeover of the omicron variant. The binding energies and the decomposition of these energies found in this study are expected to aid in the development of neutralizing agents.

Introduction

The coronavirus that causes the present coronavirus‐19 disease (Covid‐19) pandemic is the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2).[ , ] The coronavirus family, including the severe acute respiratory syndrome coronavirus 1 (SARS‐CoV‐1), is known to start the infection process by attaching to the receptors on the host cell.[ ] The receptor binding domains (RBDs) of coronavirus spike proteins bind to human angiotensin‐converting enzyme 2 (ACE2). Understanding the molecular level binding processes of the spike RBDs to the human ACE2 is obviously critical for the development of vaccinations and treatments.[ , ] There are clear differences in the spike proteins of the SARS‐CoV‐1 and SARS‐CoV‐2.[ ] In addition, a large number of variants of the SARS‐CoV‐2 also appeared since the start of the Covid‐19 pandemic and these variants are spotted by the alterations in their spike proteins.[ , , , ] The current variants of concern (VoCs) are the alpha (lineage B.1.1.7), beta (lineage B1.351), gamma (lineage P.1), delta (lineage B.1.617.2) and omicron (lineage B.1.1.529) variants according to the World Health Organization (WHO).[ ] Since the spike RBD–human ACE2 binding is the first and the critical step of the infection process, it is important to understand the binding energies and the related properties of the spike RBD–human ACE2 complexes. The characteristic of this binding mechanism is a factor affecting the transmission rate and it can be presumed that if the binding of the spike RBD of a specific variant to human ACE2 is more favorable compared to other variants, that specific variant has an increased potential to take over. Therefore, the investigation of the spike RBD–human ACE2 complexes from the binding energy point of view has the potential to provide important information for the comparison of these variants, which is the evaluated in this work.

In order to compute binding energies of the spike RBD–human ACE2 complexes, the atomic coordinates of these structures are needed. X‐ray crystallography (XRC) and cryogenic electron microscopy (cryo‐EM) are crucial for determining the molecular structure of biological molecules. In the literature, these approaches have been used to determine the atomic structures of spike protein–human ACE2 complexes[ ] and spike proteins.[ ] Despite the fact that cryo‐EM and XRC offer atomic coordinate data, binding energies and related values must be calculated using density functional theory (DFT) or other methods. The atomic coordinate data of spike RBD–human ACE2 complexes have been described in the literature,[ , ] with the prefusion conformation of spike RBDs being more favorable for attachment to human ACE2.[ ]

Various methods have been used in the literature to explore the binding characteristics of spike–ACE2 complexes.[ , ] Computational methods such as molecular dynamics (MD)[ ] and ab initio methods[ ] have been utilized to determine the binding characteristics of these compounds. For example, using classical MD simulations, a recent study estimated the binding energies of SARS‐CoV‐1–human ACE2 as −10.81 kcal mol−1 and SARS‐CoV‐2–human ACE2 as −12.86 kcal mol−1.[ ] They also calculated the binding energies of the alpha and beta mutant SARS‐CoV‐2–human ACE2 complexes as −14.66 and −13.52 kcal mol−1, respectively.[ ] Another article demonstrated the interaction of the SARS‐CoV‐1 spike protein and human ACE2 molecule using DFT simulations by fragmenting the spike protein and ACE2 molecules to make the simulations manageable.[ ] Depending on the basis set, binding energies ranging from −340.46 to −404.26 kcal mol−1 are obtained.[ ] Using a similar division technique with an interaction limit of 4.5 Å, the same research group recently examined the binding properties between SARS‐CoV‐2 and human ACE2 and obtained a binding affinity between −45.02 kcal mol−1 and −98.60 kcal mol−1.[ ] In another study, the Gibbs binding energy of the B.1.1.7 mutated spike protein to the ACE2 is found to be much higher than the original strain.[ ] There are also several studies regarding the effect of the spike protein mutations on the binding properties of SARS‐CoV‐2 in which MD simulations are employed. For example, the attachment of the SARS‐CoV‐2 variants to antibodies are investigated using MD computations.[ , ] Similarly, there are studies regarding the computation of the binding affinity between the mutated spike proteins and human ACE2 using MD methods.[ , ] There are also artificial intelligence (AI) assisted MD simulations for the prediction of the binding energy differences among mutations.[ ] In addition, molecular docking studies provide predictions on the binding properties of inhibitors attaching to SARS‐CoV‐2 spike protein.[ , ] In another study, the thermodynamic integration method in conjunction with MD is used for the computation of transmission ability of mutated SARS‐CoV‐2.[ ] Coarse‐grained model calculations are also used to gather information on the differences of the binding properties of SARS‐CoV and SARS‐CoV‐2 to the human ACE2.[ ] Anchor‐locker binding mechanism of the spike protein to the human ACE2 is analyzed in another study providing MD‐based binding energy values.[ ] The binding energy computations of the UK‐variant and the human ACE2 and miniprotein drug candidates are also performed using molecular mechanics Poisson‐Boltzmann surface area based simulations.[ ] The thermodynamics of SARS‐CoV‐2 spike–ACE2 recognition using MD‐based calculations are also studied in the literature.[ ] In another work, a combined analysis method employing MD, machine learning and free‐energy perturbation components is used to predict the binding mechanism differences between SARS‐CoV and SARS‐CoV‐2.[ ] In addition, enhanced sampling simulations are also employed to explain the attachment of SARS‐CoV‐2 spike protein to human ACE2 for exposing the binding pathways.[ ]

In this work, we have investigated the binding properties and binding energies of the spike RBDs of the original SARS‐CoV‐2 strain, delta and omicron variants to the human ACE2 considering that delta and omicron variants are on the VoC list of both WHO and U.S. Centers for Disease of Control and Prevention.[ , ] All of the atoms at the 15 Å distance to the spike–ACE2 interface are considered in the DFT computations without any division approach for the inclusion of longer‐range effects. The van der Waals dispersion and basis set superposition error (BSSE) corrections are also utilized for the accurate calculation of the binding energies. It is worth mentioning that the use of DFT in molecular analysis of the properties of SARS‐CoV‐2 is expected to give accurate results as shown in the literature.[ , ] Our DFT results indicate that SARS‐CoV‐2 spike RBD–human ACE2 binding reactions are spontaneous and that the binding energy of the spike RBD of the omicron variant has the highest absolute value which constitute a factor for the trend of the higher transmission rate of this variant at the molecular level. The roots of the difference of the obtained binding energies are further investigated using the energy decompositions to provide a better understanding of the binding mechanism.

Experimental Section

In this study, the total energies of the RBDs of the spike proteins, binding sections of ACE2 protein and the spike–ACE2 complexes are obtained using DFT. The binding energies are calculated as the difference of the bound and unbound states as shown in Equation (1)[ ]

The atomic coordinate data of the original strain, delta variant and the omicron variant are taken from the protein data bank having the PDB IDs as 6M0J,[ ] 7V8B,[ ] and 7T9L,[ ] respectively. All of these data are produced utilizing cryo‐EM method; therefore, they do not contain hydrogen atoms. Hence, hydrogens are added in the Chimera software.[ ] All of the atoms in the 15 Å neighborhood of the spike RBD–ACE2 interface are selected for the DFT computations without fragmentation to include longer‐range interactions. The selection of the 15 Å proximity is imposed by the amount of the computer memory available on the simulation server. The structures considered in the DFT simulations are shown in Figure .

Figure 1

Spike protein RBD – ACE2 interfaces for the structures of a) the original strain, b) the delta variant, and c) the omicron variant.

In order to select all of the atoms in the 15 Å proximity of the spike–ACE2 interface, the automatic intersurf operation built in the Chimera software, which enables to extract the interface surface to create a valid interaction map, is utilized.[ ] For the 15 Å neighborhood, the number of atoms in the spike RBDs, ACE2 sections, and the spike RBD–ACE2 complexes are given in Table . As it can be seen from Table 1, the number of atoms are large for DFT computations as the simulation complexity and cost increase quadratically by the number of atoms. Despite this computational complexity, the structures are not fragmented but taken as a whole for the inclusion of the longer‐range electrostatic and exchange‐correlation interactions in the DFT simulations.

Table 1

The number of atoms in the interface model of spike RBDs, ACE2 sections and spike RBD–ACE2 complexes considered in DFT computations

SARS‐CoV‐2 variant	Number of atoms in the spike protein	Number of atoms in the ACE2 molecule	Number of atoms in spike–ACE2 complex
Original strain (PDB ID: 6M0J)	1202	2061	3263
Delta variant (PDB ID: 7V8B)	1236	2197	3433
Omicron variant (PDB ID: 7T9L)	1201	2126	3327

All of the structures are geometrically optimized using the limited‐memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) method in DFT with a step size limit of 0.02 Å at 37 °C, which is set in the linear combination of atomic orbitals (LCAO) calculator.[ ] The maximum forces after the optimization step are observed as 0.096, 0.097, 0.090 eV Å−1 for the spike RBD, ACE2 section and spike RBD–ACE2 complex of the original strain; 0.079, 0.083, 0.082 eV Å−1 for the spike RBD, ACE2 section and spike RBD–ACE2 complex of the delta variant; 0.082, 0.079, 0.087 eV Å−1 for the spike RBD, ACE2 section and spike RBD–ACE2 complex of the omicron variant, respectively. These force values indicate that these optimized structures are suitable for the actual DFT simulations for the calculation of total energies.

A commercial DFT software package, QuantumATK from Synopsys, is utilized for the DFT simulations.[ , ] Basis set superposition error (BSSE) corrections[ ] and van der Waals dispersion interactions[ ] are included in the computations for increased accuracy. Using DFT, the many‐body electron structure is expressed utilizing the one electron Kohn–Sham Hamiltonian as shown in Equation (2)[ , ]

In Equation (2), ℏ is the reduced Planck's constant, m is the electron mass, and V eff is the effective potential including the potential created by other electrons and the external field which can be formulated as in Equation (3)

In Equation (3), V H is the Hartree potential which denote the electrostatic interactions of electrons, V xc is the exchange‐correlation potential caused from quantum mechanical interactions and V ext is the possible external potential imposed on the system.[ , ] On the other hand, the many‐body electron density, n(r), is given as the superposition of the occupied states as shown in Equation (4)

In Equation (4), f FD is the Fermi‐Dirac probability distribution function and ψ(r) is an occupied eigenstate. In DFT, the physical parameters of the system, including the total energies, are calculated as a functional of the many‐body electron density. The total energy of the system is expressed using the energy terms which are the functionals of the electron density as in Equation (5)[ ]

In Equation (5), T[n(r)] is the kinetic energy of the electrons, E XC[n(r)] is the exchange‐correlation energy, E H[n(r)] is the Hartree energy and E ext[n(r)] is the interaction energy of electrons due to possible external potentials. The energy values computed in our study include dispersion corrections proposed by Grimme which enable to model the van der Waals interactions as given in Equation (6)[ , ]

In Equation (6), E total is the total energy computed using Equation (4), E disp is the dispersion correction term and E total_disp_corrected is the energy value that includes van der Waals interactions. On the other hand, the basis set superposition error (BSSE) corrections have to be utilized during the calculation of the total energies of the spike RBD–ACE2 complexes since BSSE occurs due to the overlap of the localized basis sets in close proximity systems as depicted in Figure .[ , , , , ]

Figure 2

Concept of the basis set overlap during the simulation of the Molecule A–Molecule B complex.[ ]

The BSSE errors can be minimized using a counterpoise correction term (E CC) as shown in Equation (7)

The E CC term is computed using AB basis sets and taking the reciprocal molecule as ghost atoms as shown in Equation (8)

In Equation (8), denotes the DFT total energy of the Molecule A in the AB basis where the atoms of Molecule B are taken as the ghost atoms.[ ] It is worth noting that the computation of the BSSE correction term E CC requires the DFT calculations of , , , and ; therefore, BSSE corrected energy (E BSSE_corrected) calculations require five different DFT loops including the calculation of the uncorrected energy, E DFT.[ ] Considering the relative higher number of atoms in the spike–ACE2 complexes as shown in Table 1, the BSSE corrected DFT energy computations require considerable amount of computational power, which is provided by the institution of the authors.

Results and Discussion

The DFT simulations are performed using the QuantumATK software on workstations since the DFT simulations of the considered structures require high amount of memory due to the number of atoms. Double zeta polarized LCAO basis sets are used.[ , ] In addition, generalized gradient approximation (GGA) type functionals whose value is a functional of both the local and gradient of the electron density as given in Equation (9) are used[ ]

The same mesh cut‐off energy having the value of 200 Ry with an iteration limit of 10−4 eV are utilized for all of the DFT computations to have a proper comparison of the binding energies. The total energies computed using DFT with the van der Waals and BSSE corrections and the binding energies, which are calculated utilizing these energy values in Equation (1), are shown in Table . The binding energy value of the spike–ACE2 of the original strain is computed as −4.76 eV (−109.77 kcal mol−1) in our study while the binding energy was reported to be in the range of −45.02 and −98.60 kcal mol−1 employing DFT‐based calculations in the literature depending on the basis set utilized.[ ] As it can be seen from these values, the order of the binding energy obtained in our study for the spike–ACE2 of the original strain is close to the previously reported values obtained using DFT and the difference could be originated from the basis sets employed in DFT computations.

Table 2

BSSE and van der Waals dispersion corrected total energies obtained from the DFT computations and the calculated binding energies

SARS‐CoV‐2 variant	Total energy of the spike RBD [eV]	Total energy of the ACE2 section [eV]	Total energy of the spike RBD–ACE2 complex [eV]	Binding energy
Original strain	−152302.38	−264951.85	−417258.99	−4.76 eV (−109.77 kcal mol−1)
Delta variant	−155404.89	−280759.16	−436170.73	−6.68 eV (−154.04 kcal mol−1)
Omicron variant	−149543.26	−272584.70	−422139.73	−11.77 eV (−271.42 kcal mol−1)

The following results are deduced from the binding energy values shown in Table 2: i) The binding energies are all negative meaning that spike–ACE2 attachment reactions are spontaneous as expected, ii) the absolute values of the binding energies of the delta and omicron variants are higher than that of the original strain indicating a molecular level reason of the taking over of these variants compared to the original strain, iii) The binding energy of the spike RBD of the omicron variant–ACE2 have the highest absolute energy meaning that the binding of the spike RBD of the omicron variant to the ACE2 is much more favorable than the attachment of the spike RBDs of the original strain and the delta variant. The higher binding affinity of the omicron variant to the ACE2 compared to the original strain as found in our study is consistent with the results of the MD and docking based computational studies existing in the literature.[ , , , , ] On the other hand, experimental studies report differing results. The surface plasmon resonance (SPR) based experimental studies also support our findings that the binding affinity of the omicron variant to the ACE2 is higher than that of the original strain[ , ], while a noncompetitive enzyme‐linked immunosorbent assay (ELISA) based experimental study reported a comparable binding affinity for the omicron variant and the original strain.[ ] It is worth noting that the binding affinity of the omicron variant to the ACE2 is reported to be 2.4 times higher than the binding affinity of the original strain to the ACE2 in one of the SPR‐based experimental studies[ ] while this ratio is (−11.77 eV)/(−4.76 eV) = 2.47 in our DFT based study.

The reasons of the difference of the binding energies are further investigated by interpreting the changes in the energy components. The total energies computed using DFT in QuantumATK can be expressed as in Equation (10) by merging Equations (5)–(7) and then incorporating the entropy term E which denotes the contribution of entropy because of the used occupation function[ , ]

The change of each energy component is obtained using the energy decomposition reported by QuantumATK as given in Tables 3, 4, 5 for the original strain, delta and omicron variants, respectively. Energy differences denote the binding energy components obtained by the difference of the bound and unbound states as mentioned before. It is worth noting that the external energy component E ext[n(r)] is zero for all of the structures therefore not given in these tables.

Table 3

Energy components and binding energy decompositions for the original strain

Energy component	Spike RBD [eV]	ACE2 [eV]	Spike RBD–ACE2 [eV]	Bound–unbound energy difference (binding energy component)
Kinetic energy (T)	105 783.99	183 192.98	288 971.81	−5.16 eV (−118.99 kcal mol−1)
Exchange‐correlation (E xc)	−44809.44	−77846.08	−122655.58	−0.06 eV (−1.39 kcal mol−1)
Electrostatic energy (E H)	−213224.04	−370201.79	−583428.57	−2.74 eV (−63.18 kcal mol−1)
Entropy term (E S)	−0.74	−2.49	−2.71	0.52 eV (11.99 kcal mol−1)
Grimme correction (E disp)	−52.15	−94.47	−152.08	−5.46 eV (−125.91 kcal mol−1)
BSSE correction (E CC)	0	0	8.14	8.14 eV (187.71 kcal mol−1)
Total [eV]	−152302.38	−264951.85	−417258.99	−4.76 eV (−109.77 kcal mol−1)

Table 4

Energy components and binding energy decompositions for the delta variant

Energy component	Spike RBD [eV]	ACE2 [eV]	Spike RBD–ACE2	Bound–unbound energy difference(binding energy component)
Kinetic energy (T)	108 000.40	194 112.79	302 094.23	−18.96 eV (−437.22 kcal mol−1)
Exchange‐correlation (E XC)	−45819.39	−82650.91	−128465.66	4.64 eV (107.00 kcal mol−1)
Electrostatic energy (E H)	−217533.08	−392118.09	−609647.63	3.54 eV (81.63 kcal mol−1)
Entropy term (E S)	−1.16	−2.14	−2.73	0.57 eV (13.14 kcal mol−1)
Grimme correction (E disp)	−51.66	−100.81	−154.82	−2.35 eV (−54.19 kcal mol−1)
BSSE correction (E CC)	0	0	5.88	5.88 eV (135.60 kcal mol−1)
Total [eV]	−155404.89	−280759.16	−436170.73	−6.68 eV (−154.04 kcal mol−1)

Table 5

Energy components and binding energy decompositions for the omicron variant

Energy component	Spike RBD [eV]	ACE2 [eV]	Spike RBD–ACE2 [eV]	Bound–unbound energy difference(binding energy component)
Kinetic energy (T)	103 965.87	188 679.48	292 617.16	−28.19 eV (−650.06 kcal mol−1)
Exchange‐correlation (E xc)	−44235.87	−80025.54	−124258.57	2.84 eV (65.49 kcal mol−1)
Electrostatic energy (E H)	−209218.57	−381138.27	−590348.09	8.75 eV (201.77 kcal mol−1)
Entropy term (E S)	−1.71	−2.32	−2.49	1.54 eV (35.51 kcal mol−1)
Grimme correction (E disp)	−52.98	−98.05	−155.19	−4.16 eV (−95.93 kcal mol−1)
BSSE correction (E CC)	0	0	7.45	7.45 eV (171.80 kcal mol−1)
Total [eV]	−149543.26	−272584.70	−422139.73	−11.77 eV (−271.42 kcal mol−1)

The binding energy decompositions shown in Table 3, 4, 5 are plotted in Figure for clarity. As it can be seen from Figure 3, the dominant factor for the higher binding energy of the spike RBD of the omicron variant to the ACE2 is the kinetic energy term (T). Therefore, the roots of the higher binding energy of the omicron variant compared to the delta variant and the original strain are caused mainly from the higher binding energy component of the kinetic energy and not the van der Waals dispersion energy (E disp), BSSE counterpoise correction energy (E CC) or the exchange‐correlation energy (E XC). Similarly, the reason of the higher binding energy of the delta variant compared to the original strain is also caused from the higher kinetic energy component rather than the remaining energy factors. The binding energy values of the spike RBDs of the original strain, delta and omicron variants to the ACE2 reported in Tables 2, 3, 4, 5 and Figure 3 seem to be consistent with the trend of the transmission rates observed for these variants.[ , ] Therefore, it can be argued that accurate DFT computations of the binding energies of possible future variants of SARS‐CoV‐2 to the ACE2 have the potential to aid the forecast of the trend of their transmission rates ahead of their spread.

Figure 3

The decomposition of the binding energies of the spike RBD–ACE2 complexes of the original strain, delta and omicron variants (T: kinetic energy, EXC: exchange‐correlation energy, EH: electrostatic energy, ES: entropy term, Edisp: Grimme correction term, ECC: BSSE correction term).

Conclusion

In this work, the binding energies of the spike RBDs of the original strain, delta and omicron variants to the human ACE2 are obtained using DFT simulations. The investigated variants are chosen according to the variants of concern lists of the WHO and US CDC. The DZP basis sets with the GGA type exchange‐correlation functionals and a higher cut‐off energy of 200 Ry are used in DFT simulations to obtain accurate results. The van der Waals dispersion interactions and basis set superposition error (BSSE) corrections are also included in the computations for increased accuracy. All of the atoms in the 15 Å neighborhood of the spike RBD–ACE2 interface are included in the computations without the fragmentation of the structures to include longer‐range interactions. This method yields the simulated spike–ACE2 complexes to include 3263, 3433, and 3327 atoms, which are considerably large structures for the DFT computations considering today's computational standards. The limited interface volume considered in the DFT computations can be expanded by the advance of the computational capabilities. The DFT results show that the binding energies of the spike RBDs of the original strain, delta and omicron variants to the human ACE2 have the values of −4.76, −6.68, and −11.77 eV, respectively. These values show a relation with the trend of the transmission rates of these variants. In order to further investigate the binding energies; the kinetic, electrostatic, van der Waals dispersion, entropy, exchange‐correlation and BSSE correction components of the binding energies are given. These binding energy decompositions show that higher binding energies of the spike RBDs of the omicron and delta variants to the ACE2 are caused from the elevated kinetic energy components. Considering the possible relation of the obtained binding energies of the investigated variants and their transmission rate trends, it can be argued that accurate DFT computations of the binding energies of possible future variants may give clues on the trends of their transmission rates before their spread.

Conflict of Interest

The authors declare no conflict of interest.