Molecular dynamics simulations of bio-active phosphate-based glass surfaces

Classical molecular dynamics (MD) simulations were used to study the structural changes in the surfaces of bio- compatiblephosphateglasseswithcompositions(P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55 − x (x=30,35,40)toevaluatetheir effect on the solubility of the material. Direct comparison of the data for the three compositions highlighted the critical role that an enhancement in Na + concentration plays in the hydrolysis of the material, which is respon- sibleforthereleaseofnetworkcomponentsintosolution.Thecalculationsalsocon ﬁ rmthatthemostsolublema-terial is (P 2 O 5 ) 0.45 (CaO) 0.30 (Na 2 O) 0.25 , has the lowest calcium coordination number, thereby causing fewer cross links to phosphate chains. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Phosphate-based glasses (PBGs) in the CaO-Na 2 O-P 2 O 5 system have a number of biomedical applications [1,2]. Their ionic compositions mimic the mineral constituent of hard tissue, and they can therefore act as a template to repair soft tissue [3]. Based on their controlled ion release capabilities, they can also temporarily replace hard tissues whereas metal-doped PBGs are used in veterinary medicine to treat trace element deficiencies [1,4,5]. PBGs can also be employed to deliver drugs, from simple antibacterial ions such as Ag or Cu [1,6,7] or radioisotopes [8] to potentially more complex organic molecules for chemotherapy applications [9].
The success in medicinal utilization of this unique group of materials is closely related to their ability to dissolve completely in vivo in aqueous media, where the solubility rate can be tuned by altering the Ca 2+ /Na + ratio. It is known that the solubility of these ternary phosphate-based bioglasses decreases with an increase in the Ca 2+ /Na + ratio in the composition [10,11]. Much of our recent computational work has been dedicated to improving our understanding of glass dissolution processes and the structural and compositional features that affect them [12][13][14]. For example, previous theoretical studies of the bulk glass have shown that the mentioned change in solubility is caused by Ca 2+ binding more phosphate chains than Na + , thereby satisfying its bonding preference for non-bridging oxygen neighbours at the expense of sodium [12]. Also Ca 2+ binds to more PO 4 3 − tetrahedra and has a lower concentration of intra-tetrahedral phosphate bonds than Na + , where both features help to enhance the glass's durability as more Ca 2+ is included [12]. The same authors proposed that the clustering of modifier cations, which affect the solubility of these bioactive glasses by strengthening the glass network, should also be considered [14].
Thus far, theoretical studies on PBGs have been limited to the bulk material, but further detailed investigations must directly focus on the interface between the glass and the physiological environment, considering that it is at the surface that the bioactive processes occur. Moreover, among the wide range of factors affecting the solubility of the ternary phosphate-based bioglasses are the ion exchange between the surface and environment and variations in the glass surface structure. Classical molecular dynamics simulations (MD) with an accurate force field [13][14][15][16][17][18][19] have proven to be a highly suitable computational tool to obtain atomistic models of melt-derived glasses. They can provide a precise representation of the bonding and coordination environments found in the PBG surfaces, thus enabling the generation of larger statistical samples [15] than would be available from ab initio methods.
In this work we have therefore 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 = 30,35,40) and their influence on the material's bioactivity and solubility.

Structural models
The composition of the glasses and the notation used in this work are shown in Table 1 and the protocol for obtaining the models is described in details in Ainsworth et al. [1]. Starting from the volume-optimized bulk supercell, two slab geometries with 3D periodicity, but with the top and bottom faces exposed to vacuum (Fig. 1), were obtained by elongating the c vector of the simulation cell by~30 Å, producing four surface samples by exposing different sections of the bulk glass.

Evaluation of interatomic potential energies
Our simulations are based on the Born model of solids [21], where ions are assumed to interact via long-range Coulomb forces and additional short-range forces given by simple parametric functions which represent 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 via the model by Dick and Overhauser [22], 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, shown in Table 2, is the formal-charge, polarizable force field to describe inter-atomic forces, which has previously been shown to accurately reproduce the structures and properties of glasses in the system P 2 O 5 -CaO-Na 2 O, compared to experiment and ab initio methods. [1,10,12].

