The eﬀect of water on the binding of glycosaminoglycan saccharides to hydroxyapatite surfaces: a molecular dynamics study †

Classical molecular dynamics (MD) simulations have been employed to study the interaction of the saccharides glucuronic acid (GlcA) and N -acetylgalactosamine (GalNAc) with the (0001) and (01 % 10) surfaces of the mineral hydroxyapatite (HAP). GlcA and GalNAc are the two constituent monosaccharides of the glycosaminoglycan chondroitin sulfate, which is commonly found in bone and cartilage and has been implicated in the modulation of the hydroxyapatite biomineralization process. MD simulations of the mineral surfaces and the saccharides in the presence of solvent water allowed the calculation of the adsorption energies of the saccharides on the HAP surfaces. The calculations show that GalNAc interacts with HAP principally through the sulfate and the carbonyl of acetyl amine groups, whereas the GlcA interacts primarily through the carboxylate functional groups. The mode and strength of the interaction depends on the orientation of the saccharide with respect to the surface and the level of disruption of the layer of water competing with the saccharide for adsorption sites on the HAP surface, suggesting that chondroitin 4-sulfate binds to the layer of solvent water rather than to HAP. of solvent. The results of the MD simulations are in a good agreement with previous DFT calculations on GlcA, 25 showing the same interactions and very similar interaction distances and distortions in the GAG structure. The same study analysed a non-sulfonated GalNAc, so our results are not directly comparable, but any trends are similar. The discrepancies between the adsorption energies of the biomolecules obtained with DFT in vacuum and those reported in this work are due to the high stability of the GAGs in water, as well as the strong aﬃnity of the hydroxyapatite surfaces for water. Both factors contribute to making the adsorption of the GAGs onto the HAP surfaces less favourable.


Introduction
Mineralised tissues such as bone and tooth enamel have a complex composite micro-architecture, in which an underlying matrix of fibrous organic material becomes strengthened by the deposition of inorganic apatite crystals. This biological apatite is similar in stoichiometry and structure to hydroxyapatite (HAP), 1 Ca 10 (PO 4 ) 6 (OH) 2 , but with crystal defects such as substitutions by carbonate, water and various metal cations. 2 The microscopic structures of these mineralised tissues are ultimately regulated by cells such as osteoblasts, ameloblasts and odontoblasts (in bone, tooth enamel and dentine, respectively), which are able to control the growth of apatite crystals to a specific size, shape, orientation, chemical stoichiometry, and spatial distribution. It is remarkable that cells achieve this level of control when they can only manipulate the apatite via an indirect approach: by secreting various biomolecules into the surrounding tissue they create chemical conditions that are favourable for the spontaneous extracellular formation of the appropriate apatite structures.
Glycosaminoglycans (GAGs) are polysaccharide molecules that are abundant in mineralised tissues, and are known to interact strongly with apatite crystals. [3][4][5] A GAG molecule is an unbranched polymer chain of two different types of monosaccharides linked together in alternation. 6 The most prevalent GAG found in mineralised tissues is chondroitin sulfate, shown in Fig. 1(a), in which the alternating monosaccharides are glucuronic acid (GlcA) and N-acetylgalactosamine (GalNAc) shown in Fig. 1(b) and (c) respectively. This GAG can be sulfated to some extent on any of its hydroxy groups, and the most commonly sulfated position is the 4-hydroxy group of the GalNAc monosaccharide. GAGs, including chondroitin sulfate, are usually found as part of a larger macromolecular complex known as a proteoglycan, in which multiple GAG molecules are covalently attached to a central protein, and each GAG chain typically comprises 40-120 monosaccharide units. 7 There are a number of ways in which the biomolecules secreted by cells can influence the mechanism of tissue mineralisation: either by providing a template for crystal nucleation, by adsorbing onto crystal surfaces and thus inhibiting crystal growth, or by sequestering dissolved mineral ions. 8 A plethora of biomolecules have been implicated as either promoters or inhibitors of tissue mineralisation, including GAG-containing proteoglycans 9 and the collagen, 10 alkaline phosphatase 11 and SIBLING proteins. 12 It is likely, however, that mineralisation is controlled by many different organic molecules simultaneously, with exquisite control of the process effected by the cell's use of system feedback loops.
It has been found that under in vitro conditions of limited calcium availability, GAGs, including chondroitin sulfate, and their proteoglycans inhibit the deposition of hydroxyapatite, perhaps because the adsorption of GAG molecules on the mineral surfaces is an obstacle to the mechanism of crystal growth. [3][4][5] The role of proteoglycans in the mechanism of apatite formation in vivo is more difficult to interpret, because proteoglycans and GAGs have multiple functions in physiological tissues, and interact with many other biomolecules and with cells. 8,9,[13][14][15] Studies in vivo and in cell culture have often produced conflicting results, with proteoglycans observed either to promote or to inhibit apatite formation. 16 Evidence for the importance of the GAG-apatite interface in vivo has come from solid-state NMR experiments of bones, 17,18 teeth, 19 plaque 20 and kidney stones: 21 in all cases 31 P- 13 C coupling showed that the apatite mineral was predominantly in contact with GAG molecules rather than with any other organic component.
In this paper we have used computational methods to investigate the interfacial interactions between GAGs and HAP. We present molecular dynamics (MD) simulations of chondroitin sulfate saccharides adsorbed onto a slab of HAP in the presence of interfacial water. Adsorption structures have been calculated for the (0001) surface of HAP, which is the thermodynamically most stable surface, and the (01% 10) surface, which is the dominant plane of apatite nanocrystals in a biological environment. 22 We have simulated the adsorption of the monosaccharides GlcA and GalNAc 4-sulfate (GalNAc-4S) independently, rather than longer chondroitin sulfate chains, because this reduces the number of possible HAP-saccharide interactions in each simulation, and simplifies the interpretation of the results. The simulated monosaccharides were methylated in two different positions, as shown in Fig. 1, in order to the reproduce the effect of their position within a long polymer GAG chain.
Classical MD simulations have been used previously to investigate the adsorption of small organic molecules at an apatite-water interface. 23,24 For these molecules it was reported that the principal adsorption interactions were the electrostatic attractions between surface calcium ions and the adsorbate oxygen and nitrogen atoms, and that an important secondary interaction was hydrogen-bonding from adsorbate polar hydrogen atoms to the surface phosphate oxygen atoms. Recently we have used DFT calculations to predict the adsorption structures of GlcA and GalNAc (not sulfated) on HAP in vacuo. 25 It was found that GalNAc interacts with HAP principally through its hydroxy and acetyl amine functional groups, and that GlcA interacts principally through its hydroxy and carboxylate functional groups. However, the inclusion of interfacial water in the present study will give us a better understanding of adsorption processes under more realistic physiological conditions, rather than in vacuo.

