Binding of SARS-CoV-2 to Cell Receptors: A Tale of Molecular Evolution.

The magnified infectious power of the SARS-CoV-2 virus compared to its precursor SARS-CoV is intimately linked to an enhanced ability in the mutated virus to find available hydrogen-bond sites in the host cells. This characteristic is acquired during virus evolution because of the selective pressure exerted at the molecular level. We pinpoint the specific residue (in the virus) to residue (in the cell) contacts during the initial recognition and binding and show that the virus⋅⋅⋅cell interaction is mainly due to an extensive network of hydrogen bonds and to a large surface of noncovalent interactions. In addition to the formal quantum characterization of bonding interactions, computation of absorption spectra for the specific virus⋅⋅⋅cell interacting residues yields significant shifts of Δλmax =47 and 66 nm in the wavelength for maximum absorption in the complex with respect to the isolated host and virus, respectively.

Introduction

At the time of the writing of this manuscript, the situation regarding the global pandemic produced by the spread of the SARS‐CoV‐2 virus (over 50 million confirmed cases and 1.25 million deaths, with no end in sight), with dire consequences in all aspects of life, from social interactions, to the overwhelming of health and economic systems, is changing fast. Because this is a critical problem, just as the rate of virus transmission on the early stages of dissemination, the number of scientific papers on the subject (mostly preprints) increases exponentially.

SARS‐CoV‐2 is an enveloped virus of the Coronaviridae family with a single‐stranded RNA genome. Figure 1 highlights the most important structural features of the virus: besides the nucleocapsid (N) proteins, the only proteins in direct contact with the genetic material (the N‐RNA core is embedded in a lipid environment), there are membrane (M), envelope (E), and spike (S) proteins. It is the spike proteins which lead to the now familiar external morphology of the virus, but more importantly, S proteins are responsible for the interactions with receptors in the host membrane (epithelial cells in humans). These S⋅⋅⋅Receptor contacts initiate the infectious cycle of the virus.

Figure 1

Structural features of the SARS‐CoV‐2 virus and of the spike (S) protein. The envelope (E) and membrane (M) proteins as well as the lipid bilayer (LB), and nucleocapside–RNA core (N‐RNA) are highlighted in the virus. The receptor binding domain (RBD), the N‐terminal domain (NTD), and Section 2 (S2) are highlighted in one protomer extracted from the trimer constituting the S protein. We also show a cell with internal organelles and with a few enzymes that act as a virus receptors.

Each spike consists of a trimer of S proteins. Individual S proteins have been divided into two clear S1, S2 sections, with S1 containing the N‐terminal domain (NTD), and the receptor binding domain (RBD), the domain ultimately responsible for the interactions with the coupling factors present in cell membranes. Coupling factors include a variety of proteins, carbohydrates, or other types of biomolecules expressed on the surface of the cell membrane and in charge of signaling and transport, among other functions. Viruses take advantage of these molecules during the infection process. It seems well established that initial virus⇌host recognition and binding are driven by S1, and that further changes in the conformation of the S2 section mediate the viral envelope fusion to the host cell membrane.

The most commonly invoked culprit (with plenty of experimental evidence ) for the reception of SARS‐CoV‐2 is the angiotensin converting enzyme 2 (ACE2). This receptor is the subject of intensive studies aiming at finding effective therapies, and is the central focus of the ongoing race to find a vaccine. In this work, we are interested in two crucial aspects of the initial virus⋅⋅⋅cell interaction problem: to pinpoint the specific residue to residue binding sites between the structurally known spike proteins of the virus and the structurally known ACE2 receptor in cell membranes, and to understand, from a fundamental, quantum perspective, the molecular factors driving the virus⋅⋅⋅cell binding. We expect this knowledge to considerably better our understanding of the problem and to hopefully contribute to a rational design of drugs and vaccines to fight the virus.

Results and Discussion

