Quantifying the adhesive strength between the SARS-CoV-2 S-proteins and human receptor and its effect in therapeutics.

The binding affinity and adhesive strength between the spike (S) glycoproteins of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), and the human angiotensin-converting enzyme 2 (ACE2) receptor is computed using molecular dynamics (MD) simulations. The calculations indicate that the binding affinity is [Formula: see text] [Formula: see text] with a maximum adhesive force of [Formula: see text] pN. Our analysis suggests that only 27 (13 in S-protein, 14 in ACE2) residues are active during the initial fusion process between the S-protein and ACE2 receptor. With these insights, we investigated the effect of possible therapeutics in the size and wrapping time of virus particles by reducing the binding energy. Our analysis indicates that this energy has to be reduced significantly, around 50% or more, to block SARS-CoV-2 particles with radius in the order of [Formula: see text] nm. Our study provides concise target residues and target binding energy reduction between S-proteins and receptors for the development of new therapeutics treatments for COVID-19 guided by computational design.

Introduction

In late 2019 a novel severe acute respiratory syndrome coronavirus (named SARS-CoV-2) was infecting people in China, causing severe pneumonia[1]. Within a few months from the first outbreak, the novel coronavirus created a global pandemic that forced the majority of the world’s population under lockdown. The novel SARS-CoV-2 virus keeps infecting and killing a large number of people around the globe. Thus, it is imperative to understand and develop therapies that can combat COVID19, the illness caused by the SARS-CoV-2 coronavirus.

SARS-CoV-2 belongs to the coronavirus genus[2] and they usually enter to host cells by attaching and fusing to the cell membrane[3]. Cell receptors diffuse across the membrane’s external surface to reach proximity with the virus’ proteins, binding together for their posterior fusion, promoting membrane bending and virus wrapping until final uptake. Coronaviruses affinity with cell receptors occurs via a transmembrane spike (S) glycoprotein forming homotrimers on the virus’ capsid[2,4]. The S-protein is made of two functional subunits (S1 and S2) responsible for fusion to the viral-receptor adhesion. Due to their critical role in SARS-CoV-2 infections, S-proteins are the common target for developing antibodies and therapeutics for COVID19.

Several works have thus far been focused on characterizing the S-protein and its trimeric structure using cryoEM techniques. For instance, Wrapp et al.[5] have provided a cryoEM structure in the prefusion conformation and have found two states, labeled as up and down, whereby the S1 subunit is exposed and retracted, respectively. The exposed region that links to the human receptor is known as the receptor-binding domain (RBD). At the same time, Lan et al.[4] have studied the RBD bounded to the ACE2 and provided a detailed description of the S1 subunit that compose the RBD and its link to the ACE2 receptor. In particular, they found that the S-protein links to the N-terminus helix of the ACE2 protein serving as an anchor point. Moreover, only a reduced number of residues, in total 20, were in close contact with the ACE2 terminus helix and even a smaller portion was within 0.4 nm from it[4]. This observation suggests that the adhesive forces arise through short-range interactions (e.g., Hydrogen bonds and salt bridges) between these two proteins. Shang et al.[6] provided a structural basis for receptor recognition of SARS-CoV-2. They found that in addition to the residues in the S1 subunit, many glycans generated links between the two proteins. Understanding the link between the S-protein and ACE2, and in particular, the RBD is key to tackling the pandemic caused by SARS-CoV-2.

The importance of the S-protein/ACE2 interface has motivated researchers to explore the phenomenon with both experimental and computational methodologies due to a growing interest in repurposing therapeutic to treat COVID19. However, testing the efficacy of these drugs is time consuming and expensive, pushing scientists to develop predictive models based on computational tools to reduce development time. For instance, Smith et al.[7] have scanned thousands of ligands with molecular dynamics simulations of the RBD, and have ranked these ligands based on their affinity. Other studies have focused their attention on quantifying the S-protein/ACE2 receptor’s formation energy, using a full trimeric model and/or a single S-protein/ACE2 receptor[8-11]. While these studies provide useful information on the compound’s formation energy, they failed in predicting realistic interaction energies that can be indirectly contrasted with experimental measurements. This shortcoming is because the adhesive interactions between the S-proteins and receptors are short-range, and it changes as a function of the separation length. Panda et al.[12] pursued a similar approach to benchmark drugs and antibodies for SARS-CoV-2. These studies focused their attention on the binding affinity of chemical compounds to reduce the formation energy between SARS-CoV-2 S-protein/ACE2 receptor. However, a quantitative evaluation of the binding affinity between the S-proteins and ACE2 receptors, the adhesive strength of this bond and the chemo-mechanical determinants controlling coronavirus uptake are still missing. This knowledge gap significantly limits the impact of the aforementioned investigations and underlines the importance of the proposed study.

