key: cord-0667222-mao64y4r
authors: Paiardi, Giulia; Richter, Stefan; Oreste, Pasqua; Urbinati, Chiara; Rusnati, Marco; Wade, Rebecca C.
title: Three-fold mechanism of inhibition of SARS-CoV-2 infection by the interaction of the spike glycoprotein with heparin
date: 2021-03-13
journal: nan
DOI: nan
sha: a4baebaa6fd0421a8fb5744ab0511916352c5965
doc_id: 667222
cord_uid: mao64y4r

Heparin has been found to have antiviral activity against SARS-CoV-2. Here, by means of sliding window docking, molecular dynamics simulations and biochemical assays, we investigate the binding mode of heparin to the virus spike glycoprotein and the molecular basis for its antiviral activity. The simulations show that heparin binds at long, mostly positively charged patches on the spike, thereby masking the basic residues of the receptor binding domain and of the S1/S2 site. Experiments corroborated the simulation results by showing that heparin inhibits the cleavage of spike by furin by binding to the basic S1/S2 site. Our results indicate that heparin exerts its antiviral activity by both direct and allosteric mechanisms. Furthermore, the simulations provide insights into how heparan sulfate proteoglycans on the host cell can facilitate viral infection. Our results will aid the rational optimization of heparin derivatives for SARS-CoV-2 antiviral therapy.

In the last year, 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) . The spike mediates the entry of the virus into human cells by binding to the angiotensin-converting enzyme 2 (ACE2) receptor (5) . 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 (5) . 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, 5, 6) . This cleavage activates the formation of the post-fusion conformation of spike, which is necessary for SARS-CoV-2 entry into host cells (5, 6) . 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 (7) .

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

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 (8, 9, 10) . To date, three putative HBDs have been identified in the spike sequence: in the receptor binding domain (RBD), at the S2' TMPRSS2 cleavage site (15) (both present also in SARS-CoV-1, and at the novel S1/S2 furin cleavage site (10) (Fig.  1A) .

While the mechanism underlying HSPG-mediated internalization and the antiviral activity of heparin are well understood for many HSPG-dependent viruses (16) , 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 instead facilitate virus infection. In particular, we consider the mechanistic consequences of the sequence features 3 that are specific to the SARS-CoV-2 spike for the effects of heparin and HSPGs on viral infection and inflammation.

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 down-subunits. 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 up-subunit. To investigate the binding of heparin to spike, we modelled 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 (17) 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 tip of the spike head.

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 (9, 10, 11, 13, 14) . 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 modelled systems. , 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 

