A bioinformatics approach to study the role of calcium phosphate properties in bone regeneration

ABSTRACT

Calcium phosphates (CaPs) are extensively used as biomaterials for bone repair and regeneration and are considered an interesting synthetic alternative for natural bone grafts. However, the bone regenerative potential of many CaP bone graft substitutes requires further improvements. The bioactivity of CaPs is regulated by a variety of properties, including their chemical composition and (surface) structural properties. In this study, we performed a comprehensive characterization of six CaP ceramics that had been previously tested in several in vitro and in vivo studies, allowing us to in silico distinguish osteoinductive and non-osteoinductive CaP ceramics. The numerical parameters obtained from this characterization were further correlated with osteogenesis-associated gene expression and pathways of osteoblastic MG63 cells cultured on these six CaP ceramics. The results showed, among others, the influence of both chemical and structural parameters on cytoskeleton- and extracellular matrix-related gene expression and osteogenesis-associated pathways such as MAPK and BMP. The bioinformatics approach presented here allows exploring the transcriptional landscape induced by CaP ceramics and gaining fundamental knowledge about cell-biomaterial interactions.

KEYWORDS:

• Calcium phosphate
• Material characterization
• Transcriptomics
• Cluster analysis
• Correlation analysis
• Regenerative medicine
• MG63

1. Introduction

Calcium phosphate (CaP)-based biomaterials are widely used in orthopedic, craniomaxillofacial, and dental surgery for the repair and regeneration of bone defects. Despite significant improvements over the past decades and examples of CaP bone graft substitutes with an in vivo performance comparable to that of autologous bone (Yuan et al., 2010), the gold standard in bone regeneration, synthetic bone graft substitutes are in general still considered inferior to the natural bone grafts. The advantages of their use are, however, evident: they are relatively inexpensive, easy to produce in virtually unlimited amounts, and can be stored for long periods. Moreover, with bone graft substitutes, complications associated with autograft harvesting and transplantation can be avoided. In view of an increasing demand for successful treatment options for challenging bone defects (Gillman & Jayasuriya, 2021; Govoni et al., 2021), it is important to further improve the clinical performance of synthetic bone graft substitutes.

The bioactivity of CaPs, i.e., their bone regenerative potential, is influenced by a large and complex set of physico-chemical and structural properties. These properties include the chemical composition, specifically calcium, phosphorus, and other bioinorganics (Habibovic & Barralet, 2011), the phase composition (Danoux, Pereira, et al., 2016), crystallinity (Lee, Zavgorodniy, Loo, & Rohanizadeh, 2012), microstructure (Deligianni, Katsala, Koutsoukos, & Missirlis, 2000), grain size (Lapczyna et al., 2014), and porosity, among others. The most direct way to assess the contribution of an individual material property to a biological response would be to produce batches of CaPs differing in just one property at a time, and then evaluate their performance either in vitro or in vivo. However, this is difficult to achieve, as the properties of CaPs are largely intertwined, i.e., a change in one of them causes changes in others, particularly when conventional synthesis and processing techniques are used (Galvan-Chacon & Habibovic, 2017). For example, doping a CaP with an additive, such as strontium, results in the replacement of some or all calcium ions by the strontium ion in the crystal lattice, which in turn alters the crystallinity, microstructure or solubility of the material (Birgani, Malhotra, van Blitterswijk, & Habibovic, 2016). It is therefore difficult to distinguish between the direct chemical effect of the added ion and the changes its addition causes in other material properties (Birgani, 2016).

The knowledge of how each individual material property affects the interaction of the material with a biological environment is, however, crucial for developing more efficient bone graft substitutes. A number of studies addressed this issue by developing novel methods to decouple the importance of chemistry and surface structure, such as coating methods to isolate the effect of surface (micro)structure (dos Santos, Farina, Soares, & Anselme, 2009), methods to replicate the surface structure of a CaP into a polymer that does not contain the ceramic (Danoux, Sun, et al., 2016; Groen et al., 2017), the controlled modification of the surface structure (Deligianni et al., 2000; Mazón, García-Bernal, Meseguer-Olmo, Cragnolini, & De Aza, 2015), or the use of composite materials containing inorganic constituents of the ceramic rather than the ceramic itself (Danoux et al., 2015). Although these studies have made an important contribution to the current understanding of the property-function relationships of CaPs, the methods they used are not sufficient to understand the role of all properties and neglect the relevant relationships between them (Engel et al., 2008).

A different approach consists of applying bioinformatics tools based on big data analysis to establish correlations between material properties and biological responses. For example, gene expression microarray technologies allow measuring thousands of protein-encoding genes at a time, providing very large and comprehensive datasets. This approach has been successfully applied on ceramic and composite materials to understand how physico-chemical stimuli influence cell behavior (Baker et al., 2014; Eyckmans et al., 2013). For example, Groen et al. linked the transcriptional profile of MG63 cells, an osteogenic cell line frequently used to assess the osteogenic potential of ceramics (Lincks et al., 1998), to the CaP material properties by culturing the cells on osteoinductive and non-osteoinductive ceramics, as well as replicas of the structure of those materials fabricated with different chemical compositions (Groen et al., 2017). The results identified genes responding to microstructure, calcium, and phosphate ion levels.

In this study, we used a big data approach to link the properties of CaP ceramics with the gene expression levels of cells cultured on them. To complement previous studies, we performed a highly comprehensive material characterization, which, combined with a machine learning approach, allowed us to study the role of individual material properties in the biological response to the materials. The six CaPs selected, comprised pairs of materials with similar starting chemical composition but sintered at different temperatures, allowing to study the effect of similar chemical properties but different microstructures. These ceramics have been used in a large number of previous in vitro and in vivo studies (Barrère-de Groot, Yuan, & de Bruijn, 2018; Davison, Yuan, de Bruijn, & Barrere-de Groot, 2012; Davison et al., 2014; Duan et al., 2017; Groen et al., 2017; Habibovic et al., 2006; Li, De Wijn, Li, Layrolle, & De Groot, 2003; Othman et al., 2019; Wang et al., 2015; Yuan, Yang, De Bruij, De Groot, & Zhang, 2001; Yuan, van Blitterswijk, de Groot, & de Bruijn, 2006; Yuan et al., 2010), facilitating additional interpretation of the results by comparing their performance with previous data.