In this work, we investigate the chemo-mechanical interaction between S-protein and ACE2 receptors, and the resulting implications on the mechanisms for virus uptake. We computed the binding affinity and adhesive force between S-protein and ACE2 receptors and analyzed the residues in contact during the bond-breaking process. Surprisingly, our results indicate that the residues in contact change as the two proteins were pulled apart, elucidating target points to develop new therapeutics. With these findings, we investigated SARS-CoV-2 uptake to predicted the effects of the binding affinity perturbations on the uptake kinematics.

Results and discussions

After performing umbrella sampling simulations on the S-proteins/ACE2 receptor configuration, we investigated the potential of mean force (PMF) evolution as a function of the pulling distance, i.e., the reaction coordinate () (see Fig. 1a for a schematic). Figure 1b shows the results for the single and full S-protein/ACE2 receptor configuration. Focusing our attention on the full trimeric protein, we observed that, initially, the evolution of the PMF shows a metastable and a global minimum between nm. These configurations were separated by a small barrier of . For nm, the PMF’s evolution shows almost a linear behavior with the reaction coordinate up to approximately nm, where the PMF reached . Thereafter, the PMF changes slopes and tends to plateau around at the end of the sampling, when nm. The maximum error in the measure of the PMF is approximately .

Figure 1

Evolution of the PMF, force, and residues with . (a) Schematic of the setup. (b,c) Evolution of the PMF and force as a function of for a single S-protein/ACE2 receptor (blue) and for the trimeric protein (red). Shaded area in (b) represents the error band for both simulations. The stars indicate the moment when detachment has happened. (d) View of the two anchor points between S-protein/ACE2 receptor when nm. (e) Interacting residues for nm and (f) nm. (g) List of residues that are active during the pulling simulation in the RBD. Letters in red, cyan background, bold and underlined correspond to residues in contact at nm, respectively. Red dots indicate residues identified experimentally by Lan et al.[4]. Residue K-417 has been omitted in the sequence for shortness.

The PMF’s change at around nm indicates that at this point, all van der Waals interactions are off between the two proteins, as shown with the change of slope in the plot. This was also confirmed by analyzing the residues in contact (below). Remarkably, the position where the change of slope happens is very close to the selected cutoff of the van der Waals interaction set up in our model. The remaining interactions appear due to long-range electrostatic forces that are in the model. We take the curve’s inflection point as the value of the adhesive strength that characterizes the link between the S-protein/ACE2 receptor. The result obtained with a single S-protein/ACE2 receptor leads to approximately the same free energy but a slightly different path. The binding affinity is estimated to be for the full trimeric model, and for the single S-protein/ACE2 receptor indicated with stars in the plot.

The force separation between the S-protein and ACE2 receptor is obtained from the relation , hence the slope of the PMF curve in Fig. 1b. The force evolution is shown in Fig.  1c for both models. We observed that the force builds up to a maximum of pN denoting the rupture force between the S-protein/ACE2 receptor bond, in the order of magnitude expected for adhesion in cells[13]. Thereafter, the force drops significantly due to lack of contact between the residues.

The binding energy between the S-proteins/ACE2 receptor can be used to compute the dissociation constant nM (see Supplementary information (SI) for calculation of the constant). Recent works have estimated the dissociation constant of the SARS-CoV-2 virus in experimental setups, obtaining values between nM to 4.674 nM[2,4]. The almost five-fold discrepancy range among previous measurements underlines the difficulty in obtaining accurate experimental data and, also, remarks the potential impact of our computational method. Considering that our model is limited to only a small portion of the real virus/cell receptor, the agreement between the simulations and experiments is remarkable and gives confidence in our computational approach.

Molecular analysis of the adhesive mechanics

We now focus our attention on the molecular interactions between S-protein and ACE2 receptor as a function of the reaction coordinate . The analysis of the interactions was performed with the full model at coordinates nm. Using the most typical cluster configurations (see “Methods”), we first obtained the interface residues between the two molecules and performed a contact analysis between that group. We disregarded all residues whose distance was more than 0.4 nm. As expected, the number of interactions decreased when increased.

