Multi-step mechanisms of early phospholipid hydrolysis and mineralisation unveiled through combined quantum chemical calculations and experimental analysis

Abstract

Phospholipids play key roles in bone formation, with phosphatidylserine (PS) reportedly inducing more rapid mineralisation than phosphatidylcholine (PC); however, the underlying mechanisms remains unclear. This study investigated PS and PC mineralisation using experimental methods and computational chemistry. The stationary points in the potential energy surfaces of the reactions were preliminarily found using a neural network potential (PreFerred Potential in Matlantis) capable of predicting the interaction energies for arbitrary combinations of atoms, and then refined through density functional theory calculations (Gaussian16, at the B3LYP/6-31G(d,p) level of theory). When hydrolysis reactions were assumed to be the initial step in the mineralisation of phospholipids, the results were consistent with empirical analysis. PS was found to be more easily hydrolised than PC, primarily owing to the presence of a labile proton in the NH₃⁺ group of serine that facilitates proton transfer, enhancing hydrolysis of PS at lower energy thresholds. Specifically, when a single phospholipid was considered, three distinct hydrolysis routes were identified: between serine (or choline) and phosphate, between glycerol and phosphate, and between an aliphatic carbon chain and the glycerol backbone. In particular, the initial steps of hydrolysis involved the formation of a pentavalent phosphate intermediate. When calculations were performed with two adjacent phospholipid molecules, the loosely bound proton (H⁺) in the NH₃⁺ group could be readily transferred either to the P–O bond linking serine to the phosphate group; or to the P–O bond connecting the phosphate to glycerol in a neighboring PS6 molecule. These findings reveal the important roles of serine NH₃⁺ in facilitating hydrolysis of PS, and provide insights for designing novel molecules to accelerate bone regeneration.

---

1. Introduction

Phospholipids are phosphate diesters, with one ester commonly linked to head groups, such as serine or choline, and the other to the glycerol backbone, which is connected to two hydrophobic fatty acid chains. These molecules have been recognised as critical elements acting as nucleation sites for mineralisation during the initial steps of bone and tooth formation.¹ Extraction studies from mineralised bone, dentin and enamel at their initial mineralisation stages have identified several key species, including phosphatidylcholine (PC), phosphatidylethanolamine, sphingomyelin, phosphatidylinositol, phosphatidylserine (PS), phosphatidic acid, and cardiolipin.² Moreover, we previously demonstrated that cell-derived plasma membrane nanofragments (PMNFs) serve as the nucleation site for bone formation in vivo.³ These PMNFs, primarily comprising PS and PC, along with other proteins and glycoproteins in the phospholipid bilayers, were shown to initiate mineralisation within 1 to 2 days.^3,4

Between the two membrane-constituent phospholipids, PS was shown to have stronger mineralisation ability than PC.⁵ PS has also been shown to facilitate in vitro mineralisation by promoting the sequestration of calcium, which subsequently would bind free phosphate ions to form calcium phosphate minerals.⁶ PS has also been shown to promote and stabilise the formation of amorphous calcium phosphate, a precursor of hydroxyapatite.⁷ Owing to these capabilities, PS has also been applied as titanium implant coating to enhance osteoblast differentiation and bone regeneration in vivo.^8–11

However, the mechanisms underlying PS mineralisation remain unclear. In particular, the exact roles of the phosphate moiety of phospholipids in binding calcium and initiating mineral formation are not fully understood. Moreover, the distinct molecular structures of hydrophilic head groups, such as the presence of carboxylic or amino functional groups, may critically influence the mineralisation behaviour of different phospholipids, although this aspect has not been thoroughly explored.

Therefore, a more detailed mechanistic understanding of phospholipid-mediated mineralisation, including the molecular interactions between phospholipid moieties and calcium ions and the hydrolytic behaviour of phospholipids, is essential for the rational design of novel bone formation-inducing molecules.

Computational chemistry has emerged as a powerful and rapidly advancing field, driven by the development of high-speed computing systems, advanced algorithms capable of solving complex mathematical models, not only for inorganic systems^12,13 but also for biomolecules.^14,15 This approach allows for the elucidation of the physicochemical properties and reaction mechanisms of known materials and offers predictive insights for the design of novel materials. Thus, computational chemistry now plays a central role not only in explaining existing processes but also in guiding the synthesis of new functional materials. In this context, computational approaches are expected to significantly reduce the time, cost, and risks associated with material development, which traditionally relies on lengthy and resource-intensive experimental procedures.

In this study, we aimed to integrate empirical and computational analyses to compare the mineralisation of the two major phospholipid constituents of the cell membrane, and to elucidate the underlying mechanisms governing phospholipid-mediated mineralisation, with a particular focus on hydrolysis, which we hypothesise to be a key contributor to the overall reaction mechanism.

2. Experimental

2.1. Mineralisation assay

L-α-Phosphatidylserine (PS) and L-α-phosphatidylcholine (PC, L-α-lecithin) were purchased from Sigma-Aldrich (Saint Luis, MI, USA). Chloroform (FUJIFILM Wako Pure Chemical Corporation, Osaka, Japan) was used to dissolve all chemicals.

Five microliters of a 100 mM solution containing each phospholipid were dispensed onto a glass-bottom dish (Matsunami Glass Ind, Osaka, Japan) and left to dry. The phospholipids were then covered with 350 mL of 2 mM CaCl₂ solution and incubated in cell culture incubators at 37 °C, 5% CO₂, and 98% humidity. After 2 days, the mineralised phospholipids were collected in 1.5 mL tubes, centrifugally-washed three times with ultrapure water (Milli-Q) and dissolved in 0.1% HCl (FUJIFILM Wako Pure Chemical Corporation) for quantitative analysis using a simultaneous multi-element analysis atomic absorption spectrophotometer (AAS, Z-9000, Hitachi High-Technologies, Tokyo, Japan), available at the Central Research Laboratory, Okayama University Medical School. Standard solutions were prepared from a 1 ppm CaCl₂ stock solution (FUJIFILM Wako Pure Chemical Corporation) diluted in a 0.1% HCl solution.