Methods
The adsorption of GAG saccharides on HAP surfaces was modelled using interatomic potential-based molecular dynamics simulations in order to calculate the structures and energies of these systems. MD trajectories were calculated using the DL_POLY_2 code, 26 with an integration algorithm based on the Verlet leapfrog scheme. 27 The simulations were performed at constant volume (using cell parameters obtained from constant pressure calculations in the bulk) and constant temperature (310 K) using the Nosé-Hoover thermostat 28,29 with 0.5 ps for the thermostat and barostat relaxation times.
Following conventional MD simulation procedures, the system was described by a simulation cell with three-dimensional periodic boundary conditions, containing a slab of HAP with a surface on either side. The thickness of the slab was greater than 20 Å for each surface studied, and the slab was separated by a gap of approximately 53 Å from its images in the next cell. Under these conditions the system was run for 1200 ps with 400 ps of equilibration. Then the gap between the HAP slabs was filled with water molecules, leading to a simulation cell containing approximately 3000 atoms in total. Now the simulation times were increased to 1600 ps with 400 ps of equilibration. The Shake algorithm was also employed to constrain intramolecular hydrogen-bonds.

Interatomic potential model
A combination of interatomic potential models was employed to describe the interactions of this multi-component system. We have used the Glycam06 forcefield for the saccharide molecules and for non-bonded saccharide-water interactions. The Glycam06 forcefield is an established parameter set optimised specifically for carbohydrate molecules, 30 in which bond lengths and angles are described by harmonic potentials, dihedral angles by sinusoidal potentials, and non-bonded interactions by Lennard-Jones potentials. Intramolecular 1-4 (and greater) electrostatic and Van der Waals interactions are non-zero and are unscaled. The water molecules were described by the TIP3P model, to ensure compatibility with the parameterisation of the saccharide in the Glycam06 forcefield, which includes TIP3P water. The mineral slab was described by the interatomic potential model developed by de Leeuw for modelling apatite crystals. 31 In this forcefield, phosphate and hydroxy group bonds are parameterised as the sum of a Morse potential and a coulombic potential, phosphate bond angles by a harmonic potential, and non-bonded interactions by Buckingham potentials (see Table S1, ESI †). This forcefield makes use of a shell model, 32 in which each oxygen anion in the phosphate and hydroxy groups consists of both a core and a massless shell connected by a spring, in order to model the atom's electronic polarisability. The forcefield has been validated by comparison with DFT calculations, 31 and produces consistent defect and dissolution energies. [33][34][35] HAP-water interactions and HAP-GAG interactions were parameterised using non-bonded potentials derived in previous MD studies, 23,24 which are compatible with the current HAP shell model but do not consider pH effects (see Section 2.2). Specifically, the interfacial interactions were described by Buckingham potentials and Lennard-Jones potentials, and the repulsive portion of the potential was scaled for consistency with the charges on the adsorbate atoms, using a method originally described by Schröder et al., 36 and subsequently used in a number of surface adsorption studies. 23,24,31,37,38 2.2 Surface structures HAP (Fig. 2) is found experimentally to have a P6 3 /m space group, for which the hexagonal unit cell parameters are a = b = 9.432 Å, c = 6.881 Å, a = b = 901 and g = 1201. 1 The mineral contains tetrahedral phosphate groups, hydroxy groups, and calcium cations. There are two different Ca 2+ sites in the HAP crystals: 'columnar' Ca 2+ ions (Ca I ), which are arranged in singleion columns parallel to the c axis and are bound by phosphate tetrahedra, and 'axial' Ca 2+ ions (Ca II ), which form hexagonal columns, at the centre of which are the hydroxy groups.
The crystallographic structure of HAP contains columns of hydroxy groups in the hexagonal channels that are oriented parallel to the c axis. The OH groups are oriented either hydrogen "up" or hydrogen "down" in the channels, where all "up" and "down" OH À positions have occupancy of 0.5. 1 In an MD simulation the occupancy of a site is necessarily restricted to either unity or zero, so in our simulated slabs of HAP the hydroxy groups were placed in only half of the possible positions, such that within each column they all pointed in the same direction. The unit cell was then doubled in the b direction, so that neighbouring columns alternated in the orientation of the hydroxy groups. Previous DFT calculations have identified this arrangement of hydroxy groups as the lowest energy configuration of bulk HAP, 39 and it is similar to the synthetic monoclinic HAP material, where all hydroxy positions are fully ordered in a similar configuration as modelled here.
In calcium deficient environments (i.e. bones, teeth), the apatite obtained is nonstoichiometric, whereas the surface phosphate ions are usually protonated, [40][41][42][43][44][45] with the protonation rate of the ions typically ranging from PO 3 -OH at pH B10.5 to PO 2 -(OH) 2 for pH B5.5. 46 We have used three different HAP slabs for the MD simulations of the adsorption of GAGs, which are shown in Fig. 3. We have considered the low index (0001) and (01% 10) crystal surfaces, which are both present in the experimental morphology of HAP, where for the (01% 10) surface we have considered two different possible terminating planes.
The (01% 10) termination shown in Fig. 3b contains columnar Ca atoms in the surface, with the hydroxy channels remaining intact below the surface, whereas the terminating plane in Fig. 3c cuts through the hydroxy channels and leaves the columnar Ca atoms below the surface. The surface in Fig. 3b has been identified previously as the most stable (01% 10) termination, using computational methods, 47 but the surface in Fig. 3c has been identified as the principal (01% 10) surface of fluorapatite using high-resolution specular X-ray reflectivity data. 48 We note that both of these possible terminations are likely to be present transiently as a crystal grows in the (01% 10) direction.
These three stoichiometric surface structures allow us to compare our results with previous theoretical results 23,24,38,49 while also keeping the Ca/PO 4 ratio, which is important to investigate experimental findings regarding the high Ca 2+ binding capacity of chondroitin 4-sulfate. 50