First, we analyze the configuration with minimum free energy in our simulations, corresponding to nm, shown in Fig. 1d. The S-protein anchors from two locations; namely, the ends of the N-terminus helix in the ACE2 receptor and it could reach the helix on top if it, as shown in Fig. 1d. We observed that the interacting residues in the S-protein were located between positions 417-505 of the sequence, namely residues K-417, Y-449, Q-474, N-487, Y-489, G-496, T-500, G-502, Y-505, as shown in Fig. 1g—using a one-letter sequence—with red letters (see also Figure SI1, Table SI1 and Video SI1 in the SI). These residues linked to residues Q-24, D-30, E-37, Y-41, Q-42, Y-83, K-353, G-354, D-355, R-357 in the ACE2 receptor. When nm, ten residues were active—five were the same—(denoted with a cyan background in Figure 1g) with a graphical representation in Figure SI1. When nm, we observed interactions between the terminus helix and the one on top, as shown in Figure 1e involving Q-493, N-487, and Q-474 in the S-protein (bold letters in Fig. 1g), and Q-24, H-34, Y-83 in the ACE2 receptor. These residues were the most persistent ones, generating stronger links than other residues through the bond-breaking simulation. Thus, these residues can be targeted in new therapeutics strategies in COVID19. For nm, we found that only two residues were interacting, namely E-484 and Q-24 in the S-protein and ACE2 receptor, respectively (underlined in Fig. 1g). Figure 1f shows the links between these residues (see SI Video SI2). For nm, no contacts were found.

From the 27 residues in the RBD, our simulations indicate that 13 were active during the pulling simulation. In particular, we identified 13 unique residues that were active in the S-protein, and 14 in the ACE2 (see Table SI1). These residues are the same to the ones identified by Lan et al.[4] in their cryoEM analysis, with the exception of G-446, N-501 and Q-493. However, we did identify G-502 and Q-498, which are very close to the previously mentioned residues. This remarkable agreement gives confidence to our approach and simulations.

Uptake modeling and the effect of therapeutics

We now analyze the effect of the binding affinity in the endocytosis of the virus in cells. We recur to the chemo-mechanical model developed by Gao et al.[14]. The model considers the bending energy of the cell membrane, the release of chemical energy during the fusion of S-protein and receptor, the configurational entropy, and the ratio between receptor and S-protein density (see “Methods”). Also, the model needs specific parameters that are characteristic for each virus. We discuss these parameters—which are summarized in Table 1—for SARS-CoV-2 in “Methods”.

Table 1

Parameters used to model the endocytosis process in an infinite membrane. The reference temperature was taken as K. R is the radius of SARS-CoV-2 particles, is the density of spike proteins, D is the diffusivity of receptors, B is the bending modulus of lipid bilayer, is the binding affinity, and is the ration between receptor and S-protein density.

R [nm]	ξ_S[μm^-2]	D[μm^2·s^-1]	B[kCal·mol^-1]	e_RS[kCal]	ξ~=ξ_Rξ_S
30-70	2930	0.01	12.3	12.6	10^-1-10^-4

We found that the model predicts a minimum radius of nm with an optimum of nm at which the uptake time is minimum (around s, see Figure SI2) for . Particles below cannot be wrapped because the uptake is not energetically allowed. For smaller values of , we found that the minimum radius increases, in particular, for a the minimum and optimal radius are nm, and nm, respectively with a minimum wrapping time s, see Figure SI2. These trends are in agreement with Gao’s findings[14]. The predicted minimum radius of nm is in close agreement with experimental observations of SARS-CoV-2 particle size[1,15] where the minimum particle radius was 30 nm. This constraint could be implicit in the molecular architecture of the virus. However, our predictions suggest that viral particles that are smaller than nm cannot be uptaken, hence preventing their reproduction inside the host cell. Moreover, the predicted optimal radius of nm is close to the average particle size measured experimentally of nm[1,15].

Current attempts to treat COVID19 aim to repurposing therapeutics drugs and antibodies to bind between S-proteins and ACE2 receptors, thereby reducing [7,12]. Next, we provide an estimation of the reduction of needed to increment above the radius of SARS-CoV-2 particles (i.e., to stop particles from being uptaken). To this end, we modified the S-protein/receptor binding affinity , where is a reduction factor giving no affinity for and full affinity when . Thus, k represents the effectiveness of the treatment in reducing the binding affinity between S-proteins and ACE2 receptors.