For mineralisation experiments performed in different pHs, the same protocol described above was used, except that the 2 mM CaCl₂ solution was freshly adjusted to pH 6.5, 7.5 or 8.5, and immediately used to induce mineralisation of PS and PC for 1 h. The experiment was performed for 1 h to ensure stable pH conditions, as the pH tends to decrease as mineralization progresses.

2.2. Field emission scanning electron microscopy (FE-SEM) and energy-dispersive X-ray spectroscopy (EDS) analysis

For FE-SEM observation, minerals formed from phospholipid mineralisation were placed onto an aluminium holder, left to dry, and subjected to osmium coating (Neoc-STB, Meiwafosis, Tokyo, Japan) at an electrical discharge current of 10 mA and a vacuum of 10 Pa for 20 s. The samples were then observed using an SEM (S-4800, Hitachi High-Technologies), operated at 5 kV and 10 mA and equipped with an EDS detector (Apollo XV, EDAX, Mahwah, NJ, USA).

2.3. Computational methodology

The analysis of reaction dynamics requires the identification of stationary points on the relevant potential energy surfaces. The hydrolysis of phosphate compounds proceeds through mechanisms that are highly dependent on the surrounding environment.¹⁶

In particular, phosphate diesters have been reported to undergo hydrolysis extremely slowly¹⁷ and typically through the formation of a pentacoordinated intermediate, known as pentaoxyphosphorane, which has a trigonal bypiramidal geometry.^18,19 During hydrolysis, nucleophiles enter and depart from the phosphorane intermediate through an apical position, whereas electron-rich oxygen atoms lie in the equatorial position.¹⁸ The transition state leading to the formation and decomposition of the pentaoxyphosphorane intermediate is often difficult to locate and frequently requires exhaustive scanning of the relevant bond lengths. The identification of these transition states can be computationally demanding, often requiring several days, using conventional quantum chemical software.

To accelerate the search process for transition states, we employed Matlantis,²⁰ a recently developed software platform powered by deep learning, specifically designed to enhance material discovery and reaction pathway exploration. The software utilises the PreFerred Potential (PFP),^21,22 which is a neural network potential capable of predicting the interaction energies for arbitrary combinations of atoms. The use of PFP enables the optimisation of structures (including transition states), scanning of bonds and/or angles, and tracing the steepest path from the transition states to the products without the need of solving the Schrödinger equation. Following the initial localisation of the stationary points in each elementary reaction using Matlantis, all stationary structures were re-optimised and their vibration frequencies were calculated using quantum-chemical methods implemented in Gaussian16²³ at the B3LYP/6-31G(d,p) level of theory.

B3LYP is a hybrid density-functional that includes the Becke's three-parameter exchange functional²⁴ and the Lee–Yang–Parr correlation functional.²⁵ Owing to its favourable balance between computational efficiency and accuracy, B3LYP has been used extensively in the analysis of various types of reactions, including those involving biomolecules.²⁶ The 6-31G(d,p) basis set provides a combination of 6 Gaussian functions to approximate the atomic inner electron orbitals, while the valence electron orbitals are split into two, approximated by 3 and 1 Gaussian functions; additional polarization functions are included to account for angular flexibility and longer-range interactions.²⁷ All quantum-chemical calculations were performed assuming that the reactions occurred in aqueous solution using the self-consistent reaction field approach based on Tomasi's polarised continuum model (PCM).²⁸

For a more accurate evaluation of the potential energies, single-point energy calculations were carried out using the APFD functional of Austin, Frisch and Petersson with consideration of dispersion effects²⁹ in combination with the 6-311+G(2d,p) basis set. This larger basis set provides high-quality descriptions of the elements in the first three rows of the periodic table.³⁰ Accordingly, APFD/6-311+G(2d,p) single-point energies with solvent-model density³¹ correction were obtained for the B3LYP/6-31G(d,p) optimised structures. The Gibbs free energy of each stationary point was then calculated as the sum of the APFD/6-311+G(2d,p) energy and the free-energy correction obtained at the B3LYP/6-31G(d,p) level.

The actual phospholipid molecules used in the experiments contain an average of 18 carbon atoms in their aliphatic chains. However, because the aliphatic chains are not expected to possess high chemical reactivity, and to ensure computational feasibility, we considered phospholipids with truncated aliphatic chains containing only 6 carbon atoms. These simplified modes, referred to as PS6 and PC6, represent phosphatidylserine and phosphatidylcholine, respectively. The structures of PS6 and PC6 in water, optimized at the B3LYP/6-31G(d,p) level of theory, are shown in Fig. 1.

Fig. 1 (A) and (B) Molecular structures of (A) phosphatidylserine (PS) and (B) phosphatidylcholine (PC). HC1 and HC2 indicate the hydrocarbon chains. (C) and (D) The structures of phospholipids with 6 carbons (C, PS6; D, PC6) in water, optimized at the B3LYP/6-31G(d,p) level of theory.

2.4. Statistical analysis

Statistical analyses were performed using unpaired Student's t-test or one-way ANOVA, with Tukey post-hoc tests when appropriate. The StatView or GraphPad Prism 9 software was used for all analyses. The level of significance was set to p ≤ 0.05.

3. Results and discussion

3.1. Experimental results