Adsorption calculations
The monosaccharides GlcA and GalNAc-4S were placed near the HAP slabs in a number of different positions and orientations with respect to the surface. We added one saccharide molecule per surface simulation cell, which was large enough to avoid interactions between it and the saccharide molecules in neighbouring repeat cells.
For each monosaccharide and for each HAP slab, we also simulated the system where the saccharide was fully solvated, i.e. the saccharide was situated in the solution midway between the upper and lower surfaces of two slabs from neighbouring repeat cells. The adsorption energy of a saccharide on a particular surface, E ads , is given by: where E sys,ads is the average configurational energy of the system (including HAP slab, saccharide and water) when the saccharide is adsorbed at the surface, and E sys,solv is the average configurational energy of the equivalent system when the saccharide is fully solvated, i.e. midway between the HAP slabs. A negative adsorption energy indicates that adsorption of the monosaccharide at the surface is thermodynamically favourable with respect to the surface and saccharide being separated in the aqueous environment. All simulations for a particular saccharide-surface combination used the same number of water molecules, so the system energies could be compared directly. The hydration energy of the surfaces is given by: where E surf,solv is the average configurational energy of the hydrated surface, E surf is the energy of the dry surface, E H 2 O is the energy of one liquid water molecule and n is the number of solvent water molecules. The values are divided by the surface area (A) to obtain a value per unit area to allow for direct comparison between different surfaces.

Results and discussion
In order to obtain a plausible initial orientation of the monosaccharides over the three HAP surfaces, while reducing the number of possibilities and consequently the computational cost of the more realistic calculations in aqueous environment, we have first carried out simulations of both surface and GAG systems in a vacuum. The final configurations of the GAGs at the surfaces in the absence of solvent were used as starting points for a more complex scenario in the presence of water.

Surfaces in vacuum
In the absence of solvent the three surfaces evolve in various ways. The first and, to a lesser extent, the second layers of atoms in the surface slabs readjust their positions to re-balance the forces, which have become unbalanced upon the creation of the surfaces from the bulk material. We observe rotations of the phosphate groups and displacements of the topmost calcium ions to regain the maximum possible coordination. The P-O bond distances remain B1.55 Å with the Ca-O distances ranging from 2.15 to 2.88 Å. The (0001) surface is the only surface with hydroxy groups and both types of Ca 2+ in the first layer of the slab (Fig. 3a). The columnar calcium ions (Ca I ) form one extra bond, to a total of six bonds (out of nine in the bulk) with a phosphate oxygen (OPO 3 ), while the axial calcium ions (Ca II ) also form an additional bond with the oxygen of the hydroxy group, which has moved closer to the cations. In this termination of the surface we do not observe the OH groups moving more into the bulk of the material, as was reported after relaxation for a different termination. 24 The flat (01% 10) surface ( Fig. 3b) also has the two types of calcium in the top layer of the slab and they both form an additional bond with the phosphate oxygen. This surface has more calcium in the top layer than the (0001), and therefore shows the greatest rotations of the phosphate groups to enhance their coordination. We also observe that some of the columnar calcium ions move into the bulk to retain the coordination with the phosphate group that has rotated to coordinate to Ca II . In the crenellated (01% 10) surface the phosphate groups experience little rotation to form two extra bonds between their oxygen and the axial calcium ions in the trench of the surface.

