The binding of heparin to spike glycoprotein inhibits SARS-CoV-2 infection by three mechanisms.

Heparin, a naturally occurring glycosaminoglycan, has been found to have antiviral activity against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the causative virus of COVID-19. To elucidate the mechanistic basis for the antiviral activity of heparin, we investigated the binding of heparin to the SARS-CoV-2 spike glycoprotein by means of sliding window docking, molecular dynamics simulations, and biochemical assays. Our simulations show that heparin binds at long, positively charged patches on the spike glycoprotein, thereby masking basic residues of both the receptor-binding domain (RBD) and the multifunctional S1/S2 site. Biochemical experiments corroborated the simulation results, showing that heparin inhibits the furin-mediated cleavage of spike by binding to the S1/S2 site. Our simulations showed that heparin can act on the hinge region responsible for motion of the RBD between the inactive closed and active open conformations of the spike glycoprotein. In simulations of the closed spike homotrimer, heparin binds the RBD and the N-terminal domain of two adjacent spike subunits and hinders opening. In simulations of open spike conformations, heparin induces stabilization of the hinge region and a change in RBD motion. Our results indicate that heparin can inhibit SARS-CoV-2 infection by three mechanisms: by allosterically hindering binding to the host cell receptor, by directly competing with binding to host heparan sulfate proteoglycan coreceptors, and by preventing spike cleavage by furin. Furthermore, these simulations provide insights into how host heparan sulfate proteoglycans can facilitate viral infection. Our results will aid the rational optimization of heparin derivatives for SARS-CoV-2 antiviral therapy.

In the last year and a half, the COVID-19 pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has adversely affected the lives of all people around the world. SARS-CoV-2 is a lipid-enveloped positive-sense RNA virus belonging to the Coronaviridae family (1, 2, 3). Among the SARS-CoV-2 proteins, the spike S glycoprotein (spike) is highly conserved in the Coronaviridae family (76% and 96% amino acid residue sequence similarity with SARS-CoV-1 and BatCoV-RaTG13, respectively (1)). It possesses 22 N-linked glycans that likely contribute to protein stability and immune evasion of the virus (4, 5, 6). The spike mediates the entry of the virus into human cells by binding to the angiotensin-converting enzyme 2 (ACE2) receptor (7). The prefusion spike is exposed on the virion surface as a homotrimer. Each spike subunit is composed of two domains, S1 and S2, connected by the S1/S2 junction, and involved in virus attachment and fusion, respectively (7). The S1/S2 junction is a novel feature of the SARS-CoV-2 spike that consists of a multibasic sequence that is cleaved by the host furin protease (1, 3, 7, 8). This cleavage activates the formation of the postfusion conformation of spike, which is necessary for SARS-CoV-2 entry into host cells (7, 8). Intriguingly, the S1/S2 junction is part of a recently identified superantigen motif on the SARS-CoV-2 spike that has been proposed to interact with T-cell receptors and thereby elicit a hyperinflammatory response (9).

Heparan sulfate proteoglycans (HSPGs) on the surface of human host cells are composed of a core protein that harbors multiple polysulfated glycosaminoglycan chains with a similar structure to heparin (10). Experimental evidence suggests an essential role for HSPGs as coreceptors that, by binding spike, favor SARS-CoV-2 attachment to human cells (11, 12, 13). Although heparin is currently used to treat COVID-19 patients because of its strong anticoagulant activity (14), it has been shown to exert antiviral activity, primarily in its unfractionated state (up to 15 kDa corresponding to ∼50 monosaccharides), likely due to its ability to compete with HSPGs for binding to the spike (11, 12, 13, 15, 16).

The binding of spike to heparin or HSPGs is mediated by clusters of positively charged amino acid residues in the spike (hereafter called heparin-binding domains, HBDs) and negatively charged sulfate groups on the polysaccharide chains of heparin or the HSPGs (10, 11, 12). To date, three putative HBDs have been identified in the spike sequence: in the receptor-binding domain (RBD), at the S2′ TMPRSS2 cleavage site (17) (both present also in SARS-CoV-1 and at the novel S1/S2 furin cleavage site (12)) (Fig. 1A).

Figure 1

Modeling of spike–heparin interactions.A, schematic diagram of the sequence of the SARS-CoV-2 spike (S) glycoprotein head (residues 1–1140), which is composed of S1 and S2 subunits. The boxes along the sequence (Uniprot P0DTC2) show the positions of the main protein domains: N-terminal domain (16–291), receptor-binding domain (RBD) (330–530), S1/S2 site (681–686), fusion peptide (788–806), S2′ domain (810–816), central helix (987–1034), and connection domain (1080–1135). Three putative heparin-binding domains (HBDs) are indicated by blue boxes with their corresponding sequences. The S1/S2-HBD is unique to SARS-CoV-2. B, side views of the models of the spike head homotrimer before the MD simulations in closed (left) and open (right) conformations with three heparin chains bound. Only one complete heparin chain can be seen from this view. The SA, SB, and SC subunits are shown as surfaces in blue, pink, and magenta, respectively. The dark blue surface patches correspond to the putative HBDs. N-glycans covalently attached to the spike are shown in line representation, colored according to the subunit to which they are attached. The 31mer heparin chains that span from the S1/S2-HBD to the RBD-HBD are shown in stick representation, colored by element with cyan carbons. SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