The five modelled 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 3 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-400 ns in all simulations, as shown by the root mean square deviation (RMSD) relative to each starting structure ( Fig. 2B and Fig. 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 ( Fig. 2B and Fig. 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 , Tab. S1-S3). Accordingly, the RMSD of heparin indicates an induced fit along the trajectory ( Fig. 2B and Fig. 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) .

Overall, the H-bond and interaction fingerprint analysis (18) 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 ( Fig. 2A, Fig. S3 , Tab. 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) ( Fig. 2A, Fig.S3 , Tab. S4 and Tab. 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 , Tab.S1 and Tab.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 up-subunit during the simulations (Fig 2A, Fig. S3 , Tab. 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.

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 S4 and Movies S1-S3. • Heparin masks the spike S1/S2 site To assess the ability of heparin to mask the novel furin cleavage site (5, 6) , which has also been identified as a novel putative HBD (10) and as part of a superantigenic sequence (7) 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 ( Fig. 3 and Fig. S5 , Tab. S4-S7). In the initial modelled 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 , Tab. 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 ( Fig. 3B and Fig. S5 ). 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 up-subunit, but with the heparin chains also maintaining a lower shielding level at the S1/S2-HBDs of the down-subunits. This difference could be due to a more favourable arrangement of the basic patches for heparin binding between the up-and down-subunits compared to that between two down-subunits.

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 Nglycans and heparin ( Fig. 3C and Tab. 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 down-subunits is lower than for the same site of the up-subunit. Moreover, the decomposition shows that the Nglycans 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 to previous coronavirus strains. Moreover, heparin could exert its antiviral activity by masking this site, preventing cleavage by furin and hindering the activation of the post-fusion conformation of the spike glycoprotein. 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.

To confirm the hypothesis from our simulations, we measured the effect of heparin on the cleavage of a spike fragment containing the S1/S2 site. As shown in Fig. 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 pre-incubated 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 2M NaCl wash, which is known to disrupt heparin/protein interactions (19, 20) , 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 (13, 14) . 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 (30kDa, ~150 monosaccharides). As shown in Fig. 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.

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 the formation of the post-fusion 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 (7). 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 substrateimmobilized spike fragment was pre-incubated with heparin (100 M) for 10 min., washed three times with PBS or PBS containing 2 M NaCl (salt) and then exposed to furin. Single measurements are reported in plot and in Tab. S8. (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 Material and Methods section. Each value is the mean ± standard deviation of three to ten repetitions. *p<0.001, one-way ANOVA. The list of the single measurements for each experiment condition are reported in Tab. S8-S9.

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) (21) .

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 (Tab. 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 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 cluster representatives. Both the analyses show that the heparin chains do not significantly shield the residues of the RBm (Fig. 4C-D, Fig. S8 and Tab. S10). On one hand, these data suggest that the heparin chains act indirectly on these domains through an induced fit mechanism as described above. On the other hand, these data are consistent with the ability of spike to simultaneously bind HSPGs and ACE2 on the host cell surface (9, 11) .

Finally, to obtain further insights into the effect of heparin on the open RBD, we performed essential dynamics (ED) analysis ( Fig. 4E and Fig. S9 ). 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 up-subunit described by the first eigenvector (Fig. S9) . Despite different starting conformations and independent sampling, the up-subunits 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 a gating effect on host-cell receptor binding. 

Experiments have shown that HSPGs are indispensable for SARS-CoV-2 infection (9, 11) , and that heparin binds to the SARS-CoV-2 spike glycoprotein, exerting an antiviral effect (9, 10, 11, 13, 14) . Importantly, it has been demonstrated that unfractionated heparin has a 150-fold higher antiviral effect against SARS-CoV-2 than low molecular weight heparin (LMWH) (13) . Until now, 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 (9, 10, 14) . However, to understand the antiviral effect of unfractionated heparin and the role of HSPGs, 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, 5) . Therefore, considering all available experimental data to build high-quality initial models, we first employed our incremental docking and sliding window method (17) 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 co-receptors.

We identified three key mechanisms by which heparin exerts its antiviral activity (see Figure 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 (7); (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 sulfationdependent manner and that this binding was effectively due to the binding and masking of the S1/S2 basic site by heparin.

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 up-subunit 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 (22) . 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, our models and simulations suggest one direct and two allosteric mechanisms for the antiviral activity of heparin against SARS-CoV-2 and provide insights into how HSPGs can facilitate viral infection. They thus provide a basis for the rational optimization of therapeutic heparin derivatives against SARS-CoV-2. (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 ACE2independent SARS-CoV-2 internalization mediated by HSPGs binding to the closed spike conformation.

Modelling 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 (23) . 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/) (24). Standard protonation states were used. The APBS electrostatics plugin in Pymol (25) was used to compute the electrostatic potential surface calculated with PDB2PQR (26) at neutral pH using PROPKA (25) and the AMBER ff14SB force (28) . 31mer heparin chains spanning from the S1/S2 multibasic site to the RBD-HBD were modelled using the incremental docking and sliding window method developed by Bugatti and co-workers (17) using Autodock 4.2 (29) . Five model systems were generated: the closed spike with zero, 1 or 3 heparin chains, and the open spike with zero or 3 heparin chains.

The Amber20 package (30) was used to carry out the simulations. Parameters for the spike were assigned with the ff14SB (31) and GLYCAM-06j (32) force fields. Heparin was parameterized following the method published by Bugatti et al. (17) . All the glycoprotein models were placed in a periodic cubic water box using the TIP3P water model (33) 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 150mM.

Each system was energy minimized in 14 consecutive minimisation 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 A 2 on all the atoms of the systems excluding waters, counterions and hydrogens, with a cut-off for non-bonded interactions of 8 Å. The systems were then subjected to two consecutive steps of heating, each of 100000 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 (34) and 2 fs time step was used. The systems were then equilibrated at 310 K in 4 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, 4 independent replica production runs following the same protocol as for the equilibration, were carried out starting from restart files chosen randomly from the last 5ns of equilibration. During the MD simulations, a cutoff of 8 Å for the evaluation of short-range non-bonded 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 analysed using CPPTRAJ from AmberTools20 (30) and molecular graphics analysis was performed using Visual Molecular Dynamics (VMD) (35) .

-Cluster analysis of the structures was carried out using CPPTRAJ (30) 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 (30) 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 ( Fig. 2A, Fig.S3 , Tab.S1 and S2).

-Interaction fingerprint analysis (IFP) was performed using the MD-IFP python scripts (18) (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) (Tab.S3, IFP-model02, IFP-model03, IFP-model05).

-Root mean square fluctuations (RMSF) were calculated using CPPTRAJ (30) for all Calpha 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 (30) for all Calpha atoms of the individual spike subunits -SA, SB, SC -and for all the carbon, oxygen, sulfate and nitrogen atoms of the heparin chains ( Fig. 2B and Fig.S1 ). The RMSDs of the hinge regions were calculated for the C-alpha atoms of residues 527-529 for SA, SB, and SC separately (Fig. 4A and Fig.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 (30) and for the most representative clusters using NACCESS (36) . 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-685 were considered (Fig.3B-C, Fig.S4, Tab.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 co-workers (21) for the representative clusters ( Fig.4C-D, Fig.S7 and Tab.S5).

-Dihedral Principal Component Analysis (dPCA) was performed using CPPTRAJ (30) . The dihedral covariance matrix and the projection were calculated for the backbone phi/psi angles of residues 527-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 (Fig. 4B and Fig. S6 ).

-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 (30) . 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 (Fig. 4E and  Fig.S8 ). Reagents. Human recombinant furin was obtained from OriGene Technologies Inc. Rockville, MD, USA. 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., Milan, Italy) purified to remove contaminants according to established methods (37) . The purity was higher than 95% as assessed by electrophoresis in acidic buffer (38) , uronic acid quantitative determination (39) , and high-performance liquid chromatograph analysis (37) . The 13 C NMR spectrum measured following Casu et al. (40) 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 in (40, 41) . 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 (42) . To assess their integrity, the polysaccharides were re-analyzed 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 TM SARS-Cov-2 targeted proprotein convertase inhibitor screening fast kit" (Epigentek, County Blvd, Farmingdale NY). 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 Tab. S8-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-convertaseinhibitor-screening-fast-kit-p-84596.html). was supported by Erasmus+, an EMBO short-term fellowship (STF_8594) and The Guido Berlucchi foundation young researchers mobility program. Additional information Preprint server. aXriv -https://arxiv.org/abs/2103.07722 Data availability statement. All simulation trajectories are available on the BioExcel COVID-19 platform: https://bioexcel-cv19.bsc.es/#/ with the identifiers from MCV1900217 to MCV1900236. Supplementary information. Supplementary figures and tables, captions for movies M1-M3 and associated references (PDF). Enlargements of IFP-model02, IFP-model03 and IFP-model05 (tif). Movies S1, S2 andS3 (MP4).

In the following contents, a simplified nomenclature is used for the simulated systems in which they are referred to as Model01-05 and contain the following components:

Spike in closed conformation

Closed spike + 1 heparin chain

Closed spike + 3 heparin chains

Open spike

Open spike + 3 heparin chains Dashed lines indicate the H-bond interactions between the glycans and heparin. As SPR and circular dichroism spectroscopy experiments were done using unfractionated heparins (13.5-15kDa -~48 monosaccharides), we cannot exclude that longer heparin chains could maintain the interaction with these residues (1). Our simulations suggest that N-glycans may exert a shielding effect that results in a non-specific, partial and transient detachment of heparin from the spike.

Exposure of the S1/S2 site during the MD simulations. Solvent Accessible Surface Area (SASA) (Å2) of the S1/S2 site versus time (µs) for the four replica MD trajectories of the five simulated systems. The SASA of the S1/S2 site (residues 682-685) in each spike subunit -SA, SB and SC -is shown in blue, pink and magenta, respectively, and was computed using CPPTRAJ (2). 

Tab. S1 Average root mean square fluctuation (RMSF) (Å) calculated for the C-alpha atoms of the RBD (residues 330-530). The average fluctuation is calculated based on the average fluctuation of the single residues in the RBD. In model02, the single heparin chains directly bind to the RBD of monomer SC. The grey boxes represent the subunits interacting with heparin. Tab. S2 Average root mean square fluctuation (RMSF) (Å) calculated for the C-alpha atoms of the hinge region (residues 527-529). The average fluctuation is calculated based on the average fluctuation of the single residues in the hinge region. In model02, the single heparin chain allosterically acts on the hinge region of monomer SC. The grey boxes represent the subunits interacting with heparin. Tab. S3 Average root mean square fluctuation (Å) calculated for the C-alpha atoms of the S1/S2 site (residues 682-685). The average fluctuation is calculated based on the fluctuation of the single residues in the S1/S2 site. In model02, the single heparin chain binds the S1/S2 site of monomer SA. The grey boxes represent the subunits interacting with heparin. 

The data are collected for the four replicas using the H-bond analysis implemented in CPPTRAJ (2) . Each heparin chain interacts with two adjacent subunits, SA-SB, SB-SC or SC-SA, and the spike residues involved in the interaction are listed separately as residues of the first subunit and of the second subunit. The yellow boxes represent the Hbonds that were stable in more than half the simulations with an occupancy > 50% in the single trajectory considering the 3 heparin chains on each subunit. 

In Model02, one heparin chain interacts with residues of SA and SB subunits. In Models 03 and 05, each of the three heparin chains interacts with residues of SA-SB or SB-SC or SC-SA. Histograms show the interaction in the first 10 frames, last 10 frames and in all the frames in blue, orange and with lines, respectively. Blue and red boxes include the NTD-RBD and S1/S2 domains, respectively. To read the residue name and number, see the attached files IFP-model02, IFP-model03, IFP-model05. These files report the Amber numbering as follows. subunit SA residues 2-1124, subunit SB residues 1289-2411, subunit SC residue number 2574-3696 (fasta numbering for all the subunits; residues 13-1139).

Tab. S7 Solvent Accessible Surface Area (SASA) (Å2) of the S1/S2 site (residues 682-685) calculated for the most representative cluster for each trajectory using NACCESS (4) with a probe radius of 1.4 Å. Tab. S8 Furin cleavage assay. Single measurements obtained for the various experimental conditions used to investigate the ability of heparin to inhibit furin cleavage by binding at the spike S1/S2 site. Rep 

Movie S1. Spike glycoprotein in a closed conformation with 1 heparin chain bound.

First the structure is shown, focusing on the heparin binding site, then the MD trajectory for replica 1 (corresponding to the figures in the main manuscript).

Movie S2. Spike glycoprotein in a closed conformation with 3 heparin chains bound.

The structure and how the three heparin chains bind is shown followed by the MD trajectory for replica 1 (corresponding to the figures in the main manuscript).

The structure and how the three heparin chains bind is shown, focusing on the binding to the open RBD, followed by the MD trajectory for replica 1 (corresponding to the figures in the main manuscript).

Phylogenetic Analysis and Structural Modeling of SARS-CoV-2 Spike Protein Reveals an Evolutionary Distinct and Proteolytically Sensitive Activation Loop

A new coronavirus associated with human respiratory disease in China

The proximal origin of SARS-CoV-2

Site-specific glycan analysis of the SARS-CoV-2 spike

Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

A Multibasic Cleavage Site in the Spike Protein of SARS-CoV-2 Is Essential for Infection of Human Lung Cells

Superantigenic character of an insert unique to SARS-CoV-2 spike supported by skewed TCR repertoire in patients with hyperinflammation

Demystifying heparan sulfate-protein interactions

SARS-CoV-2 Infection Depends on Cellular Heparan Sulfate and ACE2

Characterization of heparin and severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) spike glycoprotein binding interactions

Heparan sulfate proteoglycans as attachment factor for SARS-CoV-2

Heparin as a therapy for COVID-19: current evidence and future possibilities

Unfractionated heparin inhibits live wild type SARS-CoV-2 cell infectivity at therapeutically relevant concentrations

Heparin Inhibits Cellular Invasion by SARS-CoV-2: Structural Dependence of the Interaction of the Spike S1 Receptor-Binding Domain with Heparin

TMPRSS2 and furin are both essential for proteolytic activation of SARS-CoV-2 in human airway cells

Heparan Sulfate Proteoglycans and Viral Attachment: True Receptors or Adaptation Bias? Viruses

Heparin and heparan sulfate proteoglycans promote HIV-1 p17 matrix protein oligomerization: computational, biochemical and biological implications

A Workflow for Exploring Ligand Dissociation from a Macromolecule: Efficient Random Acceleration Molecular Dynamics Simulation and Interaction Fingerprints Analysis of Ligand Trajectories

Hepatitis C virus and hepatitis B virus bind to heparin: purification of largely IgG-free virions from infected plasma by heparin chromatography

Characterization of a peptide containing the major heparin binding domain of human hepatic lipase

Structure of the SARS-CoV-2 spike receptorbinding domain bound to the ACE2 receptor

SARS-CoV-2 Entry Related Viral and Host Genetic Variations: Implications on COVID-19 Severity, Immune Escape, and Infectivity

Cryo-EM structure of the SARS coronavirus spike glycoprotein in complex with its host cell receptor ACE2

Electrostatics of nanosystems: application to microtubules and the ribosome

PDB2PQR: an automated pipeline for the setup, execution, and analysis of Poisson-Boltzmann electrostatics calculations

Very fast empirical prediction and rationalization of protein pKa values

Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

GLYCAM06: a generalizable biomolecular force field. Carbohydrates

Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

VMD: visual molecular dynamics

NACCESS V2.1.1 -Atomic Solvent Accessible Area Calculations. 1992-6

Further Purification of Heparin Reduces Its Bleeding Effects in the Mesenteric Vessels of Rats

Glycosaminoglycans from pig duodenum

A modified uronic acid carbazole reaction

Retention of antilipemic activity by periodate-oxidized non-anticoagulant heparins

Structure of the antithrombin-binding site in heparin

Semi-synthetic heparinoids

Heparin Inhibits Cellular Invasion by SARS-CoV-2: Structural Dependence of the Interaction of the Spike S1 Receptor-Binding Domain with Heparin

A Workflow for Exploring Ligand Dissociation from a Macromolecule: Efficient Random Acceleration Molecular Dynamics Simulation and Interaction Fingerprints Analysis of Ligand Trajectories

NACCESS V2.1.1 -Atomic Solvent Accessible Area Calculations. 1992-6

Structure of the SARS-CoV-2 spike receptorbinding domain bound to the ACE2 receptor

We acknowledge PRACE for awarding us access to Marconi100 based in Italy at CINECA (Project COVID19-54). The technical support of Alessandro Grottesi from CINECA (Italy) and Filippo Spiga from NVIDIA is gratefully acknowledged. G.P., S.R. and R.C.W. thank the Klaus Tschira Foundation and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation -Project number: 458623378 to R.C.W.) for support. G.P.

 (2) . Each heparin chain interacts with two adjacent subunits, SA-SB, SB-SC or SC-SA, and the spike residues involved in the interaction are listed separately as residues of the first subunit and of the second subunit. The yellow boxes represent the H-bonds that were stable in more than half the simulations with an occupancy > 50% in the single trajectory considering the 3 heparin chains on each subunit. The H-bonds for Model02 (data not shown) are not reported because they are very similar to these results.