Monosaccharides in vacuum
The glycosaminoglycans are polar molecules and highly soluble. We have first simulated the two monosaccharides in the absence of solvent to analyse the conformational changes that might occur once the molecule interacts with the surfaces. Under these conditions we observe that the methylated N-acetylgalactosamine-4-sulfate (GalNAc-4S) remains in the chair conformation, with the acetyl-amine group in an equatorial position and both methoxy groups in equatorial positions. These methoxy groups represent the oxygen-linked attachments to the next saccharides in the GAG polymer (see Fig. 1), whereas the sulfate group is in an axial position forming intramolecular hydrogen bonds between its O and the H of the hydroxy functional group. Methylated glucuronic acid (GlcA), is a smaller molecule and without solvent the conformation of the molecule is more variable than GalNAc-4S. The methoxy and hydroxy groups alternate between equatorial and axial positions, whereas the carboxylate also moves between the equatorial to a pseudo-axial position to form hydrogen-bonds between its oxygen and the hydrogen of hydroxy group.

Adsorption of GlcA on the surfaces in vacuum
The saccharide was placed onto the three surfaces in twelve different initial orientations following chemical intuition. Under these conditions the GlcA molecule often changed its conformation in order to direct all of its functional groups towards the surfaces, even though there is an energy penalty (173 kJ mol À1 ) associated with this change in conformation, which is, however, outweighed by the energetically favourable interactions with the surfaces. In all positions, the carboxylate group is directed towards the surface, with examples of the interactions on the (01% 10) surfaces shown in Fig. 4a and b. Fig. 4 shows that sometimes the carboxylic oxygen is bridging two axial Ca 2+ or one columnar and one axial Ca 2+ on the flat (01% 10) surface (Fig. 4a). The interactions of the carboxylate groups with the Ca 2+ of HAP has been reported previously in the literature, [51][52][53] It binds more strongly to the columnar calcium ions than to the axial calcium ions in the surface, because the columnar calcium ions are more closely spaced in the plane (3.48 Å), and the carboxylate can hence bridge more easily between two columnar calcium ions. For the axial calcium atoms (6.95 Å spacing in the surface) the carboxylate can only bind to one Ca 2+ at a time, which is less energetically favourable. On the (0001) the molecule lies flat on the surface in all cases, reaching the minimum energy value when the molecule adsorbs as shown in Fig. 4c.

Adsorption of GalNAc-4S on the HAP surfaces in vacuum
In the same way as GlcA, the saccharide was placed in twelve different initial orientations on the three surfaces of interest.
In all final configurations, the sulfate group is oriented towards the surface (Fig. 5).
On the (01% 10) surface it can form bridges between two calcium ions, all at distances similar to those shown in Fig. 5b. The oxygen of the carbonyl group is also often oriented towards surface calcium ( Fig. 5a and b).
In the case of the crenellated (01% 10) plane, the sulfate group is initially located in the ''trench'' where the hydroxy group is also interacting with the surface. In all configurations the molecule lies flat on the surface. On the (0001) surface, the hydroxy group is closest to the surface (Fig. 5c) and the sulfate is only coordinating more than one calcium ion in two configurations, whereas the molecule lies flat on the surface in four of the configurations.

Monosaccharides in water
In aqueous environment the GalNAc-4S stays in the low energy chair conformation throughout the MD simulation (Fig. S4a, ESI †) as described in the absence of solvent. However, there are fewer intramolecular hydrogen-bonds between the H of the hydroxy functional group and the O of the sulfate, with all the polar functional groups generally interacting with the surrounding water molecules. The sulfate group is negatively charged and it polarises water (more than the other functional groups), at any time forming multiple water-sulfate hydrogen bonds. There are always typically 2-3 water molecules within 2.5 Å of the sulfate oxygen, but the other polar functional groups are also observed to form hydrogen bonds to water.
GlcA remains in the low-energy chair conformation (Fig. S4b, ESI †), in the presence of solvent. The carboxylate hydroxy and methoxy groups are all in equatorial positions. Intramolecular non-bonded interactions are absent: the polar functional groups generally interact with the surrounding water molecules, rather than with each other. The carboxylate group is negatively charged,  and it therefore attracts the water molecule's hydrogen atoms. There are always typically 1-2 water molecules within 2.5 Å of the carboxylate oxygen and all other polar functional groups are also observed to form hydrogen bonds to water.