See the Computational Methods for details of our calculations. Our data shows that the RBD(S)⋅⋅⋅ACE2 complex reached well defined persistent equilibrium states long before the 600 ns of the molecular dynamics (MD) simulation time are consumed in each of the three replicas. This stability is especially encouraging in the interaction region as clearly shown in the highlighted areas of the bottom panels in Figure 2. We obtained an interaction energy ΔG int=−419.91 kcal/mol, which is in excellent agreement with calculations reported in closely related systems.

Figure 2

Equilibration of the molecular dynamics simulations. Structural fluctuations of the receptor‐binding domain of the spike protein (left) and of the ACE2 receptor (right) in the RBD(S)⋅⋅⋅ACE2 complex. Root mean square deviations of atomic positions (RMSD) plots, averaged over three independent MD replicas are shown at the top. Plots of the root mean square fluctuations (RMSF) of the residues are shown at the bottom. The RBD(S)⋅⋅⋅ACE2 interaction region is highlighted.

Virus⋅⋅⋅cell contacts

In all cases, only hydrogen bonds (HBs) were found as responsible for explicit virus⋅⋅⋅cell pair‐wise interactions. Naturally, this does not mean that other weak, long range cumulative interactions are ruled out. The number of hydrogen bonds (Figure 3) fluctuates around 5 during the intermediate to late stages of the simulation of the RBD(S)⋅⋅⋅ACE2 complex. Fittingly, as is characteristic of equilibrated systems, the number of HBs is quite stabilized for the later steps of the simulations. Table 1 lists all individual binary contacts in the form of hydrogen bonds between residues in the RBD(S)⋅⋅⋅Receptor complex found in our MD simulations with an arbitrary threshold average of 15 % occupancy on the triplicate runs.

Figure 3

Hydrogen bonding during the MD simulations. The number of fragment‐to‐fragment hydrogen bonds in the RBD(S)⋅⋅⋅ACE2 complex is averaged over three independent replicas as the MD progresses.

Table 1

Properties of virus⋅⋅⋅cell hydrogen bonds. Persistent hydrogen bonds in the RBD(S)⋅⋅⋅receptor complex exceeding the 15 % occupancy threshold during the entire 600 ns MD simulation of each replica. The arrows state the directionality of the donor→acceptor interaction in the corresponding hydrogen bond according to the classical electrostatic Xδ−−Hδ+→Yδ− description. All occupancies averaged over the three MD replicas. WBI are the Wiberg bond indices. The archetypal hydrogen bond in the water dimer is included for comparison purposes. See the specialized literature[ , , , , , ] for the formalism on how NBO and QTAIM descriptors are related to bonding.

RBD(S)  ACE2	Occ	d	−Ed→a(2)	WBI	102 ρ(r c)	|ν(r_c)||G(r_c)|	102 H(r_c)ρ(rc)
	[%]	[Å]	[kcal/mol]		[a.u.]		[a.u.]
Gly502→Lys353	38.71	1.94	6.06	0.04	2.51	0.92	6.8
Lys417→Asp30	35.88	1.55	23.20	0.13	6.46	1.16	−14.0
Asn487←Tyr83	34.65	2.11	3.28	0.02	1.54	0.91	7.4
Thr500→Asp355	30.32	1.87	8.08	0.04	2.85	0.90	8.0
Gln493→Glu35	28.44	2.21	3.68	0.02	1.46	0.98	1.1
Tyr505→Glu37	16.92	1.90	8.30	0.04	2.72	0.91	6.7
H−O−H→OH2		1.98	7.09	0.01	2.30	0.89	10.0

Note that this procedure intends to extract a representative sample from the simulations, thus, there are considerably more contacts not explicitly shown because they have lower occupancies or, because while having high occupancies, are not persistent in the three replicas. These HB contacts, whose bonding interactions are dissected below, are responsible for the attachment of the virus to epithelial cells in humans, initiating the infection process.

Quantum interactions