3.1.1. Mineralisation of phospholipids. Comparative analysis using atomic absorption spectroscopy (AAS) of phospholipid mineralisation, that is, after incubation in a 2 mM calcium solution for 2 days, showed that the amount of calcium in the minerals deposited after PS mineralisation was significantly larger than that corresponding to PC (Fig. 2A). SEM analysis (Fig. 2B) provided evidence of the distinct morphologies of the mineral structures of PS (up) and PC (down). The PS-induced minerals exhibited large, clustered, and irregularly shaped particles, indicating more extensive mineral deposition likely by continuous Ca recruitment over the surface. In contrast, PC-induced mineral clusters appeared to be slightly smaller and spherical. Note that the Ca/P ratio estimated by the EDS results was of 0.3 and 0.35 respectively for PS and PC, suggesting a very early, calcium-poor and amorphous mineral phase.⁵ These results suggest that PS significantly exhibit higher mineralisation capacity than PC and highlight the crucial role of the phospholipid head group chemistry in modulating mineralisation.

Fig. 2 Comparative analysis of mineralisation capacity of PS and PC in 2 mM CaCl₂ solution. (A) The bar graph shows calcium content measured by AAS after 2 days of PS or PC mineralisation. PS exhibited significantly higher calcium deposition than PC, indicating enhanced mineralisation capacity. **p < 0.01, unpaired t-test (N = 3). (B) SEM images of mineral structures of PS (upper panel) and PC (lower panel). Inset shows the Ca/P ratio calculated based on EDS. (C) Quantitative analysis of calcium content after PS and PC incubation for 1 hour in 2 mM CaCl₂ solutions adjusted to different pH values. *p < 0.05, **p < 0.01, one-way ANOVA, Tukey test (N = 3).

Mineralisation of PS and PC was further evaluated under different pH conditions (pH 6.5, 7.5, and 8.5) in a 2 mM CaCl₂ solution. PS consistently exhibited significantly higher calcium deposition than PC at all tested pH values (Fig. 2C), confirming its superior mineralisation capacity. A modest yet significant pH dependence was observed for PS, with mineralisation being enhanced at higher pH, whereas PC showed substantially lower mineralisation with limited responsiveness to changes in pH. The enhanced PS mineralisation under more alkaline conditions could be associated with protonation or deprotonation effects involving PS functional groups, such as the phosphate or serine head group, which could facilitate Ca²⁺ binding to phosphate and stabilise early mineral-associated complexes. To clarify these possible protonation/deprotonation processes in the PS molecule, computational analyses were performed.

3.2. Computational results

PS6 and PC6 were employed to facilitate computational analyses. In order to verify the influence of the length of the aliphatic carbon chains on the calculated energetics, we compared the activation free energies for the first step of the hydrolysis (namely, the approach of two water molecules to the phosphate group) of PS18/PS6 and PC18/PC6. As explained later, this step constitute the rate-determining step of the hydrolysis mechanism. PS18 and PC18 contain 18 carbon atoms in the aliphatic chains, and represent more realistic models of native PS and PC. The calculated activation free energies (as described in the previous section) can be seen in the second row of Table 1. The activation free energies for PS6 and PS18 differ by less than 0.8 kcal mol⁻¹, while PC6 shows an activation energy that is 2.1 kcal mol⁻¹ higher than that of PC18. Despite this difference, the overall effect of chain length on the activation barrier is relatively small, supporting the conclusion that PS6 and PC6 provide sufficiently accurate models.

Table 1 Activation free energies of the first step of hydrolysis for various reactants and different levels of theory

Level of theory	Reactants	PS6	PC6	PS18	PC18
B3LYP/6-31g(d,p)	+2H2O	36.0	38.8	35.2	36.7
ωB97X-D/6-31g(d,p)	+2H2O	34.6	38.7
B3LYP/6-31g(d,p)	+Ca + 4H2O	37.3	36.9
B3LYP/6-31g(d,p)	+Ca + 6H2O	39.5	37.3

---

Next, we evaluated the suitability of the B3LYP/6-31G(d,p) level of theory. We compared the activation free energies for the first step of hydrolysis when ωB97X-D, a functional of the electron density that includes dispersion interactions and could yield more accurate results,³² was used in the optimization of transition states and reactants in PS6 + 2 H₂O and PS6 + 2 H₂O. A comparison between the B3LYP/6-31G(d,p) and ωB97X-D levels of theory (second and third rows of Table 1) shows that ωB97X-D reduces the activation free energy for PS6 by 1.4 kcal mol⁻¹, while for PC6, the corresponding values remain essentially unchanged, with a difference of less than 0.1 kcal mol⁻¹. Given the small magnitude of these differences, we therefore concluded that B3LYP/6-31G(d,p) does provide reliable results.

In addition, the explicit inclusion of additional water molecules is expected to improve the accuracy of the calculations. We initially performed the calculations with more water molecules surrounding the two water molecules that approach the phosphate group. However, the potential energy of the system was not sensitive to the positions of the extra molecules, and in some steps along the reaction path they were highly mobile, preventing reliable optimization of their positions. Consequently, only the two approaching water molecules were retained in the final models.

When calcium ions were included, additional water molecules could be added around the calcium ion and their positions could be successfully optimised. The last two rows of Table 1 show the results obtained when two or four additional water molecules were explicitly coordinated to the calcium ion, in addition to the two water molecules approaching the phosphate group. Under these conditions, there was an increase in the activation energy of 2.2 kcal mol⁻¹ for PS6 and of 0.4 kcal mol⁻¹ for PC6. As such changes remain small, we thus assumed that inclusion of two water molecules around the calcium ion is sufficient to provide an accurate estimation of the activation free energy for the hydrolysis reaction.

The cartesian coordinates of the atoms in each of the reactants, transition states and products of the reactions described below can be found in the Supplementary Information that accompanies this article.

3.2.1. Hydrolysis of phosphatidylserine (PS6). For conciseness, we refer to the phosphate group as pho, the serine moiety as ser, and the remaining glycerol-associated group, including the hydrophobic chain, as gly.

The hydrolysis of phosphate diesters around pho can proceed either through the cleavage of P–O bonds or the adjacent C–O bonds. While C–O bonds generally hydrolyse more rapidly than P–O bonds,¹⁷ the hydrolysis of C–O bonds directly linked to pho (indicated by blue arrows in Fig. 3) in phosphate diesters, such as those considered in this work, is known to be difficult because the attack on the bond is sterically hindered.³³ Therefore, only the C–O bonds in gly (indicated by the green arrows in Fig. 3) are relevant in the present context.

