Identification of a Pentasaccharide Lead Compound with High Affinity to the SARS-CoV-2 Spike Protein via In Silico Screening

The spike (S) protein on the surface of the SARS-CoV-2 virus is critical to mediate fusion with the host cell membrane through interaction with angiotensin-converting enzyme 2 (ACE2). Additionally, heparan sulfate (HS) on the host cell surface acts as an attachment factor to facilitate the binding of the S receptor binding domain (RBD) to the ACE2 receptor. Aiming at interfering with the HS-RBD interaction to protect against SARS-CoV-2 infection, we have established a pentasaccharide library composed of 14,112 compounds covering the possible sulfate substitutions on the three sugar units (GlcA, IdoA, and GlcN) of HS. The library was used for virtual screening against RBD domains of SARS-CoV-2. Molecular modeling was carried out to evaluate the potential antiviral properties of the top-hit pentasaccharide focusing on the interactive regions around the interface of RBD-HS-ACE2. The lead pentasaccharide with the highest affinity for RBD was analyzed via drug-likeness calculations, showing better predicted druggable profiles than those currently reported for RBD-binding HS mimetics. The results provide significant information for the development of HS-mimetics as anti-SARS-CoV-2 agents.


Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of coronavirus disease 2019 (COVID- 19), which was declared a global pandemic by the World Health Organization (WHO) in 2020.Although the pandemic situation is under control at present, whether new variants will emerge is unpredictable, thus COVID-19 still poses a great threat to human health [1].SARS-CoV-2 is a member of the betacoronaviruses [2], and infects human cells via a trimeric S (spike) protein on its surface binding to the angiotensin-converting enzyme 2 (ACE2) receptor [3] through the receptor-binding domain (RBD) in an 'open' conformation [4].Heparan sulfate (HS) has been identified [5] as a coreceptor priming the S protein for ACE2 interaction [6], which then triggers the fusion of the viral membrane with the host cell membrane [7].Mutations within the S protein of SARS-CoV-2 emerged rapidly, which is assumed to alter virus-host cell interactions [8][9][10].Vaccines are critical tools to reduce the severity of virus-induced diseases; however, vaccine effectiveness against SARS-CoV-2 is threatened by the emergence of viral variants [11].Therefore, it is necessary to develop drugs targeting viral infection, including the inhibition of viral entry, virus replication and release, and virus-induced inflammation.Since HS plays important roles in virus-cell receptor interactions, identifying HS-mimetics as potential inhibitors for targeting S protein-HS-ACE2 interactions is a valid approach to attenuate SARS-CoV-2 infection [12][13][14].Here, we have established an unprecedented pentasaccharide library containing 14,112 species.Through virtual screening of the binding between RBD and the pentasaccharides, we have identified a potent lead structure with a high affinity for the RBD.Drug-likeness calculations showed better the properties of the pentasaccharide than those of the reported oligosaccharide compounds.Collectively, our results provide insights into the potential to develop HS-mimetics as therapeutic agents which prevent SARS-CoV-2 binding to target cells.

Identification of the RBD-Binding Pentasaccharide
The trimeric SARS-CoV-2 S protein engages human ACE2 with a single RBD extending from the trimer in an "open" active conformation (Figure 1) and mediates the entry of virions into target cells [15].There is a cluster of positively charged amino acid residues adjacent to the ACE2 binding site in the open RBD conformation, serving as a binding site for HS [14].To identify oligosaccharides (HS-mimetics) binding to the RBD with high affinity, we established an unprecedented virtual library of pentasaccharides, using the backbone sequences of five sugar units of glucuronic acid, iduronic acid, and glucosamine, with sulfate substitution in the possible positions, resulting in a library of 14,112 different structures.The pentasaccharide library was virtually screened for a binding to the S protein with open RBD conformation, and the top-hit pentasaccharide (designated AD08043, Figure 2) is GlcNS-β-(1→4)-GlcA-α-(1→4)-GlcNAc3S-β-(1→4)-IdoA-α-(1→4)-GlcNAc.This pentasaccharide contains one IdoA unit that increases its conformational flexibility [16], and a 3-O-sulfated glucosamine unit that is a signature of the antithrombin binding epitope of heparin and the herpes virus binding epitope of HS [17,18].Moreover, we found that the top ten pentasaccharides virtually screened, including AD08043, all exhibited a poor degree of sulfation (Table S1).
inflammation.Since HS plays important roles in virus-cell receptor interactions, id ing HS-mimetics as potential inhibitors for targeting S protein-HS-ACE2 interactio valid approach to attenuate SARS-CoV-2 infection [12][13][14].Here, we have establis unprecedented pentasaccharide library containing 14,112 species.Through virtual s ing of the binding between RBD and the pentasaccharides, we have identified a lead structure with a high affinity for the RBD.Drug-likeness calculations showed the properties of the pentasaccharide than those of the reported oligosaccharid pounds.Collectively, our results provide insights into the potential to develop HS-m ics as therapeutic agents which prevent SARS-CoV-2 binding to target cells.