The comprehensive characterization of the physico-chemical and microstructural properties of the six materials resulted in a total of 131 numerical parameters. The characterization techniques included X-ray diffraction analysis (XRD) with Rietveld analysis, inductively-coupled plasma mass spectrometry (ICP-MS) ion quantification to determine both the materials composition and its ions exchange dynamics with cell culture medium, scanning-electron microscopy (SEM), mercury intrusion, porosimetry, confocal laser profilometry and micro-computed tomography (microCT). Further clustering and correlation analyses were then performed using a transcriptomics dataset of MG63 cells cultured on the six materials obtained in a previous study (Groen et al., 2017) to identify links between materials parameters and osteogenesis-associated gene expression and pathways.

The approach presented here allows exploring the transcriptional landscape induced by ceramics, identifying which material properties are linked with particular genes or gene pathways, and represents a valuable tool for gaining fundamental knowledge about cell-biomaterial interactions in the context of bone regeneration.

2. Materials and methods

2.1. Materials characterization

2.1.1. Ceramic particles

The six CaP ceramics used in this study were supplied by Kuros Biosciences BV in the shape of particles with a size of 1–2 mm. The materials were produced with three different chemical compositions: HA, TCP, and BCP. For each chemistry, two materials were produced with different microstructure: small/large grain size, low/high microporosity, and low/high specific surface area.

2.1.2. Chemical characterization

The ceramic particles were ground in an agate mortar, reduced to fine powder, loaded into a plastic container, and analyzed using X-ray diffraction (XRD, Bruker D2 Phaser) with a step size of 2θ 0.01°, one second per step and one rotation per second, in the interval 2θ 5–80°. The data was analyzed using Profex v.4.1.0 to perform a Rietveld refinement of the crystal structure (Doebelin & Kleeberg, 2015). The powders were also analyzed by Fourier Transform Infrared Spectroscopy (FTIR) in the attenuated total reflection mode (ATR). The spectra were analyzed using SpectraGryph 1.2 (Spectroscopy Ninja).

For the elemental composition analysis, 20 milligrams of particles of each material were dissolved in 2 mL of 60% ultrapure nitric acid (HNO₃, 60% ultrapure for trace analysis, Supelco, Merck) overnight, followed by the addition of 8 mL of 1% HNO₃. Afterwards, a further 20x dilution was performed in a matrix of 1% HNO₃ and gallium standard (Ga standard for ICP-MS, VWR) to achieve a final Ga concentration of 100 ppb, which was used as an internal standard for the quantification. A standard curve was prepared with 7 points using a multi-elements standard (Thermo Fisher Scientific) and Ca and P standards (standard solution, 1000 mg/L, Aristar), in order to detect concentrations of Ca and P between 0.1 ppm and 200 ppm, and the other elements contained in the Multi elements solution (Na, Mg, S, K, Fe, Ga, and Sr) in the range 0.1 to 5 ppm.

To analyse the ion exchange between the ceramics and the cell culture medium, an identical experimental approach was applied as for the transcriptomics data experiments (Groen et al., 2017). In brief, 2–3 particles (20 mg) were placed in triplicate per condition and time point in 48 well-plates and 0.4 ml cell culture medium (α-MEM, 10% fetal bovine serum, 0.2 mM ascorbic acid, 2 mM L-glutamine, 100 U/ml penicillin, 100 µg/ml streptomycin) was added to each well. The plates were incubated in a cell culture incubator (5% CO₂, humidity control), the medium was refreshed every two days, and collected and frozen on day 2, 4, 7, and 14. Then, the samples were prepared for the analysis by dissolving 100 µL of cell culture medium in 0.9 mL 60% nitric acid (HNO₃, 60% ultrapure for trace analysis, Supelco, Merck) overnight, followed by the addition of 9 mL 1% HNO₃ containing Ga standard (Ga standard for ICP-MS, VWR) to obtain a final Ga concentration of 100 ppb, which was used as an internal standard in all samples. The standard 6-point curve was prepared in the same matrix as the samples and used to quantify Ca, P and Mg in the range of 0.1 ppm to 2 ppm.

2.1.3. Structural characterization

The surface roughness parameters of the particles were characterized with a confocal laser scanning optical profilometer (Keyence VKX-200, Keyence, Japan) using a 150x objective. For each material, three different areas of 6897.4 µm² were imaged on three different particles (n = 3). The profilometry data was processed using the Keyence Multifile Analyzer software, obtaining multiple parameters recommended in ISO 25178 to quantify standard roughness (description in Supplementary Table 1). As reported previously by Groen et al. (Groen et al., 2017), the microstructure was also observed by an SEM (XL30, Philips) in the secondary electron mode. The grain size of the ceramics was determined by measuring the lateral distance of 60 different grains from three independent SEM images using Image J. Additionally, high-resolution SEM images were acquired (Teneo, FEI) at 5 kV voltage and 2000, 10,000, and 30,000 times magnification, using an Everhart Thornley Detector to collect secondary electrons and a Trinity Detector to collect backscattered electrons.

The ceramic particles were further analyzed by µCT (SkyScan 1272 11Mp scanner, Bruker Corporation, USA) with cone-beam geometry and a 4032x2688 CCD detector. The images were acquired using a 0.25 mm aluminum filter and software-correcting for alignment, thermal drift, beam hardening, and ring artifacts, using a voxel size of 2.25³ µm³. The image dataset was reconstructed using the NRecon v.1.6.10.4 software and analyzed using CTAn (both Bruker Corporation). The VOI was carefully selected to exclude the container in which the particles were imaged. Using the 3 D analysis tool of CTAn, multiple structural properties were quantified, including the total porosity, percentage of closed pores, and trabecular thickness.