Figure 2 shows the relation between and k for various . The green region identifies all particles radii for which virus uptake is always permitted independently of any reduction in binding affinity (k). The red region identifies particles whose radii are not permitted to be uptaken due to excessive bending energy. Between these regions, virus uptake is controlled by . For instance, if the binding affinity between S-proteins and ACE2 receptors is reduced by 30% (k=0.7), the minimum radius of particles that can be uptaken increases by 21.5% for ( nm), and by 32% for ( nm, blue triangles in Figure 2). Another important aspect is that the wrapping time required for final uptake increases about 50 to 100% for these cases, respectively (see Figure SI3). This increment in the time needed for final uptake might be critical in some cases since it might give the immune system extra time to combat the infection, thus illustrating the effect of the treatment.

Figure 2

Effect of reduced binding affinity (factor) in the size of particles that can be wrapped to the cell membrane. Four receptor densities were used, namely (solid line), (dashed line), (dashed-dotted line), (dotted line). Prohibited particles are in the red zone, while favorable particles are shown with in green. In between these zones, particle will be blocked depending on the relative receptor/S-proteins density .

If we reduce the binding affinity by 50% (), the minimum radius increases by 47% when ( nm), and by 84% when ( nm, blue stars in Fig. 2). The wrapping time also increases by 150% to 425% with respect to the case when (see Figure SI3) and allowing additional time to stop infection in cells. Finally, our work provides a rough estimation of how much binding affinity has to be reduced to provide effective treatments. Our analysis indicates that a 50% reduction in might stop the adhesion of SARS-CoV-2 particles with nm. Moreover, we predict that smaller particles ( nm) are more suitable to be blocked in comparison with larger ones, based on bending energy analysis.

Conclusions

Our study reports the first attempt (to the author’s best knowledge) to evaluate the binding affinity and bond-breaking force between SARS-CoV-2 spike proteins and ACE2 receptors via computational analysis using an all-atom MD model. Our estimation of gives a dissociation constant of nM, which is in close agreement with experimental measurements ( nM). Furthermore, our analysis shows that only a reduced fraction (about 13) of the residues in the RBD interact during the protein adhesion. These residues play a critical role in the adhesion of the S-protein/ACE2 receptors, and can be used as a target for therapeutic strategies to prevent virus uptake in human cells. Finally, using a mechanistic model to analyze virus uptake, we concluded that S-proteins/ACE2 receptor’s binding has to be reduced at least by 50% to block the uptake of SARS-CoV-2 particles having radius between 30 to 70 nm. However, according to our analysis, the effectiveness of these strategies is strongly dependent on particle size and receptor density. Hence, such treatments might be more effective in blocking only a portion of the SARS-CoV-2 particles leaving others unaffected.

Methodology

Molecular dynamics simulations, Umbrella sampling and characterization

Molecular dynamics simulations were performed with the GROMACS software[16-18]. The molecular geometry was taken from different sources, including the protein data base models 6LZG for a single S-protein[19], and the PDB file 6VYB full trimeric model[2]. In addition—since we started our simulations before the PDB models were available—we used the model provided by Smith et al.[7] using the sequences available online (NCBI Reference: YP009724390.1) for the SARS-CoV S-protein’s crystal structure and the ACE2 receptor was generated using the PDB 2AJF file. While the actual numbers vary slightly, the trends are the same regardless the geometry. The model was loaded into GROMACS, where it was solvated in water using the TIP3P model to achieve a density of approximately . In order to allow for sufficient space for the pulling simulations, we generated computational cells with more than 1 nm between the proteins and end of the cells, and sufficient space on top to perform the pulling simulations. The biggest cell size has dimensions of (nm). After adding the solvent, the system had a non-zero charge and sodium (Na) ions were added as needed to equilibrate in all samples. All interatomic forces were computed with the CHARMM force-field[20]. The biggest cell had 648,265 atoms including proteins and solvent.

The solvated system was initially subjected to an energy minimization using a non-linear conjugate gradient. The forces were minimized with a convergence criterion of 1000 . After the system was relaxed, it was subjected to an NPT ensemble for 1 ns using an initial temperature of K imposed with a Berendsen thermostat[21]. Pressure was controlled to 1 bar using a Berendsen barostat. The timestep was set to fs. For all simulations carried in this work, short-range interactions were treated with a smooth force-switch cutoff of nm, and long-range electrostatics were treated using the Particle-Mesh-Ewald (PME) formalism, implemented in GROMACS[22]. Hydrogen-bonds were restrained with the LINCS algorithm[23].