Molecular dynamics simulations
The MD runs were performed in the constant volume and temperature (NVT) canonical ensemble at 300 K using the DL_POLY code (version 2.20) [23], where all atoms were allowed to move during the simulations. An Evans thermostat [24] was used 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 [25,26]. The simulation time for each production run was 1 ns, and 50 ps of equilibration.

Species distribution in the surfaces
We have considered all the atoms within 10 Å from the interface to be surface species, in agreement with previous definitions [27]. Significant reordering of the species at the interface, occurs in all three glass compositions. Figs. 2 and 3 show example topographies observed during the simulation, including some orthophosphate ions (Fig. 2b), under-coordinated phosphate units (P 3C ) (Fig. 3a), and the formation of multiple-ring structures, some of them formed by three (3 M) or four (4 M) phosphate units (Figs. 2a, c, 3b). Small 3 M rings were not reported in previous theoretical studies of the bulk material [1,10], although they have been found before in experiment [28] and have been reported to be a very significant factor in the nucleation of apatite on silicate-based bioglass [29][30][31]. Many of the ring structures observed here are directly exposed at the interface and disconnected from the main phosphate network.
As occurs in silicate-based bioglasses [18], some of the highlighted structures as well as the modifier concentrations affect the glass dissolution process. The distribution of the species near the interface can be represented using the z-profiles in Fig. 4, for the three glass concentrations studied in this work, where we have defined, rather arbitrarily, the interface as the plane passing through the centre of the topmost oxygen atoms. The z-profiles represent the portion of atom x (where x = Ca, Na, O, P), that can be found at height z of the surface. There are clear compositional changes in the interface region compared to the bulk-like region (z = 0).
All surfaces in N20 and N15 experience an increase in sodium concentration at z N 10 Å, where in some cases the Na + are found as the topmost species in the surface, similar to what has been reported for the silicate-based glasses [18]. In the case of N15 at z N 10 Å, the fraction of sodium equals the calcium value in surfaces with less sodium content.  For N25 surfaces, the increase is only significant at the interface z N 13 Å. In the PBGs, however, there is no drop in the silicon counterpart, i.e. phosphorous content, but instead the calcium content is less in N25 and N20. This effect is not surprising, if we consider that in these two systems there is more Na + than Ca 2+ in the bulk glass. Table 3 confirms that Na + is occupying more surface sites than Ca 2+ in N25 and N20, whereas in N15 where the bulk composition of Ca 2+ is bigger than Na + , the latter only becomes more prominent close to the interface.
The creation of the surfaces also causes the exposure of many nonbridging oxygens (NBOs), several of them projecting out of the interface. Fig. 5 shows that many NBOs are products of the creation of the surface. The three glass surfaces, which contain similar amounts of oxygen, express higher fractions of NBOs in the top layer compared to the bulk-like region, leading to similar NBO concentrations for N25 and N20, while N15 exhibits more exposed NBOs (Table 3). However, the behaviour of the three glass compositions is slightly different in the layers near the interface. N25 exhibits a minor drop in NBOs compared to the bulk, before attaining the maximum, while N20 shows an increment of NBOs for all z N 10.