The porosity of the particles was also analyzed by mercury intrusion porosimetry (Groen et al., 2017). Briefly, the pore size distribution, microporosity (percentage of pores smaller than 10 μm), and specific surface area were determined using mercury intrusion porosimetry (Quantachrome Instruments, Pore Master). Additional parameters concerning the shape of the pore distribution curve were quantified, and a fit to a pseudo-Voigt curve (convolution of a Lorentzian distribution and a Gaussian distribution) was calculated. Density measurements were performed to determine the true (skeletal) density of the materials using a helium pyconometer (Accu Pyc II 1340 gas pyconometer, Micrometrics) and analyzed with Pore Master software.

2.2. Cell culture and transcriptomics analysis

The dataset describing the transcriptomic profile of MG63 cells cultured on the different ceramics was obtained by Groen et al. in a previous study (Groen et al., 2017). Briefly, human osteosarcoma MG63 cells were seeded at a density of 200,000 cells in wells containing 150 µL particles, and cultured for 48 hours, using cells cultured on tissue culture polystyrene as a control. After cell culture, the total RNA was isolated, cDNA was synthesized, and hybridized on a microarray using Illumina HT-12 v4 expression Beadchips.Streptavidin Cy-3 was added to develop a fluorescent signal. The arrays were scanned using an Illumina microarray reader, and the raw intensity values were corrected for background using BeadStudio (Illumina).

2.3. Materials properties dataset processing and PCA analysis

The original dataset containing all materials parameters obtained from the materials characterization had dimensions of 6x161, where 6 is the number of different materials and 161 is the number of material parameters. Statistically significant differences between material parameters were established using a one-way ANOVA. An overview of statistical differences between material parameters can be found in the open access data repository. Material parameters that represented technical variations were excluded from further processing. Parameters that showed a strong correlation (>0.99) with another parameter were further excluded from the dataset. Incomplete material parameters (e.g., missing values) were also removed from the dataset. The final dataset obtained had dimensions of 6x76. PCA analysis was performed using R version 4.03. The hierarchical cluster plot was calculated in R version 4.03 by the complete linkage method and a Euclidian distance metric on the scaled data. Data visualization was achieved through Graphpad Prism version 9.3.1.

2.4. Gene dataset processing and material-gene correlation analysis

Heat map clustering of gene expression levels across the different materials was accomplished by selecting the top 500 genes demonstrating the highest variability across the materials. Variability is defined as the difference between a given value and the mean value. Processing and visualization were accomplished in R version 4.03. Visualization of the heat map was further processed using Graphpad Prism version 9.3.1.

For the gene-material correlation analysis, we used the same final dataset used for the PCA analysis. We further included genes belonging to the top 20% with the highest variability across the materials. Spearman correlation between genes and material properties were computed in R version 4.03 and visualized through R package ggplot version 3.1.1. Genes that demonstrated a correlation higher than 0.7 or lower than −0.7 and a p-value lower than 0.05 were considered significant. Gene ontology (GO) pathway analysis of each significant gene was performed in R version 4.03.

3. Results and discussion

3.1. Materials characterization

The main objective of the comprehensive materials characterization was to obtain a large set of numerical parameters describing the six CaP ceramics (HA1150, HA1250, TCP1050, TCP1100, BCP1150, BCP1300). The ceramics can be grouped into three categories based on their chemical composition (HAs, TCP1050, and BCPs) and per chemical composition category into two groups with different structural and topographical properties, resulting from differences in sintering conditions. Finally, the ceramics can be grouped into an osteoinductive (HA1150, TCP1050, BCP1150) and non-osteoinductive (HA1250, TCP1100, BCP1300) category, as shown in previous studies (Davison et al., 2014; Maazouz, Chizzola, Dobelin, & Bohner, 2021; Yuan et al., 2006; 2010). The complete list of parameters and a description of the least known ones is provided in Supplementary Table 1.

3.1.1. Chemical characterization

XRD analysis was performed to describe the crystalline composition and structure of the CaP ceramics. Figure 1a shows the XRD patterns of the ceramic powders, indicating the peaks corresponding to each of the phases present in the CaP materials, namely HA, β-TCP, and periclase. The Rietveld analysis details are provided in Supplementary Figure 1. For all conditions, the χ² parameter, which quantifies the difference between the experimental data and data predicted by the model, was below the acceptance limit defined as χ² < 1.5 (the closer the χ2 value is to 1, the better the fit). The results of the analysis are shown in Table 1 and represented in Figure 1b, c. During the data processing, we observed that the refinement did not give results below the acceptance limit, i.e., the model could not adequately describe the pattern, unless the possibility of Ca substitution by magnesium (Mg) was introduced into the code for refinement. The percentage of Mg substitution was calculated as the average for both CaP phases, i.e., over the total number of atoms present in the molecule without impurities. Mg is commonly introduced in the preparation of CaP, as it stabilizes the β-TCP phase and avoids the formation of α-TCP (Gibson & Bonfield, 2002; Kannan, Ventura, & Ferreira, 2007).

Figure 1. XRD analysis and Rietveld refinement. a) XRD patterns of the six CaP ceramics, with the peaks corresponding to the three main phases (HA, β-TCP, Periclase) indicated by arrows of different colors. b) Phase composition quantified by Rietveld refinement demonstrates the chemical composition differences between HA, TCP and BCP ceramics. c) Magnesium substitution percentage in the crystal lattice, obtained by Rietveld analysis. The error bars represent the estimated standard deviation in the quantification given by the model.

Table 1. Phase composition of the different materials.

