Screening coronavirus and human proteins for sialic acid binding sites using a docking approach

1.   Introduction

Interaction of sialic acids as the end caps of glycans from host glycoproteins with membrane proteins of the pathogen are considered as an initial event of the viral infectivity cycle of many viruses [1]–[3], e.g. human parainfluenza virus type 2 (HPF3) [4] influenza viruses [5] and also SARS coronaviruses [6]. Cleavage of the sialic acids from the glycoproteins [7], blocking sialic acids directly by carbohydrate binding agents (CBAs) [8],[9], or developing competitive blockers of sialic acids binding sites on the target protein [10],[11] are seen as potential routes for antiviral therapy. In this respect, identification of sialic acid binding sites and providing estimated binding energies are a key feature to support drug development.

Spike (S) glycoprotein from coronavirus is the key protein to enter the host cell, including human coronaviruses causing the common cold, as well as SARS-CoV, MERS-CoV, and also SARS-CoV-2 [12]. The protein is a homotrimeric class I fusion protein [13] built of 1,300 amino acids adopting a rod like shape of about 10 nm length. From the 22 glycosylation sites per protomer of the S protein trimer of SARs-CoV-2, about half of the sites contain fucose and 28 % mannose as end caps of the glycans, while 15 % of the sites are found to contain at least one sialic acid residue [14]. Unlike SARS-CoV S protein, S protein from SARS-CoV-2 needs to be pre-cleaved into two subunits, S1 and S2, a process which is called priming [15],[16]. The priming is done in the host cell by furin protein in the ERGIC (endoplasmic-reticulum-to-Golgi intermediate compartment) [17]. SARS-CoV and HCoV-OC43 S proteins are cleaved by transmembrane serine protease 2 (TMPRSS2) at the site of the host cell plasma membrane [17]. Out of the four domains, A to D, domains A and B are involved in binding to sialic acid and the human host cell receptor angiotensin-converting enzyme 2 (ACE2), respectively, as the initial step of host cell invasion [18],[19]. Upon binding to receptor, the entry process is then completed by the activity of the transmembrane S2 subunit initiating the fusion process [13]. Although SARS-CoV does not have the ability to bind sialic acid, experimental findings indicate a binding site for 9-O-Ac-Me-Sia within domain A of the HCoV-OC43 S protein [20],[21]. Fast on and off rates of 9-O-acatylated sialic acid supports the idea of binding via avidity.

Some evidences propose SARS-CoV-2 S protein can interact with the cell-surface via glycans containing either heparan sulfates [22] or sialic acids [23], which is the initial stage of binding. Computational experiments suggest binding of sialic acids and derivatives to a homologous site at the domain A of SARS-CoV-2 [24].

Furin is a serine endoprotease recognizing a R-X-X-R/K sequence motif and cleaves proproteins to activate them [25]. It is a transmembrane protein with three putative glycosylation sites especially active in the secretory pathway of cells, in particular the trans-Golgi network (TGN). The enzymatic reaction is catalyzed via a catalytic triad of specifically oriented amino acids such as serine, histidine and aspartate [26],[27]. Regarding the invasion of epithelia cells by SARS-CoV-2, activity of furin within the ERGIC of the infected cell by priming the spike protein of SARS-CoV-2 is seen as an important step to make the S-protein susceptible to interaction with the host cell receptor ACE2 [28]. Furin is seen as a potential target for drug development due to its essential role in the life cycle of the virus [25].

The zinc metallopeptidase ACE2 is a receptor glycoprotein specially expressed in epithelia cells and involved in regulating hypertension [29]. Especially in lungs epithelia cells it allows for an entry passage for the virus [30]. Dimeric ACE2 has a high affinity binding site for the S protein of SARS-CoV-2 which makes it the prime reason for the infectivity of the virus [28]. Binding of the S protein to ACE2 readies TMPRSS2 to S protein priming [31]. ACE2 hosts seven glycosylation sites per protomer which also contain sialic acids [32]. Experiments in which sialic acids are removed from the glycans show rather an enhanced binding affinity to the S protein of SARS-CoV-2 than a weakening of the binding indicating a minor role of the sialic acids upon binding with the S protein.