Identification of the RBD-Binding Pentasaccharide
The trimeric SARS-CoV-2 S protein engages human ACE2 with a single RBD e ing from the trimer in an "open" active conformation (Figure 1) and mediates the e virions into target cells [15].There is a cluster of positively charged amino acid re adjacent to the ACE2 binding site in the open RBD conformation, serving as a bindi for HS [14].To identify oligosaccharides (HS-mimetics) binding to the RBD with h finity, we established an unprecedented virtual library of pentasaccharides, usi backbone sequences of five sugar units of glucuronic acid, iduronic acid, and glucos with sulfate substitution in the possible positions, resulting in a library of 14,112 di structures.The pentasaccharide library was virtually screened for a binding to the tein with open RBD conformation, and the top-hit pentasaccharide (designated AD Figure 2 This pentasaccharide contains one IdoA unit that increases its conformational fle [16], and a 3-O-sulfated glucosamine unit that is a signature of the antithrombin b epitope of heparin and the herpes virus binding epitope of HS [17,18].Moreov found that the top ten pentasaccharides virtually screened, including AD08043, all ited a poor degree of sulfation (Table S1).

Analysis of the Molecular Interactions via Molecular Modelling
Next, we analyzed binding affinity of the selected pentasaccharide AD08043 with RBD in comparison with other oligosaccharides that have shown a binding to RBD (Figure 3), including the pentosan polysulfate monomer (PPS) [20,21] and (IdoA2S-GlcNS6S) 4 [6] via the molecular docking of SARS-CoV-2 RBD (Table 1).The results indicate that AD08043 binds more tightly with SARS-CoV-2 RBD (with a binding energy of −7.149 kcal/mol) than (IdoA2S-GlcNS6S) 4 and PPS (with binding energies of −6.797 and −6.247 kcal/mol, respectively).

Analysis of the Molecular Interactions via Molecular Modelling
Next, we analyzed binding affinity of the selected pentasaccharide AD08043 with RBD in comparison with other oligosaccharides that have shown a binding to RBD (Figure 3), including the pentosan polysulfate monomer (PPS) [20,21] and (IdoA2S-GlcNS6S)4 [6] via the molecular docking of SARS-CoV-2 RBD (Table 1).The results indicate that AD08043 binds more tightly with SARS-CoV-2 RBD (with a binding energy of −7.149 kcal/mol) than (IdoA2S-GlcNS6S)4 and PPS (with binding energies of −6.797 and −6.247 kcal/mol, respectively).To generate further insights into the intermolecular interactions and stability of the complexes between the pentasaccharide and RBD protein, we examined the weak nonbond forces such as hydrogen bonds and hydrophobic contacts for each of the amino acids involved in the binding (Table 2).Hydrogen bonding has a diverse role in ligand-receptor interactions.The favorability of the ligand-receptor hydrogen bonds depends upon the total energy change involved in the formation and breaking of all these hydrogen bonds.Likewise, the ligand-receptor hydrogen bonds are considered key contributors to binding due to their specificity.AD08043 shows seven hydrogen bonds involving amino acid residues of Arg355 (3.18 Å and 3.06 Å), Arg466 (3.14 Å, 3.11 Å, and 3.18 Å), Tyr396 (2.94 Å and 3.14 Å), and Leu517 (3.11 Å).In comparison, the (IdoA2S-GlcNS6S)4 oligosaccharide shows a maximum of 24 interactions, including thirteen hydrogen bonds with binding interface amino acid residues Arg355 (3.11 Å), Arg357 (3.19 Å, 2.95 Å, 3.06 Å, 3.00 Å, 2.96 Å, and 3.24 Å), and Arg466 (3.14 Å and 3.18 Å), as well as eleven hydrophobic contacts with the binding interface amino acid residues Trp353 and Asn354.PPS shows six hydrogen bonds with amino acid residues of Arg355 (2.93 Å and 3.14 Å), Arg357 (3.03 Å and 3.35 Å), and Arg466 (3.28 Å and 3.03 Å) (all of which are located in the binding interface),  To generate further insights into the intermolecular interactions and stability of the complexes between the pentasaccharide and RBD protein, we examined the weak nonbond forces such as hydrogen bonds and hydrophobic contacts for each of the amino acids involved in the binding (Table 2).Hydrogen bonding has a diverse role in ligandreceptor interactions.The favorability of the ligand-receptor hydrogen bonds depends upon the total energy change involved in the formation and breaking of all these hydrogen bonds.Likewise, the ligand-receptor hydrogen bonds are considered key contributors to binding due to their specificity.AD08043 shows seven hydrogen bonds involving amino acid residues of Arg355 (3.18 Å and 3.06 Å), Arg466 (3.14 Å, 3.11 Å, and 3.18 Å), Tyr396 (2.94 Å and 3.14 Å), and Leu517 (3.11 Å).In comparison, the (IdoA2S-GlcNS6S) 4 oligosaccharide shows a maximum of 24 interactions, including thirteen hydrogen bonds with binding interface amino acid residues Arg355 (3.11 Å), Arg357 (3.19 Å, 2.95 Å, 3.06 Å, 3.00 Å, 2.96 Å, and 3.24 Å), and Arg466 (3.14 Å and 3.18 Å), as well as eleven hydrophobic contacts with the binding interface amino acid residues Trp353 and Asn354.PPS shows six hydrogen bonds with amino acid residues of Arg355 (2.93 Å and 3.14 Å), Arg357 (3.03 Å and 3.35 Å), and Arg466 (3.28 Å and 3.03 Å) (all of which are located in the binding interface), and eleven hydrophobic contacts.The results of the molecular modeling also showed that AD08043 was located on the binding interface of RBD and ACE2 (Figure 4).and eleven hydrophobic contacts.The results of the molecular modeling also showed that AD08043 was located on the binding interface of RBD and ACE2 (Figure 4).
Table 2. Detailed information on the non-bond forces between SARS-CoV-2 RBD and the oligosaccharides.Interestingly, although the number and length of the hydrogen bonds of AD08043 is less than PPS and (IdoA2S-GlcNS6S)4, the binding affinity of AD08043 to RBD is stronger.This may be ascribed to the 3-O-sulfation substitution.

Ligand Name
Furthermore, the binding affinity of AD08043 with the RBD of several SARS-CoV-2 variants was also analyzed (Table S2).The results demonstrated a high-affinity binding of AD08043 to the RBD of SARS-CoV-2 variants (Table 3).During the interaction with AD08043, the SARS-CoV-2 B.1.1.7 (alpha) variant's RBD shows six hydrogen bonds with the binding interface amino acid residues of Arg355 (3.18 Å, 3.33 Å, and 3.23 Å) and Arg466 (2.84 Å and 2.95 Å) as well as nine hydrophobic contacts.The SARS-CoV-2 B.1.351(beta) variant's RBD shows six hydrogen bonds with the binding interface amino acid Interestingly, although the number and length of the hydrogen bonds of AD08043 is less than PPS and (IdoA2S-GlcNS6S) 4 , the binding affinity of AD08043 to RBD is stronger.This may be ascribed to the 3-O-sulfation substitution.
The SARS-CoV-2 P.1(gamma) variant's RBD shows two hydrogen bonds with the binding interface amino acid residues Asn448 (2.73 Å) and Arg509 (2.96 Å), as well as eleven hydrophobic contacts with the binding interface amino acid residues Arg346, Asn450, and Tyr451.The SARS-CoV-2 B.1.617.2 (delta) variant's RBD shows six hydrogen bonds with the binding interface amino acid residues Arg346 (2.80 Å) as well as nine hydrophobic contacts with the binding interface amino acid residues Ala348, Ala352, Asn354, and Tyr451.The SARS-CoV-2 B.1.1.529(omicron BA.1) variant's RBD shows six hydrogen bonds with the binding interface amino acid residues Arg346 (2.81 Å), Ser349 (3.06 Å), Asn354 (2.96 Å), and Asn450 (3.01 Å), as well as eight hydrophobic contacts with the binding interface amino acid residues Phe347, Ala348, and Ala352.Exceptionally, the B.1.351(beta) variant's RBD has shorter hydrogen bond lengths and more hydrophobic contacts.It should be pointed out that the interactions between AD08043 and the RBDs of the different variants of SARS-CoV-2 are located around the interface of the RBD-ACE2 (Figure 5).To provide a more comprehensive analysis, we further analyzed the binding affinity and intermolecular interactions of the top ten pentasaccharides with the RBD.The results demonstrated that the other pentasaccharides bind with a slightly lower affinity than AD08043.Additionally, the top ten pentasaccharides all bind to the RBD at the RBD-HS-ACE2 binding interface, demonstrating a comparable binding pattern with the AD08043-RBD complex (Figure S1 and Table S3).(   the MD simulations and docking.Furthermore, to explore the SARS-CoV-2 RBD-AD08043 binding mechanism during conformation changes, we intercepted a conformation every 10 ns and calculated the electrostatic potential map of the RBD.The results revealed an extended electropositive surface with dimensions and turns/loops consistent with an HSbinding site identified earlier [4,14] (Figure 6), e.g., around the binding interface between RBD and ACE2.AD08043, in almost all of the different time points (except for 10 ns, colored orange in Figure 6), was located around the RBD electropositive surfaces.Next, we conducted MD simulations to verify the reliability of the model constructed from docking.Quantitatively, the RMSD values of the structures from the MD and docking are small (Figure S2), which indicates that the structural variation is minor between the MD simulations and docking.Furthermore, to explore the SARS-CoV-2 RBD-AD08043 binding mechanism during conformation changes, we intercepted a conformation every 10 ns and calculated the electrostatic potential map of the RBD.The results revealed an extended electropositive surface with dimensions and turns/loops consistent with an HS-binding site identified earlier [4,14] (Figure 6), e.g., around the binding interface between RBD and ACE2.AD08043, in almost all of the different time points (except for 10 ns, colored orange in Figure 6), was located around the RBD electropositive surfaces.
To understand the molecular interaction mechanisms, configurations of the SARS-CoV-2 RBD-AD08043 complex were analyzed every 10 ns during the MD simulation (80 ns) (Table 4).The molecular interactions shown are those mediated by hydrogen bonds and hydrophobic contacts that play a major role in forming the SARS-CoV-2 RBD-AD08043 complexes.At 0 ns of the MD simulation, AD08043 shows a maximum of 13 interactions including five hydrogen bonds with the binding interface amino acid residues Arg355 (3.07 Å and 2.65 Å) and Arg466 (2.77 Å and 2.83 Å), as well as eight hydrophobic contacts with the binding interface amino acid residues Trp353.At 10 ns, AD08043 shows two hydrogen bonds due to the initial unstable stage of the MD simulation which is consistent with the electrostatic convergence interaction result when AD08043 is positioned on the electronegative surface of the RBD.At 20 ns, AD08043 shows five hydrogen bonds with the binding interface amino acid residues Arg357 (2.60 Å, 2.81 Å, and 2.88 Å).At 30 ns, AD08043 shows four hydrogen bonds with the binding interface amino acid residues Arg357 (2.67 Å and 2.76 Å).However, there is no non-bond interaction at 40 ns, showing that the complex is in a relatively unstable state, which is consistent with the RMSD results, as the RMSD values from the MD and docked structures fluctuate considerably around this time point.At 50 ns, AD08043 shows one hydrogen bond as well as two hydrophobic contacts with the binding interface amino acid residues Arg355 and Arg357.At 60 ns, AD08043 shows five hydrogen bonds with the binding interface amino acid residues Arg357 (2.81 Å, 2.96 Å, and 3. To understand the molecular interaction mechanisms, configurations of the SARS-CoV-2 RBD-AD08043 complex were analyzed every 10 ns during the MD simulation (80 ns) (Table 4).The molecular interactions shown are those mediated by hydrogen bonds and hydrophobic contacts that play a major role in forming the SARS-CoV-2 RBD-AD08043 complexes.At 0 ns of the MD simulation, AD08043 shows a maximum of 13 interactions including five hydrogen bonds with the binding interface amino acid residues Arg355 (3.07 Å and 2.65 Å) and Arg466 (2.77 Å and 2.83 Å), as well as eight hydrophobic contacts with the binding interface amino acid residues Trp353.At 10 ns, AD08043 shows two hydrogen bonds due to the initial unstable stage of the MD simulation which is consistent with the electrostatic convergence interaction result when AD08043 is positioned on the electronegative surface of the RBD.At 20 ns, AD08043 shows five hydrogen bonds with the binding interface amino acid residues Arg357 (2.60 Å, 2.81 Å, and 2.88 Å).At 30 ns, AD08043 shows four hydrogen bonds with the binding interface amino acid residues Arg357 (2.67 Å and 2.76 Å).However, there is no non-bond interaction at 40 ns, showing that the complex is in a relatively unstable state, which is consistent with the RMSD results, as the RMSD values from the MD and docked structures fluctuate considerably around this time point.At 50 ns, AD08043 shows one hydrogen bond as well as two hydrophobic contacts with the binding interface amino acid residues Arg355 and Arg357.At 60 ns, AD08043 shows five hydrogen bonds with the binding interface amino acid residues Arg357 (2.81 Å, 2.96 Å, and 3.28 Å) as well as two hydrophobic contacts.At 70 ns, AD08043 shows two hydrogen bonds with the binding interface amino acid residues Asn354 (3.10 Å) and Arg357 (3.04 Å), as well as eight hydrophobic contacts with the binding interface amino acid residues Arg346, Phe347, Ala348, Arg355, Lys356.At 80 ns, AD08043 shows three hydrogen bonds with the binding interface amino acid residues Arg346 (2.92 Å), Interestingly, at the late MD time points (70 ns and 80 ns), all RBD-HS-ACE2 interfacial amino acid residues are involved in non-bond interactions.Aside from the unstable time points of 10 ns and 40 ns, the measured dynamic interaction involved Arg357, which is likely a key amino acid for maintaining structural stability during the conformational changes of the complex.Notably, the interactions between AD08043 and SARS-CoV-2 RBD in the process of the MD simulation mainly involve the sulfate groups (3OS, NS), carboxyl group on the GlcA unit, and the hydroxyl groups on AD08043 skeletons (Figures S4-S12).Further, 3-O-sulfation on the GlcNAc unit seems to be critical for stabilization of the RBD-AD08043 complex, since it is directly involved in all the interactions during the entire period of MD analysis after stabilization.We have also noticed that the initial phase of the interaction involved amino acids both located at the interaction interface and outside of the interface; however, towards the end of the interaction, only interface amino acids are involved (Figure 7).By the end of the MD simulation, AD08043 is completely located on the RBD-HS-ACE2 binding interface.We also performed point mutations to re-assess the effects of key amino acids on the HS-RBD interaction.As expected, the point mutations of the amino acids involved in hydrogen bond formation, through their conversion to alanine residues, lead to the loss of RBD-HS-ACE2 binding interface interactions (Figure S13).

Frontier Molecular Orbital Analysis of the Ligands
Frontier molecular orbitals (FMOs) are usually referred to as the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO).The structural properties of the ligands used for developing the docking interactions are important in explaining the density distribution pattern on frontier molecular orbitals [22].The energetics of FMOs are crucial to the reactivities and chemical structures of ligand molecules.Differences in HOMO-LUMO energy is of importance in determining the reactivity and stability of the molecules.To explore the chemical stability of the pentasaccharide, we performed FMO analysis in comparison with the reported oligosaccharides.The energetic results show that the ∆Eg values of the HOMO-LUMO gap of AD08043, PPS, and (IdoA2S-GlcNS6S)4 are 6.572, 6.283, and 6.092 eV, indicating that the structure of AD08043 is more stable (Figure 8 and Table 5).To confirm the affinity of oligosaccharides towards the RBDs of the variants of SARS-CoV-2 in the ligand-protein complexes, binding free energy calculations were carried out.The calculated results of the Poisson-Boltzmann surface area (MMPBSA) binding free energy are represented in Table S4.The overall trend of the binding free energy values of oligosaccharides-RBD complexes, in decreasing order, is given as PPS < (IdoA2S-GlcNS6S) 4 < AD08043, respectively.The complexes AD08043 and RBD demonstrated the highest binding free energy values compared to other oligosaccharide complexes.Moreover, comparing the binding affinity of AD08043 with the RBD of the various SARS-CoV-2 variants under investigation revealed that the delta and omicron variants had a comparatively lower affinity, which is consistent with the docking results.

Frontier Molecular Orbital Analysis of the Ligands
Frontier molecular orbitals (FMOs) are usually referred to as the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO).The structural properties of the ligands used for developing the docking interactions are important in explaining the density distribution pattern on frontier molecular orbitals [22].The energetics of FMOs are crucial to the reactivities and chemical structures of ligand molecules.Differences in HOMO-LUMO energy is of importance in determining the reactivity and stability of the molecules.To explore the chemical stability of the pentasaccharide, we performed FMO analysis in comparison with the reported oligosaccharides.The energetic results show that the ∆E g values of the HOMO-LUMO gap of AD08043, PPS, and (IdoA2S-GlcNS6S) 4 are 6.572, 6.283, and 6.092 eV, indicating that the structure of AD08043 is more stable (Figure 8 and Table 5).
Figure 7. Detailed information on the interactions between SARS-CoV-2 RBD and AD08043 at 80 ns during the MD simulation period.Amino acids that produce polar and non-polar interactions with the ligand are colored in red and yellow, respectively.

Frontier Molecular Orbital Analysis of the Ligands
Frontier molecular orbitals (FMOs) are usually referred to as the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO).The structural properties of the ligands used for developing the docking interactions are important in explaining the density distribution pattern on frontier molecular orbitals [22].The energetics of FMOs are crucial to the reactivities and chemical structures of ligand molecules.Differences in HOMO-LUMO energy is of importance in determining the reactivity and stability of the molecules.To explore the chemical stability of the pentasaccharide, we performed FMO analysis in comparison with the reported oligosaccharides.The energetic results show that the ∆Eg values of the HOMO-LUMO gap of AD08043, PPS, and (IdoA2S-GlcNS6S)4 are 6.572, 6.283, and 6.092 eV, indicating that the structure of AD08043 is more stable (Figure 8 and Table 5).

