**1. Introduction**

Lubrication in natural joints is a complex multiscale process that involves interactions between constituents of articular cartilage and synovial fluid [1–3]. Although two main mechanisms—hydration repulsion and molecular synergies—have been found, the atomistic details on phenomena are still studied due to the lack of knowledge of how the cooperation between them results in facilitated lubrication [4–8]. Synovial fluid consists of up to 80% of water, macromolecular and small-molecular components.

The most recognizable components are: albumin, hyaluronan, phospholipids, *γ*globulin and lubricin [9]. Human serum albumin (HSA) deserves special attention due to its binding and transporting properties of various compounds (fatty acids, ions: K+, Na+, Mg2+, Ca2<sup>+</sup> and many other molecules [10,11]). It was demonstrated in [12,13] that albumin and *γ*-globulin play an important role in lubrication. However, their influence on lubricating properties starts to be vital when taking into account their cooperation with other synovial fluid components. Locally positively charged sites of albumin favor

**Citation:** Sionkowski, P.; Bełdowski, P.; Kruszewska, N.; Weber, P.; Marciniak, B.; Domino, K. Effect of Ion and Binding Site on the Conformation of Chosen Glycosaminoglycans at the Albumin Surface. *Entropy* **2022**, *24*, 811. https://doi.org/10.3390/e24060811

Academic Editor: Antonio M. Scarfone

Received: 11 May 2022 Accepted: 8 June 2022 Published: 10 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

interactions with the ionized carboxylic and sulfate groups in glycosaminoglycans (GAGs). Even though both macromolecules have a global negative charge under physiological conditions, it was shown that HSA could bind to negatively charged surfaces [14].

This study presents the analysis of the interactions between HSA and hyaluronic acid (HA)/ chondroitin-6-sulfate (CS6) in terms of the influence on conformations. As the topic is extensive, we are skipping here the conformation changes within albumin, focusing only on the GAGs. Such analysis of interaction between HSA and HA has been performed recently in [15]. HSA consists of a single chain of 585 amino acids, incorporating three homologous domains (I, II, and III) [16]. Domain I consists of residues 5–197, domain II includes residues 198–382, and domain III is formed from residues 383–569. Each domain comprises two subdomains termed A and B (IA; residues 5–107, IB; residues 108–197, IIA; residues 198–296, IIB; residues 297–382, IIIA; residues 383–494, IIIB; residues 495–569), these are depicted in Figure 1. Some of the domains are more prone to binding to GAGs than others; however, the binding map can alter under some conditions of various diseases [17,18]. The binding mechanism is mainly due to ionic, hydrogen bonding and hydrophobic interactions [17].

**Figure 1.** Structure of (**a**) HSA-HA, (**b**) HSA-CS6 complexes for highest affinity result in CaCl2 solution (solution is transparent on the picture). HSA is depicted as ribbon (bottom parts of picture), and its domains are colored as follows: IA-pink, IB-violet, IIA-light green, IIB-green, IIIA-light blue, and IIIB-blue. HA and CS6 are depicted as ball-stick (top parts of picture). Light blue atoms represent carbon, dark blue nitrogen, red oxygen, green sulfur and white hydrogen. Snapshots was taken using YASARA software after 100 ns MD simulations [19].

GAGs are large complex carbohydrates. Depending on the monosaccharide types and the glycosidic bonds, GAGs can be divided into four groups: (i) hyaluronic acid, (ii) chondroitin sulfate, and dermatan sulfate (iii) heparan sulfate and heparin, and (iv) keratan sulfate. First, let us underline that HA is only a non-sulfated GAG. It is vital because sulfate groups in the GAG is one of the most crucial factors influencing the interaction map between a protein and the GAG [20]. Thus, two different examples of GAGs have been picked, which are essential components of synovial fluid, to study their conformational entropy while they are in the vicinity of the albumin: CS6 and HA. CS6 was shown to be an excellent material for bone regeneration as it is the main constituent of glycosaminoglycan in cartilage. CS6 is involved in bone homeostasis and in the coordination of osteoblastic cell attachment. Kim et al. investigated the role of chondroitin sulfate's negative charge on the binding of cations (e.g., calcium and phosphate) and showed that the hydroxyapatite crystal formation was enhanced, accelerating osteogenesis.