Fig. 3 Schematic representation of the PS6 molecule. The blue arrows indicate C–O bonds adjacent to P–O bonds, whereas green arrows point to C–O bonds linking the aliphatic chains to the glycerol backbone.

As described in the Methods section, we first employed Matlantis software to analyse the reaction dynamics. After optimising the reactant structure, we scanned the distance between atoms relevant to the target reaction. More specifically, the distance between the O atom of the approaching water molecule and the P atom of PS6. This approach is required because transition states are often difficult to locate by intuitive guessing alone. In addition, relaxed scans of the atomic distances using the Gaussian method are extremely time-consuming.

Fig. 4 shows a portion of the optimised structure of PS6 (a PS molecule with six carbon atoms in the aliphatic chains) in an aqueous solution, highlighting the atoms of ser and pho. This structure was refined using Gaussian16 based on the preliminary optimisation obtained from Matlantis. The slight differences between the structures obtained by Gaussian16 and Matlantis can be attributed to the phase assumptions, that is, Matlantis assumes gas-phase reactions, whereas Gaussian16 can be adjusted for analysis under aqueous conditions (PCM approximation). Hereafter, because the structures obtained using Matlantis consistently resembled those obtained using Gaussian16, we report only the results refined using Gaussian16.

Fig. 4 Optimised structure of a PS6 molecule in water. The structure is the same as that shown in Fig. 1 but here, for clarity, the atoms belonging to the glycerol group and the carbon chains are omitted (and denoted as “glycerol”). The group linked to the phosphate group is serine; one H atom of serine appears to be loosely linked to the N atom.

First, we investigated the hydrolysis of the C–O bonds in gly. In a PS6 molecule, which contains two aliphatic chains, two C–O bonds in gly can undergo hydrolysis. The calculated activation free energies for the cleavage of the first and second C–O bonds were 39.8 and 39.7 kcal mol⁻¹, respectively.

Next, we investigated the reactions leading to hydrolysis of the bonds around the phosphate group of a single PS6 molecule in water. The reactions proceeded in four steps and followed two distinct routes, denoted as R1 and R2.

In R1, hydrolysis begins with the cleavage of the P–O bond linked to ser, whereas in R2, the reaction starts with the cleavage of the P–O bond linked to gly. In phospholipid mineralisation experiments, the carbon chains linked to gly are expected to adhere to the surface of the glass-bottom dish, which is hydrophobic, according to the manufacturer's information. Consequently, if R1 proceeds, only ser is expected to be released into the aqueous solution, whereas if R2 proceeds, the entire pho–ser unit is released into the solution.

The first step in both R1 and R2 involves the approach of two water molecules to the P atom of pho, resulting in the association of an OH⁻ group with the P atom, which is coordinated to five oxygen atoms, and a proton (H⁺) to one of the O atoms of pho. In this step, one water molecule donates the OH⁻ group to the phosphorus atom, whereas the second water molecule facilitates proton exchange, that is, it donates a proton to one of the O atoms in pho and simultaneously accepts a proton from the first water molecule.

Water molecules can approach the P atom from four different directions: (1) opposite to ser, (2) opposite to gly, (3) opposite to the P–O bond close to the NH₃⁺ moiety, or (4) opposite to the remaining P–O bond (see Fig. 4). However, Westheimer's rules¹⁸ predict that, when ser/gly are the leaving groups, the attacking H₂O should approach from the side opposite to the leaving group.

Although Matlantis assumed a gas-phase environment, the software successfully located associative transition states. Fig. 5 depicts the transition state of the first step of the route, refined using Gaussian16: the OH association resulting from two water molecules approaching P from the side opposite to ser. At the transition state, the distance between P and the approaching O (2.00 Å) is longer than the P–O bond to ser (1.74 Å), showing that ser does not leave in this step, and thus confirming the formation of the pentacoordinated intermediate. In the optimised structure of the intermediate, the P–O bond to ser elongates (1.79 Å) and becomes longer than the P–O just formed (1.77 Å), indicating progression toward bond cleavage.

Fig. 5 Transition state of the first step in R1: OH association resulting from approach of two water molecules from the side opposite to serine. The OH group at the right side of P associates at equinox position. The dashed P–O is the bond being formed, and the dashed O–H bonds are those being broken.

In the second step of R1, the hydrolysis of ser proceeds through proton transfer from the NH₃ group of ser to the O atom of the P–O bond (the transition state is shown in Fig. 6). As shown in Fig. 4, the H atom in the NH₃ moiety of ser is weakly bound to the N atom and can easily migrate towards the P–O bond, thereby promoting the release of ser from PS6.

Fig. 6 Transition state of the second step of the hydrolysis of serine: proton transfer from NH₃⁺ to the O atom of the breaking P–O bond (dashed line) leading to hydrolysis of the serine group.

Following the release of ser into the solution, the third step may proceed: a water molecule associates with the P atom in a way analogous to the first step, but now from the side opposite to gly. Again, this results in the formation of a pentavalent phosphate intermediate, where the newly added OH group and the P–O bond to gly are at the equinox positions, and the other three O atoms lie in the same plane at the equatorial positions. Although this step proceeded via a mechanism similar to that of the first step, the activation energy was slightly lower (33.3 kcal mol⁻¹). Eventually, in the fourth step, the P–O bond linked to gly is hydrolysed by proton transfer via another approaching water molecule (Fig. 7), leading to the formation of H₃PO₄ (phosphoric acid).

Fig. 7 Transition state of the fourth step: hydrolysis of the glycerol group by reaction of the pentavalent intermediate with an approaching water molecule (coming from below, left side). The dashed P–O (vertical) and O–H are the bonds being broken.