Evaluation of the HS-Analogues for Druggability via ADMET Studies
Having observed the high affinity binding between AD08043 and RBD at the interaction interface of ACE2-RBD, it should be interesting to explore the feasibility of developing AD08043 into a therapeutic drug to prevent SARS-CoV-2 from binding to ACE2.Analysis via ADMET (Table 6) revealed that AD08043 possesses a high TPSA as well as aqueous solubility due to its polar groups.The parameter of the Plasma Protein Binding (PPB) being in the range of <90% indicates a relatively high therapeutic response [23].The results also indicate that AD08043 has an appropriate in vivo volume distribution, falling within

Discussion
HS acts as a co-receptor priming SARS-CoV-2 RBD for ACE2 interaction, thus assisting the virus' attachment and entry.Previous studies have shown that HS-mimetics display potent anti-SARS-CoV-2 activity, and the concept of interfering with HS-virus interactions using HS-mimetics for prophylactic as well as therapeutic purposes against virus is drawing attention [5,12,13,21,25].Inspired by the previous findings, we employed a molecular modelling approach, aiming to find HS-mimetics binding to the SARS-CoV-2 RBD with high affinity.Virtual screening of an in-house built pentasaccharide library for their binding to the SARS-CoV-2 RBD identified one top-ranked structure, GlcNS-GlcA-GlcNAc3S-IdoA-GlcNAc, known as AD08043.This unique structure has a relatively low degree of sulfation, but displayed the highest binding affinity to the SARS-CoV-2 RBD, which reveals that the binding affinity is closely related not only to the number and site of negatively charged groups in the oligosaccharide sequence, but also to the changes in RBD protein conformation and the site of key basic amino acids during the interaction.Notably, the pentasaccharide has a 3-O-sulfate substitution that is present in the antithrombin-binding sequence of heparin and the herpes simplex virus binding domain of HS.
MD simulations were performed to determine the fluctuation and conformational changes of AD08043 binding to the RBD under physiological conditions.The result demonstrates that AD08043 at different time points (except for 10 ns) is located around the RBD electropositive surfaces, which is consistent with the mechanism of negatively charged HS-mimetics binding to positively charged amino acid patches in the target protein.In particular, Arg357 repeatedly occurred during this process, suggesting a key role of the amino acid in the conformational stability of the AD08043-SARS-CoV-2 RBD interaction.Additionally, the amino acid residues involved in each period of the MD simulation illustrate that there are non-interface amino acids that interact before the MD simulation (during docking), whereas the amino acid residues that interact at the end of the MD simulation (70 ns, 80 ns) all belong to the interface amino acids, and AD08043 is completely located on the RBD-ACE2 binding interface.This is consistent with the results obtained from our earlier study [26].This finding may represent a general mechanism of oligosaccharides interacting with amino acids at the binding interface of the RBD.
Comparing the binding affinity and intermolecular interaction of AD08043 with the RBDs of several SARS-CoV-2 variants of concern revealed the relatively lower affinity of the delta and omicron variants.One explanation is that the S protein mutations increased the proportion of hydrophobic amino acids, such as leucine and phenylalanine, which may affect the binding affinity.Nevertheless, the interactions between AD08043 and different the RBDs of the variants of SARS-CoV-2 are all located near the interface of RBD-ACE2, indicating that AD08043 may have a broader anti-SARS-CoV-2 potential.The higher stability of AD08043 than PPS and (IdoA2S-GlcNS6S) 4 , and the favorable profile of its druggability further highlights the advantages of AD08043.