A summary of the quantum descriptors for the virus⋅⋅⋅cell interactions is listed in Table 1, the corresponding pictures are shown in Figures 4 and 5. Without exception, despite being weak organic acids, residue to residue hydrogen bonds are stronger than the archetypal HB in the water dimer, this is seen in the larger binding energies, smaller distances, orbital interaction energies, of comparable magnitudes, larger bond indices, and in the properties of the bond critical points. The LYS417→ASP30 is an exceptional case because it corresponds to an unusually strong, highly ionic [H2N−H]+⋅⋅⋅−[O2C] interaction, characterized as a salt bridge in a previous work.

Figure 4

Explicit virus⋅⋅⋅cell interactions. Top: interaction between the SARS‐CoV‐2 RBD(S) and the ACE2 receptor (in blue). The snapshot was randomly extracted from the late stages of one of the three MD replicas. All persistent virus⋅⋅⋅cell hydrogen bonds listed in Table 1 are explicitly highlighted in the frames (right). The noncovalent[ , ] virus⋅⋅⋅cell interaction surface is explicitly shown in green, including the water molecules. Note that both fragments contain glycosylated glycoproteins and that all fine structure is accounted for during the calculations; however, the glycans are not shown in this picture for clarity. Bottom: NBO donor→acceptor interactions responsible for the persistent hydrogen bonds.

Figure 5

Calculated QM/MM absorption spectra of the dimers and of the virus cell complex. Calculated QM/MM absorption spectra for the monomers (blue: amino acids belonging to the virus, red: amino acids belonging to the receptor) and for the bound structures (black). λ max, the wavelength [nm] of maximum absorption for the bound structure is explicitly included. All binding energies are in kcal/mol. The convoluted absorption spectrum for the complex at the bottom shows Δλ max=47 and 66 nm shifts from the isolated host and virus, respectively. Level: B3LYP/aug‐cc‐pVDZ/AMBER.

−

102

[%]

All hydrogen bonds are well‐characterized long‐range interactions. This is clearly seen in 1) There is never a formal σ orbital between the fragments, on the contrary, orbital interactions are always of the nO→ or nO→ form. In other words, in the NBO picture (Figure 4), all explicit virus⋅⋅⋅cell contacts are stabilized by charge transfer from one lone pair in an oxygen atom in the donor residue to an antibonding orbital in the acceptor residue. 2) All properties of the bond critical points support the same picture: small ρ(r c), small bond orders, virial ratios smaller than 1, and positive bond degree parameters. Again, the nO→ in LYS417→ASP30 is the exception, with all the calculated descriptors indicating a highly ionic contact. The noncovalent interactions (NCI) calculation for the interaction region uncovers a large discontinuous noncovalent wall separating RBD(S) from the receptor. Therefore, we characterize the virus⋅⋅⋅cell binding as due to a large number of noncovalent contacts between the two proteins, enhanced by the water molecules, acting in conjunction with the specific residue to residue hydrogen bonds.

Electronic absorption spectroscopy

We concentrate on simulating the spectra of the amino acids involved to investigate whether measurable changes in their spectral response occur upon binding, while residues that are not involved in the interaction can safely be viewed as changing very little as the two structures connect. For this purpose we propose a QM/MM approach where the amino acid pairs involved in the binding are treated quantum mechanically by means of density functional theory (DFT), while the rest of the protein environment is modelled classically, through the use of the AMBER force field. In this way, the electronic structure of the QM portion is influenced by its environment by means of an electrostatic embedding paradigm, where fixed charges are assigned to the MM atoms and directly affect the QM density and computed electronic excitations. Figure 5 shows that in all six cases, binding events set off drastic changes in the spectral response of the system, explicitly seen in red‐shifted absorptions. The red shift of the absorption bands is clearly visible in the convoluted spectrum and therefore provides unequivocal evidence of virus⋅⋅⋅cell bonding. The results listed here are quite encouraging and constitute an initial step that will hopefully motivate the design of experimental protocols to detect virus infection. However, it is clear that a number of details need to be worked out before practical applications can be devised. In particular, the potential interference of signals arising from functional groups in the same region of λ max, and the ability of the dimer model to accurately mimic physiological environments, should be addressed.

Virus Mutation as a Molecular Evolution Problem