Hydrated HAP surfaces
We first investigated the flat (01% 10) surface. A plot of the distances between the oxygen atoms of the water molecules and the surface plane, considered as the topmost calcium atoms, is shown as open circles in Fig. 6a. During the MD simulations, a group of water molecules forms a monolayer on the surface (Fig. 6c), primarily interacting via their oxygen atoms with the surface calcium ions at an average distance of 2.52 Å. Hydrogen bonds are also formed between the oxygen atoms of the phosphate groups and the hydrogens of the water molecules, but few of them are at distances of less than B2 Å. The water molecules do not form a regular adsorption pattern on the surface and a frequent exchange with molecules in the bulk solution is observed, which we can express as an exchange frequency per surface Ca ion. Here B11 times water molecules coordinated to Ca 2+ exchange with bulk water during 1.2 ns of data collection time, i.e. a frequency of B9.2 s À9 . The hydration energy of the surface is À0.73 J m À2 . As already mentioned the (01% 10) surface with a crenellated termination (Fig. 3c) has been identified as the principal (01% 10) surface of fluorapatite using high-resolution specular X-ray reflectivity data. 48 This surface termination cuts through the channels of hydroxy groups. The surface contains trenches and ridges which run parallel to the c direction; in this case (open triangles in Fig. 6a) the water molecules form a more regular arrangement on the surface (Fig. 6d). In particular, the surface trenches are filled with columns of water molecules with a specific orientation, in agreement with X-ray reflectivity data, 48,54 with an average distance from calcium to oxygen of the water of 2.57 Å. In this surface the water molecules in the monolayer prefer to interact more with the surface than with each other forming more hydrogen bonds between the water and the phosphate ions. The frequency of exchanges of water molecules in the monolayer with the bulk is 26.7 s À9 , but the water molecules in the trench have an exchange frequency of less than 1 s À9 . The hydration energy of the surface is À0.81 J m À2 , i.e. the average water-HAP interactions are stronger in the crenellated termination than on the flat (01% 10) termination and the water appears to be more ordered on this termination, where order here is defined as the geometrical arrangement of the water molecules close to the surface.
The (0001) surface is the thermodynamically most stable surface in both fluorapatite 38 55 The water molecules (red stars in Fig. 5a) also interact with surface hydroxy and phosphate oxygen atoms via their hydrogen atoms, with almost all the phosphate ions involved in the interactions at distances of less than 2 Å, in a good agreement with the 1.8 Å found by Pareek et al. 54 and the 1.5 Å by Corno et al. 55 The water molecules closest to the surface do not form a regular adsorption pattern on the surface (Fig. 6b), and they are barely exchanging with the bulk solution (B1.7 s À9 ). This surface has a smaller surface area than the (01% 10) and the hydration energy is À1.36 J m À2 , indicating that the water is more strongly attracted to this surface than to either termination of the (01% 10) surface, which agrees with the much lower frequency of water exchange. The gap (B1.5 Å in Fig. 6a) between the first and second layers of water molecules on the (0001) surface is due to the stronger interactions with surface ions and also the existence of an extensive network of hydrogen-bonded interactions between the water molecules within the layer nearest the surface, at the expense of interactions with the water molecules in the second layer. Although a more stable surface may be expected to be less reactive towards water, when the geometry and composition of the surface allows strong interactions between water and surface, in addition to extensive hydrogen-bonding between the molecules within the water layer, the formation of a regular water layer is often found to be an especially favourable process. A similar situation was observed in previous work on the a-quartz (0001) surface, 56 but in that case the gap was even more pronounced (B2.5 Å) due the highly ordered nature of the first monolayer of water.