In order to compute the binding affinity between the S-proteins/ACE2 receptor, we used a combination of pulling simulations with umbrella sampling. The initial configuration of the S-proteins/ACE2 was subjected to a pulling simulation to generate the necessary configurations to perform an umbrella sampling. The pulling simulation was performed with an optimized spring constant of . This optimized spring constant was obtained by performing multiple umbrella samplings on the single S-protein-ACE2 receptor with spring constants in the range of , and optimizing the value using a quadratic fitting. Values obtained for different spring constants are shown in the SI. The pulling rate was set to along the z-direction and sampling simulations were run for a total time of ns. The pulling direction was set such that the S-protein and ACE2 receptor were pulled apart from each other. The configurations generated along the pulling simulations were systematically used to generate trajectories for the umbrella sampling, described below.

Umbrella simulations were performed for configurations separated about nm from the reference configuration. Each configuration was constrained with a spring constant of and run for ns. These simulations provided enough sampling to obtain the potential of mean force along the reaction path to display the evolution of the free energy of the system. The PMF was then estimated using the Weighted Histogram Analysis Method (WHAM)[24] using a bootstrap analysis to estimate the uncertainty in the PMF. We used 100 different measures, using 200 binning spaces along the reaction coordinate . We performed several simulations with the backbone of the proteins fixed and without any fix conditions. We found that simulations without fixing the backbone produce excessive elasticity in the proteins and lead to higher free energies.

Analysis of the atomistic configurations was performed with GROMACS cluster analysis tool. We scanned the configurations with a root mean squared displacement between the range of nm[25]. The cluster analysis yielded between three to six cluster for the analyzed configurations. In all cases shown, the most populated cluster was used when analyzing the configurations. The configurations were then analyzed with the software Pymol.

Mechanistic model of endocytosis

Gao et al.[14] developed a model considering a spherical particle being attached to an infinite membrane. The fusion of the particle with the membrane is driven due a release of the binding energy-computed above—when the S-protein and receptors are linked. It is assumed that the membrane has an equilibrium concentration of receptors, , and when the particle attaches, this concentration changes with time, e.g., , s being the arc-length. In particular, when a particle is attached to the membrane, the density of receptors, , matches to the density of spike proteins, , and far away tends to the equilibrium concentration. Considering the bending energy of a lipid bilayer—characterized through its bending modulus B and curvature , and the binding energy between S-proteins and receptors, () one can write down the following free-energy for the endocytosis process asIn Eq. (1), is the thermodynamic factor with and T denoting the Boltzmann’s constant and the absolute temperature, respectively. By requiring that the rate of free energy reduction gained in the wrapping process exactly balance the rate of energy dissipation consuming during the transport, Gao et al. found that there exists an optimal wrapping radius of the particles and a minimum radius below the particle cannot be wrapped. The wrapping time can be found aswhere R is the radius of the particle, is the speed factor (), and D is the diffusivity of the receptors in the membrane. The speed factor is found by solving the rate equationwithThe parameter defines the ratio between equilibrium receptor density in the membrane and the S-porteins in the virus particle. In the above expression, is the exponent integral function defined asThe model predicts a minimum radius for spherical particle given byand an optimal radius that is determined numerically. Gao et al. determined the wrapping time as a function of the particle radius as well as the optimal particle radius numerically.

Parameters determination

SARS-CoV-2 virus particles of sizes between nm have been reported[1,15]. This range indicates a wide range of particles, with an average size of nm. Microscopic images indicate that around spike proteins in the circumference of the virus. An elementary analysis indicates that the density of spike proteins in the novel SARS-CoV-2 virus must vary between (about spike proteins in the surface) when the average radius is taken. The computed density values are in close agreement to other coronaviruses[26]. The bending stiffness of lipid bilayers ranges between ()[27-29]. Here, following previous works, we adopt an average value of 12.3 (20 )[14,30].

The density of ACE2 receptors on the cell membrane, at thermodynamic equilibrium, is from which we compute the dimensionless ratio . Given the limited data available on the density of receptors on human cells is difficult to estimate, in particular because receptor density varies across cell types and depends on the specific receptor. Moreover, no specific data is available (to the author’s knowledge) about the density of ACE2 in epithelial cells in human lungs, the target of SARS-CoV-2.

Chen et al.[31] measured a density of 480–640 for receptor of various species, while Damioli et al.[32] measured a density of for VEGFR2 receptors. Based on these measurements, we estimated to vary between and adopt these values in our simulations. The receptor diffusivity was taken as an average value for most cells[14,30,32].

Supplementary file1

Supplementary file2

Supplementary file3