The basics of molecular and viral evolution are well established and are worth summarizing in the context of this manuscript. Both processes are number games, that is, although the probabilities of individual specific mutations occurring are quite small (rates of mutation for single nucleotide sites in viruses are typically 10−6, literally one in a million), they occur nonetheless because of the random nature of amino acid sequencing errors proper of the replication of genetic material and because of the large number of reproduction events. In the absence of artificial or sexual selection, out of the large number of mutations that do occur, some are neutral, some are harmful, and some are beneficial to the evolving entity as a function of the environment demands. Due to the evolution pressure imposed by the environment, those individuals with beneficial mutations, which start with very small populations, have a larger probability of surviving their generation and of producing descendants into the next generation after yet another round of replication, thus progressively increasing their population until eventually, after several generations, the traits brought about by those changes become the dominant characteristic. There are usually diverse ways to survive the evolution pressure, leading to speciation. Because of the very fast reproduction rates, environment driven virus mutations help them quickly develop resistance against external agents designed to fight their infections. This ability of viruses to quickly mutate is among the most serious problems faced when developing vaccines and therapies, and is particularly true for SARS‐CoV‐2.

We argue that this view of evolution as driven by environment induced molecular responses at the virus/biomolecules scales helps explain many aspects of evolution that are difficult to rationalize otherwise, namely, 1) in most cases evolution is a highly localized process, 2) because of the large number of mutation possibilities, which occur no matter how small the individual probabilities, an increase in entropy of the universe is the ultimate factor driving evolution, 3) evolution is a deterministic process driven by cumulative random changes and 4) in that sense, life itself is a deterministic process that only requires a large increase in the entropy of the universe (in other words, a long time), such that it will emerge in local environments capable of sustaining it. See for example early arguments by Schrödinger stating the apparent macroscopic stability is due to the microscopic chaos resulting from random events, and invoking a net entropy gain by the universe due to the continual energy transformation despite the heavy entropy investment in maintaining highly organized living organisms. More recently, England has discussed the statistical physics of self‐replication.

In the context of this work, the previous discussion of molecular and virus evolution leads us to hypothesize that one of the key factors in the molecular evolution problem faced by the precursor SARS‐CoV virus on its way to mutating into SARS‐CoV‐2 was solved by favoring those changes in RBD(S) that lead to an improved ability to locate available sites for hydrogen bonding in the host cell, ability that is further enhanced by the slightly basic pH found in physiological environments. This improved hydrogen bonding capabilities may be achieved in a number of ways, for example, incorporating amino acids with more acidic protons, or incorporating larger amino acids whose hydrogen bonding regions are simply closer to the receptor, among others. We support the need for improved hydrogen bonding as the selective pressure in virus mutation hypothesis in the following evidence:

Table 1 shows that fragment to fragment contacts are all in the form of hydrogen bonds. For the specific case of the SARS‐CoV‐2 virus, in five of the six identified contacts, including the most persistent HBs, the residues in RBD(S) act as donors to the corresponding hydrogen bond

Besides a small sheet and a small helix (Figure 4), there is no secondary structure in RBD(S), thus, the receptor binding domain of the spike protein has a high structural flexibility which allows the virus to probe for available hydrogen bonding sites in the receptor, which in contrast has well defined secondary and tertiary structures in the interaction region

We obtained from the GenBank the sequences of amino acids for the precursor SARS‐CoV (ID AFR58742.1) and for the mutated SARS‐CoV‐2 (ID QHD43416.1) viruses. We compare below only the 196 amino acids in the RBD(S) and highlight in red the receptor binding motif (RBM). We also underline amino acid substitution in the mutated virus.

From these sequences, we point out that

There are a total of 49 substitutions (≈25 %) in the RBD(S) of SARS‐CoV‐2

Most of the mutations occur in the actual interaction region: the rate of substitution in the binding motif is quite higher: 36 substitutions that amount to ≈50 % in RBM