In R2, the P–O bond linked to gly is first cleaved. The first step involves OH association with the P atom, similar to the first step of R1, but the water molecules now approach from the side opposite to gly. The second step of R2 entails proton transfer mediated by an intermediate water molecule through a transition state whose structure is similar to that shown in Fig. 5. However, in R2, the proton to be transferred must escape the original position, between NH₂ and the P–O bond (Fig. 4), and rotate by approximately 105° before being transferred to the approaching water molecule, leading to cleavage of the P–O bond linked to gly. Because the pentavalent intermediate is unstable, it may revert to water and PS6 rather than proceeding to the second step, hindering further progression along R2. The third step of R2 involves an OH association from the side opposite to ser, and the fourth step is the cleavage of the P–O bond linked to ser via proton transfer from NH₃.

The lines marked with circles in Fig. 8 show a comparison of the free-energy variations through R1 and R2 for PS6. These reaction dynamics are consistent with previous reports.^16–18 The activation free energies of the first step in R1 and R2 are 36.0 and 41.3 kcal mol⁻¹, respectively. The transition state energy of the second step in R2 is higher than that of the first step and requires reorientation of the H atom. The energy of the transition state of the fourth step of R2 was also high, which further disfavoured this route. The apparent discontinuity between the first and second hydrolyses arises because the byproducts from the first hydrolysis, that is, ser (R1) or gly (R2), were omitted in subsequent calculations. Instead, the energy of ser or gly in the solution was added separately to the total energy of the system during subsequent steps. Notably, the activation energy of the first step in R2 is comparable to that of the hydrolysis of the C–O bond between the aliphatic carbon chain and the glycerol backbone. These results suggest that, when a single PS6 molecule in solution is considered, hydrolysis is preferentially initiated by the cleavage of ser, followed by the hydrolysis of gly. Eventually, phosphoric acid and ser are released into the solution while the carbon chain–glycerol moiety remains adsorbed onto the glass surface.

Fig. 8 Energy variation through the four steps in the hydrolysis of PS6 and PC6. R1 is the route where the bond to serine or choline is hydrolysed first, and R2 is the route where the bond to glycerol is hydrolysed first.

3.2.2. Hydrolysis of phosphatidylcholine (PC6). First, the hydrolysis of the C–O bonds in gly of the PC6 molecule was investigated. The calculated activation energies for the cleavage of the first and second C–O bonds in PC6 were 32.8 and 36.4 kcal mol⁻¹, respectively, smaller than in the case of PS6 (which were approximately 40 kcal mol⁻¹). This result indicates that the aliphatic chains are more easily hydrolysed from PC6 than from PS6.

The lines marked with triangles in Fig. 8 show the energy landscape of the stationary states involved in the reaction dynamics during the hydrolysis of a single PC6 molecule in water. These four steps are similar to those described for PS6. For conciseness, we use cho for choline. Notably, the second and fourth steps of both R1 and R2 require a ∼180° rotation of the H atom linked to a P–O bond of pho from its state at the end of the first and third steps. Although the energy required for this H atom rotation is relatively small (less than 10 kcal mol⁻¹), the need for this conformational change implies that the hydrolysis reactions require more than four steps. In addition, rotation must occur after formation of the unstable pentavalent P intermediate.

Considering the activation energies, both R1 or R2 are eligible routes. Nevertheless, although the energy of the third transition state in R2 was higher than the corresponding value in R1, the activation energy (i.e., the energy difference between the transition state and the preceding structure) was smaller in R2 than in R1.

From the data shown in Fig. 8, it can be predicted that the hydrolysis of PC6 proceeds through both R1 and R2, with R2 being slightly preferred, because the energy of TS4 is higher in R1. However, because the hydrolysis of cho and gly from PC6 proceeds with activation energies comparable to those of the hydrolysis of the C–O bonds linking the aliphatic chains to the glycerol group, the carbon chains may be hydrolysed first, releasing the moiety cho–pho–glycerol to the aqueous solution.

3.2.3. Hydrolysis of PS6 in the presence of Ca²⁺ ions. When Ca²⁺ ions are present in large amounts in aqueous solution, they may link to one or two O atoms bonded to the central P atom of PS6, owing to electrostatic attraction. We first optimised a structure in which the Ca²⁺ ion was close to PS6; the distance between the Ca²⁺ ion and the P atom was 3.71 Å. We calculated the change in potential energy when Ca²⁺ ions moved away from the PS6 molecule through a relaxed scan to evaluate the binding affinity between Ca²⁺ ions and PS6. The energy increased gradually by approximately 14 kcal mol⁻¹ as the ion-PS6 distance increased to 5.2 Å, after which the energy became nearly constant, indicating a weak interaction beyond this point. Therefore, the subsequent hydrolysis reactions of PS6 in the presence of Ca²⁺ ions were calculated by assuming a structure in which a Ca²⁺ ion was linked to one of the O atoms of the pho group. Two additional water molecules were explicitly considered for stabilizing the ion.

In the presence of Ca²⁺ ions, the activation free energy for the hydrolysis of the first C–O bond in gly of the PS6 molecule increased to 41.3 kcal mol⁻¹.

The corresponding reaction pathways and energy profiles for the hydrolysis of ser and gly are shown by the lines marked with circles in Fig. 9. Notably, the activation energies were smaller than those for the hydrolysis of the C–O bonds in gly.

Fig. 9 Energy variation through the steps in the hydrolysis of PS6 and PC6 in the presence of Ca²⁺ ions. R1 is the route where the bond to serine or choline is hydrolysed first, whereas R2 is the route where the bond to glycerol is hydrolysed first.

As in the case without Ca²⁺, the hydrolysis proceeds through multiple steps and two different routes; however, now the activation energy required for hydrolysis of the P–O bond linked to gly (R2) is markedly reduced.