Material	HA (wt.%)	β-TCP (wt.%)	Periclase (wt.%)	Mg substitution (%)
HA 1150	99.72 ± 0.09	–	0.3 ± 0.09	3.4 ± 0.9
HA 1250	99.5 ± 0.2	0.5 ± 0.2	–	1.7 ± 0.7
TCP 1050	11.5 ± 0.4	88.5 ± 0.4	–	3.3 ± 0.2
TCP 1100	1.5 ± 0.6	98.5 ± 0.6	–	2.7 ± 0.5
BCP 1150	86.4 ± 0.9	13.5 ± 0.9	<0.08	3.5 ± 0.2
BCP 1300	87.3 ± 0.3	12.7 ± 0.3	<0.08	3.43 ± 0.05

As expected, the two HA ceramics were largely phase-pure, with small percentages of periclase (Mg oxide) due to the presence of Mg and a small percentage of β-TCP in HA1250. TCP1050 contained 11.5 wt.% HA, whereas in TCP1100, the HA percentage was about 1.5 wt.%, due to the higher sintering temperature. Both biphasic ceramics contained about 87 wt.% HA, 13 wt.% TCP, and a negligible amount of periclase (Sheikh et al., 2015; Tamimi et al., 2012).

Other parameters obtained from the refinement and included in the model were the lattice parameters, the crystallite size, and the microstrain.

ICP-MS analysis was applied to quantify the proportion of the different elements (Ca, P and Mg) in the ceramic particles upon complete dissolution. Moreover, ICP-MS analysis of the cell culture medium was performed upon immersion of the ceramic particles for 2, 4, 7, and 14 days in the cell culture medium. The quantification of the Ca and P content in the particles (Figure 2a, b) and the resulting Ca/P and (Ca + Mg)/P atomic ratio (Figure 2c, d) showed that the materials were not completely stoichiometric. Our results for the (Ca + Mg)/P ratio are in accordance with those obtained by Maazouz et al. in a recent study in which the same materials were analyzed (Maazouz et al., 2021). The ICP-MS results confirmed the presence of Mg in the CaPs in similar amounts as suggested by the XRD analysis (Figure 2e).

Figure 2. ICP-MS analysis on CaP ceramics and culture medium. a–e) Ca, P, and Mg ion content of the six CaP ceramics. f–h) Ca, P, and Mg ion content of the cell culture medium. The error bars represent the standard deviation of 3 independent measurements.

The elemental analysis of the cell culture medium upon ceramic immersion (Figure 2f,g,h) showed significant differences between the control, i.e., cell culture medium not containing a material, and the different materials, as well as among the different ceramics. Specifically, within the category of the ceramics with similar chemistry, those sintered at lower temperature had lower Ca levels in the culture media than the chemically similar material sintered at higher temperatures and even lower than in the control well (Figure 2f). A similar trend was observed for P, although in this case, the effect was less profound (Figure 2g). It was observed that the uptake of Ca and P was higher in the materials sintered at lower temperatures than in chemically similar materials sintered at higher temperatures. This uptake of Ca and P was probably due to a higher porosity/smaller grain size obtained at lower sintering temperatures, which translated into a higher specific surface area.

Various CaP ceramics have been reported to decrease the Ca content of the cell culture medium, which itself is saturated towards HA (Birgani et al., 2016; Othman et al., 2019). The observation was hypothesized to be either due to the precipitation of a new CaP layer on the surface of the ceramic, or due to the formation of calcium carbonate or calcium hydroxide (Maazouz et al., 2020; 2021).

Regarding the Mg content of the cell culture medium, the results suggested a burst release of Mg from HA1150, with a concentration at day 2 being about twice as high as that of the control and other ceramics (Figure 2h). HA1150 was the ceramic with the highest percentage of periclase among the six tested materials. No significant differences between the control and the other five ceramics were observed in terms of Mg content of the medium at any time point.

3.1.2. Structural characterization

Besides chemical composition, structural properties, including surface topography and properties of porous structure, have been shown to affect the biological response to CaP ceramics. First, SEM imaging was used to obtain qualitative information about the CaP surface topography and structure, including pores, grains, and other features. The SEM micrographs (Figure 3a) exhibited differences in surface microstructure among the different ceramics. The materials sintered at lower temperatures (left column) contained more micropores and were comprised of smaller grain, while the materials sintered at higher temperatures (right column) had scarce microporosity and larger grain size.

Figure 3. Structural analysis of CaP ceramics. a) High-resolution SEM micrographs illustrating the difference between osteoinductive (left column) and non-osteoinductive ceramics (right column). b) Pore size distribution obtained using mercury intrusion porosimetry analysis. Osteoinductive materials demonstrate to contain pores of multiple sizes. Pores were not detected for the non-osteoinductive materials except for TCP1100, which contains pores of a larger size. c) Grain size quantification demostrates low values for osteoinductive and high values for non-osteoinductive ceramics. d) Microporosity across the ceramics. e–j) Selected profilometric parameters demonstrate a larger diversity across the ceramics. k–m) Selected microCT parameters allowing a description of the 3D structures across the ceramics. The error bars represent the standard deviation of 3 independent measurements.

We confirmed these visual observations through mercury intrusion porosimetry, which allowed quantifying the pore size distribution (Figure 3b). The grain size quantification of the materials further confirmed the observations and showed that BCP1300 had the largest grain size among the six ceramics, with the lowest values found for HA1150, TCP1050, and BCP1150 (Figure 3c). The pore size distribution curves (Figure 3b) obtained from mercury intrusion porosimetry were also distinct for the different ceramics. It is important to emphasize that the mercury intrusion porosimetry quantification only included pores with sizes below 15 µm, and therefore the total porosity of the ceramics is plausibly higher than the microporosity values (Figure 3d). Taken together, these results showed different degrees of porosity among different ceramics, with HA1250 almost devoid of pores and BCP1300 having very low porosity. It is also important to note that the difference in porosity between TCP1050 and TCP1100, although statistically significant, was smaller than was observed for the HA and BCP pairs. Importantly, the ceramics with lower porosity are the non-osteoinductive ones, and porosity has been reported to play an essential role in the osteoinductive potential of CaP ceramics (Bohner & Miron, 2019).

