International Journal of Molecular Sciences Article Influence of Association on Binding of Disaccharides to YKL-39 and hHyal-1 Enzymes Agnieszka Krzemińska 1 , José-Emilio Sánchez-Aparicio 2 , Jean-Didier Maréchal 2, * , Agata Paneth 3 and Piotr Paneth 1, * 1 International Center for Research on Innovative Biobased Materials (ICRI-BioM)-International Research Agenda, Lodz University of Technology, Żeromskiego 116, 90-924 Lodz, Poland; agnieszka.krzeminska-kowalska@p.lodz.pl 2 Insilichem, Departament de Química, Facultat de Ciències, Universitat Autònoma de Barcelona, 08193 Cerdanyola del Vallès, Barcelona, Spain; joseemilio.sanchez@uab.cat 3 Department of Organic Chemistry, Medical University of Lublin, Chodźki 4a, 20-093 Lublin, Poland; agata.paneth@umlub.pl * Correspondence: jeandidier.marechal@uab.cat (J.-D.M.); piotr.paneth@p.lodz.pl (P.P.) Abstract: Disaccharide complexes have been shown experimentally to be useful for drug delivery or as an antifouling surface biofilm, and are promising drug-encapsulation and delivery candidates. Although such complexes are intended for medical applications, to date no studies at the molecu- lar level have been devoted to the influence of complexation on the enzymatic decomposition of polysaccharides. A theoretical approach to this problem has been hampered by the lack of a suitable computational tool for binding such non-covalent complexes to enzymes. Herein, we combine quantum-mechanical calculations of disaccharides complexes with a nonstandard docking GaudiMM engine that can perform such a task. Our results on four different complexes show that they are mostly stabilized by electrostatic interactions and hydrogen bonds. This strong non-covalent stabilization Citation: Krzemińska, A.; demonstrates the studied complexes are some excellent candidates for self-assembly smart materials, Sánchez-Aparicio, J.-E.; Maréchal, useful for drug encapsulation and delivery. Their advantage lies also in their biocompatible and J.-D.; Paneth, A.; Paneth, P. Influence biodegradable character. of Association on Binding of Disaccharides to YKL-39 and hHyal-1 Keywords: disaccharides; complexation; docking; DFT; GaudiMM Enzymes. Int. J. Mol. Sci. 2022, 23, 7705. https://doi.org/10.3390/ ijms23147705 Academic Editor: Marcel Popa 1. Introduction Advances in molecular-biology techniques and biotechnological processes have led Received: 15 May 2022 to a plethora of therapeutical biomolecules with a high molecular weight that includes Accepted: 9 July 2022 peptides, recombinant proteins, enzymes, monoclonal antibodies, antibody-drug conju- Published: 12 July 2022 gates, hormones, cytokines, blood factors, and nucleic acids [1]. Many of them require Publisher’s Note: MDPI stays neutral careful handling and storage conditions to ensure their quality. This is the case, for example, with regard to jurisdictional claims in of thermal proteins with instability that requires storage and transport within a limited published maps and institutional affil- range of temperature to expand their lifetime. Humidity, photo exposure, oxygen, and iations. mechanical agitation during transportation are other stress-related factors that impact stability. These therapeutics are also challenging in their delivery and bioavailability after oral intake. Most issues arise from (i) low drug solubility in physiological media; (ii) low drug stability in the physiological environment (i.e., drug degradation or denaturation); Copyright: © 2022 by the authors. (iii) a substantial reduction in drug permeability across biological barriers such as skin, Licensee MDPI, Basel, Switzerland. This article is an open access article mucosal membranes, and cell membranes; (iv) incapability of the drug to passively reach distributed under the terms and its specific molecular target; and (v) the drug interacting with intracellular biomolecules conditions of the Creative Commons before reaching its target. Currently, one approach consists of administrating these drugs Attribution (CC BY) license (https:// through parenteral injections, although novel strategies are under vivid investigation and creativecommons.org/licenses/by/ aim at enhancing stability, solubility, and bioavailability, while targeting specific tissues, 4.0/). cell populations, and/or the relevant subcellular structures. Int. J. Mol. Sci. 2022, 23, 7705. https://doi.org/10.3390/ijms23147705 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2022, 23, 7705 2 of 12 One solution to the problems described above is to use polymeric nano- and mi- crocapsules. Typically, capsular systems are made of a core–shell structure in which the inner core and the outer shell can be separately loaded with drug molecules or can be equipped with bio-functional ligands. As a smart material, such capsules are equipped with stimuli-responsive controlled-release functions for targeted and/or controlled release. These capsules are frequently constructed on saccharide motifs [2]. In this regard, efforts have been directed toward studying polysaccharides with repeating units that bear an electrolyte group (polyelectrolytes), such as chitosan (CHI) or furcellaran (FUR), or gly- cosaminoglycans such as hyaluronic acid (HA), cellobiose (CB), or heparine (HP). Such oppositely charged polysaccharides form complexes by non-covalent bonding (in particular, electrostatic, hydrogen bonds, and/or hydrophobic forces) that create three-dimensional networks that are being applied to drug- and cell-delivery carriers, as tissue adhesives or scaffolds for tissue engineering [3]. For example, positively charged CHI and negatively charged polysaccharides form complexes that can be used as a functional hemostat or as a wound dressing [3]. CHI/FUR was successfully applied as the nano-capsule carrier and delivery system of doxorubicin [4]. How these multimeric constructs impact the enzymatic decomposition in cells and how/if they lead to inhibition or activation processes have not been addressed yet. Despite a remarkable amount of experimental work performed on polysaccharides, there is very little information that can be found in the literature on the molecular description of their influence on the enzymes involved in their metabolism. In lieu of experimental data on how this complexation changes enzymatic catalysis, we decided to resort to computational tools to study these phenomena. Gaining insight into these processes from a modeling perspective does not lack chal- lenges though. First, no crystallographic structures have ever been reported in which multiple saccharides interacting between them are bound to a protein. This is, however, a prerequisite for starting the study we are aiming for. Second, the strength of the inter- actions between multiple saccharides could not be easily scored with standard force field approaches, so one needs some non-biased energetic methods such as quantum-mechanical ones (e.g., DFT) to estimate multi-saccharides affinities. Third, to get the initial models of complexes between multi-saccharides and protein receptors, protein–ligand dockings able to deal with multiple non-covalently bound ligands are necessary. Unfortunately, only a few programs allow multiple simultaneous bindings [5], although, to our knowledge, in all of them the ligands are considered as separated entities with no pre-existing non-covalent interaction networks. One way to tackle this problem is to apply GaudiMM [6], a multi- objective genetic algorithm that allows for guiding the conformational search towards complex solutions via additional molecular descriptors than the sole scoring function. For example, for the simultaneous docking of two ligands inside a protein host, one can consider the docking scoring function of the ligand complex together with additional eval- uations of intermolecular (i.e., complex-protein) hydrogen-bond networks. Therefore, the combination of GaudiMM with other approaches such as quantum-mechanical calculations offers a unique framework to pursue our goal: learning how the complexation of different disaccharides influences their ability to bind to the enzymes responsible for their metabolic decomposition. In this study, we have selected four complexes of disaccharides as models. This selection is based on (1) a variety of interacting charged groups, (2) the different net charges of the complexes, and (3) a different number of hydrogen bonds that exist in the complexes. In particular, we have studied the CHI/FUR complex binding to human cartilage chitinase 3-like Protein 2 (YKL-39), an enzyme that hydrolyzes the glycosidic bond in chitin. Bearing in mind that other polyelectrolytes, such as HA, HP, or CB, assembled in strongly bounded complexes can form nano-capsules, we have studied the interaction of HA/CHI, HA/CB, and HP/CB with lysosomal human hyaluronidase-1 (hHyal-1). This enzyme is responsible for HA hydrolysis and can be easily tracked experimentally since it Int. J. Mol. Sci. 2022, 23, 7705 3 of 12 can be detected in plasma and urine [7–9]. The basic properties of individual disaccharides and their complexes studied herein are summarized in Table 1. Table 1. Net charges, number of charged and hydroxyl groups, and number of intermolecular hydrogen bonds of the individual disaccharides and their complexes used in this study. Net # of NH3 + # of OH # of SO3 − # of COO− Net # of Disaccharide Complex Charge Groups Groups Groups Groups Charge H-Bonds HP −4 0 4 3 1 HP/CB −4 6 FUR −1 0 4 1 0 HA/CB −1 2 HA −1 0 5 0 1 HA/CHI +1 3 CB 0 0 8 0 0 CHI/FUR +1 4 CHI 2 2 6 0 0 2. Results and Discussion 2.1. Characterization of Interactions in Polyelectrolyte Complexes Interactions and strength of the four disaccharides complexes (HA/CHI, HA/CB, HP/CB, and CHI/FUR) have been characterized at the quantum level. The three-dimensional struc- tures of the complexes are stabilized by non-covalent interactions, which can be strong or weak depending on changes in their charge distribution upon complexation [10]. Usually, it is a challenging task to evaluate the strength of the polyelectrolyte complexation energy because it depends on many environmental factors such as pH, ionic strength of the solu- tion, the salt content, type of polyelectrolytes, or the presence of other charged entities [11]. In the literature, there are not many theoretical results that describe the strength of the polyelectrolyte complexes. Key hydrogen bond interactions identified in our studies are presented in Figure 1. Three inter-ligand hydrogen bonds in the HA/CHI complex have been identified: one between the amine-hydroxyl pair, one between the two hydroxyl groups, and one between the amine-carboxyl pair. In the HA/CB complex, two inter-ligand hydrogen bonds were identified: one between the carboxyl and hydroxyl groups and the other between the hydroxyl group and the oxygen from the amide group. The HP/CB complex has the most favorable hydrogen bond interaction network, due to the six hydro- gen bonds between the hydroxyl groups from CB and the sulfate groups from HP, and an additional two hydrogen bonds between the hydroxyl groups of CB and the one carboxyl group of HP. Finally, the CHI/FUR complex contains four inter-ligand hydrogen bonds: one between the oxygen atom from the galactose ring of FUR and an amine group of CHI, one hydroxyl-pair hydrogen bond, one hydrogen bond between an amine group of CHI and the sulfate group of FUR, and one hydrogen bond between the hydroxyl group of CHI and the sulfate group of FUR. The interaction energy for complex formation (Table 2) has been decomposed into the electrostatic, exchange, induction, and dispersion components. The HA/CB com- plex presents the lower formation energy (−43.9 kcal·mol−1 ), while the HP/CB complex presents the strongest one (−116.0 kcal·mol−1 ). Noticeably, those complexes are, respec- tively, with the lowest and highest amount of hydrogen bonds, therefore supporting the importance of such interactions in the stabilization of the complex. The other two com- plexes, HA/CHI and CHI/FUR, have intermediate energies of −109.2 kcal·mol−1 and −71.3 kcal·mol−1 , respectively. Electrostatic attraction plays a dominant role in HP/CB for- mation and consists of 38.3% of the overall complexation energy, which exceeds the effects reflected in second-ordered terms (dispersion and induction). In the other three complexes, the electrostatic is about 1/3 of the total energy. Furthermore, the electrostatic attraction is larger for the CHI/FUR complex than the sum of the induction and dispersion terms, while for HA/CHI and HA/CB the electrostatic attraction has a similar value to second-ordered terms. Inductive attraction forces are equally important for HA/CHI (16.3%) and HP/CB (16.5%) formation, a little bit less for the HA/CB (14.3%) complex, and the lowest for CHI/FUR (12.8%). Surprisingly, a large dispersion contribution to the total complexation energy was computed for the HA/CB complex (18.4%) in contrast to the other complexes, and HA/CB the electrostatic attraction has a similar value to second-ordered terms. Induc- tive attraction forces are equally important for HA/CHI (16.3%) and HP/CB (16.5%) for- Int. J. Mol. Sci. 2022, 23, 7705 mation, a little bit less for the HA/CB (14.3%) complex, and the lowest for CHI/FUR 4 of 12 (12.8%). Surprisingly, a large dispersion contribution to the total complexation energy was computed for the HA/CB complex (18.4%) in contrast to the other complexes, where dis- persion is only about where dispersion 10%.about is only For complexes characterized 10%. For complexes by strong attraction characterized by strongforces, it is attraction not surprising that the repulsion exchange must be proportionately large, which forces, it is not surprising that the repulsion exchange must be proportionately large, which is ob- served for the HA/CHI (45.4%) and CHI/FUR (46.7%) complexes. Relatively is observed for the HA/CHI (45.4%) and CHI/FUR (46.7%) complexes. Relatively small small repul- sion was computed repulsion was computedfor thefor HA/CB (34.6%) the HA/CB and HP/CB (34.6%) (34.1%)(34.1%) and HP/CB complexes, perhaps complexes, due to perhaps the duelack of lack to the the specific ionic nature of the specific of CB. The ionic nature strongest of CB. electrostatic The strongest attraction electrostatic and weak- attraction and est exchange weakest repulsion exchange result result repulsion in the in strongest interaction the strongest energy energy interaction obtained for the HP/CB obtained for the complex. HP/CB complex. Figure 1. Optimized Figure 1. Optimized 3-dimensional 3-dimensional structures structures of of polyelectrolyte polyelectrolyte complexes complexesofofHA/CHI, HA/CHI, HA/CB, HA/CB, HP/CB, HP/CB,andandCHI/FUR CHI/FURwith withtheir theirkey keyhydrogen hydrogenbond bondinteractions interactionsmarked markedininlight lightblue. blue. Table 2. The Table 2. The contributing contributing components components toto the complexation energy the complexation energy of of complexes complexes computed computed by using by using −−11 and as a percentage of the total the SAPT0/def2TZVP theory level. Values are reported in kcal·mol and as a percentage of the total the SAPT0/def2TZVP theory level. Values are reported in kcal · mol SAPT0 energy. SAPT0 energy. Energy Component Energy HA/CHI HA/CB HP/CB CHI/FUR HA/CHI HA/CB HP/CB CHI/FUR electrostatic Component −322.6 27.0% −46.5 32.6% −139.7 38.3% −324.9 29.9% exchange electrostatic 541.8 −322.6 45.4% 27.0% 49.3 −46.5 34.6% 32.6% 124.1 −139.7 34.1% 38.3% 508.1 −324.9 46.7% 29.9% induction exchange −193.9 541.8 16.3% 45.4% −20.5 49.3 14.3% 34.6% −60.2 124.1 16.5% 34.1% −139.7 508.1 12.8% 46.7% dispersion induction −134.5 −193.9 11.3% 16.3% −−26.3 20.5 18.4% 14.3% −40.3 −60.2 11.1% 16.5% −114.8 −139.7 10.6% 12.8% dispersion −134.5 11.3% −26.3 18.4% −40.3 11.1% −114.8 10.6% total energy total energy −109.2 −109.2 - - −43.9 −43.9 -- −116.0 −116.0 -- −71.3 −71.3 -- 2.2. Binding of Isolated Ligands and CHI/FUR Complex to YKL-39 According to Suginta and co-workers [12], chitooligosaccharide (GlcNAc2 ) binds to YKL-39 with hydrogen bonds and hydrophobic interactions, with Trp360, in the center of the chitin-binding cleft, playing a prominent role. It is believed that an arriving chitin chain would first attract the center of the binding site. The amino acids located in the center of the binding cleft make a highly ordered net of interactions, where Leu105 and Phe301 Int. J. Mol. Sci. 2022, 23, 7705 5 of 12 interact with one subsite of the ligand, and Tyr104, Asp213, Tyr269, and Trp360 stabilize the other part of the ligand. The binding affinity of GlcNAc2 to YKL-39, equal to −5.1 kcal·mol−1 , was established by isothermal-titration calorimetry [12]. We used the crystal structure 4P8V deposited in the PDB by these authors for computing the binding affinity with PRODIGY. The obtained value of −5.7 kcal·mol−1 is in very good agreement with the experimental value, indicating the good performance of the PRODIGY methodology, which, in contrast to the free-energy-perturbation (FEP) and thermodynamics-integration (TI) methods, is less computationally time-consuming. Although using FEP or TI would probably lead to a more accurate binding-affinity value, the result obtained from PRODIGY should serve our purpose, since we want to compare several binding affinities for different ligands and systems, and we do not focus on absolute values. Results of the computed binding affinities of the GlcNAc2 , FUR, CHI, and CHI/FUR complexes bound to the YKL-39 enzyme, along with a representation of the key interactions, are presented in Figure 2 and Table 3. In the case of the GlcNAc2 ligand, one of its subsites interacts through hydrogen bonds with Tyr104 (3.17 Å) and Asp213 (2.59 Å) and additional van der Waals interactions with Leu210 (4.08 Å) and Tyr269 (3.70 Å), plus it is stabilized by Trp360 (3.45 Å). The other subsite of the ligand interacts via van der Waals interactions with Leu105 (3.02 Å) and Phe301 (3.91 Å). Interestingly, both CHI and FUR ligands present a higher affinity when bound to YKL-39 (−7.2 kcal·mol−1 in both cases) than their native ligand GlcNAc2 . The CHI ligand has three hydrogen bonds with Tyr104 (3.61 Å), Asp213 (2.40 Å), and Tyr269 (3.32 Å) and one hydrophobic interaction with Trp360 (3.85 Å), whereas the FUR ligand has three hydrogen bonds with Asp213 (2.97 Å), Tyr104 (3.22 Å), and Phe301 (3.30 Å) and one hydrophobic interaction with Leu105 (3.72 Å). Finally, the highest binding affinity for the YKL-39 enzyme was observed for the CHI/FUR complex (−8.3 kcal·mol−1 ). This complex has four hydrogen bonds with Leu105 (3.72 Å), Asp213 (2.89 Å), Phe301 (3.27 Å), and Trp360 (2.84 Å) and two hydrophobic interactions with Tyr104 (3.77 Å) and Phe301 (3.27 Å). It can be observed that, in the case of the CHI/FUR complex, the cationic CHI creates more interactions with the binding center, while the anionic FUR mostly stabilizes the structure of the complex and is held by the central residue Trp360. As is shown in Figure 2, the CHI and native GlcNAc2 ligands occupy the same space in the YKL-39 binding site. The sulfate group of the FUR ligand hinders its deeper binding into the YKL-39 binding site, yet it provides more attraction with cavity amino acids. Finally, the most space-consuming CHI/FUR complex binds shallowly, but its charged character, derived from different functional groups, provides the most attractive interactions with the YKL-39 residues. Table 3. Binding affinity (BA) in kcal·mol−1 computed for CHI/FUR, CHI, FUR, and GlcNAc2 bound to YKL-39 and key distances in Å between YKL-39 amino acids and ligands. Hydrogen bonds are marked in bold. BA Tyr104 Leu105 Leu210 Asp213 Tyr269 Phe301 Trp360 GlcNAc2 −5.7 3.17 3.02 4.08 2.59 3.70 3.91 3.45 CHI −7.2 3.61 4.55 4.21 2.40 3.32 4.86 3.85 FUR −7.2 3.22 3.72 6.81 2.97 4.84 3.30 4.02 CHI/FUR −8.3 3.77 2.65 6.23 2.89 4.29 3.27 2.84 2.3. Binding of Isolated Ligands and Complexes HA/CHI, HA/CB, and HP/CB to hHyal-1 We examined the binding to hHyal-1 of the set of ligands CHI, HA, CB, and HP and the complexes HA/CHI, HA/CB, and HP/CB. Results are presented in Figures 3–5, and Tables 4–6. In all cases, the binding of individual polysaccharides as well as their complexes are dominated by hydrogen bond interactions with the hHYal-1 binding center. Int. J. Int. J. Mol. Mol. Sci. 2022, 23, Sci. 2022, x FOR PEER REVIEW 23, 7705 66 of of 12 12 Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 7 of 12 is almost equal to that of the other ligands (−6.4 kcal·mol−1); even though it does not pre- sent any ionic character, its eight alcohol (OH) groups can create hydrogen bonds with the residues, Asp129 (3.04 Å), Glu131 (1.87 Å), Tyr202 (2.08 Å), Tyr247 (3.54 Å), and Trp321 (3.19 Å), of hHyal-1. Multiple hydroxyl groups in saccharide ligands tend to find interactions with adjacent Trp321, for HA, HP, CHI, and CB although Trp321 may also create π-stacking or van der Waals hydrophobic attraction. Non-covalently stabilized complexes bind more strongly than single disaccharides to hHyal-1. They can create more attractive interactions, as observed in the binding affinity equal to −7.2 kcal·mol−1 for HA/CHI-hHYal-1 and HA/CB-hHYal-1 or the −7.9 kcal·mol−1 computed for HP/CB-hHYal-1. From a structural point of view, most interactions are pre- served in the disaccharide complexes, with Asp129 being further than Glu131. In those complexes, hydrogen bond interactions are preponderant. The anionic-cationic HA/CHI complex interacts with Asp129 (3.03 Å), Glu131 (2.52 Å), Tyr202 (2.47 Å), Tyr247 (2.35 Å) via hydrogen bond interactions, and it is stabilized by two π-stacking interactions with Trp321 (3.46 Å) and Trp324 (5.27 Å). Table 4. Binding affinity (BA) computed for HA/CHI, HA, and CHI bound to hHyal-1, and repre- sentation of key distances between hHyal-1 amino acids and ligands. Values are reported in kcal·mol−1 and Å for BA and distances, respectively. Hydrogen bonds are marked in bold. BA Ser76 Asp129 Glu131 Tyr202 Tyr247 Trp321 Trp324 HA −6.3 7.93 2.20 2.05 3.02 2.99 3.53 3.46 CHI −6.3 Figure 2. Ribbon 3.19 2. Ribbon representation 3.75 of 2.91 YKL-39 with 8.15 docked CHI/FUR, CHI, FUR, 6.16 CHI, FUR, and GlcNAc 2 saccha- 3.41and GlcNAc5.30 Figure representation of YKL-39 with docked CHI/FUR, 2 saccha- HA/CHI −7.2 rides, respectively. 4.71 respectively. Key amino acids 3.03Key amino 2.52 that create 2.47interactions with ligands are marked 2.35 ligands are3.46 in black. 5.27 rides, acids that create interactions with marked in black. 2.3. Binding of Isolated Ligands and Complexes HA/CHI, HA/CB, and HP/CB to hHyal-1 We examined the binding to hHyal-1 of the set of ligands CHI, HA, CB, and HP and the complexes HA/CHI, HA/CB, and HP/CB. Results are presented in Figures 3–5 and Tables 4–6. In all cases, the binding of individual polysaccharides as well as their com- plexes are dominated by hydrogen bond interactions with the hHYal-1 binding center. The anionic HA and HP ligands and the cationic CHI ligand exhibit the same binding affinity (−6.3 kcal·mol−1). They all interact with Asp129 (average donor–acceptor distance 3.14 Å), Glu131 (average donor–acceptor distance 2.39 Å), and Trp321 (average donor– acceptor distance 3.21 Å) via hydrogen bonds. In all three cases, Glu131 is closer to the saccharides than Asp129, as it serves as a proton donor for the hydroxyl group for hydrol- ysis. Asp129 is also a key mechanistic residue, since it neutralizes the positive charge on the nitrogen and facilitates transition-state stabilization. It may stabilize Glu131 via proton sharing3. and Figure thusrepresentation Ribbon can occupy aofposition farther the hHyal-1 withfrom theHA/CHI, docked ligand. These HA, andtwoCHI residues have saccharides, Figure 3. Ribbon representation of the hHyal-1 with docked HA/CHI, HA, and CHI saccharides, from left to right, the greatest respectively. impact Keybinding on the final amino acids that value. affinity create interactions withthe Interestingly, ligands are marked low energy in struc- from left to right, respectively. Key amino acids that create interactions with ligands are marked black. ture of HA displays a hydrogen bond network with Tyr202 (donor–acceptor distance 3.02 in black. Å) and Tyr247 (2.99 Å). As shown in Figure 3, the arrangement of those residues is strongly related to Ser245, which does not interact directly with HA but is crucial for the stabilization of the interaction between Tyr202 and Tyr247 observed in the X-ray. A highly sulfated HP glycosaminoglycan required more space, which results in shallower binding into the cleft, with still strong interactions with more external residues such as Tyr75 (2.54 Å) and Lys144 (4.00 Å). Cationic CHI besides Asp129 (3.75 Å), Glu131 (2.91 Å), and Trp321 (3.41 Å), attracts Ser76 (3.19 Å) with hydrogen bonds as well. The binding affinity of CB Int. J. Mol. Sci. 2022, 23, 7705 Figure 3. Ribbon representation of the hHyal-1 with docked HA/CHI, HA, and CHI saccharides, 7 of 12 from left to right, respectively. Key amino acids that create interactions with ligands are marked in black. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 8 of 12 Table 5. Binding affinity (BA) in kcal·mol−1 computed for HA-CB, HA, and CB bound to hHyal-1 and key distances in Å between Hyal-1 amino acids and ligands. Hydrogen bonds are marked in bold. BA Asp129 Glu131 Tyr202 Tyr247 Asp292 Trp321 Trp324 Glu325 HA −6.3 2.20 205 3.02 2.99 2.96 3.53 3.46 6.37 Figure 4. Ribbon representation of the hHyal-1 with docked HA-CB, HA, and CB saccharides, re- CB 4. Ribbon Figure −6.4 representation 3.04 1.87 2.08 of the hHyal-1 3.54docked with 2.84 HA-CB,3.19 HA, and2.60 4.29 CB saccharides, spectively. Key amino acids that create interactions with ligands are marked in black. HA/CB −7.2 respectively. 2.87 Key amino 2.60 acids that 2.54 3.27 create interactions 3.86 are marked with ligands 4.15 in black. 3.68 3.15 Figure 5. Ribbon representation of hHyal-1 with docked HP-CB, HP, and CB saccharides, respec- Figure 5. Ribbon representation of hHyal-1 with docked HP-CB, HP, and CB saccharides, respectively. tively. Key amino acids that create interactions with ligands are marked in black. Key amino acids that create interactions with ligands are marked in black. Table 4. Table 6. Binding Binding affinity affinity(BA) (BA)computed in kcal·molfor −1 computed for HP-CB, HP, and CB bound to hHyal-1 and HA/CHI, HA, and CHI bound to hHyal-1, and represen- key distances in Å between Hyal-1 amino acids and ligands. Hydrogen bonds are marked in bold. tation of key distances between hHyal-1 amino acids and ligands. Values are reported in kcal·mol−1 and Å for BA BAand Tyr75 Asp129 distances, Glu131 respectively. Lys144bonds Hydrogen Tyr202 Tyr247 inTyr286 are marked bold. Asp292 Trp321 HP −6.3 2.54 3.46 2.20 4.00 5.79 6.32 8.56 5.01 3.23 BA Ser76 Asp129 Glu131 Tyr202 Tyr247 Trp321 Trp324 CB −6.4 4.04 3.04 1.87 5.30 2.08 3.54 8.58 2.84 3.19 HA −6.3 7.93 HP/CB −7.9 2.20 3.39 2.67 2.05 1.91 3.61 3.02 3.09 2.99 3.72 2.52 3.53 3.03 2.91 3.46 CHI −6.3 3.19 3.75 2.91 8.15 6.16 3.41 5.30 HA/CHI −7.2 4.71 3.03 2.52 2.47 2.35 3.46 5.27 3. Computational Methods In this study, we employed a two-step computational protocol. First, we carried out Table 5. Binding affinity (BA) in kcal·mol−1 computed for HA-CB, HA, and CB bound to hHyal-1 and quantum-mechanical calculations, which allowed us to ascertain the structures and sta- key distances in Å between Hyal-1 amino acids and ligands. Hydrogen bonds are marked in bold. bility of the individual disaccharide ligands and their polyelectrolyte complexes. The BA three-dimensional Asp129 Glu131 structures Tyr202 resulting Tyr247 fromAsp292 the first step were then Trp321 used in optimized Trp324 Glu325 docking calculations to determine their binding affinity and to analyze the key interac- HA −6.3 2.20 2.05 3.02 2.99 2.96 3.53 3.46 6.37 CB −6.4 tions inside 3.04 YKL-39 or2.08 1.87 hHyal-1 enzymes. 3.54 2.84 3.19 2.60 4.29 HA/CB −7.2 2.87 2.60 2.54 3.27 3.86 4.15 3.68 3.15 3.1. Quantum-Mechanical Calculations TableSingle disaccharide 6. Binding affinity (BA)ligands and−1their in kcal·mol polyelectrolyte computed complexes for HP-CB, HP, were to and CB bound optimized at hHyal-1 and the distances key ωB97x-Din[13]/def2TZVP Å between Hyal-1 [14]amino levelacids of theory using Gaussian09 and ligands. [15].are Hydrogen bonds Subsequently, the marked in bold. symmetry-adapted perturbation theory (SAPT) [16] method was employed to study the BA Tyr75 non-covalent Asp129 interactions Glu131 between Lys144 the saccharides Tyr202 that form Tyr247 a polyelectrolyte Tyr286 Asp292 complex. Trp321 HP −6.3 2.54 This method allows for the decomposition of the overall complex interaction energy, ac- 3.46 2.20 4.00 5.79 6.32 8.56(10) 5.01 3.23 (10) (20) CB −6.4 4.04 cording 3.04to Equation 1.87 (1), into5.30 its electrostatic, 2.08 , exchange, 𝐸𝐸𝑒𝑒𝑒𝑒𝑠𝑠𝑠𝑠3.54 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒ℎ , induction, 8.58 2.84 𝐸𝐸𝑖𝑖𝑖𝑖𝑖𝑖,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 3.19 + 3.39 𝐸𝐸 (20) 2.67 (2) (20) 3.09 (20) 3.72 𝑒𝑒𝑒𝑒𝑒𝑒ℎ−𝑖𝑖𝑖𝑖𝑖𝑖,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 + 𝛿𝛿𝐻𝐻𝐻𝐻 , and dispersion, 𝐸𝐸𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 + 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒ℎ−𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , components. HP/CB −7.9 1.91 3.61 2.52 3.03 2.91 (10) (10) (20) (20) (20) (20) (2) 𝐸𝐸𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆0 = 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 + 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒ℎ + 𝐸𝐸𝑖𝑖𝑖𝑖𝑖𝑖,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 + 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒ℎ−𝑖𝑖𝑖𝑖𝑖𝑖,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 + 𝐸𝐸𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 + 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒ℎ−𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 + 𝛿𝛿𝐻𝐻𝐻𝐻 (1) (2) The subscript resp describes orbital-relaxation effects, while the term corre- 𝛿𝛿𝐻𝐻𝐻𝐻 sponds to the higher-order induction effects obtained from the Hartree–Fock (HF) inter- (10) action energy. The electrostatic ( 𝐸𝐸 ) is a Coulombic charge-clouds-interactions term, Int. J. Mol. Sci. 2022, 23, 7705 8 of 12 The anionic HA and HP ligands and the cationic CHI ligand exhibit the same binding affinity (−6.3 kcal·mol−1 ). They all interact with Asp129 (average donor–acceptor distance 3.14 Å), Glu131 (average donor–acceptor distance 2.39 Å), and Trp321 (average donor– acceptor distance 3.21 Å) via hydrogen bonds. In all three cases, Glu131 is closer to the saccharides than Asp129, as it serves as a proton donor for the hydroxyl group for hydrolysis. Asp129 is also a key mechanistic residue, since it neutralizes the positive charge on the nitrogen and facilitates transition-state stabilization. It may stabilize Glu131 via proton sharing and thus can occupy a position farther from the ligand. These two residues have the greatest impact on the final binding affinity value. Interestingly, the low energy structure of HA displays a hydrogen bond network with Tyr202 (donor–acceptor distance 3.02 Å) and Tyr247 (2.99 Å). As shown in Figure 3, the arrangement of those residues is strongly related to Ser245, which does not interact directly with HA but is crucial for the stabilization of the interaction between Tyr202 and Tyr247 observed in the X-ray. A highly sulfated HP glycosaminoglycan required more space, which results in shallower binding into the cleft, with still strong interactions with more external residues such as Tyr75 (2.54 Å) and Lys144 (4.00 Å). Cationic CHI besides Asp129 (3.75 Å), Glu131 (2.91 Å), and Trp321 (3.41 Å), attracts Ser76 (3.19 Å) with hydrogen bonds as well. The binding affinity of CB is almost equal to that of the other ligands (−6.4 kcal·mol−1 ); even though it does not present any ionic character, its eight alcohol (OH) groups can create hydrogen bonds with the residues, Asp129 (3.04 Å), Glu131 (1.87 Å), Tyr202 (2.08 Å), Tyr247 (3.54 Å), and Trp321 (3.19 Å), of hHyal-1. Multiple hydroxyl groups in saccharide ligands tend to find interactions with adjacent Trp321, for HA, HP, CHI, and CB although Trp321 may also create π-stacking or van der Waals hydrophobic attraction. Non-covalently stabilized complexes bind more strongly than single disaccharides to hHyal-1. They can create more attractive interactions, as observed in the binding affinity equal to −7.2 kcal·mol−1 for HA/CHI-hHYal-1 and HA/CB-hHYal-1 or the −7.9 kcal·mol−1 computed for HP/CB-hHYal-1. From a structural point of view, most interactions are pre- served in the disaccharide complexes, with Asp129 being further than Glu131. In those complexes, hydrogen bond interactions are preponderant. The anionic-cationic HA/CHI complex interacts with Asp129 (3.03 Å), Glu131 (2.52 Å), Tyr202 (2.47 Å), Tyr247 (2.35 Å) via hydrogen bond interactions, and it is stabilized by two π-stacking interactions with Trp321 (3.46 Å) and Trp324 (5.27 Å). 3. Computational Methods In this study, we employed a two-step computational protocol. First, we carried out quantum-mechanical calculations, which allowed us to ascertain the structures and stability of the individual disaccharide ligands and their polyelectrolyte complexes. The three-dimensional structures resulting from the first step were then used in optimized docking calculations to determine their binding affinity and to analyze the key interactions inside YKL-39 or hHyal-1 enzymes. 3.1. Quantum-Mechanical Calculations Single disaccharide ligands and their polyelectrolyte complexes were optimized at the ωB97x-D [13]/def2TZVP [14] level of theory using Gaussian09 [15]. Subsequently, the symmetry-adapted perturbation theory (SAPT) [16] method was employed to study the non-covalent interactions between the saccharides that form a polyelectrolyte complex. This method allows for the decomposition of the overall complex interaction energy, according to (10) (10) (20) (20) Equation (1), into its electrostatic, Eelst , exchange, Eexch , induction, Eind,resp + Eexch−ind,resp + (2) (20) (20) δHF , and dispersion, Edisp + Eexch−disp , components. (10) (10) (20) (20) (20) (20) (2) ESAPT0 = Eelst + Eexch + Eind,resp + Eexch−ind,resp + Edisp + Eexch−disp + δHF (1) Int. J. Mol. Sci. 2022, 23, 7705 9 of 12 (2) The subscript resp describes orbital-relaxation effects, while the δHF term corresponds to the higher-order induction effects obtained from the Hartree–Fock (HF) interaction en- (10) ergy. The electrostatic ( Eelst ) is a Coulombic charge-clouds-interactions term, the exchange (10) (20) (20) (2) interactions (Eexch ) is a repulsive force, the induction ( Eind,resp + Eexch−ind,resp + δHF ) cov- (20) (20) ers both polarization and charge-transfer effects, while the dispersion (Edisp + Eexch−disp ) is an attractive interaction from electron–electron dynamic correlation. The SAPT com- putations can be very accurate, however, the convergence heavily depends on the size of the system. To study the polyelectrolyte-complex interactions, the SAPT0 [17]/def2TZVP theory level has been used as implemented in the Psi4 [18] program. 3.2. Protein-Systems Setup and Interaction Analysis The initial coordinates of YKL-39 and hHyal-1 were taken from the Protein Data Bank, using PDB entries 4P8V [12] and 2PE4 [19], respectively. Protein backbone and side-chain atoms were checked with ChimeraX [20] and Discovery Studio (version 21.1.0.20298, Das- sault Systèmes, Paris, France), while hydrogen atoms were added with the tLeaP module of the Amber program [21]. The GAFF [22] force field was used to prepare the parameters for ligands and complexes. For the case of YKL-39, its native ligand chitooligosaccharide GlcNAc2 , FUR, CHI, and CHI/FUR complexes were studied; whereas for the hHyal-1 case, CHI, HA, HP, CB, and the complexes HA/CHI, HA/CB, and HP/CB were studied. In all cases, the ligands or polyelectrolyte complexes were docked with the GaudiMM software into their corresponding protein (i.e., YKL-39 or hHyal-1). Key distances and position analysis were evaluated with ChimeraX, DiscoveryStudio, the protein–ligand interaction profiler (PLIP) [23], and LigPlot+ [24]. For the resulting complexes, key-interaction analysis and binding-affinities calculations were performed. The thresholds for key interactions were carefully selected to keep the values as restrictive as possible but also to account for a broad range of interactions: the maximum cut-off for the determination of binding-site atoms was 8.5 Å, the maximum distance of carbon atoms for a hydrophobic interaction was 4.0 Å, the maximum distance between acceptor and donor in hydrogens bonds was 4.1 Å, the minimum angle at the hydrogen-bond donor was 100◦ , the maximum distance between ring centers for stacking was 7.5 Å, and the maximum deviation from an optimum angle for stacking was 30◦ [25,26]. The affinities of binding values were computed with the PROtein binDIng enerGY prediction program (PRODIGY) [27]. The Gibbs free energy of binding (∆Gbind ) is described through the thermodynamic term of Gibbs free energy: ∆Gbind = Gbound − Gfree = RT lnK (2) where Gbound defines the energy of the bound form of the complex, Gfree corresponds to the sum of the energies of the protein and ligand separately, R is the gas constant, T is the absolute temperature, and K is the equilibrium constant. 3.3. Docking Protocol Most of the current molecular docking protocols are unable to dock ligands that are non-covalently bonded (e.g., hydrogen-bonded) complexes. GaudiMM (Genetic Al- gorithms with Unrestricted Descriptors for Intuitive Molecular Modeling) is a versatile platform [6] that turned out to be an excellent environment for non-standard-docking tasks. It comprises a multi-objective genetic algorithm that can optimize all the needed variables (objectives) at once and behaves as an evolutionary algorithm for complex molecular prob- lems such as protein–ligand docking. At first, GaudiMM generates an initial population of potential solutions called individuals, which are further evaluated with several objectives. Each of them is assigned a fitness value to find the best results. The best structures are selected by the algorithm and propagated to the next generation. During the next iterations, the initial population of structures evolves to end up with the best set of protein–ligand complexes. For this study, four objectives have been selected: ‘Clashes’ that evaluates steric Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 10 of 12 Int. J. Mol. Sci. 2022, 23, 7705 10 of 12 and GAFF-based parametrization [22] for the other atoms of the system, and the ‘HBonds’ hindrance that descriptor and accounts for protein–ligand helps to eliminate hydrogen-bond close contacts stabilizations. between ligand and proteinThe latest atoms, is the the most relevant for the prospect of this work because it helps to displace AutoDock ‘Vina [28] scoring function’ that is widely used in the literature, an energetic eval-the GaudiMM conformational uation based on exploration the Amber99SB towards solutions forcefield with good for standard intermolecular protein residues and hydrogen GAFF-basedbond interactions. Out of all the structures received in the docking process, parametrization [22] for the other atoms of the system, and the ‘HBonds’ descriptor that the best possible position had accounts been chosen. The for protein–ligand parameters presented hydrogen-bond such as stabilizations. TheClashes, latest isVina the score, energy, most relevant andthe for HBonds were prospect of taken into because this work consideration it helpsto to select the set displace theofGaudiMM best structures. Protocols conformational used for running exploration GaudiMM towards solutionsare provided with in the Supplementary good intermolecular hydrogen Material. bond interactions. Out The best structures were then compared with the position of all the structures received in the docking process, the best possible of the native ligand position hadinbeen the case of 4P8V-bound structures, with particular attention given to the chosen. The parameters presented such as Clashes, Vina score, energy, and HBonds were position of the gly- cosidicinto taken bond around the amino consideration acids to select theofset theofactive center (see Figure best structures. Protocols6). For usedligands bound for running to 2PE4, where GaudiMM there is no are provided in good structure in the the Supplementary Protein Data Bank with a bound ligand, Material. the comparison was made The best structures with were the compared then position ofwiththe HA the tetrasaccharide from the position of the native published ligand in the studies, case in which thestructures, of 4P8V-bound hydrolysiswithreaction of theattention particular peptide bond givenwas investigated to the position of[29]. the glyco- sidic Since bondGaudiMM around thesimulations amino acids inof thethe present activeform do(see center not allow Figurefor6).aFor dynamical ligands change bound of the to 2PE4,protonation where therestates, is noneither of the ligands good structure in thenor of theData Protein protein Bank along withthe search,ligand, a bound proto- nation the states corresponding comparison was made with to the the position optimal of pH theconditions for both proteins HA tetrasaccharide from the have been published used. in which the hydrolysis reaction of the peptide bond was investigated [29]. studies, Figure 6. Figure 6. Comparison ComparisonofofCHI/FUR, CHI/FUR, CHI, and CHI, FUR and FURligand positions ligand obtained positions by docking obtained withwith by docking the native the GlcNAc native 2 molecule in complex with the 4P8V structure. The glycosidic bond overlap is high- GlcNAc molecule in complex with the 4P8V structure. The glycosidic bond overlap 2 lighted in the black box—Cartesian coordinates of atoms in the box are provided in the Supplemen- is highlighted in the black box—Cartesian coordinates of atoms in the box are provided in the tary Materials. Supplementary Materials. 4. Conclusions Since GaudiMM simulations in the present form do not allow for a dynamical change of theThis work combines protonation quantum-mechanical states, neither calculations of the ligands nor withalong of the protein atypical the protein–ligand- search, protona- tion states corresponding to the optimal pH conditions for both proteins haveof docking experiments to investigate the molecular nature of the interaction theused. been com- plexes of disaccharides, as models for polysaccharides, in different key targets. First, quantum-mechanical calculations were used to show that the disaccharides are 4. Conclusions capable Thisofwork forming very quantum-mechanical combines strong complexes stabilized calculationsby with electrostatics and hydrogen atypical protein–ligand- bonds. Among docking the complexes experiments studied to investigate here, CHI/FUR the molecular natureand HA/CHI of the haveof interaction been the shown ex- complexes perimentally of to be disaccharides, asuseful modelsforfordrug delivery or asinan polysaccharides, antifouling different surface biofilm, and our key targets. calculations suggest that two other First, quantum-mechanical complexes,were calculations HA/CBusedand HP/CB, to show arethe that promising drug- disaccharides encapsulation are and delivery capable of forming candidates. very strong Our results complexes show stabilized bythat the CHI/FUR, electrostatics HA/CHI, and hydrogen HA/CB,Among bonds. and HP/CB complexesstudied the complexes are mostly stabilized here, CHI/FUR byand electrostatic HA/CHIinteractions, which have been shown amount to about experimentally 30% to be of the useful fortotal drugcomplexation delivery or as energy. A further an antifouling attraction, surface biofilm,approxi- and our mately 30% of calculations the complexation suggest that two other energy, comesHA/CB complexes, from the andinduction HP/CB,and dispersion drug- are promising com- ponents, which cover polarization, charge transfer, and electron–electron dynamic-corre- encapsulation and delivery candidates. Our results show that the CHI/FUR, HA/CHI, lation effects. HA/CB, The remaining and HP/CB complexescontribution are mostlytostabilized the total bycomplexation electrostaticenergy is the effects interactions, which responsible amount for repulsion. to about 30% of theThe totalelectrostatic complexation stabilization energy. A dominates over the further attraction, contribution approximately 30% of the complexation energy, comes from the induction and dispersion components, which cover polarization, charge transfer, and electron–electron dynamic-correlation effects. The remaining contribution to the total complexation energy is the effects responsible for Int. J. Mol. Sci. 2022, 23, 7705 11 of 12 repulsion. The electrostatic stabilization dominates over the contribution from the hydrogen bonds. Strong non-covalent stabilization makes these complexes excellent candidates for self-assembly smart materials for drug encapsulation and delivery. Their advantage lies also in their biocompatible and biodegradable character [30]. Given that the studied complexes are intended for medical applications, it is worth looking at their enzymatic decomposition in cells, which has not been studied at the molecular level to date. The first stage of evaluating their metabolic consequences is their binding to degrading enzymes and describing the key interactions. This requires specific tools for docking, which include their non-covalent interactions inside the complex stabilization and their non-covalent interactions outside of the complex. Herein, for this nonstandard docking task, we used the GaudiMM engine equipped with a special HBond objective, which drives the conformational exploration towards solutions with hydrogen bonding interactions between the ligand complex and the protein. The separate ligands FUR and CHI and their CHI/FUR complex were docked into YKL-39, which is a specific biomarker for osteoarthritis. Both FUR and CHI interact with the YKL-39 enzyme more strongly than the native GlcNAc2 ligand and even more strongly than the separate ligands FUR and CHI bind to their CHI/FUR complex. This complex requires more space than the separate ligands, but at the same time it can create more attractive interactions due to its larger surface. Moreover, its dual anionic/cationic nature allows for finding the optimal orientation inside the enzymatic cavity that leads to the complex of lower energy. The same situation has been observed for the HA/CHI, HA/CB, and HP/CB complexes that were docked to hHYal-1, which is an enzyme involved in hyaluronic acid biodegradation. The complexes, show higher affinities toward the enzyme active site, than the individual ligands. This observation is missing for both the experimental and computational studies, so it requires more focus in the future. Furthermore, our calculations indicate that it is worth considering HP as a potential anionic saccharide to build up the polyelectrolyte complexes, as its sulfate groups create a very stable complex even with neutral cellobiose. This hypothesis is also worth testing more intensively and thoroughly. Supplementary Materials: The supporting information can be downloaded at: https://www.mdpi. com/1422-0067/23/14/7705/s1. Author Contributions: Conceptualization and methodology, P.P. and J.-D.M.; calculations and data curation, A.K. and J.-E.S.-A.; original draft preparation, A.K. and A.P. All authors have read and agreed to the published version of the manuscript. Funding: This research was conducted as a part of the International Research Agendas PLUS program of the Foundation for Polish Science, co-financed by the European Union under the European Regional Development Fund (MAB PLUS/2019/11). Access to the computational facilities of PLGrid at Cyfronet, Kraków, is gratefully acknowledged. J.-D.M. and J.-E.S.-A. thank the Spanish Ministry of Science and Innovation for financial support (PID-2020-116861GB-I00). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no competing financial interests. References 1. Greish, K.; Fang, J.; Inutsuka, T.; Nagamitsu, A.; Maeda, H. Macromolecular therapeutics: Advantages and prospects with special emphasis on solid tumour targeting. Clin. Pharmacokinet. 2003, 42, 1089–1105. [CrossRef] [PubMed] 2. Pitakchatwong, C.; Schlegel, I.; Landfester, K.; Crespy, D.; Chirachanchai, S. Chitosan nanocapsules for pH-triggered dual release based on corrosion inhibitors as model study. Part. Part. Syst. Charact. 2018, 35, 1800086. [CrossRef] 3. Ishihara, M.; Kishimoto, S.; Nakamura, S.; Sato, Y.; Hattori, H. Polyelectrolyte complexes of natural polymers and their biomedical applications. Polymers 2019, 11, 672. [CrossRef] [PubMed] Int. J. Mol. Sci. 2022, 23, 7705 12 of 12 4. Milosavljevic, V.; Jamroz, E.; Gagic, M.; Haddad, Y.; Michalkova, H.; Balkova, R.; Tesarova, B.; Moulick, A.; Heger, Z.; Richtera, L.; et al. Encapsulation of doxorubicin in furcellaran/chitosan nanocapsules by layer-by-layer technique for selectively controlled drug delivery. Biomacromolecules 2020, 21, 418–434. [CrossRef] 5. Eberhardt, J.; Santos-Martins, D.; Tillack, A.F.; Forli, S. AutoDock Vina 1.2.0: New docking methods, expanded force field, and python bindings. J. Chem. Inf. Model. 2021, 61, 3891–3898. [CrossRef] 6. Pedregal, J.R.-G.; Sciortino, G.; Guasp, J.; Municoy, M.; Maréchal, J.-D. GaudiMM: A modular multi-objective platform for molecular modeling. J. Comput. Chem. 2017, 38, 2118–2126. [CrossRef] 7. Jedrzejas, M.J.; Stern, R. Structures of vertebrate hyaluronidases and their unique enzymatic mechanism of hydrolysis. Proteins 2005, 61, 227–238. [CrossRef] 8. Davies, G.J.; Gloster, T.M.; Henrissat, B. Recent structural insights into the expanding world of carbohydrate-active enzymes. Curr. Opin. Struct. Biol. 2005, 15, 637–645. [CrossRef] 9. Rye, C.S.; Withers, S.G. Glycosidase mechanisms. Curr. Opin. Chem. Biol. 2000, 4, 573–580. [CrossRef] 10. Rathee, V.S.; Sidky, H.; Sikora, B.J.; Whitmer, J.K. Role of associative charging in the entropy–energy balance of polyelectrolyte complexes. J. Am. Chem. Soc. 2018, 140, 15319–15328. [CrossRef] 11. Laaser, J.E.; McGovern, M.; Jiang, Y.; Lohmann, E.; Reineke, T.M.; Morse, D.C.; Dorfman, K.D.; Lodge, T.P. Equilibration of micelle–polyelectrolyte complexes: Mechanistic differences between static and annealed charge distributions. J. Phys. Chem. B 2017, 121, 4631–4641. [CrossRef] 12. Ranok, A.; Wongsantichon, J.; Robinson, R.C.; Suginta, W. Structural and thermodynamic insights into chitooligosaccharide binding to human cartilage chitinase 3-like protein 2 (CHI3L2 or YKL-39). J. Biol. Chem. 2015, 290, 2617–2629. [CrossRef] 13. Chai, J.D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. [CrossRef] 14. Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. [CrossRef] 15. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, USA, 2009. 16. Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes. Chem. Rev. 1994, 94, 1887–1930. [CrossRef] 17. Parker, T.M.; Burns, L.A.; Parrish, R.M.; Ryno, A.G.; Sherrill, C.D. Levels of symmetry adapted perturbation theory (SAPT). I. Efficiency and performance for interaction energies. J. Chem. Phys. 2014, 140, 094106. [CrossRef] 18. Turney, J.M.; Simmonett, A.C.; Parrish, R.M.; Hohenstein, E.G.; Evangelista, F.; Fermann, J.T.; Mintz, B.J.; Burns, L.A.; Wilke, J.J.; Abrams, M.L.; et al. Psi4: An open-source ab initio electronic structure program. WIREs Comput. Mol. Sci. 2012, 2, 556. [CrossRef] 19. Chao, K.L.; Herzberg, O. Structure of human hyaluronidase 1, a hyaluronan hydrolyzing enzyme involved in tumor growth and angiogenesis. Biochemistry 2007, 46, 6911–6920. [CrossRef] 20. Goddard, T.D.; Huang, C.C.; Meng, E.C.; Pettersen, E.F.; Couch, G.S.; Morris, J.H.; Ferrin, T.E. UCSF ChimeraX: Meeting modern challenges in visualization and analysis. Protein Sci. 2018, 27, 14–25. [CrossRef] 21. Case, D.H.M.; Belfon, I.Y.; Ben-Shalom, S.R.; Brozell, D.S.; Cerutti, T.E.; Cheatham, G.A., III; Cisneros, V.W.D.; Cruzeiro, T.A.; Darden, R.E.; Duke, G.; et al. Amber 2020; University of California: San Francisco, CA, USA, 2021. 22. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. [CrossRef] 23. Salentin, S.; Schreiber, S.; Haupt, V.J.; Adasme, M.F.; Schroeder, M. PLIP: Fully automated protein–ligand interaction profiler. Nucleic Acids Res. 2015, 43, W443–W447. [CrossRef] 24. Laskowski, R.A.; Swindells, M.B. LigPlot+: Multiple ligand-protein interaction diagrams for drug discovery. J. Chem. Inf. Model. 2011, 51, 2778–2786. [CrossRef] 25. McGaughey, G.B.; Gagne, M.; Rapppe, A.K. pi-Stacking interactions. Alive and well in proteins. J. Biol. Chem. 1998, 273, 15458–15463. 26. Hubbard, R.E.; Haider, K. Hydrogen Bonds in Proteins: Role and Strength. In Encyclopedia of Life Sciences; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2010. 27. Xue, L.C.; Rodrigues, J.P.; Kastritis, P.L.; Bonvin, A.M.; Vangone, A. PRODIGY: A web server for predicting the binding affinity of protein-protein complexes. Bioinformatics 2016, 32, 3676–3678. [CrossRef] 28. Trott, O.; Olson, A.J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2010, 31, 455–461. [CrossRef] 29. Zhang, L.; Bharadwaj, A.G.; Casper, A.; Barkley, J.; Barycki, J.J.; Simpson, M.A. Hyaluronidase activity of human hyal1 requires active site acidic and tyrosine residues. J. Biol. Chem. 2009, 284, 9433–9442. [CrossRef] 30. Yadav, S.; Sharma, A.K.; Kumar, P. Nanoscale self-assembly for therapeutic delivery. Front. Bioeng. Biotechnol. 2020, 8, 127. [CrossRef]