Estimating the overall acidity of large biomolecules is an open problem, there are in fact many available algorithms to calculate isoelectric points of peptides and proteins.[ , , , , ] Here, we take a pragmatic approach to determine relative acidities between the precursor SARS‐CoV and the mutated SARS‐CoV‐2 viruses. We took averages of the experimental isoelectric points (IEP) for the 49 amino acids involved in the mutation, that is, we calculated the average IPE in the replaced amino acids in RBD(S) and found 6.50 and 6.44 for the precursor and for the mutated viruses, respectively. We also calculated the same averages for the interaction motifs only and obtained 6.44 and 5.98 over the 36 mutations. Thus the mutated virus is collectively considerably more acidic in the interaction region, which improves its ability to donate protons to hydrogen bonds.

Summary and Conclusions

The most important contributions from this work may be summarized as follows:

We pinpoint the specific residue (in the virus) to residue (in the cell) interactions during the initial virus⋅⋅⋅cell binding.

We characterize the virus⋅⋅⋅cell molecular attachment as the result of a large number of noncovalent contacts between the receptor binding domain in the spike protein of the virus and the angiotensin converting enzyme 2 in the receptor expressed in epithelial human cells. This large surface of noncovalent interactions is enhanced by water molecules located in the interfragment region and acts in conjunction with an extended network of hydrogen bonds.

The need for improved hydrogen bonding is the selective pressure in virus mutation. Thus, one of the key factors in the molecular evolution problem faced by the precursor SARS‐CoV virus on its way to mutating into SARS‐CoV‐2 was solved by favoring those changes in RBD(S) that lead to an improved ability to locate available sites for hydrogen bonding in the host cell.

Mutated SARS‐CoV‐2 is more contagious than the precursor SAR‐CoV because it is better at finding hydrogen bonding sites in the receptor cell. More specifically, virus mutations have produced SARS‐CoV‐2, a virus that is more acidic in the interaction region than SARS‐CoV.

The calculated QM/MM absorption spectra show that as a result of virus⋅⋅⋅cell binding, absorption bands are shifted and spectral intensities are quenched, thus suggesting that significant changes occur at the electronic level as a consequence of binding.

Computational Methods

The starting point of our calculations was the complex between the receptor binding domain of the S protein, RBD(S), and the ACE2 receptor. Cartesian coordinates for the RBD(S)⋅⋅⋅ACE2 complex were taken from the protein data bank (PDB ID: 6LZG ) and then treated with CHARMM‐GUI, the graphical user interface of CHARMM[ , ] to include missing hydrogen atoms at pH 7.4, to ensure that all glycans are included, and to construct the force field. The entire system was enclosed by a truncated octahedral box such that the smallest atom⋅⋅⋅wall distance was set to 9 Å, then the available volume in the box was filled with TIP3P water molecules. NaCl molecules were added until a physiological 0.154 M concentration was attained and, finally, counterions were added to restore charge neutrality. This procedure lead to a system comprising a total of 171161 atoms, with 52 735 water molecules, 9 585 atoms in the receptor, and 3 034 atoms in RBD(S).

The system was subjected to a steepest descent energy minimization in order to correct for potential inconsistencies in atom coordinates that may arise during the procedure of randomly filling the available space with water, NaCl, and counterions. Once minimized, we ran triplicate all‐atom MD simulations under the conditions summarized in Table 2 and described next. First, there were three equilibration steps lasting a total of 0.625 ns with 1 fs time intervals, during which the structural constraints were progressively relaxed until finally being totally lifted. These structural constraints were imposed by harmonic constants that prevent deformation of the backbone (k bb), side chain (k sc), and dihedrals (k d). Then, the system underwent a production step lasting 600 ns with time intervals of 2 fs. For all MD runs, the Lennard–Jones potential was softened starting at 0.8 nm until eventually vanishing at 1.0 nm. Also, the cutoff radius for electrostatic interactions was set to 1.0 nm. All these simulations were conducted by using the CHARMM36m forcefield[ , ] as implemented in GROMACS 2019.4 at 310.15 K and 1 bar.

Table 2