Virtual Screening
For the virtual screening validation test, the crystal structure of the SARS-CoV-2 RBD (PDB ID: 6VSB) was used.The protein was protonated at a physiological pH and was then converted to .pdbqt from the .pdbformat using AutoDock MGLTools [28].The pentasaccharides in the library were prepared using the built-in script 'prepare_ligand4.py' in MGLTools to assign polar hydrogen atoms and convert the format from .pdb to .pdbqt.The binding interface of the SARS-CoV-2 RBD as the target of the screening was predicted via Schrodinger SiteMap [29].In the rescoring procedure, the following amino acids at the binding interface were allowed to be flexible: Arg346, Phe347, Ala348, Ser349, Ala352, Trp353, Asn354, Arg355, Lys356, Arg357, Lys444, Asn448, Asn450, Tyr451, Arg466, and Arg509, which is also consistent with previous work [4].The box size was set to 60.00 Å × 60.00 Å × 60.00 Å with Glide [30,31] as the receptor-based virtual screening program.

Virtual Screening
For the virtual screening validation test, the crystal structure of the SARS-CoV-2 RBD (PDB ID: 6VSB) was used.The protein was protonated at a physiological pH and was then converted to .pdbqt from the .pdbformat using AutoDock MGLTools [28].The pentasaccharides in the library were prepared using the built-in script 'prepare_ligand4.py' in MGLTools to assign polar hydrogen atoms and convert the format from .pdb to .pdbqt.The binding interface of the SARS-CoV-2 RBD as the target of the screening was predicted via Schrodinger SiteMap [29].In the rescoring procedure, the following amino acids at the binding interface were allowed to be flexible: Arg346, Phe347, Ala348, Ser349, Ala352, Trp353, Asn354, Arg355, Lys356, Arg357, Lys444, Asn448, Asn450, Tyr451, Arg466, and Arg509, which is also consistent with previous work [4].The box size was set to 60.00 Å × 60.00 Å × 60.00 Å with Glide [30,31] as the receptor-based virtual screening program.

