Journal Non-Crystalline Solids A molecular dynamics study of the e ﬀ ect of water di ﬀ usion into bio-active phosphate-based glass surfaces on their dissolution behaviour

Classical molecular dynamics (MD) simulations were used to study the e ﬀ ects of water on the structural re- laxation of the surfaces of bio-active phosphate-based glasses with compositions (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55-x ( x = 0.30, 0.35 and 0.40). Direct comparison of the data for the three compositions showed that surfaces with x = 0.30 experienced the highest calcium di ﬀ usion, as well as highest sodium concentration at the glass/water interface, con ﬁ rming these systems as the most soluble of the three compositions studied. Our results also show the importance of surface hydroxylation in the simulation of these types of bio-glasses, which causes di ﬀ erential relaxation of the surfaces, leading to changes in network polymerization that modulate the di ﬀ usion of water and modi ﬁ ers into and out of the glasses, with direct impact on dissolution. Data curation, Formal Data Investigation,


Introduction
The phosphate-based glasses (PBGs) in the CaO-Na 2 O-P 2 O 5 system have ionic compositions which are similar to the mineral constituent of hard tissues. As such, they are viable for bio-medical purposes [1,2], e.g. in tissue repair or replacement [3,4]. PBGs are also capable of releasing their component ions in a controlled manner, making them suitable for treating trace metal deficiencies in animals [5][6][7], and controlled drug delivery, e.g. the slow release of antibacterial ions such as silver or copper [5,8,9], radio-isotopes [10], or potentially more complex organic molecules for chemotherapy [11].
The success of medicinal utilization is closely related to the ability of the materials to dissolve completely in vivo in aqueous media, where the solubility rate can be regulated by changing the Ca 2+ /Na + ratio. For example, it is known that the solubility of these ternary phosphatebased bio-glasses decreases with an increase in the Ca 2+ /Na + compositional ratio [12]. Recent computational work has been dedicated to improving our understanding of glass dissolution and the structural and compositional features that affect the proces [13][14][15]. Thus far, however, theoretical studies on PBGs have been limited to the bulk material. Our recent study [16] of the surfaces of three PBG compositions in vacuum highlighted the migration of Na + and Ca 2+ to the region closer to the interface, but further detailed investigations must consider the influence of the solvent and the pH on the relaxation of the surfaces, on any ion migration, and on the solubility of these ternary phosphatebased bio-glasses. For example, experimental studies of silicate glasses have reported that water diffusion into the glass surface plays a key role in the stability of the material [17][18][19][20], which we therefore need to consider in PBGs as well.
Classical molecular dynamics simulations (MD) with accurate force fields [21] have been demonstrated to be a highly suitable computational tool to obtain atomistic models of melt-derived glasses. They can also provide a precise representation of the bonding and coordination environments found in the PBG surfaces, thus enabling the generation of larger statistical samples [22] than would be available from ab initio methods. In this work we have used MD to explore structural changes in the glass surfaces of compositions (P 2 O 5 ) 0. 45 (CaO) x (Na 2 O) 0.55-x (x = 0.30, 0.35 and 0.40) and the effects of the solvent on the bioactivity and solubility of the materials.

Structural models
The compositions of the glasses and the notation used in this work are shown in Fig. 1, and the protocol for obtaining the models uses a melting temperature of 2500 K and a cooling rate of 5.50 K/ps, is described in detail in Ainsworth et al. [5]. Starting from the volume-optimized bulk supercell, two slab geometries with 3D periodicity, but with the top and bottom faces exposed to vacuum, were obtained by elongating twice one vector of the simulation cell, in two different directions, by~30 Å, thereby producing four completely different surface samples by exposing different sections of the bulk glass, thus ensuring that the surfaces represented as much as possible a cross section of the most important surface features, including oxygen sharing by phosphates The vacuum region was then filled with water molecules at a density of 1.27 g/cm 3 , which is still higher than the experimental water density, but that is due to the shell water model employed in the simulations. The full configuration results in a total of 4731 species in the PBG surface-water system. Hydroxylated surfaces were created by protonating explicitly all the non-bridging oxygens (NBOs) found in the upper-most layer of the slabs. The use of a reactive water potential [23,24] could spontaneously produce a fully dissociated water sample at the interface levels, similar to that described in the sections below. However, due to the intrinsically limited timescales of the MD, the formation of a fully dissociated sample, reflecting experimental findings, can be hampered by the kinetic barriers of the water dissociation process. Therefore, we have chosen the explicit hydroxylation technique as a suitable starting approach to this complex behaviour. In order to preserve the charge neutrality in the system and simulate experimental conditions, we have associated the OH − groups, resulting from the water dissociation to produce the protons, with the modifier cations in the top layer of the surface, in accordance with one of the hydrolysis mechanisms proposed from ab initio-MD simulations for silicate-based bioglasses [25][26][27]. The selection of the OH − -coordinated modifier type (Na + or Ca 2+ ) and its location in the surface was performed randomly. The number of these OH − species varies from surface to surface, covering a range of possible pH conditions (but not exceeding the relevant values for the PBC application of 6<pH<8). These relatively mobile OH − can migrate along the simulation box and interact with the phosphate chains, in some cases breaking the P-O-P bridges as described in Section 3.2 and figures therein. Both types of systems (hydrated by molecularly adsorbed water and hydroxylated) are investigated in our simulations.

Interatomic potential models
Our simulations are based on the Born model of solids [28], where ions are assumed to interact via long-range Coulomb forces and additional short-range forces, which are described by simple parametric functions representing electron-mediated interactions, e.g. Pauli repulsions and Van der Waals dispersion attractions between neighbouring electron charge clouds. The electronic polarizability of the ions is also accounted for, using the model by Dick and Overhauser [29], in which each polarizable ion is represented by a core and a massless shell, connected by a spring. The spring constant and the distribution of the ion charge between the core and shell determine the polarizability of the ion. In our case only the oxygen anion was considered polarizable, while cations were modelled as rigid ions.
The potential model used here for the PBGs, derived by Ainsworth et al. [5], is a formal-charge, polarizable force field to describe the inter-atomic forces, which has previously been shown to accurately reproduce the structures and mechanical properties important for the accurate description of thermodynamical and transport properties of two phases of crystalline P 2 O 5 . It has further been employed to simulate glasses in the system P 2 O 5 -CaO-Na 2 O which showed good agreement with structural properties derived from experiment and ab initio methods [10,25]. The water potential is based on the polarizable water model originally derived by de Leeuw and Parker [30], with later modifications introduced by Kerisit and Parker [31], whereas the Baram-Parker potential is used to represent the OH groups [32]. The additional Buckingham parameters for the interaction between the phosphorus and the oxygen of the hydroxy group (P-Oh pair in table S1) were derived and successfully used before in the simulations of bioactive glasses. [33,34] The full set of parameters is listed in Table S1 (Supplementary information, SI).

Molecular dynamics simulations
The MD simulations were performed in the constant volume and temperature (NVT) canonical ensemble at 300 K using the DL_POLY code (version 2.20) [35], where all atoms were allowed to move during the simulations. An Evans thermostat [36] with velocity-verlet integration algorithm [37] was applied and the timestep between successive integrations of the Newtonian equation of motion was set to 0.1 fs. The motion of oxygen shells was treated by assigning them a small mass (0.2 a.u.) to simulate the fast adaptation of the electronic polarization to the change in ionic positions [38,39]. The simulation time for each run was 1 ns, with 50 ps in which the system reaches stable NVT condition.

Validation of the force field for the water-glass interaction
We have used the same protocol mentioned above (Section 2.1) to create a smaller sample of the N25 glass (Table 1) containing 270 atoms and a box size of length 15.86 Å. This small model is computationally feasible to simulate using ab initio MD (AIMD). This glass surface was also created as described in Section 2.1, and the system was then simulated using classical MD (Section 2.3) with the force field described in Section 2.2. We then extracted the final configuration and removed all solvent, bar one randomly selected water molecule interacting with the material. The energy of this system with a single water molecule was evaluated using both classical MD and the ab initio code CP2K [40], the latter with the parameters described in Ainsworth et al. [41]. The parameters employed in the energy evaluation were the same as used for the MD simulation of the larger and fully hydrated model, which is the subject of the rest of the paper. The enthalpy of hydration was then calculated using Eq. (6) The classical force field selected in this work is not able to simulate the full dissociation of the water molecules on the glass surface, but the AIMD simulations did not observe spontaneous dissociation either. Despite the approximations in both techniques, the difference in hydration enthalpies between the two methods is within 9% (−86.8 kJ/ mol for classical and −77.0 kJ/mol for ab initio) and the chosen force fields hence provides a good initial description of the PBG-water interface.
Moreover, we have extended our calculations with other PBG compositions in [41], using both the force field and AIMD and calculated the diffusion coefficients for the different ions in the glasses. The values are given in the SI. The agreement between the AIMD values and the classical force field MD is acceptably close and the relative trends amongst the atomic species is in good agreement. These results show that the force field not only matches high level AIMD results in the description of atomic diffusion in the melt.

Degree of polymerization
The degree of surface polymerization is estimated for each glass concentration by calculating the number of non-bridging oxygens per network-forming cation via the following equation [20], where Nc is the number of cations in the glass (Ca 2+ , Na + and P 5+ ), M is the number of network modifier cations (Ca 2+ and Na + ) and n is the charge of a cation in the top two layers of the slabs.

Modifiers clustering on the surfaces
A statistical measure of the inclination of cations of species x to form clusters is the ratio R x-x of the total number of other x atoms observed in its coordination sphere (N obs ) and the corresponding number N hom which would result from a homogeneous distribution, where all the ions x uniformly occupy the available volume [10].
where CN x-x is the x-x coordination number calculated for the distance r c , the minimum value obtained for their radial distribution function (r c =5 Å). N x is the number of x ions at the surface and V s is the volume of the surface section of the slab.

Coordination of modifiers to water
The coordination of the modifier cations to the water was estimated using the radial distribution function (RDF) M-OW in the surface slab (where M = Ca 2+ or Na + , OW is the water molecule's oxygen atom). The value of the integration of the first peak n(r) with r = 3.25 Å, is then multiplied by the mean number of cations that have real 'contact' with the water (M slab /M s ).
where n(r) = 4πg(r)ρr 2 dr, M s is the number of cations at the surface (z > 10 Å in Figs. 2,3) and M slab is the total number of cations in the simulation slab.

Diffusion coefficients of the modifiers
DL_POLY provides a tool to calculate the mean square displacement (MSD) of the different species. The MSD can be defined as: where R i is the position vector of atom i, and N is the total number of atoms i.
We have modified the code to consider only the z-coordinate of the position vector to obtain the MSD in the direction perpendicular to the surface and the corresponding diffusion coefficient as the slope of: where C is a constant and only the values of the last 15 ps of the 1 ns total simulation time are considered to obtain the most linear section of the MSD.

Enthalpy of hydration
For each surface slab the hydration energy is given by: where E surf,solv is the average configurational energy of the hydrated surface, E surf is the average configurational energy of the dry surface, E H2O is the energy of one liquid water molecule and n is the total number of solvent water molecules. The values are then divided by N (z), the number of water molecules with the centre of mass less than or equal to the height z of the surface. In this case, this is equivalent to the distance of the first hydration shell obtained from the radial distribution functions Ca-OW, Na-OW and O-HW (see figure S1) i.e., in direct contact with either surface of the slab.

Solvent accessible surface area
The solvent accessible surface area was computed over the full trajectory of the simulation and the complete slab, using Mdtraj [42] module in the Anaconda Python distribution. The Shrake and Rupley [43] algorithm was employed with the parameters probe-radius=1.4 Å, number of spherical points to be considered was 960.

Results and discussion
We have considered all the atoms within 10 Å from the interface to be surface species, in agreement with previous definitions [44]. The 10 Å distance definition was established from the starting geometry, considering as the interface the plane formed by the topmost oxygen atoms; this reference was kept constant in time. In nature, and under experimental conditions, immediately after contact with water the top phosphate groups of the glass will experience a process of hydrolysis and be rapidly hydroxylated by protonation of the external non-bridging oxygens (NBOs). Here, we will consider two scenarios, i.e. both pristine glass surfaces and hydroxylated surfaces interacting with water, analysed below. We started the simulations considering just molecular water interacting with the pristine glass surfaces.

Species distribution in the non-hydroxylated surfaces
A good representation of the distribution of the species in the surfaces can be obtained from the z-density profiles (Fig. 2), from the glass bulk (z = 0) to the glass surfaces at the water interface (z > 0). Significant relaxation of the species at the interface occurs in all three glass compositions, when water is present, although the behaviour is different for each glass composition. In N25 surfaces, the amount of both modifiers decreases at z > 12 Å. In N20 surfaces, the amount of Na + remains constant while the amount of Ca 2+ decreases at z > 10 Å. In N15 surfaces, the amount of both modifiers remains constant at z > 10 Å. This behaviour contrasts with what was previously observed for dry surfaces [16], at the vacuum interface, where all surfaces experienced a significant increase in sodium at z > 10 Å.
In both scenarios (dry and hydrated as cut, from here on just "as cut"), Na + is the modifier cation that prevails in the top layers of the glass at the water interface (Table 1) in all three systems studied. However, the levels of Na + at the interface decrease with the increase its content in the glass. This behaviour can be explained by the degree of polymerization of the glass in those layers of the slab (NBO/Nc), which under these conditions is similar for N25 and N15, but the compactness of the mesh in the presence of water makes it difficult for Na + to migrate in N25.
Additionally, the presence of larger network cavities at the interface facilitate the full solvation of some Na + and Ca 2+ that, in the case of the N25 glass composition, leave the interface region and move into the solution, causing a drop in the concentrations at z > 12 Å. Furthermore, the water molecules that infiltrate the surfaces coordinate internal cations that do not migrate further towards the interface, but Na + diffuses faster than Ca 2+ (see 3.4.3), thus contributing to a major Na + population at the interface.

Distribution of species after hydroxylation of the surface
After the surface hydroxylation, the density of network formers (P 2 O 5 ) at the surface in N15 systems is higher and the cavities are not as large compared to those described before in 3.1. One consequence is a lower penetration of the water molecules into these surfaces. Fig. 3 also reveals more water molecules in the deepest regions of one of the N25 surfaces than in N15 surfaces, in which only 20 water molecules are found at z ≤ 10 Å.
We also note a major presence of Na + at the interface when z > 15 Å, with few Na + escaping into the solution. Hydroxy ions do not enter the N15 surfaces as deeply as they do in other surfaces, as they can coordinate more easily to Ca 2+ , due to the higher Ca 2+ concentrations at the interfaces of these systems.
As was observed in the "as cut" surfaces, Ca 2+ migration towards the water interface is less than for Na + , reflected by the increase in Na + at the interface of the surfaces. We can confirm from Table 1 that the Na + are prevalent at the interfaces, and, perhaps not surprisingly, the levels of Na + are higher in the systems with high Na + concentration in the bulk. Also, the (NBO/Nc) values show a less polymerized surface in N25, facilitating the migration of Na + towards the surface and the water into the glass. It is known that a quick release of Na + will favour Na + /H + (solvent) exchange, leading to high alkalinity, whereas the increasing amount of OH − induces the breaking of PeOeP bonds [45,46].

Surface network
In our previous study of the surfaces in the absence of water [16], the surface region was denser than the bulk glass and multiple ring structures were observed, some of them formed by three-(3 M) or fourmembered (4 M) phosphate units. As a result of the availability of water molecules to satisfy dangling bonds, the number of small ring structures has decreased as result of surface solvation to a negligible quantity, which is more pronounced if hydrolysis is considered in the model. This trend is independent of the composition of the glass.
The disruption of the phosphate network at the interface, due to the cleavage and the presence of the solvent, produces some cavities and trenches that will vary in size depending on the glass composition and the consideration of the hydrolysis of the phosphate ions. A deeper analysis can be done following the Q n distributions, where Q n denotes a phosphorus atom sharing some oxygen with n others, were calculated using only the top 10 Å of the slabs of the three glasses. For all three compositions, at surface level the phosphorus atoms prefer to share two oxygen atoms with two other P atoms, which behaviour is similar to that observed in the bulk in both theoretical and experimental studies [5,47]. There is an evident increase in Q 1 and decrease in Q 2 species when surface hydroxylation is considered, compared to the molecularly hydrated surfaces (Table 2). Furthermore, in the hydroxylated surfaces there is a small drop in Q 1 with an increase in Ca 2+ content in the glass (N25 to N20), in agreement with experimental findings in the bulk glasses [1] and with previous simulations of the bulk [5], thus showing the strong link between the concentration and ratio of the cation modifiers, and the alteration of the phosphate network [5,13,15,47]. We observe an increase in non-attached phosphate groups in the hydroxylated surfaces, confirming the importance of surface hydrolysis in the dissolution process of the glasses. As we were able to observe during the simulation of N25, where a calcium vacancy caused a mobile hydroxy ion bound to a nearby calcium to approach the phosphorus and finally bind to it, which leads to the breaking of the PeOeP and the detachment of the group (see Figure S3) The data in Table 2 show that the hydroxylated surfaces with composition N25 present the largest shift to lower values of n in Q n , indicating a lower network connectivity (NC).
Lower values of NC are closely related to the bioactivity and an increase in the solubility of the glass, but, as stated in previous studies [13], the value of NC is only associated with the number of phosphates and not with the Ca 2+ /Na + ratio in the PBG composition. It is therefore relevant to analyse the interactions of the cation modifiers in the surfaces.

Cation clustering at the surfaces
Within the definition of Eq. (2), a value of R > 1 denotes aggregation of the cations at the interface. Sodium cations tend to aggregate more strongly than calcium in all systems. Ca 2+ aggregation only occurs in the surfaces and increases upon hydroxylation, while Na + aggregation in the hydrated surfaces only increases noticeably in the N20 system and in the hydroxylated N20 and, to a lesser extent, N25 systems (see Table 3).
The clustering peaks in the N20 compositions. It could be argued that such behaviour is linked to the amount of Q 2 and Q 1 phosphates at surface level (Table 2), which in N20 is very similar to the N25 composition, indicating that more NBOs are available in this composition to interact with the modifiers.
We also note that Ca 2+ tend to aggregate least in the N25 surfaces, which is consistent with their lower concentration. Moreover, considering that there is more Na + in N25 systems and it is indeed more mobile (see also 3.4.3), more Na + will be available to favour the dissolution process in systems with N25 composition, as shown in Table 1 for the hydroxylated surfaces, thus confirming the relevance of considering hydroxylation in these simulations.

Cation coordination
The averaged coordination numbers (CN) of both Ca 2+ and Na + modifiers with the network oxygens, in the surface region during the simulation, are summarized in Table 4. The addition of water increases the availability of NBOs in the interface region, boosting the CN for Ca 2+ at the surfaces, but, similar to the behaviour in the bulk [5], [CN Ca-NBO/Ca-BO ] > [CN Na-NBO/Na-BO ] for the surfaces of the three compositions. Ca 2+ has a higher field strength than Na 2+ and will bond to more NBOs, thereby cross-linking more phosphates and thus strengthening the glass structure [13]. If surface hydroxylation is neglected, the results suggest a weaker network in N15, than in N20, considering its lower CN for Ca 2+ , while hydroxylated N25 surfaces have the lowest Ca 2+ coordination numbers and the trend for the CN for Ca 2+ is as expected N25<N20<N15 The latter result is in better agreement with experimental findings showing that higher Ca 2+ concentrations leads to lower solubilities, i.e., N25 glass concentrations are more soluble than N15 systems [13,48], again indicating the importance of hydroxylation in PBG surface models.
We observe that the coordination numbers of the calcium to water (1.76-1.43) drop after hydroxylation in both surface slabs with N15 composition, which is in part a consequence of the previously described decrease in water diffusion into the bulk of the glass. Table 5 summarises the values obtained via Eq. (3) for the hydroxylated surfaces of N25 and N15 compositions. The numbers estimate the amount of water molecules interacting with each cation in the top 10 Å of the surface slab. As expected, we observe that the modifier cations in the top layer   of the surfaces are coordinated to a greater number of water molecules than the more internal cations. The Ca 2+ are more hydrated than Na + , which effect is more pronounced in the N25 surfaces (Table 5). Ca 2+ also has a larger coordination shell than sodium and Ca 2+ are more strongly coordinated to the NBOs of the network. Na + are more mobile (see Section 3.4.3) and will dissolve faster and more easily than Ca 2+ if more water molecules approach them, e.g. in the case of the N25 surfaces. Table 5 also shows the greater propensity of Ca 2+ for water coordination in N25, suggesting that the ions are more liable to dissolve than in the N15 systems. In N15 we observe the lowest coordination to solvent, indicating that those Ca 2+ are both less accessible to water and better satisfied in their coordination to the phosphate chain, making N15 surfaces more resistant to dissolution.
The numbers in Table 4 do not provide information regarding the cohesion of the network, although such information is available if we compute how many of the oxygen coordinated to Ca 2+ belong to different phosphate chains. Compared to the bulk glasses, the number of phosphate chains bonded to each modifier (Table 6) have dropped in the hydroxylated surfaces to 2.8-3.0 chains bonded to Ca 2+ (cf. 2.5-2.9 in the molecularly hydrated surfaces) and 2.4-2.5 chains bonded to Na + (cf. 2.5-2.6 in the molecularly hydrated surfaces), whereas in the bulk 3.9 − 4.0 chains are bonded to Ca 2+ and 3.2 − 3.3 chains are bonded to Na + [13]. However, the difference between the N25 and N15 systems is not significant enough to justify the suggestion of a weaker N15 network, where the Ca 2+ in those surfaces, is on average bound to 2.8 phosphate chains.
These results support a behaviour where a lower Ca 2+ content in the PBG will lead to weaker structures and higher solubility rates, i.e., systems with N25 concentration should be more soluble.

Cation diffusion
Next, we have computed the diffusion coefficients of the modifier cations in the direction perpendicular to the surface (z in Figs. 2, 3) for the surface slabs (Table 7) with N25 and N15 concentrations, (see the same numbers for a full bulk simulation in Table S5). The values of the diffusion coefficients are more reliable at higher temperatures, but the values in the tables (T = 300 K) will provide a good qualitative analysis between the two concentrations. The presence of solvent increases the self-diffusion of both Na + and Ca 2+ in the direction perpendicular to the newly created interfaces compared to the values obtained in the bulk. If surface hydroxylation is neglected, we note that Na + diffusion is faster than Ca 2+ , in all systems, but also that both modifier cations diffuse faster in N25, where Ca 2+ can diffuse up to two times faster than in N15. The largest values for the diffusion coefficient of Na + (6.7 × 10 −12 m 2 /s) and Ca 2+ (4.6 × 10 −12 m 2 /s) in one of the surface slabs of N25 is related to larger trenches occurring in the surfaces, thereby enhancing the movement of the modifiers through the solvent column. These results confirm the indications already discussed above, i.e., that the strength of the network is more related to the ability of Ca 2+ to coordinate the chains' oxygens as it is indeed the more static modifier, particularly in N15 systems. The trend in the diffusion coefficients is also in good agreement with the shown in Table S2 and previous theoretical [49,50] and experimental findings [51][52][53] in silicate-based bio-glasses; a higher Na + content leads to higher Ca 2+ diffusivity. In the slabs with hydroxylated surfaces, Na + also diffuse faster than Ca 2+ , in this case by an order of magnitude. The diffusion coefficients for Ca 2+ in the direction perpendicular to the surface in N25 are doubled compared to the values in N15, while Na + values are similar for all systems. The diffusion of the modifier cations is generally reduced after hydroxylation of the surfaces. In the case of Ca 2+ in N25 the decline is by two orders of magnitude; the hydroxy bonds reduce the network size, as seeing in Table 6, and mentioned the previous section, a Ca 2+ is shared by more phosphate chains, thereby lowering the attractive electrostatic interactions between Ca 2+ and the P 2 O 5 mesh, in agreement with experimental observations in silicate glasses [54].
In addition, mobile OH − that reach deeper layers of the surface coordinate the cations, thereby stopping their diffusion towards the interface, which behaviour mostly affects the N25 surfaces. Interestingly, in N15 surfaces the deepest OH − are part of the network species and their oxygen is not coordinated to a cation. Besides, the OH − approaching the surface also contribute to reducing the dynamics of the cations in that region, mainly calcium. This effect is accentuated most in the N15 surfaces with highest Ca 2+ concentrations at the interface. Although diffusion in the direction perpendicular to the surface is reduced by hydroxylation of the surfaces for all concentrations studied here, Ca 2+ remains more static in N15 surfaces, thus providing additional strength to the phosphate network and making dissolution more difficult in these systems.

Solvent accessible surface area
The solvent accessible surface area (SASA) was estimated probing the area around the atoms in every slab with a spherical molecule of radius 1.4 Å (approx. the water molecule) Fig. 4 shows the variation of the SASA in time for the different slabs and in both conditions of the simulation.
We can see that when the surfaces are hydrated as cut all slabs but two have similar SASA, particularly at the end of our simulation times. The two systems falling of the trend are one in the N25 composition and the second one in N20. The latter one going slightly higher than the group and the former staying all the time with lower values of SASA. This suggests little difference in solvation amongst the group of slabs (eight surfaces) disregarding of the glass composition. We need to highlight that the lowest values in one slab of N25 is due to its two surfaces that are less rough and more hollowed than the rest, at the interface, as we have mentioned above in the text and we can see Figure  (S4), showing snapshots of one surface taken at the end of the simulation, although the characteristic is sustained during the simulation as the trend in Fig. 4 confirms. In this surface slab, there are less network formers sticking out of the surface and the interfacial chains are not as long as they are in other systems. One composition N20 with one slab with a surface with high roughness, i.e., with phosphate chains sticking out at the interface, and because of the Ca 2+ concentrations and the  clustering of these modifiers observed at the interface, the SASA values are the highest (see Figure S6, for a snapshot). It is also good to note that systems with N20 composition undergo the greatest number of transformations along the simulation time, due to the already mentioned agglomeration of the cations, but mainly the high number of network fragments in direct contact with water, due to orientation towards the solution.
When the hydroxylation is considered all the systems experience and increase in the SASA, consistent with the migration of the cations to towards the interface, they get even more indistinguishable amongst the three compositions, the trenches and cavities observed before (without hydroxylation) are filled with migrated cations and "free" OH − , but the hydroxylation is not causing a relaxation of the N25 surfaces that leads to the full closure of the big trenches observed "as cut". In these two surfaces the water can still infiltrate deeper and increase the dissolution rate of the system. Again, the slab with N20 composition displays the highest SASA, this system the water molecules do not infiltrate as deep as they do in more hollow interfaces.

Conclusions
We have carried out a comprehensive molecular dynamics study of a few phosphate-based glass surfaces, with three different CaO and Na 2 O compositions, to link their atomic-level structure to their solubility and bioactivity. We have initiated the study of the role of the water and its diffusion into the glass surface in the dissolution process of pure phosphate bio-glasses, considering both molecular and dissociative hydration of a number of PBG surfaces. In future, this study could be extended to include additional surface compositions to improve the statistics and insight into the glass behaviour at the water interface.
Our results of the molecularly hydrated surfaces are in agreement with our previous findings for dry surfaces [16] and other studies of the bulk glasses [13,15], indicating that Ca 2+ are responsible for strengthening the network, cross-linking various chains through coordination to NBOs. Ca 2+ dissolves more easily from the N25 (lowest Ca 2+ composition) surfaces, where the coordination numbers are lowest. Moreover, these surfaces experience the highest sodium concentrations at the interface. These factors are the strongest contributors to a weaker network in the N25 surfaces and therefore facilitators of the dissolution process. Like silicate glasses, the network dissolution has a direct impact on bioactivity. The dissolved species will nucleate others, such as calcium phosphates, thereby encouraging, for example, the formation of the bone mineral apatite Furthermore, this work emphasises the necessity to consider both the solvent and hydroxylation of the surfaces to study the dissolution of bio-active phosphate-based glasses. Our results underline the impact of solvent on the migration of the modifiers to the surface, linked to the different relaxation of the species at the interfaces, while interacting with the water molecules. Further relaxation is experienced with hydroxylated phosphate groups, whereas the presence of mobile OH − ions particularly affects the diffusion of water molecules into the surface. Also affected are the diffusion coefficients of the modifiers in the direction perpendicular to the surface plane, as well as their coordination to the phosphate chain. Here, a reactive water potential combined with longer simulation times or enhanced-sampling techniques may overcome the water dissociative barriers and would be useful to study the full process of hydrolysis of the PBGs.
Our simulations were also able to show the detachment of phosphates in a N25 surface, in good agreement with previous studies, indicating that the increasing amount of OH − induces the breaking of PeOeP bonds [45,46]. Thus, a variation of the pH of the solution can strongly affect the solubility of these systems, although, in line with the previous literature, our findings indicate that only the molecular water is mobile enough to diffuse into the bulk glass and cause further hydrolysis [18].
Finally, also in agreement with experimental findings [48], our simulations confirm N25 as the most soluble glass composition of the three samples analysed here. Table 7 Diffusion coefficients (D z ) at 300 K of the modifier cations in the direction perpendicular to the surfaces of the systems with compositions N25 and N15.