Adsorption of the GlcA on the hydrated surfaces
The saccharide was placed onto the surfaces in five different orientations (labelled from 1 to 5 in Table 1), corresponding to  Table 1 shows that GlcA adsorbs favourably onto the (0001) surface for two configurations (2 and 3). The adsorption energies on the (01% 10) surfaces are positive for all configurations, suggesting unfavourable adsorption of the saccharide on to these surfaces.
Although the values of adsorption energy are less positive in the case of the flat terminated (01% 10) surface than for the crenellated termination and for the (0001) surface. In all cases the GAG and the solvent compete for the adsorption sites on the surfaces and generally the saccharides interact with the mineral surfaces via a monolayer of interfacial water, but the few exceptions to this general behaviour, corresponding to the configurations with the least positive or negative adsorption energies in Table 1 are discussed in the sections below.
3.7.1 Adsorption of GlcA on the flat (01% 10) hydrated surface. Table 1 shows that for all configurations the adsorption energy of GlcA onto the flat (01% 10) surface in the presence of water is unfavourable. All configurations considered show interactions by oxygen of the carboxylate to calcium of the surface (Fig. S7, ESI †), but all the initial bridges and hydroxy interactions are no longer present because GlcA does not lie as flat on the surface as it does in the absence of solvent, with only the carboxylate and occasionally the hydroxy groups being located at the level of the first layer of water. In the system with the least positive adsorption energy, the strongest interactions occur between the carboxylate oxygen atoms and axial calcium ions (Fig. 7) while the hydroxy groups form hydrogen-bonds.
The GAG has undergone a conformational change and the carboxylate group is now in the axial position, not in the most energetically favoured equatorial position. The GlcA forms hydrogenbonds with the surface in all five configurations studied (Fig. S8, ESI †), but these interactions by themselves are not enough to stabilize the adsorption of the saccharide.
A combination of a strong and long-lasting interaction through the carboxylate group, the formation of hydrogen-bonds and the relative position of the molecule with respect to the surface plane (discussed further in Section 3.9) are the main factors contributing to the relatively small positive energy of configuration (2).
3.7.2 Adsorption of GlcA on the crenellated (01% 10) hydrated surface. The crenellated (01% 10) surface has a higher hydration energy than the flat surface, indicating greater surface reactivity, but all GlcA adsorptions are again energetically unfavourable (Table 1). Similar to the flat termination, in all configurations interactions occur between the oxygen of the carboxylate and the calcium in the surface, with the carboxylate group piercing the first hydration layer and sharing the trench with water molecules. In two configurations, the adsorbate is lying flat on the surface at the level of the water molecules in the first layer of hydration. Only in three configurations are hydrogen-bonds between hydroxy groups and phosphates also observed (Fig. S9, ESI †). The energetically most favourable configuration (3) is obtained when the organic molecule is not very deep in the surface trench, as indicated by the first peak shifted along the r axis of the RDF Ca 2+ -O carboxylate plot (Fig. S10, ESI †), but with the carboxylate still able to form a bridge between two calcium atoms in the bottom of the trench (Fig. 8).
It is noted that configuration 3 does not form hydrogenbonds to the surface as found for the most strongly bound configuration on the flat (01% 10) surface. GlcA in configuration 3 is adsorbed on the crenellated surface, with the molecule pinned to the surface in a vertical position in which the hydroxy groups are more exposed to the solution, hence interacting more with the water molecules, which also stabilises this configuration.
When the organic molecule is lying along the trench in the surface, all functional groups are able to interact strongly with the surface atoms. The interactions between the oxygen atoms of the saccharide functional groups and the calcium ions help to lower some of the positive adsorption energies compared to the flat (01% 10) surface.
3.7.3 Adsorption of GlcA on the hydrated (0001) surface. GlcA adsorbs on the (0001) surface, without causing considerable distortions to the top surface atom arrangements, as expected for the very close-packed (0001) surface. The hydroxy groups of the surface dissolve into the water in four of the simulations, similar to what was observed in previous works. 31,39 In four configurations, the GlcA is incorporated within the first layer of adsorption,  where in one (4) the carboxylate group is the only group able to interact with the calcium ions in the surface. Table 1 shows two adsorption configurations with favourable adsorption energies. In the most stable configuration, the saccharide is adsorbed on the surface where it interacts through its carboxylate and two hydroxy groups with surface oxygen atoms (Fig. 9). The other stable adsorption presents similar interaction distances. Fig. 8 also shows that the GlcA has experienced a change in conformation, where the carboxylate has rotated into the axial position to interact with the surface. It has been reported in previous work that binding to material surfaces may induce these kinds of conformational changes in organic molecules 25 to enhance adsorption.
The intensity of the RDF peak (Fig. S11, ESI †) is less than for the (01% 10) surfaces, indicating a weaker interaction, because the oxalate group is not able to form bridges with the calcium ions. For all the configurations in which GlcA remains in the first layer of adsorption, it forms hydrogen-bonds with the surface.
As occurs on the (01% 10) surfaces, the combination of the various interactions lowers the positive adsorption energies, reaching negative values with a change in conformation of the organic molecule. However, when the GAG remains in its most stable conformation in solution the adsorption is even weaker than on the (01% 10) surfaces, for two reasons: The interactions between surface and carboxylate are not particularly strong, because it cannot interact with more than one calcium ion at a time, and the relative position of the organic molecule with respect to the surface in certain cases removes the solvent from the monolayer, as is discussed in more detail in Section 3.9.