In this study spatial distribution and estimated binding energies of sialic acids to either S protein, ACE2 or furin as negatively charged monomers are investigated. The data are compared to the binding properties of other sugar molecules present in glycans. Available crystal structures of the three proteins are taken as targets. The results are discussed in terms of their implication of the viral infectivity cycle. The questions addressed are (i) whether one of the proteins has a preferred binding site for sialic acids over the other proteins and (ii) whether the estimated binding energies of sialic acids are preferential over the energies of the sugar molecules.

Single sugar molecules are taken as a model system for probing the mode of interaction of the sialic acids with the proteins. Visible inspection is used to screen the poses of the sugar molecules for proper orientation so that they could be linked to the glycan chain of the protein via a α2-6 linkage. In such a pose the linking hydroxy moiety is facing the environment rather than being oriented towards the protein. The interactions of the sialic acid are compared with those of the other sugar molecules which exist in the glycans and are identified as potential end caps.

2.   Experimental section

The protein structures were taken from the Protein Data Bank (www.rcsb.org): domains A of human coronavirus OC43 (HCoV-OC43) spike (S)-proteins (PDB ID: 6NZK) including ligand (l) 9-O-AC-Me-Sia (Sp-l), HCoV-OC43 (PDB ID: 6OHW) with no ligand (Sp), S protein of SARS-CoV-2 (PDB ID: 6VSB) as the original S-protein of CoV-2, (Sp2^or), as well as the peptidase domain of human proteins angiotensin-converting enzyme 2 (ACE2) (PDB ID: 6LZG), and serine-protease furin (PDB ID: 6HZA), which includes in its holo form the synthetic peptide RRKR-Amba (furin-h). Furin-h has also been used in its apo form (furin-a).

The protein structures were either used (i) as experimental structures from the protein data bank without any minimization protocol (mp-0), (ii) structures with the side chains minimized but backbone atoms restraint (mp-1), or (iii) fully minimized structures with both, side chain and backbone atoms minimized (mp-2). Short minimization (using Molecular Operation Environment (MOE) suit, www.chemcomp.com) was done by applying steepest descent, as well as consequent conjugated gradient and truncated-Newton calculations using the Amber10 force field.