The analysis of the interaction between HSA and HA or CS6 in the presence of various species such as water and ions is a meaningful task; as such, interaction is closely related to synovial systems' unique properties [4]. The use of molecular modeling allows for researchers to evaluate the influence of various factors, such as the presence of ions and solvation on the properties of proteins, including their ability to bind ligands [21]. Following this, the present study aims to evaluate the effect of Na+, Mg2+, Ca2<sup>+</sup> cations on the affinity of this two specific GAGs to HSA using: firstly docking, and secondly molecular dynamics (MD) methods. The studies are focused on the description of changes in conformation of the GAGs in the vicinity of the HSA.

GAGs are important complexes that participate in many biological processes through the regulation of specific proteins. Hence, their secondary structure and stability are very important to study. Both of above-mentioned properties can be well quantified by conformational entropy. As conformational entropy, we understand the Shannon entropy computed for bivariate histograms of chosen pairs of dihedral angles (see, Section 2.1) [22]. Variations of such entropy are considered to measure important properties of the biochemical processes [23,24]. In this research, we also refer to the informative interpretation of entropy. Basically, one can compute this entropy by studying the motion of the molecule. There are various methods dedicated to characterize these motions (e.g., NMR relaxation methods [25], AFM-unfolding [26], and neutron spectroscopy [27]). Theoretical studies based on classical mechanics approach are a reasonable alternative but computationally demanding [28,29]. The intermediate approach, we follow, is the computation of conformational entropy from all-atom MD simulations, see [23,30,31].

In more detail, conformation description of the GAGs relays on the analysis of their structures bound to HSA domains in aqueous ionic solution. This analysis is carried out to check whether there are any differences in the conformation of the glycosidic linkages between each oligosaccharide monomer of the GAG, when the kind of the ion is changed in aqueous solution. The linkages are investigated basing on specific dihedral angles. In the present paper, conformational entropy is computed from the frequency distribution of those angles' values. We anticipate that those angles determine important characteristics, such as shape and stiffness. As the conformational entropy is calculated from the distribution of the angles, it is expected to be a crucial feature (enclosing the quantitative description in one numerical value).

#### **2. Materials and Methods**

We have performed all-atom simulations of the two model biosystems (one is HSA with HA and the other is HSA with CS6) in aqueous ionic solution. First, a molecular docking procedure has been executed to obtain preliminary information on the stability of the structure and to find the most energetically optimal places where each GAG attaches to the HSA. Next, energetically best-docked structures (sorted from the strongest connection to the weakest connected), with added water solution of chosen ions (Na+, Mg2+, Ca2<sup>+</sup> and Cl−), have been subjected to MD simulations.

Chemical structures of HA and CS6 were obtained from Pubchem and modified to obtain chains of around 8 kDa (24 units). This modification relied on connecting units of selected GAGs until polymers of desired length were obtained. To acquire the most stable complexes, we docked GAG ligand (HA or CS6) to HSA using the VINA method [32] with their default parameters and point-charge force field [33] initially assigned according to the AMBER14 force field [34] (the HA molecule was parametrized by applying the GLYCAM06 force field [35]). Then, we damped the system to mimic the less polar Gasteiger charges used to optimize the AutoDock scoring function. The simulation was done with the YASARA molecular modeling program [19]. In each case (HA and CS6), 10 of the best distinctive complexes which differs in the position of GAG docked to HSA (best complexes means the complexes of the highest energy of binding and RMSD of complexes from the best binded complex values, which were computed by VINA) with −10 kcal/mol free energy of binding prepared for MD simulation.

MD simulations of HSA (PDB code: 1e78) with GAG have been run with YASARA software. Optimization of the hydrogen bonding network was included in the setup to increase the solute stability and a pK*a* prediction (based on a Henderson–Hasselbalch equation) to fine-tune the protonation states of the protein residues at the given pH = 7.4 [36,37]. Optimization was done in three main steps: (a) pKa prediction was included to consider the influence of the pH on the hydrogen bonding network, (b) nonstandard amino acids and ligands were fully accounted for with the help of a chemical knowledge library in SMILES