Molecular Dynamics
To examine whether the fluctuation and conformational changes of the model constructed from protein-ligand docking are reliable, a molecular dynamics (MD) simulation of the binding between the RBD of SARS-CoV-2 S protein and the selected pentasaccharide was performed.The AMBER14 force field was utilized for the MD simulation, as implemented in the YASARA program.For the long-range coulomb forces beyond the 8 Å cutoff, the MD simulation used periodic boundary conditions and the particle-mesh Ewald method.NaCl was used at a concentration of 0.9%, and the density of HOH was 0.997 g/mL in the MD cell.No restraints were applied during the MD simulation and the settings employed in the second equilibration dynamics were used.The energies and coordinates were saved every 100 ps with a total simulation length of 80 ns at constant temperature (298 K) and uncontrolled pressure in an NVT ensemble.The structural stability of the RBD-pentasaccharide complex was examined by analyzing the average values of potential energy with root mean square deviation (RMSD) throughout the trajectory.The RMSD profiles of the MD structures (Figures S2 and S3) show that the variation of the RMSD values tends to be stable, which assumes that the structures at equilibrium have been obtained and the last MD structure was chosen for the analysis.
The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) was used to compute the binding free energies of the complexes [36].The calculations were conducted using the YASARA program.The cumulative free energy of the binding is calculated as follows: where G complex is the total MMPBSA energy of the protein-ligand complex, and G protein and G ligand are the highest occupied molecular orbital energies of the separated protein and ligands.