Next, we performed confocal profilometry to quantify the ceramic surface properties. The 3D height maps again exhibited differences in microstructural features depending predominantly on the sintering temperature (Supplementary Figure 3). The root mean square height (Sq) (Figure 3e), a parameter frequently used to quantify surface roughness (Semnani, 2017; Webb et al., 2012), was higher for the non-osteoinductive ceramics, expect for the BCP ceramic, with Sq being highest for BCP1150. This is an interesting observation considering that roughness is generally considered to improve an osteogenic phenotype (Faia-Torres et al., 2014). A similar trend was observed for the peak density parameter (Spd), a parameter that represents the number of peak features per mm² (Figure 3f). The arithmetic mean curvature (Spc) is a feature that represents the arithmetic mean of the principal curvature of the peaks of the surface, i.e., a low Spc indicates that the points of contact of the surface have rounded shapes, while higher values mean that they have pointed shapes. Spc values also followed a similar trend as the Sq parameter (Figure 3g). Also, the developed interfacial area ratio (Sdr) was significant different between TCP1150 and BCP1300, and BCP1150 and BCP1300, with large yet non-significant higher levels detected for HA1250 and TCP1100 (Figure 3j). The upper material ratio (Smr1), which physically represents the areal material ratio that divides the areas that would be removed by initial abrasion with another surface from the core surface of the material (Supplementary Figure 4), was only significantly different between HA1150 and the two TCP ceramics (Figure 3h). The lower material ratio (Smr2), which represents the percentage of the measurement area that contains valley structures able to hold liquid on the surface and improve lubricity (Supplementary Figure 4), was significantly different between BCP1150 and the two TCP ceramics (Figure 3i).

MicroCT images also provided information about the surface topography and structure. In addition, since microCT is a 3D imaging technique, it allows the quantification of parameters such as the pore interconnectivity and the number of closed pores. Due to the resolution of the acquisition, 2.25³ µm³ voxel size, the microCT results are complementary to those from the mercury porosimetry, as they predominantly cover larger pores and 3D macrostructure. The surface-to-volume ratio of the non-osteoinductive ceramics was higher than their osteoinductive counterparts (Figure 3k). While this result may be considered counterintuitive, as it is to be expected that materials with higher porosity have a higher surface-to-volume ratio, this is probably due to a relatively low resolution of the technique, not taking into account the pores smaller than 2.25³ µm³. Of interest, the trabecular bone pattern factor (TBPF), a measure of bone microarchitecture, also allowed distinguishing osteoinductive and non-osteoinductive ceramics and was highest for HA1150 and lowest for BCP1300 (Figure 3l). The percentage of closed porosity, i.e., pores without connection to the exterior of the ceramic particles, was relatively low for all ceramics (Figure 3m). It is noticeable that the least osteoinductive material of the set, TCP1100, exhibited the highest percentage of closed porosity. This is relevant because the presence of surface pores has been recently suggested to be a key parameter for CaP-induced osteoinduction in vivo (Bohner & Miron, 2019). To conclude, these methods provided a detailed material properties characterization of each ceramic.

3.2. Correlations among material properties

After establishing an extensive dataset describing different properties of the six ceramics, we first aimed to investigate which of the material parameters are correlated. For this, Spearman correlations analysis was performed and plotted in a heatmap (Figure 4).

Figure 4. Heat map showing Spearman correlations among material properties of the CaP ceramics. Positive correlations between material parameters are marked in blue, while negative correlations are marked in red.

The Mg substitution percentage showed correlations with seven surface topographical parameters, namely the arithmetic mean peak curvature, the developed interfacial area ratio, the root mean square gradient, the density of peaks, the reduced peak and valley height, and the peak material volume. This is an interesting finding, as it is well known that adding bioinorganics to the CaP lattice can significantly affect the ceramic microstructure (Birgani et al., 2016; Farzadi, Bakhshi, Solati-Hashjin, Asadi-Eydivand, & Osman, 2014; Tahmasebi Birgani et al., 2016). Specifically, Mg substitution has been shown to change the size and morphology of CaP nanoparticles (Batra, Kapoor, & Sharma, 2013; Nabiyouni, Ren, & Bhaduri, 2015) and the crystallinity of amorphous CaP powders (Mehrjoo et al., 2015). Here, Mg substitution percentage showed a strong negative correlation with Spc and Sdr (r = −0.95), indicating that a higher Mg substitution was associated with more pointed shapes and a smaller surface area. The specific correlation graphs are depicted in Supplementary Figure 5A. Of interest, the addition of Mg to CaPs can also lead to the formation of periclase, which, based on our analysis, correlated with the Mg release from the ceramic, negatively at day 2 and positively at day 14. Furthermore, it was observed that the fraction of periclase positively correlated with the degree of anisotropy of the material, suggesting that the presence of periclase leads to a more anisotropic material, perhaps due to an anisotropic spatial distribution of the phases.

We also observed a strong correlation between the levels of Ca and P in the medium (Supplementary Figure 5B), which was already evident in Figure 2f, g, showing a similar trend for the content of both elements, depending on the ceramic type tested.

The grain size was shown to be inversely correlated with the trabecular pattern factor, which describes the concavity and convexity of the material, and is related to the connectivity of the material (Supplementary Figure 5C) (Hahn, Vogel, Pompesius-Kempa, & Delling, 1992).

A number of other correlations was observed between porosity, profilometry, and microCT parameters. It is plausible that due to the definitions of these parameters and their scale of quantification (i.e. submicrometer-scale in porosimetry and profilometry, micrometer and milimiter-scale in microCT), they actually may describe the same structural features.

3.3. PCA and hierarchical cluster analysis of material properties