Surface network
Other structural attributes of the surface, that are useful to analyse and that contribute to explain the different behaviour of NBOs in N20, are the distribution of the Q n species (where Q n denotes a phosphorous atom sharing an oxygen with n others). The Q n distributions for the surfaces (using only the top 10 Å of the slab) of the three glasses are shown in Table 4. For all three compositions, at surface level the phosphorus atoms prefer to share two oxygen atoms with two other P atoms, which is similar behaviour to that observed in the bulk in both theoretical and experimental studies [1,10], but there is a significant increase in the cross-linking of just two P tetrahedra. Moreover, there is an evident increase in Q 1 and decrease in Q 2 with respect to previous calculations of the bulk (Table 4), which is partially due to the surface fracture, but is also strongly linked to the composition of the cation modifiers, which are known to alter the phosphate network [1,10,12,14]. We observe the same trends in the surfaces and bulk, i.e. an increase in Q 1 and decrease in Q 2 with the increase in Ca 2+ content in the glass (N25 to N20), as reported in theoretical studies of the bulk glass [1,10], but there is a small drop in Q 1 from N20 to N15, in agreement with experimental findings in the bulk glasses [1,32].
The small variations in the surface Q n distribution among the different Ca 2+ compositions are mirrored by the absence of trends as a function of Ca 2+ /Na + ratios in the PO 4 3− chain distribution presented in Fig.   Fig. 5. z-profiles. Fraction of NBO and BO in the slab model, calculated as in Fig. 4.  6. The average length of the chains is 5.7 (N15), 5.2 (N20) and 5.3 (N25). As is observed in the bulk glass, there is only a small difference in the chain lengths among the compositions. The values are shorter than reported for the bulk, because only the first 10 Å from to the interface are considered where there is an increase in the number of chains with lengths =2, 3. Another distinctive feature in the surfaces is the increase in Q 0 due to the presence of non-connected phosphates, mainly detected as undercoordinated phosphorus (Fig. 3a) and free orthophosphates (Fig. 2b). P 3C have been reported before in melt glasses [10] and appear here during relaxation. In the case of N25 there is also a 1% increase of Q3 phosphates caused by the formation of ring structures (Figs. 2a, 3b).
The data in Table 4 show that the surface with composition N20 is the one with the largest shift to lower values of n in Q n , indicating a lower network connectivity (NC) and thus exhibiting more NBOs per total amount of oxygen.
Lower values of network connectivity are closely related to the bioactivity and increase in the solubility of the glass, but, as stated in previous studies [12], 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
Previous MD simulations of silicate glasses have underlined a link between the aggregation of modifier cations in certain areas of the glass and its stability and bioactivity [14,20].
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 [8].
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. Within the definition of Eq. (1), a value of R N 1 denotes aggregation of the cations at the interface. The values of R for the three glass surfaces are reported in Table 5 and clearly show clustering of both modifiers in the region near the interface. As expected from the z-profiles and Table 3, the effect is more pronounced for Na + than for Ca 2+ . In the case of calcium there is no difference in the clustering trends for the N25 and N20 compositions, whereas for sodium the tendency to aggregate varies among all compositions. In Fig. 4, it can also be seen that the amount of sodium increases more rapidly for z N 10 Å in N20, although in the N20 slabs the Na + remains in a region away from the interface. The difference observed among N25, N20 and N15 is more related to the initial Na + composition in the materials.

Cation coordination
We have analysed the coordination numbers (CN) of both Ca 2+ and Na + modifiers in the surface region, as collected in Table 6. Owing to the creation of the surface from the bulk, it is expected that the coordination numbers of the ions closer to the interface are lower than those found in the bulk (Fig. 7). We also found that the decrease in coordination is at the expense of Me-NBO coordination (where Me = Ca 2+ , Na + ) with the two extreme calcium compositions (N25, N15) the most significantly affected for calcium coordination, while the sodium experiences the same drop in NBO coordination compared to the bulk values for the three compositions, showing the lowest coordination number in N25.
Similarly to the behaviour in the bulk, [CN Ca-NBO/Ca-BO ] N [CN Na-NBO/ Na-BO ] for surfaces of the three compositions. Calcium has a higher field strength than sodium and will bond to more NBOs, cross-linking more phosphates and strengthening the glass structure [12].    (Table 6) indicates that the chains closer to the interface in these surfaces are structurally weaker in comparison with those of different compositions in this study. In the same way the higher values obtained for calcium coordination in N15 indicate that these surfaces are more resistant to processes like dissolution.
In order to confirm the suggestions obtained through the CN we have also computed the number of chains bonded to each modifier ( Table 7).
It is reported for the bulk glasses that, on average, a Ca 2+ is bonded to 3.9-4.0 chains, while a Na + is only bonded to 3.2-3.3 chains [12]. We note that these values have fallen in the surfaces to 2.8-3.0 chains bonded to Ca 2+ and 2.5-2.7 chains bonded to Na + , suggesting that a lower calcium content in the PBG will lead to weaker structures and higher solubility rates.
We have calculated additional significant structural features that also reflect the differences in dissolution rates and have been reported before in theoretical studies of the bulk, such as the intra-tetrahedral bonding, a phenomenon describing how many oxygens of the same phosphate tetrahedron are coordinated to a given modifier cation [12] (Fig. 8). We found that an average Ca 2+ had 0.8 intra-tetrahedral bonds at surface level (0.9-1.1 in the bulk [12]), while an average Na + had 1.5-1.6 ( Table 8) compared to 1.2-1.3 in the bulk [12].
As occurs in the bulk, the lesser tendency of calcium for intra-tetrahedral bonding is likely to contribute strongly to the higher number of phosphates found around a single cation, and hence the strengthening effect on the network. The slight increase in the intra-tetrahedral bonding for sodium at surface level is related to the higher Na + concentrations observed compared to those in the bulk. As there are more Na + near the interface, these ions have to satisfy their need of coordination through using more than one oxygen from the same phosphate group.