The structures of the sugar molecules N-acetylgalactosamine (GalNAc), N-acetylglucosamine (GlcNAc), the sialic acids N-acetylneuraminic acid (Neu5Ac), and N-Glycolylneuraminic acid (Neu5Gc) as well as the sugar molecules fucose, galactose, mannose, were taken from PubChem (pubchem.ncbi.nlm.nih.gov) and generated using the MOE suit. The sialic acid methyl 9-O-acetyl-sialic acid (9-O-Ac-Me-Sia was obtained from PDB data bank (PDB ID: 6NZK). All structures were operated through a short minimization as described above.

2.1.   Structure preparation and minimization

For docking to HCoV-OC43 (Sp-l), which included the sialic acid 9-O-AC-Me-Sia bound to it, the sialic acid was removed and only the sialic acid binding site was considered for docking. The amino acids of the protein forming the binding pocket were selected by applying a radius of 7 Å around each of the atom of the ligand. The same protein HCoV-OC43 but without ligand (Sp) was used in LeadIT to identify the pockets (number of pockets identified by LeadIT: 9 (mp-0), 10 (mp-1), 11 (mp-2)). Swiss-Model (https://swissmodel.expasy.org) was used to fill the existing gaps in Sp2^or via homology modelling with itself. This protocol generated a S-protein monomer, hither forth referred to as Sp2^h. In addition, the three subunits A, B, C of Sp2^or are used separately for docking. Only domain A (N-terminal domain) of each S protein was used for docking.

The furin structure contained a cyclic synthetic peptide RRKR-Amba. Docking was done with furin in its holo form (furin-h), including the RRKR peptide from which the Amba appendix was deleted, as well as in the absence (furin-a) of this peptide. Number of pockets identified by LeadIT were 9 for both, furin-h/furin-a (mp-0), 9 (furin-h) and 12 (furin-a) (mp-1), 9 (furin-h) and 7 (furin-a) (mp-2). Both catalytic domain and P domain were used for docking.

The ACE2 structure used was part of a complex of ACE2 being in contact with the receptor-binding domain of the SARS-CoV2. The S protein was removed and ACE2 P domain was used for docking. Number of pockets identified by LeadIT were 12 using mp-0 and mp-1 structures, and 14 using mp-2 structures.

2.2.   Docking software

LeadIT (BioSolveIT, Germany): Putative pockets (P) for ligand binding were suggested by the software and used for docking. The docking software was used in its default mode and the score values representing estimated free energy ΔG (kJ/mol) were taken for data analysis. With the HYDE routine, the respective values were corrected by a dehydration enthalpy.

Decoy molecules were obtained using DUD-E in combination with the ZINK database [33]. Seven independent runs were conducted from which each run generated 50 molecules. Overlapped ligands found in several runs were deleted and finally 300 ligands ranked by the ZINK database entry number were chosen.

2.3.   Analysis software

Receiver operating characteristic curve (ROC) and area-under the curve (ARC) calculations are done using Origin9 (OriginLab Corporation, Northampton, MA).

3.   Results

The docking approach follows the intention to mimic interaction of sialic acids of the S proteins with the human receptor proteins ACE2 and furin (Figure 1). The estimated binding energies of the S proteins correspond to those sialic acids which are located at the sites of the human proteins, while those values of the human proteins correspond to those values of the sialic acids at the sites of the S proteins. The same scheme applies to the estimated binding energies of the sugar molecules.

3.1.   Binding affinities from the average over various binding sites

A total of eight sugar molecules including 3 sialic acids and 5 sugar molecules which are commonly used as the building block of polysaccharides, are docked to various S protein structures as well as to ACE2 and furin (Figure 2). The values for each of the sugar molecules are a result of averaging over (i) the structures derived from the three minimization protocols (mp-0, mp-1, mp-2), (ii) the best binding energies in the individual pockets which are identified by the docking software of the respective protein, and (iii) the three sialic acids (Figure 3, Sia) and accordingly over the five sugar molecules (Figure 3, Su). Docking results are shown for the experimentally derived structures for Sp-l without the ligand 9-O-AC-Me-Sia at the binding site. In addition, the afore mentioned results are averaged further over the three protein structures Sp, Sp2^h and Sp2^or, as well as the two structures furin-a and furin-h.

The top 20 poses are inspected and used for data analysis. Data analysis are based on either the sialic acids and the sugar molecules ranked as number one (rank-1) or the best scored and oriented sugar molecules (oriented). The oriented sugar molecule is identified by visible inspection opting for an orientation of the respective O-sites on the sugar molecules so that the O-sites could be sterically linked to a putative glycol-chain. In case that none of the 20 poses show an adequate orientation for being selected as ‘oriented’, the finding was considered as having no-docking pose.

The binding energies are the lowest (sialic acids (−20.1 ± 1.6) kJ/mol, sugar molecules (−16.7 ± 0.7) kJ/mol for the oriented ligands docked to Sp-l compared to those values obtained to docking them to the other experimental structures (Figure 3A, blue bars, and Suppl. Table 1a). The sites at the furin-a/h reveal values as low as (−11.0 ± 0.5) kJ/mol for sialic acids and (−13.7 ± 0.3) kJ/mol for the sugar molecules. Correcting the poses for hydration penalty using the HYDE routine does not change the pattern between sialic acids and sugar molecules (Figure 3, orange bars and Suppl. Table 1a). Averaging over the rank-1 ligands, reveals a similar pattern as mentioned for the oriented sugar molecules with values calculated for Sp-l as (−20.6 ± 1.4) kJ/mol for sialic acids and (−17.3 ± 1.3) kJ/mol for the sugar molecules (Figure 3B, and Suppl. Table 1b). The respective values for furin unfold as (−13.3 ± 0.7) kJ/mol for sialic acids and (−14.4 ± 0.3) kJ/mol for the sugar molecules.

The sialic acids bind most strongly to the experimentally identified binding site followed by strong binding sites at furin proteins. When using available S-protein structures without the identified specific site estimated binding energies are indifferent with those values derived for the human proteins. There are only a few structural features for which the sialic acids bind better than the sugar molecules. In many cases the sugar molecules are at least equally, good binders as the sialic acids.

3.2.   The individual docking sites

3.2.1.   Maximum difference values between sialic acids and sugar molecules

The respective difference between the estimated binding energies of sialic acids and sugar molecules are an average taken over the largest difference identified for each pocket of the protein and consequent averaging over either the 3 sialic acids or the 5 sugar molecules. (Suppl. Table 2)

Identifying the biggest difference values on the basis of the LeadIT values for each of the proteins shows that the values for the oriented molecules are as high as (6.1 ± 0.8) kJ/mol for (Sp-l), (5.2 ± 2.1) kJ/mol for ACE2 and (6.9 ± 1.5) kJ/mol for furin-a (Suppl. Table 2). Selecting the rank-1 poses, the numbers are almost similar with (6.0 ± 1.2) kJ/mol for Sp-l, (6.2 ± 1.5) kJ/mol ACE2, and (6.5 ± 1.5) kJ/mol for furin-a.

Focusing on the HYDE values, the differences of the values for the oriented poses are as high as (14.5 ± 4.1) kJ/mol for Sp2^or, (8.9 ± 3.0) kJ/mol for ACE2, and 13.5 kJ/mol for both furin proteins furin-h (standard deviation (SDEV) ± 5.5 kJ/mol) and furin-a (STDEV ± 6.7 kJ/mol) (Suppl. Table 2). The respective HYDE values for the rank-1 poses are (15.7 ± 4.8) kJ/mol for Sp2^or, (12.8 ± 6.2) kJ/mol for ACE2 and (13.9 ± 6.7) kJ/mol for furin-a.

The difference between sialic acid and sugar binding energies can be as high as 6–7 kJ/mol (LeadIT) and 14–15 kJ/mol (HYDE). There is also a trend that the differences for the human proteins are slightly higher than the differences for the viral proteins.

3.2.2.   Best binding sites

All absolute values for both, the oriented and rank-1 poses, are in the range of −20 to −30 kJ/mol (Table 1 and Suppl. Table 3). Searching for the absolute lowest oriented values obtained from all docking poses for the individual proteins reveals that the values do not differ between the sialic acids and sugar molecules ((−23.9 ± 0.7) kJ/mol (sialic acids) versus (−23.9 ± 3.4) kJ/mol (sugar molecules)) for the S-protein structures and those of the human proteins ((−23.5 ± 2.2) kJ/mol (sialic acids) versus (−21.3 ± 0.9) kJ/mol (sugar molecules)) for LeadIT (Table 1). For the S-proteins in two cases, Sp and Sp2^h, the sugar molecules have lower values than the sialic acids. In the case of the human proteins, it is Neu5Ac which adopts the lowest value (e.g. Neu5Gc to ACE2: −24.9 kJ/mol). There is no difference between the energies of the sialic acids at the S proteins and the human proteins.

The pattern is similar if the rank-1 poses are selected (Suppl. Figure 2). The difference between sialic acids and sugar molecules of the S-protein is marginal with values of (−24.1 ± 0.5) kJ/mol (sialic acids) versus (−24.6 ± 2.8) kJ/mol (sugar molecules), but slightly enhanced when comparing the energies amongst the human proteins ((−28.2 ± 2.1) kJ/mol (sialic acids) versus (−23.1 ± 1.4) kJ/mol (sugar molecules)). Amongst the sialic acids, the estimated binding energies are lower for binding to the human proteins than for the sugar molecules.

For the S proteins Neu5Gc (−31.0 kJ/mol) scores highest for the oriented poses looking at the HYDE values (Table 1). When average over the best poses of all the S protein the binding energies are slightly in preference for the sugar molecules ((−27.3 ± 2.9) kJ/mol (sialic acids) versus (−28.3 ± 3.4) kJ/mol (sugar molecules) (p = 0.7)). The lowest values are obtained for the sialic acids 9-O-AC-Me-Sia and Neu5Gc on furin-h (−28.0 kJ/mol) followed by Neu5AC -26.0 kJ/mol at ACE2. The HYDE values of the sialic acids binding to the human proteins are not better than those values for binding to the S protein (p = 0.9). The overall pattern as described remains the same when selecting the rank-1 ligands, except that the lowest estimated binding energy for Neu5Gc (−31.0 kJ/mol) with ACE2 is followed by energy values of Neu5Gc and 9-O-AC-Me-Sia (−28.0 kJ/mol) with the furin proteins (Suppl. Table 3).

In terms of absolute estimated binding energies, the binding energies do not differ so much between the sialic acids and the sugar molecules. When collecting the data over a series of the different protein conformations, the average values are in preference for the sialic acids.

3.3.   Binding poses

The positions of the best poses of the sialic acids on the S-protein coincide with the site identified in the structure of Sp^−l (Figure 4). In case of the human proteins, some sites are close to the identified active sites such as the catalytic site of furin (e.g. orange site in Figure 4B and site in Figure 5B, III)).

The poses of the best oriented sialic acids are stabilized by 8–9 hydrogen (h) bonds (9 h-bonds for Sp-l and furin-a; 8 h-bonds for ACE2) (Figure 5A) and a slightly lower number of hydrogen bonds for the best poses identified by HYDE (6 h-bond for all proteins) (Figure 5B). Additional other hydrophilic residues such as serine and threonine as well as asparagine and glutamine residues are in close contact with the sialic acids. This pattern is also found for the rank-1 positions (10 h-bonds for Sp-l and furin-a, 11 h-bonds for ACE2 for LeadIT (Suppl. Figure 1A) and 6 h-bonds for Sp2^or and furin-h as well as 8 h-bonds for ACE2) for HYDE (Suppl. Figure 1B). HYDE poses show lower number of h-bonds than those for LeadIT which is in concert with the idea of considering dehydration upon binding as a penalty. In many cases the best sites according to the docking results of LeadIT, Neu5Gc Sp-l, and HYDE, Neu5Gc Sp2 ^or, 9-O-Ac-Me-Sia furin-h, are identical for the oriented and the rank-1 ligands.

3.4.   Specificity of binding sites

The best binding sites in terms of absolute binding energies for all three proteins are evaluated using 300 decoys. All decoys chosen for the individual sugar molecules show the highest contribution of molecules with 5–7 rotatable bonds (Suppl. Figure 3). Due to 9-O-AC-Me-Sia (8 rotatable bonds) some decoys also have up to 11 rotatable bonds. Thus, they are chosen in the range of rotatable bonds identified for the sugar molecules which is not exceeding the number of 8 rotatable bonds (Suppl. Table 5).

The AUC values for the oriented poses for LeadIT sample from the ROC plot follows Sp-l > furin-a > ACE2 (rank-1: ACE2 > Sp-l > furin-a (Suppl. Figure 2A and Suppl. Table 5). Looking at the best respective HYDE values the sequence remains the same, however the AUC numbers are high than for the LeadIT values (oriented Sp-l > furin-a > ACE2 (rank-1: ACE2 > Sp-l > furin-a) (Suppl. Table 6).

The AUC values for the oriented and ranked-1 poses for the absolute best HYDE samples are above 0.99 independent of the protein showing a fully specific binding of the sugar molecules (Suppl. Figure 2B and Suppl. Table 5). The AUC values for the corresponding LeadIT values are all below 0.1 with those for furin-h the highest of 0.07 for both oriented and ranked-1 poses compared to 0.04 (oriented) and 0.06 (rank-1) for ACE2 and 0.02 (oriented/rank-1) for Sp2^or. (Suppl. Table 6)

Whilst scoring by LeadIT suggests that the sites are less specific for the individual sugar molecules, the poses scored by HYDE suggest highly specific sites.

4.   Discussion

In this study the LeadIT is used to perform docking of the sialic acids and sugar molecules. The software fragmentizes the molecule at rotatable bonds and reassembles the molecule within the identified binding pockets. It shows good performance [34] and has been used in docking approaches of dipeptidyl peptidase-IV (DPPIV) inhibitors rich in hydroxy groups [35]. The number of rotatable bonds of the sugar molecules is in the range of 0 (fucose) to 8 (9-O-AC-Me-Sia (Suppl. Table 4) and thus, below a reported threshold of ≥ 10 which is reported to lower the performance of the software [36],[37]. No preferred binding topology are considered letting the program placing the sugar molecules without constrains like reported for docking of sugar molecules to e.g. Ca-dependent lectins [38].

4.1.   Interpretation of the results and biological interpretation

A putative binding site of a sialic acid derivative to the S protein is identified experimentally and supported by the docking data. The presented data propose, that this site is a better sialic acid binding site than a sugar binding site. No other site on the A domain of S protein seems to be able to compete with this site in terms of absolute values or whether sialic acids are preferred over sugar molecules. SARS-CoV-2 binds to human proteins, here e.g. ACE2, via their sialic acids is in preference over binding of their sugar molecules. Also, the best estimated binding energy of the sialic acids on a human protein is found for furin, making the interaction of SARS-CoV2 S protein with furin in the ERGIC [17] a moderately enhanced interaction.

Absolute binding numbers including those for HYDE are within a narrow range (28–34 kJ/mol) and with these values, corresponding to an apparent binding constant in the micro molar range, indicating that the sugar molecules are moderate binders [39]. This finding is in accordance with results from experiments in which sialic acids, when connected with the sugar molecules and interact with e.g. CD45 [40] or with sialic acid-binding immunoglobulin-type lectins (siglecs) [41],[42], are classified as week binders with binding affinities in the range of 0.1–3 mM.

While the LeadIT sites are not very specific for the individual sialic acid, the HYDE sites are very specific. Since the differences between the binding energies of the sialic acids and sugar molecules are marginal, this finding supports the idea that the accumulation of the individual binding sites for the polysaccharide chains on the target protein with its sialic-acid caps drives binding via avidity rather than affinity.

Binding of negatively charged sialic acids of host glycoproteins to the S proteins is the entry route of some coronaviruses besides the proteinaceous interaction with the receptors [43]. The S protein can be the target of the polysaccharides of the human proteins, in as much human ACE2 has 7 potential N-glycosylation sites and 3 O-glycosylation sites [44], furin harbours only three potential N-glycosylation sites [45]. On the other hand, the viral S protein of HCoV-OC43 has 23 glycosylation sites (22 for SARS-CoV-2 [46]) sites with 12 sites confirmed to be glycosylated [47],[48] and about 15 % of these sites containing sialic acids [14]. The modest good numbers for sialic acids binding to furin proposes that SARS-CoV-2 S protein will target the furins more thoroughly than the ACE2s which is proposed to increase the efficiency of furin cleavage. It is proposed that this binding of the S protein to the furin is as efficient as binding to ACE2s [17].

Sialic acid binding sites may not interfere with the active sites but are close to active sites, supporting the idea of their role as an anchoring tool locking the approaching proteins into the proper position for interaction.

5.   Conclusions

Sialic acids have a slightly preferred binding mode over the sugar molecules. An overall weak binding affinity in combination with the number of glycosylation sites on either of the proteins, proposes that binding is driven by a combination of avidity and affinity. The weighting of whether either of the scenario prevails could depend on the strength of the affinity or putative number of glycosylation sites on the proteins.

Binding of the SARS-CoV-2 S protein to furin is proposed to be due to affinity rather than avidity, based on somewhat higher binding affinities in combination with the lower number of glycosylation sites on the protein. This scenario is suggested to be reversed in the case of binding of the protein to ACE2 where the moderate strength of the interaction is compensated by avidity.