Adsorption of GalNAc-4S on the hydrated surfaces
Similarly to GlcA, the GalNAc-4S molecule was placed onto the surfaces in five different orientations (labelled from 1 to 5 in Table 2), corresponding to the fifth configurations with lowest energy values obtained for the adsorption simulations in vacuum.
The values in Table 2 indicate the adsorption of the saccharide on to the three surfaces covered in this study is thermodynamically unfavourable. For some of the configurations analysed, the molecule experiments complete desorption from the material. The adsorption energies are generally higher than those obtained for GlcA. The GalNac-4S also prefers to interact with the material through the monolayer of water and as it occurs the case of GlcA the least positive values of adsorption energies are for the (01% 10) surface with the flat termination.
In the sections below we will discuss the geometries leading to the least unfavourable adsorptions on each surface.
3.8.1 Adsorption of GalNAc-4S on the hydrated flat (01% 10) surface. The GalNAc-4S is sulfated in position 4, Fig. 1c. This sulfate group has three prominent negatively-charged oxygen atoms that can interact with HAP calcium ions, and it has more volume than the carboxylate present in GlcA. This functional group also attracts more water than the carboxylate.
In configuration 1 (Table 2), the organic molecule desorbs from the surface. The other configurations exhibit interactions between the sulfate and the surface, where in all cases the sulfate, and in one case the hydroxy group, are the only groups of the GAG within the first adsorption layer. The configuration with the   least positive adsorption energy (Fig. 10) is the only configuration forming hydrogen-bonds with the surface (Fig. S12, ESI †), which as we described in Section 3.7.1 for GlcA is a determinant factor to increase the strength of binding. Configuration 2 in Fig. 10 is showing a very similar adsorption to the one obtained for GlcA (Fig. 7). As no other functional group of GalNAC-4S is interacting with the surface, the stabilization of the adsorption of GalNAC-4S compared with the GlcA is due to the nature of the sulfate, the adsorption position of the molecule and the more frequent occurrence of hydrogen-bonds, between the hydroxy group of the organic molecule and the phosphate ions, shown in the RDF plot in Fig. S12 (ESI †). The intensity of the first peak is higher than the one obtained for a similar configuration of adsorbed GlcA (Fig. S8, ESI †), even though GalNAc-4S only has one hydroxy group while there are two in GlcA. However, we suggest that hydrogen-bonding is more transient, breaking and re-forming frequently during the simulation. 3.8.2 Adsorption of GalNAc-4S on the crenellated (01% 10) hydrated surface. During the simulation, the sulfate group of the saccharide inside the trench gradually moves upwards. In three configurations (2, 3 and 4) the GAG as a whole is lying flat on the surface, at the level of the first layer of hydration. In the other two configurations (1 and 5), only the sulfate and in one case the carbonyl groups are level with the first water layer or sharing the trench with water molecules. The adsorption energy is positive for all the five configurations (Table 2), with the weakest bindings seen in any of the systems analyzed in this study, indicating that the adsorption of GalNAc-4S on the (01% 10) surface is thermodynamically highly unfavourable.
Configuration 1 (Fig. 11) exhibits the strongest interaction between surface and adsorbate, with the sulfate group as the only functional group inside the trench, and an additional interaction between the carbonyl oxygen and the calcium ion of the surface. With the main part of the molecule outside the trench, no formation of hydrogen-bonds with the surface is observed (Fig. S13, ESI †).
As the GalNAc-4S molecule has a single hydroxy group, the presence of peaks in Fig. S13 (ESI †), also supports the fact that in the other configurations the organic molecule resides further into the trench of the surface, to allow the hydroxy group to interact with surface phosphate oxygens. However, here the several interactions occurring at the same time do not produce the same effects observed for GlcA on this surface.
The presence of the sulfate and the acetylamine groups in the molecule makes it a bigger system compared to GlcA, and the diffusion of the GalNAc-4S in the trench causes more rotations of the phosphate groups, and hence less coordination between Ca 2+ -O phosphate .
The distance between calcium ions is decreased in comparison with the pure surface, causing more strain in the surface and further instabilities. A combination of these effects due to the presence of the molecule, and the loss of interactions with the solvent (discussed in Section 3.9), causes the relative instabilities of these systems.
3.8.3 Adsorption of GalNAc-4S on the hydrated (0001) surface. In two configurations (1 and 5) the GAG desorbs from the surface ( Table 2). The strongest interaction occurs when the oxygen atoms of the carbonyl and occasionally the sulfate groups of the saccharide interact with the surface calcium ions (Fig. 12) and in this case only the interacting groups are at the level of the first hydration layer. Configuration 3 is another example of the importance of the position of the adsorbate on the surface and the formation of hydrogen-bonds between the organic and the surface, but the interactions between the surface and the sulfate are transient during the entire simulation time.   of solvent. The results of the MD simulations are in a good agreement with previous DFT calculations on GlcA, 25 showing the same interactions and very similar interaction distances and distortions in the GAG structure. The same study analysed a non-sulfonated GalNAc, so our results are not directly comparable, but any trends are similar. The discrepancies between the adsorption energies of the biomolecules obtained with DFT in vacuum and those reported in this work are due to the high stability of the GAGs in water, as well as the strong affinity of the hydroxyapatite surfaces for water. Both factors contribute to making the adsorption of the GAGs onto the HAP surfaces less favourable.