We further wanted to establish if the extensive material characterization is sufficient for distinguishing the ceramic materials from each other. For this, we employed a principal component analysis (PCA). To reduce the dimensionality of the data and improve the accuracy of the algorithm, we removed material parameters that represented technical variations and those that showed a strong correlation (>0.99) with other parameters. The final dataset included 76 material properties. Of interest, the PCA plot shows grouping of the materials, with on the one side the osteoinductive ceramics HA1150, BCP1100, and TCP1050, and on the other side the non-osteoinductive ceramics BCP1300, HA1250, and TCP1100 (Figure 5a). Here, the non-osteoinductive ceramics had more distinct properties than the osteoinductive ceramics as demonstrated by their larger distance between them. The data demonstrates that the physico-chemical set of properties we investigated allows separating ceramics that are osteoinductive and non-osteoinductive.

Figure 5. Characterization of material parameters allows identifying common characteristics in the ceramic materials. a) PCA plot of the six CaP ceramics based on the quantified materials parameters. The osteoinductive materials BCP1150, HA1150, and TCP1050 cluster together on the left side of the plot. The non-osteoinductive ceramics appear to be more distinct between each other. b) Hierarchical cluster analysis of the material properties illustrates common and distinct material characteristics across the materials. Cluster four, containing parameters of grain size and Ca/P/Mg content in the medium, allows for distinguishing osteoinductive and ­non-osteoinductive ceramics.

Next, we wanted to investigate which properties are responsible for this separation by employing a hierarchical cluster analysis (Figure 5b). Through this visualization, it becomes clear that cluster 4 allows distinguishing the osteoinductive and non-osteoinductive ceramics. This cluster contains the grain size, bone surface/volume ratio, and the Ca/P/Mg concentration in the medium, confirming the observations shown in Figures 2 and 3. Of interest, also in cluster 2 a separation between the osteoinductive and non-osteoinductive ceramics is possible except for BCP1300, which followed a similar pattern as BCP1150, HA1150, and TCP1050 concerning the surface roughness parameters.

Among the three osteoinductive materials, HA1150 appears to be the most distinct regarding its physico-chemical properties, as seen from the PCA plot. The difference between BCP1150 and TCP1050 is derived from both the structural and chemical parameters. HA1150 possesses a lower porosity (30%) than TCP1050 (44.6%) and BCP1150 (48.4%) as shown in Figure 3 and cluster 3. The release of Mg was also higher for HA1150, as shown in Figure 3, especially on day 2, and contained more periclase (cluster 7). Looking at its performance in vivo, HA 1150 induced less de novo bone formation upon implantation in the femoral muscle of dogs, and the onset of osteoinduction occurred later relative to BCP1150 (Yuan et al., 2006). These results further support the proposition that the properties of the ceramics indeed impact their osteoinductive potential.

BCP1150 and BCP1300, despite having a similar chemical composition (Figure 2b), appeared on the opposite sides of the PCA plot. The properties that make them different are both structural and chemical concerning the ion exchange with the culture media (cluster 4). Their porosity and microporosity differ strongly, with over a 10-fold difference, as is shown in Figure 4d-e and cluster 3, and is noticeable in the SEM images (Figure 4a). Moreover, the two ceramics differed in their grain sizes (Figure 4b) and bulk densities. In an earlier in vivo study, BCP1150 had superior bone formation properties compared to BCP1050 (Habibovic et al., 2006), which the authors suggested to be due to the microstructure difference between the two ceramics, possibly coupled to the ions exchange dynamics differences on the surface.

3.4. Correlations between material properties and gene expression

We further aimed to correlate material properties with gene expression levels of MG63 cells cultured on these ceramics. This will allow us to investigate which genes or pathways associated with osteogenesis are activated due to individual material parameters. To demonstrate that the ceramics induce unique gene expression profiles, we performed a hierarchical cluster analysis of the 500 most variable genes between the ceramics (Figure 6a). Four main clusters allowed separating the gene expression profiles across the ceramics. The largest differences were observed between BCP1150 and HA1250. Of interest, TCP1050 and TCP1100, despite being different in their properties and differing significantly in their osteoinductive potential, appeared close in the clustering analysis, suggesting that their influence is reflected by the differential expression of a relatively small number of genes. The two HAs, HA1150 and HA1250, are close in the cluster hierarchy, although in this case, the number of genes with similar expression levels are smaller than TCP1050 and TCP1100. BCP1150 and BCP1300 induced the largest difference in gene expression of the ceramics with similar chemical composition.

Figure 6. CaP ceramics induce distinct gene expression profiles in MG63 cells and allow a correlation analysis with material properties. a) Hierarchical cluster analysis of the top 500 genes with strongest gene expression variation across the six CaP ceramics. Ceramics with identical chemical compositions cluster together except for the BCP ceramics. b) Heatmap showing the frequency of significant correlations between CaP ceramic properties and pathways associated with osteogenesis.

Next, we correlated the material properties with gene expression of each ceramic through Spearman correlations. The genes with a p-value lower than 0.05 and a correlation value higher than 0.7 or lower than −0.7 were grouped into eight categories based on Gene Ontology (GO) terms associated with osteogenic signaling pathways (Vermeulen, Tahmasebi Birgani, & Habibovic, 2022). Such an approach was previously employed to elucidate associations between surface topographical structures and gene expression levels in hMSCs (Vasilevich et al., 2020). The GO categories included the cytoskeleton, skeletal development and osteogenesis, matrix organization and development, calcium signaling, phosphate signaling, the mitogen-activated protein kinase (MAPK) pathway, the BMP pathway, and the WNT pathway. Next, we plotted in a heatmap the number of genes from each category correlated to an individual material property of which the variability between the ceramics was sufficient to allow correlations with gene expression levels (Figure 6b).