While the mechanism underlying HSPG-mediated internalization and the antiviral activity of heparin are well understood for many HSPG-dependent viruses (18), it has only recently begun to be investigated for SARS-CoV-2. Indeed, among the various coronaviruses, the SARS-CoV-2 spike is unique in possessing three different putative HBDs that overlap with motifs that have distinct functions (RBD, TMPRSS2, and S1/S2 sites), indicating that spike–HSPG or spike–heparin interactions may have multifaceted effects on SARS-CoV-2 infection. In this context, we here report the results of molecular modeling and simulation along with biochemical assays aimed at revealing how heparin exerts its antiviral effects on SARS-CoV-2 and how HSPGs can, on the other hand, facilitate virus infection. In particular, we consider the mechanistic consequences of the sequence features that are specific to the SARS-CoV-2 spike for the effects of heparin and HSPGs on viral infection and inflammation.

Results

Long mostly basic patches on the spike glycoprotein accommodate heparin polysaccharide chains

We considered the inactive closed and active open prefusion conformations of the homotrimeric SARS-CoV-2 spike head, which consists of three subunits: SA, SB, and SC. In the inactive closed conformation, the ACE2 receptor-binding face on the RBD is not accessible in any of the three subunits, which we refer to as downsubunits. In the active open conformation, the ACE2-binding site of the RBD is accessible only in the SC subunit, which we refer to as the upsubunit. To investigate the binding of heparin to spike, we modeled five systems with the spike: (i) in the closed conformation; (ii) in the closed conformation with a single heparin chain bound; (iii) in the closed conformation with three heparin chains bound; (iv) in the open conformation; and (v) in the open conformation with three heparin chains bound.