Visualization of Protein-Ligand Complex
LigPlot analysis was used to study the protein-ligand interactions by automatically generating schematic 2D representations of the protein-ligand complexes.The results are shown in a simple and informative diagram of the intermolecular interactions and their strengths, including hydrogen bonds, hydrophobic interactions, and atom accessibilities.The electrostatic potential map of the SARS-CoV-2 S protein's RBD was generated from the crystal structure (PDB: 6VSB) and visualized using PyMoL (version 2.5.4 by Schrödinger).

Quantum Chemical Analysis
The electronic analysis of the selected pentasaccharide and other oligosaccharides was carried out using GAUSSIAN 09 [37], while the orbitals were shown using Multiwfn [38] and VMD 1.9.3 [39] software.The gradient-corrected (density functional theory) DFT with the three-parameter hybrid functional (B3) [40] for the exchange part, and the Lee-Yang-Parr (LYP) correlation function [41], was utilized to compute the molecular structure, vibrational frequencies, and energies of the optimized structures [42].Moreover, to further explain the dispersion interactions that the B3LYP function is unable to describe, B3LYP-D3 was employed [43].Meanwhile, the basis set 6-311 + G(d,p) was augmented by polarization functions on heavy atoms, polarization functions on hydrogen atoms, and diffuse functions for both hydrogen and heavy atoms [44].The optimized geometries have been used to calculate the HOMO and LUMO energy parameters in this study.∆E g = E LUMO − E HOMO (2) where ∆E g is the energy difference between the calculated excited state (LUMO) and the ground state (HOMO), E LUMO is the lowest unoccupied molecular orbital energy, and E HOMO is the highest occupied molecular orbital energy.