In contrast, when ser was hydrolysed first (R1), the free energy variation along the first two steps was nearly identical to that observed in the absence of Ca²⁺ (Fig. 8). However, when the P–O bond to gly breaks after the P–O bond to ser, the energy barrier is substantially higher than that calculated in the absence of Ca²⁺ ions. These results indicate that when Ca²⁺ ions are present, R2 occurs preferentially to R1.

3.2.4. Hydrolysis of PC6 in the presence of Ca²⁺ ions. A relaxed scan of the distance between a Ca²⁺ ion and PC6 was performed to evaluate the affinity of PC6 for Ca²⁺ ions. The results showed that the potential energy reached a minimum at a P–Ca²⁺ distance of 2.99 Å, which is considerably shorter than that observed for PS6. As the Ca²⁺ ion moved away from the P atom, the energy gradually increased, exceeding a 20 kcal mol⁻¹ rise at approximately 7 Å. The shorter P–Ca²⁺ distance and the larger energy increase upon separation of Ca²⁺ from P indicate that the Ca²⁺ ions exhibit stronger affinity with PC6 than with PS6, which is likely due to the higher dipole moments of PC6 (22.6 Debye) compared with PS6 (17.1 Debye).

The activation free energy of the hydrolysis of a C–O bond between an aliphatic carbon chain and the remainder of PC6 in the presence of Ca²⁺ ions was 35.3 kcal mol⁻¹, which was smaller than that in the case of PS6 in the presence of Ca²⁺ ions (approximately 41 kcal mol⁻¹).

The hydrolysis reactions of PC6 in the presence of Ca²⁺ ions followed the routes shown by the lines marked with triangles in Fig. 9. These results suggest that R2 likely undergoes hydrolysis first. However, it should be reiterated that the hydrolysis of the P–O bond to gly requires a 180° rotation of an H atom. In comparison with the Ca²⁺-free case (Fig. 8), the presence of Ca²⁺ ions tend to facilitate the hydrolysis process, as evidenced by a modest decrease in the activation energies. Again, hydrolysis of the aliphatic carbon chain from the glycerol backbone may proceed first, because of its lower activation energy.

3.2.5. Hydrolysis of phospholipid–Ca²⁺–phospholipid systems. In aqueous solution, phospholipids associate with Ca²⁺ ions to form networks. We designed complexes of two PS6 or PC6 molecules coordinated to a single Ca²⁺ ion to model the underlying reaction mechanisms (Fig. 10 shows a schematic view of the complex PS6–Ca²⁺–PS6). Evidence for the formation of this type of complex has been observed experimentally, at least for high concentrations of calcium ions and PC.³⁴ Such complexes capture the key chemical events that are not observed in single-molecule systems. Depending on which bond is hydrolysed first—indicated by the blue or green arrows in Fig. 10—either a single serine is released into solution (ser–pho break, blue arrows) or the entire ser1–pho1–Ca²⁺–pho2–ser2 unit remains near the surface, anchored by the opposite glycerol group (gly–pho break, green arrows).

Fig. 10 Schematic view of the complex PS6–Ca²⁺–PS6. For simplicity, only two carbon atoms are shown in the aliphatic chains instead of five. The blue arrows point to bonds linking P to serine groups and green arrows point to bonds linking P to glycerol groups.

3.2.5.1. Hydrolysis of PS6–Ca²⁺–PS6. Here, the activation energy for C–O bond hydrolysis between the aliphatic chain and core structure increased to 44.1 kcal mol⁻¹.

In the case of a single PS6 molecule, R1 hydrolysis proceeded more readily than that of R2, because the NH₃ group is positioned sufficiently close to efficiently transfer a proton to the ser–pho bond. Similarly, this preference persists in the PS6–Ca²⁺–PS6 system. Surprisingly, however, the required free energy is less than 30 kcal mol⁻¹, which is substantially lower than the energy observed for a single PS6 molecule (Fig. 8 and 9). In addition, the presence of two phosphate (pho1 and pho2), serine (ser1 and ser2), and glycerol (gly1 and gly2) groups allows for unique interactions. Notably, the NH₃ group of ser1 may transfer a proton to the bond linking gly2 to pho2 owing to the favourable spatial arrangement.

The lines marked by circles in Fig. 11 show that the hydrolysis of ser1 was slightly more favourable than that of gly1. If ser1 is hydrolysed first, the subsequent cascade can occur at pho1–gly1, pho2–ser2, or pho2–gly2. For clarity, only the most favourable case, the hydrolysis of the bond between pho2 and ser2, is depicted. Similarly, if gly1 is cleaved first, although various reactions are possible, only the most favourable, the pho2–gly2 hydrolysis, is shown.

Fig. 11 Energy variation throughout the hydrolysis steps of PS6–Ca²⁺–PS6 and PC6–Ca²⁺–PC6.

These results suggest that the difference in energy variation for the sequences of hydrolysis: ser1 → ser2 or gly1 → gly2 is not large, indicating similar probabilities for each pathway. However, after the first two steps, the energy divergence increased. Specifically, the sequence gly1 → gly2 → ser1 → ser2 may occur with a probability higher than ser1 → ser2 → gly1 → gly2. Assuming that the hydrophobic tails of the complex PS6–Ca²⁺–PS6 are attached to the glass surface, gly1 → gly2 pathway would result in the release of the hydrophilic complex ser1–pho1–Ca²⁺–pho2–ser2 to the aqueous solution.

This difference in hydrolysis sequence impacts the mineralisation location: the gly1 → gly2 → ser1 → ser2 cascade facilitates pho release and subsequent mineralisation in the solution. In contrast, the ser1 → ser2 → gly1 → gly2 pathway retains the two pho groups close to the glass surface.

Notably, hydrolysed ser, now bearing a NH₂ group, may act as a proton acceptor and facilitate reactions subsequent to ser1 → ser2, for instance. The increase in the mineralization of PS observed experimentally under higher pH values (Fig. 2C) can be interpreted as evidence supporting this effect. As hydrolysed ser residues become proton acceptors, an increase in pH (higher concentrations of OH⁻ ions) would enhance proton mobility away from the growing calcium-phosphate complex. The specific reactions undergone by hydrolysed ser residues should be further investigated.