The model of the SARS-CoV-2 spike head (Fig. 1B) was based on the SwissModel model (available at https://swissmodel.expasy.org/repository/species/2697049) to which we covalently attached 18 N-glycans per subunit (4). To identify continuous positively charged paths on the protein surface at which the long heparin chains could bind, the electrostatic potential of the spike head protein was computed for both the closed and open conformations. Analysis of the electrostatic potential suggests that a long polyanionic heparin chain can follow similar paths on the surfaces of the two conformations of the spike head that differ only in the interactions with the RBD-HBD (Fig. 1B). We used the incremental docking and sliding window method (19) to model complexes of the spike head with heparin chains of 31 monosaccharides (31mer). The docked heparin chains run along mostly basic patches, passing partially through surface grooves, from the S1/S2 basic motif (R682, R683, R685), via the channel between the N-terminal domain of the same spike subunit and the RBD-HBD of an adjacent spike subunit (R346, R355, K356, R357) to the top of the spike head. These models are consistent with the ability of heparin to bind both active open and inactive closed conformations of spike, as shown by Clausen and coworkers using the S ectodomain mutated to lock the RBDs into a closed or open conformation (11).

Due to the structural similarity between heparin and the polysaccharide chains of HSPGs, it is expected that the latter also bind along the basic paths identified on the trimeric spike head. By burying these basic regions, heparin can be expected to hinder the binding of HSPGs, reducing the amount of SARS-CoV-2 that can tether to the cell surface, thus decreasing the binding to the ACE2 receptor and hence infection. This is one of the mechanistic explanations that the models provide for the experimental data demonstrating that heparin can diminish SARS-CoV-2 infection (11, 12, 13, 15, 16). To investigate further mechanisms by which heparin can exert antiviral activity against SARS-CoV-2, we proceeded to carrying out molecular dynamics (MD) simulations of the modeled systems.

Heparin stabilizes the spike glycoprotein in the closed conformation

The five modeled systems were each simulated in aqueous solution in four replica all-atom MD simulations, each of 1 μs duration (Movies S1–S3). Representative structures obtained after clustering analysis of the spike in closed and open conformations with three heparin chains bound show that the heparin chains maintain their alignment along the long positively charged surface patches on the spike during the simulations (Fig. 2A). The spike structures reached convergence within ∼2 to 400 ns in all simulations, as shown by the root mean square deviation (RMSD) relative to each starting structure (Figs. 2B and S1). Moreover, when in complex with heparin, the spike showed an approximately 1 Å lower RMSD, indicating that the binding to heparin stabilized the homotrimer structure (Figs. 2B and S1). Further evidence of the structural stabilization was provided by the root mean square fluctuations (RMSF), which tended to lower values for the RBD, the hinge region associated with the opening of the RBD (see next section for further details), and the S1/S2 site when bound to heparin (Fig. S2, Tables S1–S3). Accordingly, the RMSD of heparin indicates an induced fit along the trajectory (Figs. 2B and S1). The complementarity of the interactions and the change in the width of the basic groove in both the closed and open spike structures bound to heparin suggest an induced fit upon binding that results in the well-defined partially grooved basic path where heparin lies (Fig. 2A).

Figure 2

Stability of the spike–heparin complexes.A, representative closed (left) and open (right) structures obtained after MD simulation of spike bound to three heparin chains are displayed as molecular surfaces with electrostatic potential mapped onto them to show the partially grooved positively charged path occupied by a heparin 31mer. Heparin is shown in stick representation colored by element with cyan carbons. The yellow dashed lines show H-bonds between the spike and heparin. The insets for the closed and open conformations highlight the H-bonding interactions between heparin and the residues in the RBD (T345, R346, N354, R355, N360) and S1/S2 (R682, R683, R685) HBDs shown in stick representation with carbons colored according to the subunits to which they are bound. B, structural convergence of the spike and heparin chains in the five simulated systems. RMSD versus time is shown for one replica MD simulation for each modeled system. Corresponding plots for all replica MD simulations are shown in Fig. S1. HBD, heparin-binding domain; MD, molecular dynamics; RBD, receptor-binding domain.

Overall, the H-bond and interaction fingerprint analysis (20) of all the MD simulations show that each heparin chain maintains stable interactions with two adjacent spike subunits in both the closed and open conformations of spike (Figs. 2A and S3, Tables S4–S6). Heparin binds through H-bonding interactions (with over 90% occupancy throughout all the MD trajectories) to the basic residues of the S1/S2-HBD in the first subunit and of the RBD-HBD of the second subunit of both the open and closed spike conformations. However, along the heparin path, additional, less specific binding regions, which differ between the open and the closed conformations, were identified in over 75% of the MD frames, and these can further stabilize the complexes.

In the closed spike model, the binding of either one or three heparin chains hinders the opening of spike by stabilizing the closed conformation through the simultaneous binding of the RBD of one subunit (residues T345, R346, N354, R355, S359, N360, N450) and the N-terminal domain of the adjacent subunit (residues N165, C166, T167, E169, V171, Q173, F220, N280, N282, T284, T286) (Figs. 2A and S3, Tables S4 and S6). In the closed conformation of spike, the RBDs are essentially strapped into the downorientation by heparin, which prevents spike activation to the open conformation. Finally, heparin binding is also stabilized by hydrogen-bond interactions adjacent to the multibasic S1/S2 site (residues N606, S686, S689, S690) (Fig. S3, Tables S1 and S3).

The simulations of spike in the open state do not show any tendency for heparin to induce closure of the spike. However, through the simultaneous binding to the up-RBD of subunit SC (residues T345, R346, N354, R355, R357, N360) and the N-terminal domain of the adjacent subunit (residues R34, T167, E169, Q173, L176, R190, H207, T208, F220, S221), heparin induces a change in the orientation of the RBD in the upsubunit during the simulations (Figs. 2A and S3, Tables S5–S6 and next section for further details). As for the closed models, some polar residues near to the multibasic S1/S2 site (N606, Y674, S686, Q689, S690) permanently interact with the heparin chain in the open state.

Overall, comparing the localization of the known emerging SARS-CoV-2 spike glycoprotein variants (Uniprot P0DTC2, on 29.09.2021) with the binding residues involved in the interaction with heparin, our data indicate that the natural variants do not directly negatively affect the binding of heparin (Fig. S4). However, we cannot exclude that specific mutations, such as P681R near to the S1/S2 basic domain, could facilitate and increase the affinity for heparin (and thus, also for HSPGs) of the spike glycoprotein.

Notably, the simulations of both open and closed forms of spike with heparin bound suggest an aspecific modulatory effect of N-glycans on the binding between spike and heparin, see Figure S5 and Movies S1–S3.

Heparin masks the spike S1/S2 site

To assess the ability of heparin to mask the novel furin cleavage site (7, 8), which has also been identified as a novel putative HBD (12) and as part of a superantigenic sequence (9) in the spike of SARS-CoV-2, H-bond formation to the site was monitored along the trajectories, and surface exposure of the S1/S2 site to the solvent was analyzed for representative structures after clustering (Figs. 3 and S6, Tables S4–S7). In the initial modeled complexes, the first monosaccharides of the heparin chain interact with the basic residues of the S1/S2-HBD. Salt links with R682, R683, and R685 are maintained along all the trajectories (>90% occupancy), indicating strong interactions of the S1/S2-HBD with heparin (Fig. 3A, Tables S4–S6). The calculated solvent-accessible surface area (SASA) of this multibasic site shows the persistence of the interactions during the simulations of the closed spike bound to heparin (Figs. 3B and S6). Heparin approximately halves the surface exposed in the closed models by binding directly to the S1/S2-HBD, and the three heparin chains simultaneously mask all three multibasic sites on the closed spike trimer. In the open conformation, the reduced SASA in the presence of heparin indicates a significant shielding effect, primarily at the S1/S2-HBD of the upsubunit, but with the heparin chains also maintaining a lower shielding level at the S1/S2-HBDs of the downsubunits. This difference could be due to a more favorable arrangement of the basic patches for heparin binding between the up- and downsubunits compared with that between two downsubunits.

Figure 3

Interaction of heparin with the spike S1/S2-HBD showing shielding of the furin cleavage site by heparin in closed and open spike conformations and inhibition of furin cleavage by heparin.A, the spike S1/S2 site in closed (left) and open (middle, right) representative conformations is shown with and without heparin. The molecular surfaces of the SC and SB subunits are shown in magenta and gray, respectively. Basic residues of the S1/S2-HBDare shown in blue and labeled. N-glycans are shown in line representation colored according to the corresponding subunit. Heparin is shown in stick representation colored by element with cyan carbons. B, SASA of the S1/S2-HBD plotted as a function of time for one trajectory for each simulated system: closed, closed with one heparin chain, closed with three heparin chains, open and open with three heparin chains. Corresponding plots for the replica trajectories are shown in Fig. S6. C, average SASA of the S1/S2-HBD computed for the cluster representatives in the MD simulations for different factors contributing to the burial of the S1/S2 site surface. Standard deviations are reported in Table S7. D–E, effect of heparin on the cleavage of spike by furin. D, wells coated with a peptide fragment containing the S1/S2 basic site of spike were: left untreated (CTRL), incubated with 13.6 kDa heparin (10 μM) alone or exposed to furin (25 ng/well) in the absence or presence of heparin (10 μM). In addition, the substrate-immobilized spike fragment was preincubated with heparin (100 μM) for 10 min, washed three times with PBS or PBS containing 2 M NaCl (salt), and then exposed to furin. E, wells coated with the spike fragment were incubated with furin in the absence or presence of 13.6 KDa unfractionated heparin (UFH), 4.0 kDa low-molecular-weight heparin (LMWH) or unsulfated K5 polysaccharide (K5) at the indicated concentrations. At the end of the incubations, spike cleavage was measured as described in the Experimental procedures section. Each optical density (O.D.) measurement was performed at 450 nm and is reported as the mean ± standard deviation of three to ten repetitions. Single measurements for each experimental condition are given in panel D (as single points •) and in Tables S8 and S9. ∗p < 0.001, one-way ANOVA. HBD, heparin-binding domain; SASA, solvent-accessible surface area.

To assess the shielding effect of heparin relative to that of spike glycosylation, the SASA of the S1/S2-HBD calculated for the representative clusters was decomposed into the area exposed without consideration of the N-glycans and heparin sugars, the area exposed accounting for the N-glycans, and the area exposed accounting for both the N-glycans and heparin (Fig. 3C and Table S7). In agreement with the previous analysis, these calculations indicate that heparin directly binds to the S1/S2-HBDs, halving the exposed surfaces in both the closed and open conformations. Again, for the open conformation, the shielding effect exerted by heparin on the S1/S2-HBDs of the two downsubunits is lower than for the same site of the upsubunit. Moreover, the decomposition shows that the N-glycans of spike make little contribution to the shielding of the multibasic sites.

In summary, when comparing the binding of one or three heparin chains to the spike homotrimer, all the data point to the ability of heparin to occupy and shield the S1/S2 site without a significant shielding contribution of the spike N-glycans. In addition, our simulations demonstrate a key role of the S1/S2 site in the binding of heparin and indicate its involvement in binding to HSPGs. The surface area analysis shows that the S1/S2-HBD is exposed to the solvent and available to interact with heparin or HSPGs, supporting the hypothesis that this HBD, which is specific to SARS-CoV-2, may contribute to the increased affinity of SARS-CoV-2 spike for heparin and HSPGs compared with previous coronavirus strains.

In conclusion, heparin could exert its antiviral activity by masking the S1/S2 site, preventing cleavage by furin, and hindering the activation of the postfusion conformation of the spike glycoprotein, which is required for the fusion with the cell membrane and the subsequent internalization of the virus. Interestingly, we observed that, in the closed spike model, a single heparin chain occupies the S1/S2-HBD of only one subunit. This suggests that, to exert its full antiviral activity, heparin needs to be administered at doses high enough to saturate all three S1/S2-HBDs of the spike trimer; otherwise, some of them will still be available for tethering to the HSPGs on the host cell surface.

The direct binding of heparin to the S1/S2 basic site prevents spike cleavage by furin

Experiments have demonstrated that the binding of heparin to spike occurs via the basic domain present in the RBD region (13, 16, 21, 22). Our simulations additionally show a direct binding of heparin to the S1/S2 multibasic site. To confirm this hypothesis arising from our simulations, we measured the effect of heparin on the cleavage of a spike fragment containing the S1/S2 basic site. As shown in Figure 3D, heparin effectively inhibits spike cleavage by furin. The effect is specific since heparin alone does not interfere with the colorimetric assays. To assess whether the inhibition of furin cleavage is dependent on the binding of heparin to the S1/S2 basic site, heparin was preincubated with a substrate-immobilized spike fragment, unbound heparin was removed by PBS wash, and then furin was added. In these experimental conditions, heparin prevents spike cleavage by furin, indicating that it physically associates with the immobilized spike fragment. On the other hand, a 2 M NaCl wash, which is known to disrupt heparin/protein interactions (23, 24), completely detached heparin from the spike fragment, restoring cleavage by furin (Fig. 3D). The capacity of heparin to bind spike and to inhibit SARS-CoV-2 infection is dependent on its length and degree of sulfation (15, 16). We therefore compared the inhibitory potency of unfractionated heparin (13.6 kDa, ∼ 46 monosaccharide units), a low-molecular-weight heparin (4.0 kDa, ∼ 13 monosaccharide units), and unsulfated K5 polysaccharide (30 kDa, ∼150 monosaccharides). As shown in Figure 3E, unfractionated heparin exhibits a higher furin inhibitory potency than the low-molecular-weight heparin, whereas the unsulfated K5 polysaccharide does not prevent furin cleavage. These results demonstrate that, for binding to the spike fragment and protecting it from furin cleavage, the length of the heparin chain is of importance and the presence of sulfated groups is necessary. Importantly, very long heparin chains can block a functional domain in the absence of a direct, specific interaction, simply by making contact with basic residues outside the functional domain itself and thereby masking it by “aspecific hindrance.” This latter possibility is here ruled out by the use of a very short spike fragment that demonstrates how the inhibition of furin cleavage by heparin is indeed due to the specific interaction of the GAG with the 681RRAR684 residues in the S1/S2 site.

These experiments corroborate the hypothesis that heparin exerts its antiviral activity in a length- and dose-dependent manner by masking the S1/S2 site on spike, thereby inhibiting furin cleavage and, consequently, the formation of the postfusion conformation. These results also suggest that, by binding to the S1/S2 site, heparin could affect the superantigenic character of the spike by preventing its binding to TCR cells (9).

Heparin allosterically affects the hinge region of the RBD and directly interacts with the basic residues of the RBD-HBD

To assess the effect of heparin binding on the RBD, we analyzed the stability of the hinge region associated with the activation of the RBD, the exposure of the up-RBD (on subunit SC) along the trajectory, and the shielding by heparin of the residues of the RBD involved in the interaction with ACE2, hereafter referred to as the receptor-binding motif (RBm) (25).

From a comparison between the crystal structures of the spike in closed and open conformations and the RMSFs in the simulations (data not shown), we identified residues 527-PKK-529 as the hinge region responsible for the conformational change that induces the opening of the spike protein. Importantly, no direct interactions were observed to occur between these residues and heparin (Tables S4–S6), prompting us to investigate possible allosteric effects induced in this region by the binding of heparin to spike. For this purpose, we calculated the RMSD of the hinge region of each subunit along the trajectory and performed dihedral principal component analysis (dPCA) for this region. As shown in Figures 4, A and B and S7 and S8, the closed spike shows structural stability and limited structural variability of the three hinge regions along all the simulations. The closed conformation with one heparin bound shows structural stabilization of the two subunits that are directly involved in the interaction with heparin. However, the SC subunit, which does not directly interact with heparin, shows a higher RMSD (Figs. 4A and S7), a higher RMSF (Fig. S2, Table S2), and more structural variability (Figs. 4B and S8), suggesting a compensatory effect due to the movement of the hinge region to maintain the closed conformation. Interestingly, RMSD plots and dPCA of the hinge region of the closed-SC subunit in the complex with three heparins bound show a restored stability of the hinge region, indicating the ability of heparin to stabilize the closed conformation of spike and hinder the opening of the RBD. For the open spike, the binding of three heparin chains results in the subunits SA and SB with down-RBDs having structurally stable hinge regions, whereas the hinge region of the SC subunit, which has the up-RBD, has an increased RMSD (Figs. 4A and S7) and altered sampling of the conformational space (Figs. 4B and S8), indicating an intriguing effect of heparin on the opening of the spike conformation, which could affect host-cell receptor binding.

Figure 4

Changes in the conformational flexibility of the spike hinge region (residues 527–529) upon heparin binding and the effect of heparin on the motion of the spike RBD and the accessibility of the receptor-binding motif (RBm).A, RMSD vs. time and (B) dPCA of the hinge region for a trajectory for each simulated system showing that heparin-binding results in structural stabilization of the hinge region of the closed spike and changes the motion of the hinge region in the open spike. The SASA of the RBm along the corresponding trajectories and in the most representative cluster for each system is reported in (C) and (D), respectively, and shows that heparin binding does not reduce the exposure of the RBm (standard deviations are reported in Table S10). E, the change in motion of the up-RBD induced by the presence of heparin is shown by the superimposition of ten conformations extracted at equal time intervals along the trajectories (from magenta to blue) and projected onto the first essential dynamics eigenvector without (left) and with (right) heparin. The RBD is shown in cartoon representation, and heparin is omitted for ease of visualization. Corresponding plots for the replica trajectories are shown in Figs. S6–S10. See text for details. dPCA, dihedral principal component analysis; RBD, receptor-binding domain; SASA, solvent-accessible surface area.

To evaluate if the induced fit promoted by heparin causes the masking of the ACE2 binding residues in the RBm, we calculated the SASA of these residues along the trajectory and their accessibility in the representative clusters. Both the analyses show that the heparin chains do not significantly shield the residues of the RBm (Figs. 4, C and D and S9, Table S10). In addition, we superimposed the ACE2-RBD crystal structure (PDBid 6M0J (25)) onto the open spike with heparin (Fig. S10). This superposition confirms that heparin does not interact with the residues of the RBm involved in the interaction with ACE2-RBD. Furthermore, preliminary molecular dynamics simulations of spike–heparin–ACE2 complexes (G.P., data not shown) indicate that the presence of heparin further stabilizes the spike–ACE2 interaction, consistent with the observed ability of heparin to speed up the formation of spike-ACE2 complexes in vitro (11). Thus, the simulations indicate that the same allosteric induced fit mechanism may enable spike to simultaneously bind HSPGs and ACE2 on the host cell surface (11, 13). On the other hand, heparin will, in molar excess, exert an antiviral effect by binding spike, sequestering the virus, and preventing it from interacting with HSPGs and hindering presentation of the RBD to ACE2 on the host cell surface.

Finally, to obtain further insights into the effect of heparin on the open RBD, we performed essential dynamics (ED) analysis (Figs. 4E and S11). The analysis on the closed conformation shows an overall stabilization of the spike by heparin without significant effects on the RBDs (data not shown). ED analysis on the open spike conformation reveals that the binding of heparin results in a different direction of the motion of the RBm loop (residues 472–489) in the upsubunit described by the first eigenvector (Fig. S10). Despite different starting conformations and independent sampling, the upsubunits consistently show this difference across all the replica simulations. These differences suggest that the presence of heparin (or HSPGs) could affect the motion of the RBD, possibly having an inhibitory gating effect on host-cell receptor binding.

Discussion

Experiments have shown that HSPGs are indispensable for SARS-CoV-2 infection (11, 12) and that heparin binds to the SARS-CoV-2 spike glycoprotein, exerting an antiviral effect (11, 12, 13, 15, 16). The binding of heparin to spike occurs with measured Kd values ranging from 40 PM to 64 nM, (12, 13, 26). Interestingly, these values indicate similar or higher affinity than the affinity of heparin for antithrombin III (Kd values ranging from 1.5 to 281 nM (27, 28)), which provides the main mechanism for heparin's anticoagulant activity. Taken together, these findings suggest that heparin exerts its antiviral activity with a potency that may be similar to that of its anticoagulant activity and that the doses of heparin currently used to treat COVID-19 patients ((14) and references therein) may be sufficient to obtain an antiviral effect. However, the routes of administration of heparin required to exert the two activities in vivo are likely to be different: heparin is administered systemically to prevent thromboembolism, while a topical administration of heparin via aerosol has been proposed to maximize its antiviral effect (29) and like classical antiviral agents (30, 31, 32), would be most effective in the first phase of the disease.

Heparin administered systemically may have further therapeutic effects than its anticoagulant and antiviral activity. The S1/S2 site masked by heparin has been identified as a superantigenic site able to overactivate the TCR (9). Therefore, heparin could have an anti-inflammatory effect, especially in patients with multisystem inflammatory syndrome. Moreover, a systemic treatment would not be expected to interfere with the action of antibodies against SARS-CoV-2 because the MD simulations show that heparin does not interact with the main regions targeted by antispike antibodies (33): it does not shield either the RBm or the RBD interface oriented toward the core domain of the adjacent monomer, and it only marginally interacts with the N-terminal domain of spike (Fig. S12).

It has been demonstrated that unfractionated heparin has a 150-fold higher antiviral effect against SARS-CoV-2 than low-molecular-weight heparin (LMWH) (15). Recent data show that long chain heparin inhibits spike–cell interactions, whereas different LMWHs were less effective (34). Furthermore, a spike mutant lacking the arginine-rich furin cleavage site had lower affinity to heparin compared with the spike glycoprotein with the conserved S1/S2 basic domain (34). Accordingly, we demonstrated that the capacity of unfractionated heparin to inhibit furin cleavage of spike is significantly higher than that of LMWH heparin, suggesting that this activity contributes to the antiviral effect of heparin and further pointing to the importance of the length of heparin in its action against SARS-CoV-2 and the role of the S1/S2 basic site of spike glycoprotein. Until now, however, molecular models to investigate the binding of heparin or HSPGs to spike have been limited to models of short heparin chains (≃6–8 monosaccharides) binding to the basic domain of the spike RBD without N-glycans (11, 12, 16) while, to fully understand the antiviral effect of unfractionated heparin and the role of HSPGs in SARS-CoV-2 infection, it is necessary to model the binding of long polyanionic chains to the spike head. This is a challenging task due to the length, variable sulfation pattern and flexibility of the polysulfated glycosaminoglycan chains, and because of the large size and flexibility of the spike head for which some regions and N-glycans are structurally poorly defined (4, 7). Therefore, considering all available experimental data to build high-quality initial models, we first employed our incremental docking and sliding window method (19) for docking 31mer heparin chains to models of the spike head glycoprotein. We then performed multiple microsecond MD simulations to refine the models of the spike head with zero, one, or three heparin chains bound to study the dynamic effects of heparin binding. Despite a total simulation length of over 20 microseconds, the sampling of the configurational space of the spike–heparin systems was inevitably incomplete. Nevertheless, these simulations provided sufficient configurational sampling to allow us to explore the main dynamic features associated with the predicted heparin-binding modes and thus, to identify mechanisms by which long heparin chains exert their antiviral activity and HSPGs can act as coreceptors.

From our results, three key mechanisms by which heparin exerts its antiviral activity can be identified (see Fig. 5): (i) heparin directly competes with HSPGs for the same binding sites on the spike head, burying the same basic surface regions, and hindering the binding to HSPGs of both the closed and open conformations of the spike head; (ii) heparin masks the S1/S2 multibasic site (unique to SARS-CoV-2), preventing the cleavage by furin and the activation of the prefusion conformation of the spike, as well as modulating the triggering of the hyperinflammatory response due to the binding to T-cell receptors of this superantigenic site but not interfering with the immune system (9); (iii) heparin masks the basic residues of the RBD and allosterically acts on the hinge region that is suggested to be key for the opening of the spike and the movement of the RBD to expose the ACE2 receptor-binding face. Notably, both the direct and the allosteric mechanisms require long heparin chains. To validate the results of our simulations, we carried out an enzymatic assay that showed that heparin inhibits cleavage of a spike fragment containing the S1/S2 site by furin in a dose-, length-, and sulfation-dependent manner and that this binding was effectively due to the binding and masking of the S1/S2 basic site by heparin.

Figure 5

Schematic diagram showing the proposed mechanisms by which heparin and HSPGs can affect SARS-CoV-2 infection of a host cell. The results of the molecular dynamics simulations and biochemical assays suggest the following mechanisms: Heparin exerts antiviral activity by hindering the binding of both closed and open conformations of spike to HSPGs. Heparin binds the S1/S2-HBD of spike in both the conformations, preventing furin cleavage (1), hinders the opening/activation of the closed spike (2), and acts allosterically on the hinge region associated with the movement of the RBD, while directly masking the RBD-HBD (3). Moreover, based on the structural similarity between heparin and HSPGs, we expect that HSPGs are able to bind both the closed (4) and open (5) conformations of spike and, in the presence of the ACE2 receptor, favor the activation of the closed spike and its interaction with the ACE2 receptor (4) and subsequent furin cleavage and SARS-CoV-2 infection. However, we cannot exclude a possible ACE2-independent SARS-CoV-2 internalization mediated by HSPGs binding to the closed spike conformation. Created with BioRender.com. ACE2, angiotensin-converting enzyme 2; HBD, heparin-binding domain; HSPG, heparan sulfate proteoglycan; RBD, receptor-binding domain; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

Based on our MD simulations with heparin and on the structural similarity between heparin and the heparan sulfate chains of HSPGs, we can infer that: (i) HSPGs bind along the same paths on the surface of spike occupied by heparin; (ii) HSPGs do not directly induce the activation of the closed spike, suggesting that their role in viral infection may be solely to increase the concentration of the virus on the host-cell surface. However, HSPGs may also mediate the activation of the closed spike by exposing it to the human ACE2 receptor on the host cell. HSPGs may also favor the formation of a ternary complex (9) with the open spike and the ACE2 receptor by binding the basic residues of the RBD of the upsubunit while not masking the RBm. Notably, both the open and closed spike models with heparin bound suggest a modulatory effect of N-glycans on the binding between spike and heparin or HSPGs. Finally, emerging SARS-CoV-2 spike glycoprotein variants have demonstrated increased infectivity (35). We cannot exclude that this capability could in part be due to mutations near to the basic domains mentioned above that increase binding to HSPGs.

In conclusion, from our models, simulations, and experiments, we identify three mechanisms for the antiviral activity of heparin against SARS-CoV-2 and obtain insights into how HSPGs can facilitate viral infection. These results provide a basis for the future development of heparin derivatives endowed with well-balanced antiviral and anticoagulant activities for the effective treatment of COVID-19 patients.

Experimental procedures

Modeling of the systems

The initial models of the SARS-CoV-2 spike head protein in closed and open conformations were taken from the SwissModel website (https://swissmodel.expasy.org/repository/species/2697049) and were based on two structures determined by cryo-electron microscopy: PDBid 6ACC (seq. identity 76.47%) and PDBid: 6ACD (seq. identity 76.47%), respectively (36). The spike models were completed by adding 18 N-glycans per subunit, covalently attached in accordance with the experimentally determined glycomic profile (4), using the Glycam web (http://glycam.org/) (37). Standard protonation states were used. The APBS electrostatics plugin in Pymol (38) was used to compute the electrostatic potential surface calculated with PDB2PQR (39) at neutral pH using PROPKA (40) and the AMBER ff14SB force (41). 31mer heparin chains spanning from the S1/S2 multibasic site to the RBD-HBD were modeled using the incremental docking and sliding window method developed by Bugatti and coworkers (19) using Autodock 4.2 (42). Five model systems were generated: the closed spike with zero, one, or three heparin chains, and the open spike with zero or three heparin chains.

All-atom molecular dynamics (MD) simulation

The Amber20 package (43) was used to carry out the simulations. Parameters for the spike were assigned with the ff14SB (41) and GLYCAM-06j (44) force fields. Heparin was parameterized following the method published by Bugatti et al. (19). All the glycoprotein models were placed in a periodic cubic water box using the TIP3P water model (45) with 10 Å between the solutes and the edges of the box. Na+ and Cl-ions were added to neutralize the systems and to immerse them in solvent with an ionic strength of 150 mM.

Each system was energy minimized in 14 consecutive minimization steps, each of 100 steps of steepest descent followed by 900 steps of conjugate gradient, with decreasing positional restraints from 1000 to 0 kcal/mol A2 on all the atoms of the systems excluding waters, counterions, and hydrogens, with a cutoff for nonbonded interactions of 8 Å. The systems were then subjected to two consecutive steps of heating, each of 100,000 steps, from 10 to 100 K and from 100 to 310 K in an NVT ensemble with a Langevin thermostat. Bonds involving hydrogen atoms were constrained with the SHAKE algorithm (46), and 2 fs time step was used. The systems were then equilibrated at 310 K in four consecutive steps of 2.5 ns each in the NPT ensemble with a Langevin thermostat with random velocities assigned at the beginning of each step. For each system simulated, four independent replica production runs following the same protocol as for the equilibration were carried out starting from restart files chosen randomly from the last 5 ns of equilibration. During the MD simulations, a cutoff of 8 Å for the evaluation of short-range nonbonded interactions was used, and the Particle Mesh Ewald method was employed for the long-range electrostatic interactions. The temperature was kept constant at 310 K with a Langevin thermostat. Coordinates were written at intervals of 100 ps.

Production simulations were carried out on the Marconi100 accelerated cluster (https://www.hpc.cineca.it/hardware/marconi100), which is based on the IBM Power9 AC922 CPU architecture with each accelerator having four Volta V100 NVIDIA GPUs. Each simulation was carried out on a single GPU

Analysis of MD simulations

MD trajectories were analyzed using CPPTRAJ from AmberTools20 (43), and molecular graphics analysis was performed using Visual Molecular Dynamics (VMD) (47).

—Cluster analysis of the structures was carried out using CPPTRAJ (43) for the last 100 ns of all the trajectories considering backbone C-alpha atoms for the protein residues and all the carbon, oxygen, sulfate, and nitrogen atoms for the heparin chains. N-glycans were excluded from the analysis. The hierarchical agglomerative (bottom-up) approach was used with a minimum distance between the clusters of 3.0 Å and the distance between clusters defined by the average distance between members of two clusters (Fig. 2A).

—Hydrogen bond (H-bond) analysis was performed using CPPTRAJ (43) along all the trajectories for frames at intervals of 10 ns (coordinated collected every 100 ps collected with a stride of 100 frames) and setting 3.5 Å as the upper distance for defining a H-bond between heavy atoms. All the atoms including the hydrogens of the systems were considered (Figs. 2A and S3, Tables S1 and S2).

—Interaction fingerprint analysis (IFP) was performed using the MD-IFP python scripts (20) (https://github.com/HITS-MCM/MD-IFP). The interactions between the spike and the heparin chains were computed along all the trajectories for frames at intervals of 10 ns (coordinates were written every 100 ps and then analyzed with a stride of 100 frames) (Table S3, IFP-model02, IFP-model03, IFP-model05).

—Root mean square fluctuations (RMSF) were calculated using CPPTRAJ (43) for all C-alpha atoms of the individual spike subunits—SA, SB, SC—and for all the carbon, oxygen, sulfate, and nitrogen atoms of the heparin chains (Fig. S2).

—Root mean square deviations (RMSD) were calculated using CPPTRAJ (43) for all C-alpha atoms of the individual spike subunits—SA, SB, SC—and for all the carbon, oxygen, sulfate, and nitrogen atoms of the heparin chains (Figs. 2B and S1). The RMSDs of the hinge regions were calculated for the C-alpha atoms of residues 527 to 529 for SA, SB, and SC separately (Figs. 4A and S5).

—Solvent-accessible surface area (SASA). Two separate SASA analyses were carried out for the SA, SB, and SC subunits separately: along the trajectory using CPPTRAJ (43) and for the most representative clusters using NACCESS (http://wolf.bms.umist.ac.uk/naccess/). In both the analyses, the van der Waals radius of the solvent probe was assigned a value of 1.4 Å. For the analysis of the S1/S2-HBD site, residues 682 to 685 were considered (Figs. 3, B and C and S4, Table S4). In the case of the receptor-binding residues, all the residues of the RBD (residues 330–530) were considered along the trajectory but only the RBm residues (K417, G446, Y449, Y453, L455, F456, A475, F486, N487, Y489, Q493, G496, Q498, T500, N501, G502, Y505) suggested by Lan and coworkers (25) for the representative clusters (Figs. 4, C and D and S7, and Table S5).

—Dihedral principal component analysis (dPCA) was performed using CPPTRAJ (43). The dihedral covariance matrix and the projection were calculated for the backbone phi/psi angles of residues 527 to 529 of the SC monomer. The first four eigenvectors and eigenvalues were extracted and the first two principal components were plotted for all of the systems. All the systems were transformed into the same principal component space to evaluate the variance across the replicas (Figs. 4B and S8).

—Essential dynamics (ED) analysis was performed using Principal Component Analysis (PCA) of the unbiased MD simulations. PCA was performed along all the trajectories individually with CPPTRAJ (43). The principal modes of motion were visualized using VMD. The first normalized eigenvectors for model04 and model05 were plotted along the trajectory, and the direction of motion was defined by visual inspection (Figs. 4E and S8).

Reagents

Human recombinant furin was obtained from OriGene Technologies Inc. Conventional heparin (13.6 kDa) was obtained from a commercial batch preparation of unfractionated sodium heparin from beef mucosa (Laboratori Derivati Organici S.p.A.) purified to remove contaminants according to established methods (48). The purity was higher than 95% as assessed by electrophoresis in acidic buffer (49), uronic acid quantitative determination (50), and high-performance liquid chromatograph analysis (48). The 13C NMR spectrum measured following Casu et al. (51) showed 78% N-sulfate glucosamine, 80% 6-O-sulfate glucosamine, and 59% 2-O-sulfate iduronic acid. Low-molecular-weight (LMW) beef mucosal heparins (4.0 kDa) were obtained by controlled nitrous acid degradation of unfractionated heparin as described (51, 52). The capsular E. coli K5 polysaccharide (30,0 kDa) has the same structure [(→4)-β-D-glucuronic acid-(1→4)-α-D-N-acetyl-glucosamine-(1→)]n as the heparin precursor N-acetyl heparosan and was prepared as described in (53). To assess their integrity, the polysaccharides were reanalyzed before the beginning of the experiments (data not shown).

Furin cleavage assay

The ability of heparin to inhibit spike cleavage at the S1/S2 site by furin was evaluated by using the colorimetric assay “CoviDrop SARS-Cov-2 targeted proprotein convertase inhibitor screening fast kit” (Epigentek). Briefly, a 2.0 kDa SARS-CoV-2 spike fragment containing the intact S1/S2 681RRAR684 HBD is tagged with polyhistidine and biotin at its N-terminus and C-terminus, respectively. The spike fragment is then immobilized onto microplate wells through histidine/Ni-NTA. The cleavage of the substrate at the S1/S2 site removes the C-terminal S2 fragment of spike that is linked to biotin, causing a decrease of the signal generated by avidin/biotin binding that is detected by an appropriate colorimetric reaction (recorded by the absorbance in a microplate spectrophotometer at a wavelength of 450 nm. See Tables S8 and S9 for values of single measurements). Furin cleavage inhibition blocks the reduction of the signal; consequently, the extent of spike cleavage is inversely proportional to the signal intensity. The assay was performed according to the manufacturer's instructions (https://www.epigentek.com/catalog/covidrop-sars-cov-targeted-proprotein-convertase-inhibitor-screening-fast-kit-p-84596.html).

Data availability

All simulation trajectories are available on the BioExcel COVID-19 platform: https://bioexcel-cv19.bsc.es/#/ with the identifiers from MCV1900217 to MCV1900236. Preprint server. aXriv - https://arxiv.org/abs/2103.07722

Supporting information

This article contains supporting information (16, 20, 25, 43, 54, 55, 56, 57).

Conflict of interest

The authors declare that they have no conflicts of interest with the contents of this article.