Theoretical Research of Interactions Between Glycosidases and Glycosaminoglycan Ligands with Molecular Docking and Molecular Dynamics Methods

Computational methods have proved to be sometimes a single tool available for investigation of glycosaminoglycan-protein interactions. A two-stage process including molecular docking with its subsequent detalization using the methods of molecular dynamics is a prospective approach to theoretical modeling of protein-glycosaminoglycan complexes. This review deals with specific features of protein-glycosaminoglycan interactions studied by computational methods, docking and scoring function algorithms, and the use of molecular dynamics results with short-time (ps and ns) changes for processes developing within much longer time scales (ranging over several orders of magnitude). The data obtained with help of computational methods contribute the disclosure of biological interaction mechanism, elaboration of enzyme activity control and grounding of rational recommendations for novel therapeutic means development of high-molecular sort. The results of molecular docking of heparanase, chondroitinlyase ABC, chondroitinase B, and hyaluronidase were shown. The approach to productive design of molecules of compounds (regulating enzyme activity for novel drug derivative obtaining) is representative. The investigations of such kind are directed on ascertainment of action mechanism of these agents in biosystems for production of high efficacy of drug preparations of enzyme nature. It is shown that the molecular dynamics method provides modeling of all degrees of freedom in a protein-ligand complex and draws special attention to protein structure flexibility as a considerable challenge in the development of molecular docking. Computational data are reviewed in the aspect of complex formation between proteins and glycosaminoglycan ligands.