3.2.5.2. Hydrolysis of PC6–Ca²⁺–PC6. The PC6–Ca²⁺–PC6 complex, modelled similarly to the PS6 system but with cho replacing ser, exhibited an activation energy for the hydrolysis of a C–O bond between an aliphatic carbon chain and the core structure of 36.2 kcal mol⁻¹, which is smaller than that calculated for PS6–Ca²⁺–PS6 (approximately 44 kcal mol⁻¹).

As shown by the lines marked with triangles in Fig. 11, a substantial decrease in free energy occurs after cho1 is hydrolysed, probably because the hydrolysed cho1 becomes positively charged, leaving an electron behind. Surprisingly, this decrease is greater than that observed for monomeric PS6–Ca²⁺ (Fig. 9).

The next step involving the hydrolysis of cho2 follows an unusual pathway: after proton transfer initiates the formation of a pentacoordinated phosphate at pho2, cho2 is cleaved without the involvement of a transition state: a proton from a water molecule is transferred to cho2 while the water molecule itself receives a proton from pho1. Such reaction dynamics allow the cho1 → cho2 sequence to proceed smoothly, leaving two electrons behind. However, no further decrease in free energy was observed after cho2 hydrolysis, possibly because the remaining complex gly1–pho1–Ca²⁺–pho2–gly2 becomes neutral.

Following the release of the two cho groups, gly1 and gly2 differ only in their geometrical orientation; thus, their hydrolysis reactions are similar to each other.

Overall, the most probable hydrolysis sequence for the PC6–Ca²⁺–PC6 complex is cho1 → cho2 → gly1 → gly2, suggesting that the two pho groups and the Ca²⁺ ion would remain close to the glass surface, at least for some time. However, aliphatic carbon chains may be cleaved first because the activation energy for the hydrolysis of the C–O bond between them and the glycerol backbone is smaller. It should be noted that the low sensitivity of PC mineralization to pH observed experimentally (Fig. 2C) is consistent with the computational results because hydrolysed cho is a positive ion and would associate with the OH⁻ ions made more abundant at higher pH, rather than facilitating proton transfer processes that promote mineralisation.

This is consistent with experimental observations: solutions containing PC appeared opaque, whereas those containing PS remained transparent. If the aliphatic carbon chains were hydrolysed first, the large hydrophilic block ([glycerol backbone]–pho–cho)₂Ca²⁺ is released into the solution, which can be expected to contribute to the solution opacity.

Furthermore, as the presence of Ca²⁺ lowered the activation energy for the hydrolysis of a single PC6 molecule, it is possible that, in complexes larger than PC6–Ca²⁺–PC6, where the ratio of Ca²⁺ ions to PC6 molecules exceeds 1 : 2, the activation energies for the hydrolysis of the relevant P–O bonds would be further reduced. Therefore, larger systems should be investigated in future.

3.2.5.3. Comparison between PS and PC. The computational results are summarised in Table 2. From the discussion above, PS can be hydrolysed more easily than PC primarily because of the labile H⁺ on the NH₃⁺ group in ser, which facilitates proton transfer to P–O bonds. This lowers the activation barriers and promotes hydrolysis. This mechanistic insight aligns with the experimental observations showing faster mineralisation of PS. Additionally, hydrolysed ser may act as a proton acceptor, enhancing subsequent reactions.

Table 2 Activation free energies (E_ac, in kcal mol⁻¹) of the most probable reactions in the phospholipid hydrolysis

System	1st bond	Eac	2nd bond	Eac
One molecule
PS6	P–O_ser	36.0	P–O_gly	32.1
	C–O@gly	39.8	C–O@gly	42.5
PS6–Ca^2+	P–O_gly	35.7	P–O_ser	30.1
	C–O@gly	41.3
PC6	P–O_cho	38.8	P–O_gly	34.5
	C–O@gly	32.8	C–O@gly	36.4
PC6–Ca^2+	P–O_gly	34.7	P–O_cho	33.7
	C–O@gly	35.3
Two molecules
PS6–Ca^2+–PS6	P–O_ser1	28.1	P–O_ser2	25.5
	P–O_gly1	31.1	P–O_gly2	24.3
	C–O@gly	44.1
PC6–Ca^2+–PC6	P–O_cho1	37.9	P–O_cho2	34.9
	P–O_cho1	37.9	P–O_gly2	36.4
	C–O@gly	36.2

---

In contrast, in systems containing PC molecules, the hydrolysis of the C–O bonds linking the aliphatic carbon chains to the glycerol backbone competes with the hydrolysis of the P–O bonds. Therefore, the aliphatic carbon chains may be hydrolysed first.

Although our findings support the notion that hydrolysis is a key initial step driving the faster mineralisation of PS, further investigation into downstream mechanisms and consideration of larger systems are needed to fully understand the differences in mineralisation efficiency among phospholipids.

4. Conclusions

In this study, we examined the differences in the mineralisation efficiencies of PS and PC through experiments and quantum chemical calculations.

Experimentally, PS has been found to mineralise more than PC. Computational analysis suggests this difference is primarily due to the NH₃⁺ group in PS, which contains a loosely bound proton that can be readily transferred to P–O bonds, either at phosphate–serine or phosphate–glycerol, thereby lowering the activation energy for hydrolysis. Furthermore, after being released into the aqueous solution, hydrolysed serine contains an NH₂ group that is highly capable of accepting protons. This capability may facilitate further reactions; therefore, the hydrolysed serine can participate in subsequent reactions. In contrast, PC is more likely to undergo hydrolysis at the C–O bonds linking the aliphatic chains to the glycerol backbone, before cleaving the P–O bonds in phosphate–choline or phosphate–glycerol.