The effect of the aqueous environment
A common issue affecting the stabilization of both saccharides at the three hydroxyapatite surfaces is the interaction of the organic molecule with the solvent molecules near the surface, which are identified in Section 3.6 and represented in dark blue in Fig. 6b-d. A similar situation was reported in previous studies of the behavior of water at quartz surfaces, when a desorbed Si(OH) 4 molecule remained near the quartz/water interface 56 in a process calculated to be endothermic. During adsorption on the flat (01% 10) and (0001) surfaces both saccharides generally displace B10% of the water molecules closest to the surface, which value rises (up to 20% for the (0001) surface), if the molecule lies flat on the surface. In the case of the crenellated (01% 10) surface, both GAGs affect significantly the order of the water adsorbed in the bottom of the trench. The displacement of the water molecules in the first hydration layer prevents the formation of hydrogen-bonds between the water and the phosphate, and there are also fewer interactions between Ca 2+ and oxygen of the water, hence the system becomes less stable than the hydrated surface in the absence of GAGs. The strongest binding of GAGs is obtained for the systems in which the first hydration layer is less affected by the adsorbate, confirming that the interactions between calcium ions and water are very important, as was also shown in a previous study. 24 To overcome the effect of the loss of water at the surface, the saccharide has to mimic the water interaction in the same area, in a similar way as described in Section 3.7.3 for configuration 3 of GlcA on the (0001) surface (Fig. 8), where less than the typical amount of water is displaced by the adsorbate and where interactions between the oxygen of the hydroxy group and the surface calcium ions replace the Ca 2+ -Ow interaction of the lost water (see Fig. S14 and S15 in the ESI †). The hydroxy groups of the saccharide also coordinate to the remaining water molecules closest to the surface. As the (0001) surface is also the surface with the largest value for its hydration energy, and is the one where the water monolayer is most pronounced (Fig. 5a), the competition between GAG and water for the adsorption sites on the surface will have the highest impact on GAG adsorption energies. The crenellated (01% 10) surface, where the water molecules are highly ordered (Fig. 6d), which ordering is broken after adsorption, will be the second surface most affected by the competition mentioned above.
Recent simulations of the adsorption of peptides considering pH effects 46 report considerable changes of adsorption sites and geometries as a result of change in the pH of the solvent, leading to a change in adsorption energies with more negative values at lower pH. In view of those findings and suggestions by Posner and Park 40,48 regarding the influence of phosphate protonation on the structure and binding strength of the water monolayer, we can speculate that at the typical pH of B6-7.4 in the kidneys, 57 the changes in the structure of the water monolayer could lessen the competition between the organic and the water in certain adsorption sites, thereby increasing the binding probabilities for the GAGs. Some proteoglycans for example, have shown an increase in affinity at low pH. 50

Conclusions
We have performed a series of molecular dynamics simulations to investigate the adsorption of two monosaccharides of a common glycosaminoglycan at three major surfaces of hydroxyapatite in aqueous environment. The results have shown that the solvent is competing with the monosaccharides for certain adsorption sites at the surfaces, but that the water molecules also solvate the organic molecules. The results are in a good agreement with previous simulations of HAP-citric acid interactions 23 and experimental solid state NMR findings identifying the presence of a layer of water between the GAG and the surface of HAP. 58 In other studies of interactions of hydroxyapatite with amino acids, such as HAP-glycine, 59 the water also competes with the biomolecule for adsorption sites. Our results also suggest that when the surface-bound water molecules form a regular arrangement, as in the case of the crenellated (01% 10) surface, the adsorption of GalNAc-4S is thermodynamically unlikely.
Previous experimental work established that in general, there is a direct relationship between the ability of a macromolecule to bind to HAP and its ability to block crystal growth 60 where the interaction occurs mainly through surface calcium ions, 61 demonstrating that the calcium binding capacity of GAGs is a significant factor in their ability to inhibit HAP growth. 50 In the present study we have shown that only those monosaccharide functional groups able to bind to Ca 2+ and causing less perturbation of the surface hydration layer are viable for adsorption. Sulfate and carboxylate groups are the desirable functional groups for adsorption to take place, with the sulfate group showing the greater affinity for calcium, also in agreement with experimental findings demonstrating the importance of the sulfate group for the adsorption capacities of chondroitin 4-sulfate. 62 The values of the adsorption energies obtained in this work support experimental reports, which considered that proteoglycans and GAGs with more sulfate groups are better adsorbates to HAP than chondroitin 4-sulfate. 3,5 The present work suggests that further stabilization of the surface-saccharide complex compared to the purely hydrated surface is due to hydroxy groups in the organic molecule replacing the missing water-water and waterÀsurface hydrogen-bonds, while other stronger interactions also take place, e.g. between the carboxylate and sulfate functional groups and the surface.
The GlcA displaces fewer surface water molecules and maintains multiple surfaceÀorganic interactions for longer periods of time on the (0001) than on the flat (01% 10) surface, leading to stronger binding on the (0001), but it is also possible that the kinetic barrier to desorption is too large for this to be observed on the timescale of the simulations. The interactions of the saccharide with the surfaces are within the energy penalty of a conformational change in the adsorbate, while the GalNAc-4S is not able to sustain long lasting interactions with (0001), but its adsorptions onto the flat (01% 10) surface show very small positive values.
However if the saccharides are forming chondroitin 4-sulfate, the conformational changes observed in the monomers, specifically those stabilizing the interactions with (0001), are energetically forbidden by the presence of the glycosidic linkage and the length of the chain. Our simulations suggest that amino acids and GAGs behave differently when adsorbing onto HAP surfaces in the presence of solvent. The amino acids or peptides manage to displace the water and bind favourably to the surfaces, 24,53 whereas the trajectories and values of the adsorption energies obtained here for the GAGs support previous experimental findings which show that disaccharides present in proteoglycans do not bind to HAP for long enough to affect crystal growth. 63 If water is present, chondroitin sulfate binds to the layer of water rather than to HAP. 58 Our results also support suggestions from experiments at pH 4 10.5, that GAGs act as templates for HAP formation, instead of moderators of crystal growth by binding to the surfaces. 58 The same experimental work showed that when HAP was formed in situ at pH B7À8 in the presence of the GAG, the polysaccharide was more closely attached to the surface, confirming the importance of the pH and therefore the inclusion of both solvent and surface protonation in future simulations to obtain more accurate adsorption energies.