Introduction
Ever-increasing complexity of biomedical studies requires a considerable increase in research budget. The broad range of objects for investigation includes glycosaminoglycans (GAG), components of endothelial glycocalyx [1,2] which together with endothelial cells forms a double protective layer on the vascular wall [3,4]. Pharmacological control over the state of glycocalyx is important clinically and diagnostically for evaluation of the of circulatory system efficacy. So far, only replacement therapy has been used for this purpose [5]. High heterogeneity of GAG concerning chain length, composition and sequence [6] hampers the study of their 3D conformation and molecular dynamics [7]. Bearing in mind repetitive structure of GAG and flexibility of GAG molecules in water, it seems reasonable to employ computational methods for their modeling which often becomes a single tool to provide a better insight into GAG interactions [8]. The approaches associated with quantum and molecular mechanics, molecular docking and dynamics, and coarse graining (for large molecules within a longer time scale) have been used for the investigation of GAG oligosaccharides [7]. Interaction of interleukin-8 with hyaluronan, chondroitin sulfate, dermatan sulfate and their sulfated derivatives was studied with use of large combination of experimental (fluorescence spectroscopy, NMR) and theoretical (molecular docking, molecular dynamics) methods [9]. The results of investigation established that the sulfation pattern determines the strength of binding. Minimal size of GAG fragment (required to achieve specific binding to interleukin-8) has been became tetrasaccharide that agrees with experimental data also. Combination of experimental and computational approaches was productive for construction of GAG models and their complexes with proteins. At present the number of these studies begins to increase. So if the structural databases contain tens of thousands of high-resolution protein structures but the number of GAG-protein complexes composes fewer hundred [10] evincing the trend to their increase [11]. It should be remember the consideration of GAG-protein complex models (from side of structural change of GAG) show the glycosaminoglycans are present predominantly in the extracellular matrix or cell surface. So they are included in endothelial glycocalyx of blood vessels. The effects of GAG consist in tensility and compressibility of tissue, proliferation and recognition of cells, receptor way for viral entry, offering of several opportunities for therapeutic reactions. Molecular composition of GAG includes repeating pair units of hexosamine and uronic acid (or galactose), that are often variably sulfated. The latter conducts to tremendous structural heterogeneity and high charge density hampering structural analysis and elucidation of binding modes. This state stipulates the expediency of computational modeling methods use for understanding the structure and function of these molecules. The hexosamine may be an N-sulfated (GlcNS) or an N-acetylated (GlcNAc) glucosamine or galactosamine (GalNS, GalNAc), variably O-sulfated at the 3, 4 and/or 6 positions. The uronic acid may be a glucuronic acid or an iduronic acid, and may be 2-O-sulfated. It should be noted the employment the Glycosaminoglycan Builder (a point-and-click structure modeling utility /on GLYCAM-Web (http://glycam.org/gag)/) for facilitation 3D structure modeling of GAG fragments [12]. The use of the Glycosaminoglycan Builder contributes the determination of sulfation, acetylation influence of polymeric GAG units on the structure of these compounds. The interaction of heparin tetrasaccharide with chemokine CCL5 is able to show the strong dependence on the pattern and extent of sulfation as well as pH [13]. Molecular docking and molecular dynamics methods were used for elucidation of effective inhibitors of hyaluronidase [14]. This approach skews the accent of investigation of protein-ligand complexes to protein component side. Our review depiction touches on enzyme-GAG complex formation. The main purpose of present research is the consideration of protein conformation changes in these interactions in order to elicit the regulation mechanism of endothelial state. Identification of GAG-binding sites on the protein molecule surface is based on GAG data combined with results of rigid-body docking or electrostatic potential of the protein (with reliance on the dynamic molecular docking data). AUTODOCK, FINDSITE [7], GAG-Dock computational method [15], FRED, Glide [16] etc. have been used for docking.

Enzyme-Glycosaminoglycan Docking
Computational methods have been increasingly used in theoretical study of protein-GAG complexes. Interestingly, GAG represents a unique biological molecule which performs its numerous functions via specific and nonspecific interactions [8]. GAG is located extra-and intracellularly, demonstrating a wide-range time and space dynamics together with size variability rare for other biomolecules. Moreover, GAG offers specific chemical space and information over several orders of magnitude greater than that of other biopolymers, although not all GAG sequences are naturally occurring [8]. Presumably, the size of GAG chain depends on composition of its links, being independent of changes in their sequence [7]. Iduronic acid (IdoA) increases the volume and rigidity of GAG chain, while glucuronic acid (GlcA) impairs its plasticity. This may account for the fact that block-copolymer GAG, such as heparan sulfate, contains microarchitectural elements for multivalent binding with growth factors and collagens [7]. It should be noted that molecular dynamic modeling with regard to large systems, e.g., biomolecules or long-time scales, generally employs calculation of interaction energy and force on the basis of classical mechanical force fields [17]. Upon rare transitions between rotameric states the molecular dynamics is modeled at time periods of 100 ns and longer or with increased approbation to ensure data convergence. Studies of interactions between GAG-cleaving enzymes (glycosidases and carbogen-oxygen glycosidases) with the use of computational methods have provided interesting results only with low molecular weight ligands due to software limitations during reckoning with high molecular large chain GAG.
Human heparanase breaks down heparan sulfate, a proteoglycan constituent of extracellular matrix and basal membrane, and releases heparin/heparan sulfate oligosaccharides. This in turn causes the release of growth factors accelerating tumor growth and metastasizing, which actualizes the search for heparanase inhibitors. A three-dimensional structure of human heparanase was developed using a homology-modeling approach to identify heparanase inhibitors and design antitumor drugs [18]. Stimulation of heparanase synthesis and release (the enzyme is produced by normal and tumor cells) can increase invasion and metastasizing. These processes can be suppressed and blocked by specific heparanase inhibitors. Various GAG inhibitors were docked with heparanase to find out which amino acid residues in the protein interact with these sulfonated saccharides [19]. Apart from electrostatic interactions with heparan-binding domains, hydrophobic interactions contribute to increased binding affinity for some GAG inhibitors. The docking revealed a large binding site extending at least two saccharide units beyond the cleavage site (toward the nonreducing end) and at least three saccharides toward the reducing end (toward heparin-binding site 2). The obtained results allow rational design of heparanase-inhibiting molecules for anticancer drug development by targeting the two heparin/heparan sulfate recognition domains [19]. Docking of roneparstat, a non-anticoagulant 100% N-acetylated and glycol-split heparin acting as a potent heparanase inhibitor, with human heparanase has demonstrated interaction of a roneparstat molecule with one of two heparin-binding domains of heparanase and interaction of the enzyme with two fragments of roneparstat (from two individual molecules) or with heparin-binding domain 1 or 2, which is consistent with different stoichiometry of enzyme-inhibitor binding [20].
Chondroitin lyase is a GAG depolymerizing lyase capable of changing GAG chemical structure, thus promoting anti-tumor activities by inhibiting angiogenesis and tumor metastases [21]. In an attempt to elucidate the mechanisms underlying chondroitinase АВС I activity, complexes between the enzyme and its substrates (chondroitin sulfate and dermatan sulfate) were constructed by molecular docking. The substrates were located inside the constructed active center ensuing a model of its crystal structure. A modified AMBER force field was used to assign the potentials for both the enzyme and the tetrasaccharide substrates. The starting structural models of the enzyme-substrate complexes were subjected to energy minimization. The major part of the protein was fixed, and only amino acids that were part of the active site subset were allowed to move during the minimization. The procedure did not distort the ring conformation of the monosaccharides. Four structurally conserved amino acids: His-501, Tyr-508, Arg-560 and Glu-653 involved in the catalytic activity of chondroitinase ABC I were identified. Low-energy enzyme-substrate complexes were obtained by systematic energy minimization methods where the active site amino acids and the substrates were allowed to move freely. Unfavorable steric contacts were removed by an initial energy minimization to obtain good starting structure for the enzyme-substrate complex. This structure was further subjected to total energy minimization to obtain local minimum energy of the structural complex [21]. Its model suggests that catalytic residues in the enzyme are positioned to cleave chondroitin sulfate more favorably than dermatan sulfate.
Dermatan sulfate is the sole GAG substrate of chondroitinase B [22]. Wide range of pH optimum for the enzyme activity hampers determination of the precise role of the active site amino acid residues. Dermatan sulfate tetrasaccharide was used for docking. The initial orientation of its structure relative to chondroitinase B was obtained by superimposing the non-reducing end of the tetrasaccharide onto the disaccharide in the co-crystal structure of chondroitinase В with a disaccharide product of dermatan sulfate degradation. All the manipulations of the structures and docking were done using docking modules of INSIGHT II and the AMBER force field modified to include carbohydrates. Optimal orientation of the tetrasaccharide substrate in the active center of the enzyme with reasonably low steric hindrance was selected for further energy minimization. Docking and energy minimization resulted in repositioning of the tetrasaccharide substrate to achieve maximum contact with the active site cleft of the enzyme. In final orientation the tetrasaccharide completely occupied the -2, -1, +1, and +2 subsites of its active site [22].
Enzyme activity of hyaluronidases is directed towards degradation of hyaluronan (a sole nonsulfated GAG of endothelial glycocalyx) and, with lower efficiency, of chondroitin and chondroitin sulfate [23]. Molecular docking of 3D model of bovine testicular hyaluronidase with chondroitin sulfate trimers (hexosaccharides) and heparin tetramers (octasaccharides) revealed eight binding sites at which free binding energy for the ligands is at least 2-fold higher than that of free ligands ( Figure 1) [24]. Conformational mobility of hyaluronidase (up to denaturation temperature) was evaluated based on its substantial domain stability. The time of dynamics was determined as described [18,24] during conformational changes [25]. If changes in energy considering temperature-dependent fluctuations reach a plateau, calculations can be stopped. If energy continues changing, calculations should be performed until new plateau is reached. For hyaluronidase docking these were fifteen 100-ps intervals [24]. New plateau was reached after 1.2 ns. Reversibility/irreversibility of conformational transitions of bovine testicular hyaluronidase determined by the return to the initial within a given time period (50-100 ps) indicates that the protein is stable (i.e., there are no heat-induced conformational oscillations). Thus, electrostatic forces of GAG ligand interaction with 3D model of hyaluronidase induce reversible and irreversible conformation transitions in the enzyme leading to its stabilization or inactivation. Effective modifications of a biocatalyst lay the foundation for production of its stabilized forms and subsequent use in medical practice. It should be noted moreover that binding of chondroitin ligands at ch6, ch3 and ch1 sites of hyaluronidase stabilized the enzyme, increasing its denaturation temperature by 10 -15°C [26].
Fine influence of disaccharides on hyaluronidase properties confirmed with presence of two C-mannosylation sites of enzyme (Thr-130 and Thr-321). The polar mannose is attached to the non-polar tryptophan residue (Thr-130), that induced the conformational changes of biocatalyst molecule with negative regulation of its secretion and attenuation of its activity [27].
The feature of chemical modification of hyaluronidase by chondroitin sulfate (low molecular mass /30-50 kDa/ and high molecular mass /120-140 kDa/) is the preservation of appreciable remain endoglycosidase activity of enzyme (68-78%) after its deep modification (degree of aminogroup modification 82-98%) [28]. The data of theoretical investigation showed that diverse modification of biocatalyst by chondroitin sulfate (in respect to its lysine residues) don't lead to substantial alteration of accessibility of hyaluronidase active site for its substrate ( Figure 2) [29].
Making 3D structure of covalent hyaluronidase-chondroitin sulfate conjugate demonstrated the achievement of blockade altogether surface lysine residues due to coupling two oligomeric chondroitin sulfate fragments (total 320-480 saccharide rings). Random topology of building the complete chondroitin sulfate polymer (one of 18! factorial) allows to obtain the 3D structure of hyaluronidase-chondroitin sulfate conjugate when the enzyme is enclosed the chondroitin sulfate coat fully almost ( Figure 3). Two sites of biocatalyst globule (without surface lysine residues) are exclusion. One of such sites is the area of enzyme active site. The obtained results are agreed with appreciable remain endoglycosidase activity of hyaluronidase after its deep modification with chondroitin sulfate of different molecular mass [29].
Thus, the theoretical analysis of 3D structure of hyaluronidase allows forecast the opportunity of different chemical modification of all extended to solvent volume the lysine residues of this protein without appreciable decrease of enzyme activity, if the reaction conditions don't contribute the interaction of masked into enzyme globule of lysine residues. At present it should be noted that polysaccharides are studied actively and proposed for extension of doctor arsenal of drugs on the base polymer-protein conjugates [30,31]. Recombinant human hyaluronidase (pegvorhyaluronidase alfa) is used successfully for clinical trials in respect to treatment of pancreas cancer [32]. Named examples demonstrate that modification of biocatalysts has been conduced their conversation to contemporary effective drugs.

Computational Approaches in the Study of Glycosaminoglycan-Protein Complex Formation
Computational methods have found wide application in the investigation of interactions between GAG and proteins. GAG structure-function relationships illustrate the specificity of complex formation between GAG and proteins, such as growth factors, antithrombin, cytokines and cell adhesion molecules [33]. The complexity of GAG-protein interactions is based on conformational flexibility and underlying sulfatation patterns of GAG, the role of metal ions, and the effect of pH on the affinity of binding. The structure of GAG interactions with proteins, GAG-binding sites, and 3D structure of GAG have been computationally modeled to understand the mechanisms underlying these processes. Successful development of GAG-based drugs is determined by the relationship between GAG structure and its interactions with proteins. Calculations of free energy of binding between heparin fragments of varied length and platelet-endothelial adhesion molecule 1 (PECAM-1) have revealed a region of low GAG-binding affinity in domains 5-6 of PECAM-1 and a region of high affinity in domains 2-3, which is consistent with experimental data and ligand-protein docking studies [34]. A conformational movement observed between domains 2 and 3 allows binding of heparin fragments of increasing size (pentasaccharides to octasaccharides) with an increasingly higher binding affinity. In general, free energy calculations show that the binding of heparin to protein surfaces is dominated by strong electrostatic interactions for longer fragments with equally important contributions from van der Waals interactions and vibrational entropy changes against a large unfavorable desolvation contribution due to high charge density of these molecules. When combining docking simulation with cluster analysis to extract adequate docking structures from the many possible output structures the AUTODOCK 3.0 was used to predict the structure of basic fibroblast growth factor (bFGF) bound to heparin [35]. The energy minimization (calculated by AMBER8 or AUTODOCK 3.0) could not be used satisfactorily by themselves to select a proper heparin-binding complex from the output structures. Their majority generated by AUTODOCK 3.0 was fairly close to each other in atom geometry and in averaged geometry and was close to that of the crystal form of the complexes. Exact prediction of heparin-binding structures of these proteins (basic fibroblast growth factor, antithrombin and annexing V) shows that the approach used in this study is effective in the docking of ligands that have a variety of conformations due to the presence of multiple rotatable bonds and charged chemical groups [34]. Two putative heparin-binding peptides have been identified near the C-and N-terminal regions of promatrilysin (matrix metalloproteinase 7, MMP-7). However, molecular modeling suggests a more extensive binding site or cradle from crossing multiple peptide filaments [36]. Binding of MMP-7 through heparin/heparan sulfate with cell surface facilitates its direct proteolytic attacks, activation of other MMP or regulation of cell surface proteins. The research emphasizes the importance of MMP anchoring to cell surface in controlling tissue remodeling. Molecular docking and molecular dynamics can be effectively used in the study of these processes.

The Study of Protein-Glycosaminoglycan Interactions
When highly purified GAG are not available for experiments the molecular docking becomes relevant and reasonable method for theoretical studies. This method is used as a computational tool in prediction and simulation of preferable orientation for GAG and/or their fragments and for protein molecules (with calculations of minimal energy for GAG binding to the protein molecule surface) [6]. Although GAG-protein binding is based predominantly on electrostatic interactions, the contribution of van der Waals and hydrophobic interactions to GAG-protein complex formation should not be ignored. It is noteworthy that GAG-protein binding can be specific and nonspecific, depending on the negative charge of GAG and positive charge of protein domains. Thus, GAG-protein interactions are characterized by different levels of specificity [6,37]. The binding can complicate further by desolvatation, entropic contribution and similar components. Two participants of binding are GAG carbohydrate chain and protein amino acid chain. GAG-protein interaction is much more difficult for elicitation than protein-protein interaction, since 3D structure of GAG polysaccharide is not as rigid as that of protein [6]. The performance conditions pose a significant restriction to docking when molecular surface of protein remains rigid, which could prevent consideration of the effects of induced-fitting on the protein binding sites [17]. The dominant driving force of denaturing deformation of bovine testicular hyaluronidase beetles the difference in potentials between the positively charged middle part of enzyme molecule and the negatively charged areas of the right and left periphery regions ( Figure 4) [38]. The charge of biocatalyst molecule (according to reckonings, pH 7.5) is + 8.5 a.u. (atomic units) and a dipole moment of 631 D with the vector headed up. Hyaluronidase should be denatured already at 320 K in the absence of the substrate or other GAG oligosaccharides. At low temperature under the effect of Na and Cl ions the difference of potentials decreases considerably and hyaluronidase molecule maintains structure even in the absence of substrate. The increase of temperature leads to molecular deformation due to an increase in ion mobility and to enzyme inactivation.
Screening of virtual ligands to predict the structure of carbohydrate-protein complexes is reasonable with the use of molecular dynamics simulation and automated ligand docking with µs time scale. Longer simulation times improve the reliability of comparison between calculated and experimental data with considerable prolongation of calculations. It is important to select simulation conditions and force fields providing the best agreement with experimental practice. Force fields for carbohydrates can be consistent with protein parameters [17]. At the current stage of development of computational methods, the data (obtained with relatively short chains of GAG oligosaccharides [6] and small-size binding sites of proteins) demonstrate higher reliability, which facilitates more rapid calculations [16]. If protein structure is far from being stable or has destabilizing contacts in some segments the modeled protein undergoes fundamental structural changes during molecular dynamics. The dependence of potential energy of protein structure on time demonstrates a gradual energy decrease which reaches a plateau after certain time period varying over the base line/value. A long plateau indicates that the system is energetically stable after the given time period and protein structure remains stable after reaching equilibrium. Long-time simulation [18] or determination of time interval after which test parameters (e.g., energy, upon its minimization) [24] stop changing are necessary to assess stability of modeled protein and to test potential structural tensions including conformational changes in side chains of the protein. Protein-ligand binding processes can occur via free energy redistribution. The understanding of cause-and-effect relationship with subsequent dynamics of the system (as function of time) is fundamentally important for complete determination of function of a complex system and progress in molecular design [39]. Data on short-time scale dynamics are implicated for long-time scale processes which are exemplified by biocatalysis. Fluctuations of adenylate cyclase (AdK) on a ps scale are related to µs and ms dynamics. ThermoAdk (hyperthermophilic enzyme homolog) and mesoAdk (mesophilic homolog) demonstrate different turnover rapidity at the same temperature. A comparative analysis of dynamics for these enzymes (which allows identification of atomic fluctuations) has revealed their critical importance for enzyme activity. The correspondence of ps and ns loop activities flexibility between mesophilic and hyperthermophilic Adk coordinates the relationship between local short-time dynamic scale and more slow global dynamics of the entire protein molecule. The conceptions of structure, binding, energy movement and entropy over these time scales are an interrelated and united with interdisciplinary approach to gaining detailed insight into the interplay occurring at a molecular level.
Obviously it is extremely difficult to achieve this. Perturbations on short time scales (fs to ps) lead to ensuing functional dynamics on multiple longer time scales (with diapason in several orders of magnitude) [39].
It should be noted that in reality an individual protein molecule (globule 3-4 nm in diameter) interacts with an oligosaccharide within ps time intervals, while elementary reactions occur within fs (synthetic) or tens fs (dissociative). Under experimental conditions numerous molecules interact at different times which accounts for the observed duration of reactions (minutes). The process occurs within ms (but not s) for long molecules (whose parts can interact independently) since the other partaker of the reaction is not selected in the incubation volume but in the sphere whose diameter is equal to the polymer length. Computer calculates simulation time (interaction time) after the test when protein molecule meets the chosen oligomer [24]. The calculated time period is measured in ps, i.e., this is the time after encounter of protein molecule with oligosaccharide which binds to the best fit binding center of protein. If binding occurs at several centers on the protein the dynamics of the protein structure is initiated with all these centers. In this case all molecules (protein, oligosaccharide (s), water and salt) are allowed to move in compliance with the temperature conditions. Dynamic equilibrium (with fluctuations) is generally reached after 50-100 ps (this is the stage of molecular dynamics at physiological temperatures). If the protein conformation changes reversibly or irreversibly after binding with the ligand the process duration increases by one or two orders of magnitude. The protein structure does not change globally while amino acid residues in the attached ligand proximity relocate to provide strong binding with relaxation times ranging from tens to hundreds ps. It should be noted the changes in distances between certain of amino acid residues in the 3D structure of free hyaluronidase (at 320 K in dependence on computation time of observation, ps) ( Figure 5) [38]. Highly appreciable the translation of Glu-105 and Arg-59 was (i.e., between the peak of enzyme chain with Glu-105 at its top and the positively charged region around Arg-59) representing an important initial stage of computed denaturing of biocatalyst. Distance between other amino acid residues of hyaluronidase (residues Asp-147 and Glu-149 responsible for catalytic activity located in enzyme active site and amino acid residues beside and around them, i.e., Asp-147 and Asn-150, Glu-149 and Lys-162, Trp-148 and Glu -149) turned out to be less labile, but relative orientation of the noted amino acid residues changed considerably ( Figure 5). These data have been showed significant distortion of hyaluronidase active site and complete stoppage of the entrance for enzyme substrate, which definitely predicts the inactivation of the biocatalyst. Simulation of complex structure for bovine testicular hyaluronidase with the chondroitin ligands detected the its stabilization at 320 and 340 K compared to free biocatalyst structure. Binding of the chondroitin ligands at sites 6, 3 and 1 of hyaluronidase ( Figure 4) increased the virtual denaturing temperature of the enzyme, considerably slowed down the process of biocatalyst structure deformation, and changed the nature of conformational translations in the enzyme molecule [38].
It is noteworthy that the initial state is determined by docking but not by dynamics. Protein and partially ligands are considered as rigid structures. The ligand molecule is artificially unevenly moved in the protein proximity (depending on the algorithm) and binds to the protein center where its binding is energetically more preferable than that to the solvent while the ligand, water and salt are capable for movement in contrast to the protein molecule which remains stationary. Calculations are repeated several times to reveal all potential binding sites. It should be noted that prior to the initiation of local dynamics docking is performed without solvent. There is no mobile solvent else, continuous medium similar to the solvent is imitated. Only after that protein molecule is surrounded with water and salt so that the total charge of the system, including sodium and chlorine ions, is equal to zero. Binding's centers on the protein are ranged according to the energy of the ligand removal. When the bound ligand is deleted from the protein the vacant site can be occupied by water and salt. The difference between the system energy and the ligand energy in the solvent with salt is calculated. A database of protein structures with different locations of ligands (ranging from the strongest to the weakest binding) and their combinations is created. The structures are by turn multiply subjected to molecular dynamics evaluation at 305-310 К.
Traditional approach with using of experimental and theoretical data of investigation has based on the fulfillment of molecular docking and molecular dynamics for interaction of protein/enzyme with GAG ligands. The subsequent experimental verification of obtained computational results implies the achievement of corresponding concordance (possibly consecutively mutually corrected) between these suites of materials. The other way is likely thoroughly and consists in the initial experimental study of protein-GAG interactions and elaboration of mechanism of their implementation with computational manners. Latter has been contributed the optimum experimental GAG reaction on protein object for breakthrough of determined and nominative aims of work.
Analysis of molecular dynamics implies that the differences between cold-adapted protein to its heat-adapted analog are most likely based on the flexibility of the external protein part (upon protein-water interaction) than on the flexibility of active center. In other words, the rigidity/flexibility of protein surface allows optimum adjustment of thermodynamic parameters. A seamless fusion of experimental and calculated data will provide a strategy for their efficient use in further studies.

Molecular Docking and Scoring Function
An active search for common criteria of the choice of docking algorithm and scoring function has been performed on the basis of binding parameters of specific protein targets. Docking programs consist of two major parts: search algorithm for conformational, rotational and translational space available for ligand molecules on the protein binding centers and an objective scoring function which minimizes the process duration [16]. This function is calculated by rough measurements of binding affinity (or receptor-ligand complementarity) and serves two purposes: the first, to differentiate multiple positions of a single ligand in the binding center of protein/receptor; the second, to determine binding affinity and to improve classification of test compounds. The crucial role of a specific category of docking estimation has determined a great number of scoring functions which were divided into 3 classes [16]. By broad application and considerable number of citations the empirical scoring function (the first class) is regarded as the most important. It approximates free binding energy and describes various types of interactions between protein/receptor and ligand. The second class of scoring function is based on force fields of molecular dynamics. Binding affinity is determined by the sum of electrostatic and van der Waals energies of interaction between protein/receptor and ligand. The third class includes scientific (knowledge-based) scoring function based on statistical analysis of experimentally (rentgenostructurally) determined protein-ligand structures. The third class scoring function is characterized by unpredictable results of initial virtual screening and difficult explicability of failures and lacks the ability to be improved directly. These properties limit the use of this scoring function [16]. Scoring function for interactions with highly charged ligands (GAG) has low or zero specificity [17]. Selection of docking algorithm specific for the target can increase the effectiveness/realization of virtual screening. However, docking algorithm cannot be selected independently on scoring function incorporated in the molecular docking software. Moreover, none of the currently available scoring functions can be applied with the same quality to all types of protein centers even at correct positions of docking for active compounds in a randomly composed library. Scoring function is more important for general control over the process than docking algorithm [16]. Modern docking tools are quite reliable for common use in pharmaceutical industry.
Virtual screening is a more accurate, effective, less expensive and rational way for drug search than conventional high-performance experimental screening. Recent studies have turned molecular docking into a tool of increasing importance for pharmaceutical research. The purpose of molecular docking is to predict the structure of the protein-ligand complex (approbation of ligand conformations in active center and their ranging with the use of scoring function). It should be noted that with six translational and rotational degrees of freedom similar to that with conformational degrees of freedom there is a great number of potential types of binding between these molecules [40]. However, calculation of all conformations is very expensive. Docking of a flexible ligand with a rigid protein/receptor has been performed for a long time and still remains the major method. It is very important to consider the flexibility of the ligand and protein/receptor, since both protein and ligand change conformation to form a complex compliant with minimal energy. However, when protein structure is considered flexible also, calculation costs increase considerably. This makes common the approach when the ratio between accuracy and calculation time is aimed at flexibility of the ligand with protein remaining rigid during docking. As a result, almost all docking software was adapted to this methodology: docking of a flexible ligand and a rigid protein/receptor. Meanwhile, intrinsic protein mobility is relevant to the process of ligand binding. Therefore, consideration of protein/receptor flexibility poses a challenge in the improvement of docking techniques [40]. Molecular dynamics allows modeling of all degrees of freedom in the ligand-protein complex implying the development of docking in the direction of flexible ligand and flexible protein/receptor. It should be remembered that molecular dynamics creates a problem of inadequate approbation and high cost of calculations which limits its use in the screening of large chemical databases [40]. At present time it should be proposed that our understanding of ligand binding mechanisms has been evolved not only course according to the description of the induced fit (when initial binding stage assumes following conformational changes for optimization of complex between high molecular and ligand components) but with selection by ligand the optimal conformation of target partner existing in multiple equilibrium conformations [41]. Computational methods deserve intense application for achievement of reconcile data due to these approaches and determination of criterions of commanding course of ligand binding. Improved accuracy and moderate/low cost of software seem as prospective ways to boost the effective use of docking. The recommendation of such kind is suggested for follow-up and future investigations on this research topic.
A machine-learning-based scoring function has been constructed [40,42]. This approach considerably improved (in comparison with classical scoring functions) both meaning of estimation (prediction of binding affinity) and sense of classification (prediction of relative hierarchy of test compounds). This is an important achievement, bearing in mind that scoring function is the major limiting factor for efficient virtual screening based on molecular docking. The results ground the suggested strategy of virtual screening using AUTODOCK, AUTODOCK VINA and Le Dock with subsequent repeat evaluation of data with the help of a machine-learning-based scoring function and averaging of the evaluations and sorting of test compounds for their further experimental evidential investigation [42].
Taken together, these data emphasize important aspects in the performance of efficient molecular docking that requires well-defined scorinng function and sufficient computer power. This becomes indispensable for using computational chemistry methods to solve the given task. A double-stage approach (molecular docking with subsequent elaboration of its data by molecular dynamics methods) is a promising tool to the modeling of protein-GAG complexes [7]. Software based on a harmony between docking algorithm and scoring function is used for this purpose. The choice of specific software is associated with the reliability of the results obtained (GAG length, size of protein binding domain, total size of protein molecule, etc.) and the development of new software. Indeed, critical examination of 5 docking tools (AUTODOCK, FRED, CDOCKER, FlexX and GOLD) with the aim of selecting the most appropriate one for docking (with rather shallow and featureless binding site on lectin) has shown that even widely used docking programs have certain limitations [43]. Optimal orientation of molecules is crucial for formatting coordinated bonds which stresses the importance of the ability of selected docking tools to reproduce optimal binding conformation of the protein interacting with oligosaccharides. More than fifty docking tools differing in algorithm and scoring function are currently available. They all require 3D structure of target protein and ligand. It should be noted as instructive corollary that future improvements of docking tools are associated with removal of current limitations, such as consideration of protein molecule in the majority of cases as a rigid structure, and partial speculation or full ignorance of the contribution of desolvation which is very important in docking of highly solvated ligands similar to oligosaccharides [43]. Qualified preparation and efficient performance of multi-aspect docking researches are impossible without organized cooperation between specialists in computational methodology and life science profile.

Conclusion
Elucidation of the mechanisms (responsible for interactions between the components of biological systems, including high-molecular-weight compounds) is the need for well-grounded recommendations concerning experimental modification of these compounds. The further development of technologies for pharmaceutical production of macromolecular drugs has been builded on naturally occurring biologically active substances this one is conducive the broad use of computational chemistry techniques. Transition from the in silico concepts to production of novel drugs is a clear tendency in modern biopharmaceutics [44]. The ponderable challenge of molecular docking progress is demonstrated as using of notion of enzyme structure flexibility due to a mode of molecular dynamics with modeling of all freedom degrees in enzyme-ligand complex. An important necessity of increasing computer power, improving scoring functions, developing new docking software, considering the effects of salts and pH on molecular interactions of biological agents creates a noticeable gradient towards acceleration in the development of computational techniques to study molecular interactions in biological systems. Rapid improvement of computer technologies generates the task of creating highly-coordinated sophisticated automated systems for diagnostics, monitoring, therapy and prevention of diseases. Development of molecular docking of enzyme with glycosaminoglycans promotes progress in this direction with harmonious verification of computational data by experimental methods.

Funding
The study was supported by Ministry of Health Care of Russian Federation and Russian Foundation for Basic Research (grant 18-015-00056).