ADMET Study
The potential druggability of the selected pentasaccharides was evaluated using the web-based platform ADMETlab2.0(https://admetmesh.scbdd.com/service/evaluation/index (accessed on 20 July 2023)).ADMETlab2.0[45] utilizes a series of high-quality prediction models trained via a multi-task graph attention framework to conveniently and efficiently implement the calculation and prediction of the physicochemical properties, pharmacokinetics, and toxicology of the ligands.Molecules were uploaded in SMILES format.The web server automatically standardizes the input SMILES strings and computes all the endpoints including physicochemical properties, medicinal chemistry friendliness, and the ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties.

Conclusions
The multidimensional biological functions of HS offer a wide range of pharmacological development potentials through targeting HS-protein interactions.Although hundreds of HS-binding proteins have been identified, knowledge of the exact structures of HS binding to a defined protein is still limited.One primary challenge in development of HS-mimetics is to identify a specific HS structure that interacts with a target protein with high affinity.A typical example is the heparin pentasaccharide binding to antithrombin [46].With the rapid progress in the illustration of protein structures, virtual screening provides a powerful tool for the preliminary identification of HS structures.As such, computer-aided study is an efficient and economical way to rationally design drug candidates, and our work provides an example, showing the potential of this approach.This established procedure may be applied to a broad-spectrum screening of oligosaccharides-protein interactions, to explore potential novel drug candidates.
Despite the high throughput and efficiency of virtual screening, the most challenging task in the development of oligo/polysaccharide drugs is production of the compounds.The successful synthesis of the anticoagulant fondaparinux [47] demonstrates the possibility of chemically synthesizing oligosaccharides; however, the complicated chemical process including protection, activation, coupling, and de-protection, requires multiple steps, often with low yield.Nevertheless, the emerging synthetic biology is expected to break the bottle-necks of organic synthesis, making the synthesis of complicated carbohydrate structures possible [48], adding another pathway to the preparation of HS-mimetics for pharmaceutical development purposes.Yet, when considering the highly complicated structures of oligo/polysaccharides, more efficient synthesis procedures will be required, with tremendous effort needed to meet the requirement of diverse oligosaccharides structures.

Figure 1 .
Figure 1.Structure of the trimeric SARS-CoV-2 S protein showing the open and closed form of RBD (PDB ID: 6VSB) [19].

Figure 1 .
Figure 1.Structure of the trimeric SARS-CoV-2 S protein showing the open and closed form of RBD (PDB ID: 6VSB) [19].

Figure 4 .
Figure 4.The binding molecular model of oligosaccharides of the HS-analogue and SARS-CoV-2 RBD.The green surface indicates the RBD-ACE2-RBD binding interface.AD08043 is colored in purple, PPS is colored in teal, and (IdoA2S-GlcNS6S)4 is colored in yellow.

Figure 4 .
Figure 4.The binding molecular model of oligosaccharides of the HS-analogue and SARS-CoV-2 RBD.The green surface indicates the RBD-ACE2-RBD binding interface.AD08043 is colored in purple, PPS is colored in teal, and (IdoA2S-GlcNS6S) 4 is colored in yellow.

Figure 7 .
Figure 7. Detailed information on the interactions between SARS-CoV-2 RBD and AD08043 at 80 ns during the MD simulation period.Amino acids that produce polar and non-polar interactions with the ligand are colored in red and yellow, respectively.

Figure 7 .
Figure 7. Detailed information on the interactions between SARS-CoV-2 RBD and AD08043 at 80 ns during the MD simulation period.Amino acids that produce polar and non-polar interactions with the ligand are colored in red and yellow, respectively.

Figure 9 .
Figure 9. Illustration of the pentasaccharide structures in the library.The carbon atom serial number is shown in the monosaccharide plane structure.

Figure 9 .
Figure 9. Illustration of the pentasaccharide structures in the library.The carbon atom serial number is shown in the monosaccharide plane structure.

Table 1 .
Results of the molecular docking between SARS-CoV-2 RBD and the ligands.

Table 2 .
Detailed information on the non-bond forces between SARS-CoV-2 RBD and the oligosaccharides.

Table 3 .
Detailed information on non-bond forces between AD08043 and RBDs of different variants SARS-CoV-2.

Table 4 .
Detailed information on the non-bond forces of SARS-CoV-2 RBD and AD08043 at different time points during the MD simulation period.