format, and (c) the use of the SCWRL algorithm allows for finding the globally optimal solution. In the case of the HSA-HA system, the complex has been immersed in one of the three aqueous 0.9% salt solutions NaCl, CaCl2 or MgCl2. In the case of HSA-CS6, it was 2% water solution of the same salts in pH = 7.0. After necessary minimization of the model system to remove clashes, the simulation was run for 100 ns using the AMBER14 force field [34] for the HSA, GLYCAM06 [35] for HA and CS6, and TIP3P for water. The van der Waals forces' cut-off distance was set to 10 Å [38]. The particle Mesh Ewald algorithm was used for computing long-range interactions (e.g., electrostatic interactions) [39]. Simulations were performed under following conditions: temperature 310 K and pressure of 1 atm (NPT ensemble) [37]. Periodic boundary conditions were applied to a box of size roughly equal to 120 Å3. Berendsen barostat and thermostat were used to maintain constant pressure and temperature (relaxation time of 1 fs) [40]. The equations of motions were integrated with multiple time steps of 1.25 fs for bonded interactions and 2.5 fs for non-bonded interactions. In the considered simulations, the time step between saved states of the system equals 100 ps. Thus, the time series for 100 ns of simulations obtained 1000 save points. Snapshots of the two complexes after 100 ns of MD simulation have been presented in Figure 1.

All analyses and computation have been performed using YASARA and in-house written data processing programs in Python 3.8 [41].

#### *2.1. Backbone Angles Determination*

The method of entropy calculation, used in this study, relies on computation of the frequency distribution of the backbone's dihedral angles (Φ,Ψ), usually presented in the form of a Ramachandran-type plot [31,42]. It provides a simplistic view of the conformation of a molecule by clustering angles (Φ,Ψ) into district regions (bear in mind that there are two sets of such pairs of angles, which will be discussed later).

The CS6 consists of glucuronic acid (GlcA) and galactosamine (sulfated at C-6 atom of galactosamine; GalNAc), while HA consists of glucuronic acid (GlcA) and acetylglucosamine (GlcNAc). Linkages between the two monosaccharides are as follows: in the case of CS6, it is [4)-*β*-GlcA-(1−→3)-*β*-GalNAc-)1−→], and for HA: [4)-*β*-GlcA-(1−→3)-*β*-GlcNAc-)1−→] (see, Figure 2). Using abbreviation G for GlcA and N for GlcNAc (in HA) or GalNAc (in CS6), we can depict the glycosaminoglycans as linear heteropolysaccharide chains consisted of repeating disaccharide units [31]. In the presented study, the chains consist of NG units repeated 24 times. Two sets of torsion angles describe the conformations around the glycosidic linkages: Φ1−<sup>4</sup> and Ψ1−<sup>4</sup> (N-G linkage), and Φ1−<sup>3</sup> and Ψ1−<sup>3</sup> (G-N linkage) [31]. Thus, in the GAG chains, there are 24 of Φ1−<sup>4</sup> and Ψ1−<sup>4</sup> angles, and 23 Φ1−<sup>3</sup> and Ψ1−3. These angles can be written as follows:

$$\begin{aligned} \Phi\_{1-4} &= \text{O5(N)-C1(N)-O1(N)-C4(G)},\\ \Psi\_{1-4} &= \text{C1(N)-O1(N)-C4(G)-C5(G)},\\ \Phi\_{1-3} &= \text{O5(G)-C1(G)-O1(G)-C3(N)},\\ \Psi\_{1-3} &= \text{C1(G)-O1(G)-C3(N)-C4(N)}.\end{aligned} \tag{1}$$

All dihedral angles in Equation (1) have been presented in Figure 2. The available conformational space of the GAGs' chains depends mainly on the two torsion angles. As mentioned before, a goal of the present study is to investigate how various features of the system (ions, HA vs. CS6) affect the frequency distribution of Φ and Ψ torsion angles. As the informative feature of above-mentioned frequency distribution, conformational entropy has been used [24,31]. Note that, in information theory, entropy is the measure of the information tied directly to the variable (univariate or multivariate) [22,43]. Using this analogy, we can analyse which pairs of dihedral angles are more or less informative for specific GAGs and ions.

**Figure 2.** Structures of (**a**) HA and (**b**) CS6 with dihedral angles: Φ1−4—red circles, Ψ1−4—blue circles (N-G linkage), Φ1−3—green circles, Ψ1−3—violet circles (G-N linkage).