The cytoskeleton-related genes had numerous correlations with structural parameters identified using the microCT analysis. These parameters include the degree of anisotropy of the material, the fractal dimension, the trabecular thickness, and the trabecular bone pattern factor. The degree of anisotropy accounts for the existence of preferential alignment along a spatial axis in the material and is negatively correlated with the expression of genes related to the cytoskeleton (Supplementary Figure 6a). This suggests that cells on isotropic surfaces have higher cytoskeleton-related gene expression levels compared to anisotropic surfaces. Previous research with micro-topographical surfaces that induced cell morphologies with a high level of anisotropy similarly associated this observation with a decreased expression of cytoskeleton-related genes (Vermeulen et al., 2020). Also the fractal dimension, a quantification of an object’s surface complexity through surface features that are repeated over different spatial scales, had negative correlation coefficients with cytoskeleton gene expression (Supplementary Figure 6b).

The correlations with the trabecular pattern factor and the trabecular thickness (Supplementary Figure 6c, d) indicated that cytoskeleton-related gene expression was enhanced on structures with thicker and more interconnected trabeculae. In addition, the ceramic grain size appeared to be correlated to the cytoskeleton-related genes (Supplementary Figure 6e), although both positive and negative correlations were observed. Previous literature has shown that smaller CaP grain sizes led to a better attachment of fibroblasts and osteoblasts to the materials surfaces (Dasgupta, Tarafder, Bandyopadhyay, & Bose, 2013; Veljović et al., 2012). The porosity parameters also showed numerous correlations with cytoskeleton genes (Supplementary Figure 6f–i), which is not surprising since porosity, pore size, and pore distribution influence cell attachment and morphology (Deligianni et al., 2000; Guadarrama Bello, Fouillen, Badia, & Nanci, 2017). Several roughness parameters, the root mean square height, the developed interfacial area ratio, the autocorrelation length, and the kurtosis, exhibited correlations with the cytoskeleton genes (Supplementary Figure 6k–n), which was expected as roughness can affect cytoskeleton formation (Anselme & Bigerelle, 2011; Faia-Torres et al., 2014; Linez-Bataillon, Monchau, Bigerelle, & Hildebrand, 2002). Finally, Ca, P, and Mg ion ceramic content and levels in the cell culture medium also influenced cytoskeleton formation, a relationship that has been widely reported (Bennett & Weeds, 1986; Murakami et al., 2010; Rosado & Sage, 2000).

The genes related to skeletal development and osteogenesis showed correlations with both chemical and structural parameters (Supplementary Figure 7). The composition of the ceramic and the inorganic content of the cell culture medium have been reported to influence osteogenesis processes (Habibovic & Barralet, 2011; Landi et al., 2008; Otsuka et al., 2008). Related to this is also the microstrain, influenced by the Mg substitution, which determines the stability of the crystalline phase, and therefore indirectly influences the inorganic content of the medium. Besides these chemical parameters, correlations for genes related to osteogenesis were found for microstructure parameters such as root mean square height (Sq) and pore size distribution parameters, which was also reported before (Dalby et al., 2007; Wilkinson et al., 2011).

The set of genes related to matrix organization and development was predominantly correlated with structural parameters (Supplementary Figure 8). To the best of our knowledge, the effects of trabecular thickness, fractal dimension connections, and connectivity on the expression of matrix-related genes have not been studied directly yet. However, it is plausible that these structural dimensions change cellular morphology, influencing the expression of matrix-related genes (Le et al., 2017). The correlation with Sq was expected as surface roughness is known to influence matrix development (Cai et al., 2020).

The Ca signaling gene set had a small number of correlations with the Ca levels in the cell culture medium (Supplementary Figure 9). In contrast, the Ca signaling was correlated mainly with structural parameters. As discussed, this may be an indirect effect, as (surface) structural properties largely determine the available surface area for interaction with the biological environment, which in turn affect the ions exchange dynamics with the environment, including the exchange of Ca. Furthermore, since cell morphology can influence Ca signaling (Ribeiro, Reece, & Putney, 1997), also here, a connection can be made with structural parameters that influence cell shape.

Only a few genes were found to correlate with the phosphate signaling, mainly affected again by the pore size (Supplementary Figure 10). For example, both Mg levels and pore connectivity showed correlations with the xenotropic and polytropic virus receptor 1 (XPR1), which encodes a protein involved in phosphate homeostasis by mediating the phosphate export from cells (Chande & Bergwitz, 2018; Mailer et al., 2021).

The MAPK pathway affects genes related to a variety of cellular processes, such as proliferation, differentiation, stress response and apoptosis, which can be modulated by the materials properties. It is further related to the phosphorylation and activation of RUNX2 through various signaling pathways such as collagen-integrin binding or BMP signaling (Franceschi & Xiao, 2003). According to our results, MAPK showed abundant correlations with structural parameters (Supplementary Figure 11a–i). Similarly, correlations with Ca, P, and Mg levels in the medium were observed (Supplementary Figure 11j–s), which confirms the findings from previous studies (Barradas et al., 2012).

The BMP pathway, an important contributor to osteogenic processes (Wang, Wang, et al., 2019), was mainly influenced by structural parameters (Supplementary Figure 12). Previous studies where a chemistry-topography decoupling strategy was used, showed a mild effect of topography on BMP-2 expression (Danoux, Sun, et al., 2016; Sun et al., 2016).

Finally, we observed correlations between WNT pathway genes and both chemical and topographical parameters (Supplementary Figure 13). Examples include the WNT5b ligand, a crucial component of the WNT pathway (Perkins, Suthon, Miranda-Carboni, & Krum, 2021), with a positive correlation existing with Smr1 and the P levels in the cell culture medium. Next to WNT5b, we also found a positive correlation between the Mg content of the cell culture media and receptor tyrosine kinases Ryk, WNT1-Inducible-Signaling Pathway Protein 1 (WISP1), and Zinc Finger RANBP2-Type Containing 1 (ZRANB1), all proteins related to WNT signaling (Green, Nusse, & van Amerongen, 2014; Tran, Hamada, Schwarz-Romond, & Bienz, 2008; Xu, Corcoran, Welsh, Pennica, & Levine, 2000). Why this correlation exist needs further investigation. Considering the importance of WNT signaling for osteogenesis, it is not surprising that CaP ceramics can influence this pathway (Wang, Pan, et al., 2019).