Conditions for the MD simulations of the RBD(S) ACE2⋅⋅⋅complex. All times intervals (0‐600.625) in nanoseconds and all constants in kJ mol−1 nm−2. Harmonic k bb, k sc, k d constants prevent deformation of the backbone, sidechains, and dihedrals, respectively.

Parameter	0–0.125	0.125–0.325	0.325–0.625	0.625–600.625
kbb	400	400	0	0
ksc	40	40	0	0
kd	4	4	0	0
Ensemble	NVT	NPT	NPT	NPT
Thermostat	v‐rescale	v‐rescale	v‐rescale	Parrinello‐Raman
Barostat	–	Berendsen	Berendsen	Nose‐Hoover

–

The RBD(S)⋅⋅⋅ACE2 interaction energy (ΔG int) was estimated by the MM/PBSA method as implemented in GROMACS. Dielectric constants were set to ϵ solute=1, ϵ solvent=80 at the simulation temperature. In short:

Here, ΔE virus⋅⋅⋅cell is the gas‐phase energy of the RBD(S)⋅⋅⋅ACE2 complex, ΔG p, ΔG np are the solvation energies due to the polar and nonpolar interactions, respectively, and T ΔS is the entropy contribution. More precisely, ΔG p was computed under the Poisson–Boltzmann model, ΔG np was estimated using the solvent accessible surface area, and the entropy term was obtained from the model of Duan and co‐workers. To finally estimate ΔG int, we took 300 points in the 140–170 ns interval of one of the MD replicas.

It has been recently shown that randomly chosen configurations from late stages of MD simulations are adequate sources to obtain deep insight into interfragment bonding. Accordingly, aiming at understanding the fundamental forces driving the attachment of RBD(S) to host cells, virus⋅⋅⋅cell bonding interactions were dissected following these steps:

Persistent residue (in the virus) to residue (in the host cell) contacts during the 600 ns of the MD simulations were identified using the VMD program with a cutoff radius of 3.5 Å.

One frame was randomly chosen from the late stages of one MD run.

We extracted all extended interacting pairs in the chosen frame, kept them in the configurations they had in the interacting system (this is more accurate to understand the virus⋅⋅⋅cell bonding interactions than reoptimizing the isolated pairs), and

Computed accurate interaction energies using highly correlated domain based local pair‐natural orbital‐coupled cluster (DLPNO‐CCSD(T)/aug‐cc‐pVDZ) single‐point energy calculations[ , ] on the dimers and in the monomers. The ORCA suite of programs, version 4.0.1.2, was used to this end

Dissected the intermolecular interactions using the tools provided by the natural bond orbitals (NBO[ , , ] as implemented in NBO7.0 ) and by the quantum theory of atoms in molecules (QTAIM[ , , ] as implemented in AIMall )

Calculated QM/MM absorption spectra for the monomers and for the dimers. All TD‐DFT calculations were carried out using the B3LYP/aug‐cc‐pVDZ model chemistry[ , , , ] (tests using the dispersion corrected B3LYP‐D3, ωB97xD, functionals yielded essentially identical results). QM/MM Electrostatic Embedding was exploited, in which only the extended dimers were considered as the quantum region, and the rest of the system as the MM region, which was modeled by the Amber force field and by assigning to atom types the same charges used in the MD runs. A large number of excited states are needed to guarantee that both the intensities and shapes of the absorption spectra are accurately reproduced. Therefore, in this work the first 20 excited states were computed at the TD‐DFT/QM‐MM level in each case. Vertical excitations were shifted by −0.7 eV to account for the systematic error due to the choice of functional. This value was chosen in order to match the experimental absorption maximum for Tyrosine.[ , ] All Spectra were then convoluted with Gaussian line shapes with full width half maximum (FWHM) of 0.6 eV. All QM/MM calculations were carried out with Gaussian16

We isolated the interaction region by including everything within a 3.0 Å radius from the last atom at the end of each amino acid (1325 atoms in total) and calculated the interfragment non covalent interaction (NCI as implemented in NCIPLOT ) surface using the promolecular densities approximation.[ , ]

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.