These findings highlight the unique reactivity of PS, which may underlie its superior mineralisation behaviour. Further studies are needed to validate the downstream effects of the serine proton-accepting capacity and to elucidate the mechanisms contributing to phospholipid-mediated mineralisation to guide the development of more efficient bone-repair-inducing molecules.

Author contributions

K. S. performed the final calculations related to the large phospholipid–Ca²⁺ ion complexes and participated directly in writing the manuscript; R. K., T. S., and Y. C. participated in the initial phase of the implementation of the experiments/calculations; K. Y. participated in the analysis/comparison of the results of experiments and calculations, and in the revision of the manuscript.; E. S. H. and N. K. were both involved in the conceptualisation, data curation, formal analysis, funding acquisition, investigation, methodology, and writing of the manuscript (E. S. H. supervised the experiments; N. K. supervised the calculations).

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI): the cartesian coordinates of the atoms in all stationary points, optimized through Gaussian16 at the B3LYP/6-31G(d,p) level of theory, along the potential energy surfaces related to the chemical reactions investigated are given in the Supporting Information file. See DOI: https://doi.org/10.1039/d5tb01744k.

Acknowledgements

The authors acknowledge the support from JSPS KAKENHI Grants (No. JP21K09963 and JP25K13073), and the Japan Science and Technology Agency, JST, FOREST Program (Grant Number JPMJFR210X). This work was also partly supported and conducted at the Joint Research Center for Environmentally Conscious Technologies in Materials Science (Projects No. 31018 and 02003, Grant No. JPMXP0618217637) at ZAIKEN, Waseda University.

References

1. J. T. Irving, Nature, 1958, 181, 704–705.
2. R. E. Prout, A. A. Odutuga and F. C. Tring, Arch. Oral Biol., 1973, 18, 373–380.
3. E. S. Hara, M. Okada, N. Nagaoka, T. Hattori, T. Kuboki, T. Nakano and T. Matsumoto, ACS Biomater. Sci. Eng., 2018, 4, 617–625.
4. E. S. Hara, M. Okada, T. Kuboki, T. Nakano and T. Matsumoto, J. Mater. Chem. B, 2018, 6, 6153–6161.
5. E. S. Hara, S. Oozawa, N. Nagaoka, M. Okada, S. Kasahara and T. Matsumoto, Mater. Adv., 2021, 2, 4423–4431.
6. J. M. Cotmore, G. Nichols, Jr. and R. E. Wuthier, Science, 1971, 172, 1339–1341.
7. R. E. Wuthier, Calcif. Tissue Res., 1975, 19, 197–210.
8. M. Bosetti, M. Santin, A. W. Lloyd, S. P. Denyer, M. Sabbatini and M. Cannas, J. Mater. Sci.: Mater. Med., 2007, 18, 611–617.
9. M. Bosetti, A. W. Lloyd, M. Santin, S. P. Denyer and M. Cannas, Biomaterials, 2005, 26, 7572–7578.
10. J. Hatakeyama, H. Anan, Y. Hatakeyama, N. Matsumoto, F. Takayama, Z. Wu, E. Matsuzaki, M. Minakami, T. Izumi and H. Nakanishi, J. Oral Sci., 2019, 61, 111–118.
11. A. Merolli, M. Bosetti, L. Giannotta, A. W. Lloyd, S. P. Denyer, W. Rhys-Williams, W. G. Love, C. Gabbi, A. Cacchioli, P. T. Leali, M. Cannas and M. Santin, J. Mater. Sci.: Mater. Med., 2006, 17, 789–794.
12. K. Noda, Y. Jagawa, A. Fuwa and N. Kunioshi, J. Phys. Chem. A, 2022, 126, 8658–8673.
13. K. Anzai, N. Kunioshi and A. Fuwa, Appl. Surf. Sci., 2017, 392, 410–417.
14. J. L. Medina-Franco and F. I. Saldívar-González, Biomolecules, 2020, 10, 1566.
15. N. Flores-Holguín, J. Frau and D. Glossman-Mitnik, Sci. Rep., 2022, 12, 506.
16. S. C. L. Kamerlin, N. H. Williams and A. Warshel, J. Org. Chem., 2008, 73, 6960–6969.
17. R. Wolfenden, C. Ridgway and G. Young, J. Am. Chem. Soc., 1998, 120, 833–834.
18. F. H. Westheimer, Acc. Chem. Res., 1968, 1, 70–78.
19. S. Mikkola, T. Lönnberg and H. Lönnberg, Beilstein J. Org. Chem., 2018, 14, 803–837.
20. Matlantis (https://matlantis.com/), software as a service style material discovery tool.
21. S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee and T. Ibuka, Nat. Commun., 2022, 13, 2991.
22. S. Takamoto, D. Okanohara, Q. J. Li and J. Li, J. Materiomics, 2023, 9, 447–454.
23. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision B.01, Gaussian, Inc., Wallingford CT, 2016.
24. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
25. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
26. H. Fahmi, J. Am. Chem. Soc., 2017, 139, 6780–6786.
27. W. J. Hehre, L. Radom, P. R. Schleyer and J. A. Pople, Ab Initio Molecular Orbital Theory, John Wiley and Sons, New York, 1986, pp. 76–82.
28. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
29. A. Austin, G. Petersson, M. J. Frisch, F. J. Dobek, G. Scalmani and K. Throssell, J. Chem. Theory Comput., 2012, 8, 4989–5007.
30. J. Foresman and A. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, Inc, Wallingford, CT, 3rd edn, 2015, p. 11.
31. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
32. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
33. G. K. Schroeder, C. Lad, P. Wyman, N. H. Williams and R. Wolfenden, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 4052–4055.
34. C. Altenbach and J. Seelig, Biochemistry, 1984, 23, 3913–3920.

---

Footnote

† The two corresponding authors contributed equally from design of the research to production of the manuscript.

---