The bioinformatics method used here proved capable of identifying numerous correlations between material properties and gene pathways, most of which are supported by the literature. This allows to select a certain material property and analyze its influence at a single gene level. For example, several strong correlations were observed between the microporosity and a range of genes belonging to various categories as described before (Table 2). Previous research associated microporosity with cell adhesion (Hoai & Nga, 2018), osteogenic differentiation, and in vivo bone formation (Duan et al., 2017). In addition, the presence of micropores has been related to changes in cytoskeleton (Hoai & Nga, 2018), ECM production (Lien, Ko, & Huang, 2009), and migration towards/through the pores (Hoai & Nga, 2018; Li et al., 2015; Mandal & Kundu, 2009; Pamula, Filova, Bacakova, Lisa, & Adamczyk, 2009). Finally, cell proliferation also seems to be affected by the microporosity of the substrate (Hoai & Nga, 2018; Hornez et al., 2007; Li et al., 2015).

Table 2. Genes correlated with the microporosity parameter with known relevant functions.

Gene/protein	p-value	Functions	References
MAPRE1/EB1	<0.001	Growth and stabilization of microtubules	Bu and Su (2001)
		Cell migration	Nehlig, Molina, Rodrigues-Ferreira, Honore, and Nahmias (2017)
TRIB2	0.004	Cell fate regulation	Fang et al. (2021)
MSRB2	<0.001	Actin filaments polymerization	Bai et al. (2020)
PSMB8	<0.001	Cell differentiation	Arimochi, Sasaki, Kitamura, and Yasutomo (2016)
TIMP1	<0.001	ECM remodeling	Prihadi, Sugandi, Siregar, Soejono, and Harahap (2018)
		MSCs proliferation and osteogenic differentiation	Liang et al. (2019)
YMLE1L	<0.001	Cell proliferation	Stiburek et al. (2012)
STAT1	0.005	Osteoclastogenesis and bone density	Place et al. (2021)
		Osteoblastogenesis	Mizuno et al. (2020)
		Bone remodeling	He et al. (2020)
		MSCs osteogenic differentiation	Xiao et al. (2019)
TRIM32	0.005	Cell migration	Kano, Miyajima, Fukuda, and Hatakeyama (2008)
VAMP3	<0.001	Cell migration	Riggs et al. (2012)
		Cell adhesion	Luftman, Hasan, Day, Hardee, and Hu (2009)
		ECM remodeling	Kean et al. (2009)

This study also showed some unexpected results, which require further investigation. For example, based on the previous literature, one would expect stronger correlations between the Ca signaling pathway, the Ca content of the ceramics and the cell culture medium (Barradas et al., 2012). Similarly, no correlations were found between phosphate signaling and the P content of the medium. This might be because we only looked at a 48-hour time point in the transriptomics study and may have missed some processes, such as a strong upregulation of BMP-2 at earlier time points. Similarly, we obtained correlations for STAT1, which is known to upregulate RUNX2, but not for RUNX2 itself, which may be expressed at a later time point. Furthermore, some other parts of the experimental setup need to be considered for future studies, such as the fact that the cell culture medium used for the analysis of the inorganic content was collected in the absence of cells, which may be different from the setting where cells are present.

Taken together, the main achievement of this study is that we have shown the potential of the bioinformatics approach to explore correlations between materials properties and transcriptomics results, thereby taking into consideration many more parameters than would be possible without such a tool. The results pointed to some interesting correlations, many of which were supported by the available literature. Nevertheless, a follow-up validation phase is essential to verify the results concerning the correlations found among materials properties and between material properties and transcriptomics. In the future, we also foresee that other ceramic or composite designs, and cell types, such as primary osteoblasts or osteoclasts, will be included in our bioinformatics approach to further elucidate material parameters relevant for bone regeneration.

4. Conclusions

In this study, we used a bioinformatics approach to explore the correlations between properties of six CaP ceramics with known osteoinductive and bone regenerative potential and gene expression of MG63 cells cultured on these materials. A comprehensive material characterization and data analysis were performed to obtain a large set of numerical parameters describing both the chemistry and the microstructure of the ceramics. The results showed that the ceramics could be divided into an osteoinductive and a non-osteoinductive group based on their material parameters. A Spearman analysis of the material properties revealed a number of correlations between chemical and structural parameters, confirming that the properties of CaPs are largely intertwined. The Pearson correlation analysis of the materials properties dataset with the transcriptomics gene expression dataset identified a large number of correlations, most of which were supported by the literature. The approach here proves the value of detailed (bio)materials characterization in combination with statistical analysis to identify the genes and gene pathways regulated by the individual material properties.

Acknowledgments

The authors acknowledge Dr. Arno Wilbers (DSM Resolve) for the help with the acquisition and reconstruction of the microCT datasets, and Dr. Nicola Döbelin (RMS Foundation) for the assistance help with the Rietveld analysis.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The raw/processed data to reproduce these findings is available in the following repository: https://doi.org/10.34894/VWMT6T

Additional information

Funding

This research has been made possible with the support of the Dutch Province of Limburg (program “Limburg INvesteert in haar Kenniseconomie/LINK”; grant nos. SAS-2014-00837 and SAS-2018-02477). SV acknowledges the financial support by the European Union Interreg Vlaanderen-Nederland project “BIOMAT on microfluidic chip” (grant no. 0433). PH gratefully acknowledges the Gravitation Program of the Netherlands Organisation for Scientific Research (NWO) (project Materials-Driven Regeneration; grant no. 024.003.013) and the NWO Vidi grant “Bone Microfactory” (grant no. 15604).