Conclusions
The glass bioactivity is closely linked to the ability of its components to disconnect from the surface [18,33], and we thus need to investigate all the different surface features that could be involved in this process. The surfaces of the three compositions are fragmented in a similar fashion, but the slightly higher values of Q 1 and Q 0 in N20 suggest that these surfaces would be more prone to attack by solvent. However, the network modifiers play a key role in the dissolution. All the surfaces of the compositions studied here experience an enrichment in sodium, but it is the higher Na + fraction in the bulk that determines the sodium content in the surfaces. Furthermore, in N20 sodium has a preference to aggregate more in a region just below the interface, and there are therefore not as many Na + ions exposed at the interface, i.e. closer to the solvent, (see concentrations in Table 3), as in other compositions. Previous studies have shown that Na + -water interactions are potential initiators of surface degradation [34], whereas it is also known that a quick release of Na + will favour Na + /H + (solvent) exchange, leading to high alkalinity where the increasing amount of OHinduces the breaking of P\ \O\ \P bonds [18,34,35]. The surfaces with higher concentrations of sodium at the interface will thus dissolve more quickly. The calcium, however, does not show a distinct tendency towards clustering in the studied range of compositions. The Ca 2+ cations are responsible for strengthening the network, cross-linking various chains through coordination to NBOs and it is in the N25 (lowest Ca 2+ composition) surfaces, where calcium coordination numbers are lower, that the cations are bonded to fewer phosphate chains. These factors are the strongest contributors to a weaker network in the N25 surfaces and therefore facilitators of the dissolution process. Similar to 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, apatite formation.
Like silicates, immediately after contact with water the top phosphate groups will be rapidly hydroxylated by protonation of most external NBOs. Table 3 suggests that the network formers in the surfaces with N25 and N15 compositions should be hydroxylated to the highest degree. This is the initial step in hydrolysis, but the difference in hydroxylation between these two concentrations may not be enough to override other factors like network strength.
Small 3 M rings are understood to be relevant to the long-term solubility of the glasses, since metal poly-phosphates are highly insoluble [36]. In this study, they are only formed at the interface z N 14, with very few cations in the first coordination sphere, and at almost all times disconnected from the main phosphate network. Experimentally it is known that, once in solution, the phosphate rings are very stable [2,37], but their very low concentrations observed in the MD, the lack of solvent in the model employed and no other information from higher levels of theory, make it difficult to predict their impact.
In summary, we have carried out a series of molecular dynamics simulations of different phosphate-based glass surfaces, with three different CaO and Na 2 O compositions, in an effort to link their structure to their aptitude towards dissolution and bioactivity. Our results are in a very good agreement with previous predictions based on bulk  Fig. 8. Na + coordinating to two oxygens of the same phosphate (intra-tetrahedral coordination highlighted in balls) Color code: Na violet, Ca green, O red, P pink. materials, highlighting N25 as the most soluble material, and they show the critical role of the enrichment of sodium at surface levels. However, following similar work on the dissolution of silicate surfaces [37], further work including explicit solvent and pH effects is needed to assess clearly the diffusion of the ions and more complex phenomena like the formation of the small rings and their role in the solubility of the PBG materials.