Full factorial design of experiment-based and response surface methodology approach for evaluating variation in uniaxial compressive mechanical properties, and biocompatibility of photocurable PEGDMA-based scaffolds

The goal of this study is to fabricate biocompatible and minimally invasive bone tissue engineering scaffolds that allow in situ photocuring and further investigate the effect on the mechanical properties of the scaffold due to the prevailing conditions around defect sites, such as the shift in pH from the physiological environment and swelling due to accumulation of fluids during inflammation. A novel approach of incorporating a general full factorial design of experiment (DOE) model to study the effect of the local environment of the tissue defect on the mechanical properties of these injectable and photocurable scaffolds has been formulated. Moreover, the cross-interaction between factors, such as pH and immersion time, was studied as an effect on the response variable. This study encompasses the fabrication and uniaxial mechanical testing of polyethylene glycol dimethacrylate (PEGDMA) scaffolds for injectable tissue engineering applications, along with the loss in weight of the scaffolds over 72 h in a varying pH environment that mimics in vivo conditions around a defect. The DOE model was constructed with three factors: the combination of PEGDMA and nano-hydroxyapatite referred to as biopolymer blend, the pH of the buffer solution used for immersing the scaffolds, and the immersion time of the scaffolds in the buffer solution. The response variables recorded were compressive modulus, compressive strength, and the weight loss of the scaffolds over 72 h of immersion in phosphate-buffered saline at respective pH. The statistical model analysis provided adequate information in explaining a strong interaction of the factors on the response variables. Further, it revealed a significant cross-interaction between the factors. The factors such as the biopolymer blend and pH of the buffer solution significantly affected the response variables, compressive modulus and strength. At the same time, the immersion time had a strong effect on the loss in weight from the scaffolds over 72 h of soaking in the buffer solution. The biocompatibility study done using a set of fluorescent dyes for these tissue scaffolds highlighted an enhancement in the pre-osteoblasts (OB-6) cell attachment over time up to day 14. The representative fluorescent images revealed an increase in cell attachment activity over time. This study has opened a new horizon in optimizing the factors represented in the DOE model for tunable PEGDMA-based injectable scaffold systems with enhanced bioactivity.


Introduction
In the field of orthopaedics, a large portion of research has been focused on developing improved bone repair techniques and interventions. While mature bone can naturally heal minor defects under normal circumstances, large traumatic, congenital, or surgical bone loss may lead to complications in the healing process [1]. As a result, bone tissue scaffolds have emerged as a popular potential technique to repair bone fractures and other osseous defects. These scaffolds are threedimensional (3D) structures that mimic the extracellular matrix and are inserted at bone defect sites as a template for bone regrowth and healing. In addition, they can promote bone healing by providing mechanical strength to the healing site, introducing osteoblasts for repair, and mimicking the natural tissue environment before slowly degrading and paving the way for tissue regeneration [2,3]. Thus, suitable tissue engineering scaffolds must possess biodegradability, biocompatibility, and adequate mechanical strength [2,4]. To this end, multiple types of scaffolds have been investigated, ranging in material composition and fabrication technique [3,[5][6][7][8].
Hydrogels are well explained as a 3D network of hydrophilic polymers that can encompass a large volume of liquid, viz. water or other physiological solutions, and display tunable properties governed by their intrinsic configurational and compositional features [9][10][11][12]. The hydrogel matrix's structural stability and other mechanical properties depend highly on the nature of crosslinking that holds the polymer network together. The polymer network conformed by the covalent linkages between polymer chains has distinct differences in terms of properties and applications as compared to polymer matrices with intertwined and entangled chains [13,14]. Thus, a crosslinking mechanism that harnesses an external light source is beneficial during in situ curing of bioactive scaffolds and drug delivery systems [15].
Among the plethora of options in hydrogels that have demonstrated good structural stability and bioactivity, polyethylene glycol (PEG) based photocurable hydrogel systems have gained immense interest, primarily due to scaffold stability, biocompatibility, and relatively more straightforward synthesis methods [16][17][18]. A suitable photoinitiator and the desired wavelength of light energy can initiate crosslinking via free-radical polymerization in pristine or composite PEG-based photocurable polymers [19][20][21][22]. Various functionalized groups have been studied with PEG, including diacrylate and dimethacrylate linkage formulated PEG gaining importance in photocurable PEG-based polymers [20,[23][24][25][26][27][28][29]. Properties of the polyethylene glycol dimethacrylate (PEGDMA) gel matrix may be adjusted by varying factors such as monomer concentration and molecular weight, as well as crosslinking density governed mainly by the photoinitiator amount and light energy applied [18,20].
The movement of polymer chains within the photocrosslinked hydrogel matrix largely depends on the liquid molecules encapsulated in the polymer matrix. With the inherent stress of the polymer chains due to their configuration and other network defects, the ability of the matrix to expand to a specific limit is varied [11,19,30,31]. Nevertheless, the hydrogel expansion or swelling degree scenario governs several vital parameters for cell activities, especially mobility, and proliferation [32,33]. Moreover, there is a requirement for adequate mechanical strength of the tissue engineering scaffold for adequate physiological acceptability and cell activity in situ [34][35][36]. The compressive modulus of the scaffold provides the metric of stiffness that significantly impacts the subsequent cell activities of attachment and proliferation [37][38][39][40]. Therefore, the scaffold's compressive modulus and subsequent swelling degree are quintessential parameters in predicting scaffold compatibility in tissue engineering applications. Due to different configurational changes occurring within the polymer network because of a change in the pH of the extracellular environment, the PEGDMA scaffold may be experiencing varying crosslinking changes within the network, which can be highlighted as a change in the mechanical property, such as compressive modulus and strength [41,42].
This study has been sculpted around the variation of the compressive mechanical properties of PEG-DMA scaffolds concerning the change in pH of the surrounding environment. The full factorial design of experiment (DOE) model has been constructed to analyze the parameters crucial to changing the compressive modulus/strength of the prepared photocrosslinked scaffolds. ASTM standard cylindrical scaffolds were used for compressive mechanical testing at different time points and pH environments. Lithium phenyl-2,4,6-trimethylbenzoylphosphinate (LAP) was used as the photoinitiator, which performs the crosslinking of the PEGDMA matrix in the blue light region, thus, avoiding the use of direct UV in an in situ environment.
Moreover, bioceramic additives may be used with the polymer material during scaffold fabrication to improve the strength, cell attachment, and osteoconductivity of synthetic polymers such as PEGDMA [43,44]. The most common of these bioceramics are calcium phosphate (CaP)-based due to the matrix of natural bone having an abundance of CaP-like minerals [45]. Hydroxyapatite (HA) in the form of CaP is found in bone and can be used as a bioceramic additive in polymeric tissue scaffolds [43,[46][47][48]. The similarities between these ceramics and the bone matrix make them an ideal material to incorporate into bone tissue scaffolds. HA alone can achieve a compressive strength near that of the natural bone [48]. In addition, CaP-based bioceramics have increased polymeric scaffold strength [43,44,46,48]. Therefore, it is possible that combining bioceramics and synthetic polymers may keep the beneficial properties of each material in the composite material while helping improve the overall strength. Thus, nano-hydroxyapatite (nHA) was incorporated in the biopolymer blend with PEGDMA and evaluated using the DOE model alongside the pristine PEDGMA scaffolds. Moreover, the regression equation generation for the DOE model was done using a response surface methodology (RSM), which is a statistical technique involving scrutiny and optimization of multiple factors within an experiment model to select the best-fit model and predict outcomes based on given parameters [49,50].

Materials
PEGDMA, with a molecular weight of 400, was purchased from Polysciences, Inc. (USA). LAP was procured from Allevi, Inc. (USA). nHA was bought from Sigma Aldrich (MO, USA). Invitrogen LIVE/DEAD Viability/Cytotoxicity Kit for mammalian cells was purchased from ThermoFisher Scientific (USA). Gibco, by Life Technologies, alpha Minimum Essential Medium (α-MEM), Fetal Bovine Serum (FBS), Phosphate-buffered Saline (PBS, pH 7.4), Dulbecco's PBS, and penicillin/streptomycin were also purchased from ThermoFisher Scientific (USA). All solutions/dilutions in water were prepared using Ultrapure (Type I) water obtained from Millipore Sigma ultrapure water system.

Fabrication of scaffolds
The scaffolds were fabricated by mixing PEGDMA with 0.5% (w/v) LAP dissolved in water in a 1:1 volume ratio. The nanofiller nHA, at a concentration of 2% (w/v), was dissolved in the biopolymer blend containing LAP and then vortexed at high speed until homogenous throughout the solution. The prepared polymer blend was injected using a 1 ml syringe into ASTM standard cylindrical molds 12 mm in length and 6 mm in diameter. Photocuring the scaffold was done using blue light for 5 min and then promptly removed from the mold for further mechanical testing.

Compressive mechanical testing
All compressive mechanical tests were carried out using an ADMET eXpert 2600 universal testing machine with an Interface S-type load cell. The scaffold samples were subjected to uniaxial compression testing at a constant loading rate of 0.005 mm s −1 .
Six (6) replicates were used from each sample/factor group for testing, with n = 108 for each response variable during the entire DOE model processing. In addition, force vs. change in length data were collected and analyzed using standard mathematical equations and approaches [51] to obtain compressive mechanical property parameters, such as modulus and strength.

DOE structure
The full factorial DOE model comprised of three (3) factors with two (2) levels for Factor 1, and 3 levels for Factors 2 and 3 each, respectively, as shown in table 1. Factor 1 was divided into two levels, pristine PEG-DMA vs. PEGDMA mixed with nHA, to study the effect of the nanofiller incorporation on the mechanical properties of the scaffold. It has also been reported that during defect healing, the bone tissue defect repair process expresses a varying pH, from acidic to alkaline [52][53][54][55]. Therefore, factor 2 was designed to incorporate three levels for pH -5, 7.4, and 9, which takes into account the effect of an acidic, normal physiological, and alkaline pH environment on the scaffold's mechanical properties. Since immersion time in a specific medium/buffer influences the water uptake and the changes in the scaffold's mechanical properties with time, factor 3 was constructed to contain three levels for the immersion time of the scaffolds in the respective medium. The immersion media was prepared using 1× PBS and subsequently adjusted for the pH per the requirement. All scaffolds were immersed in the medium at 37 • C for the desired time and then immediately tested under the universal testing machine.

In vitro cell culture and attachment
Two-dimensional mammalian cell culture of murine preosteoblast (OB-6) cell line was done using α-MEM culture medium containing 15% FBS and 1% penicillin-streptomycin. We used a murine preosteoblast (OB-6) cell line that we received from Dr Beata Lecka-Czernik at the University of Toledo for our cell culture experiments. All cell cultures and respective in vitro experiments were maintained at an ambient of 37 • C with 95% air/5% CO 2 , and the culture medium was replaced every three days. After reaching 80% cell confluence, the cells were digested and subcultured using 0.25% Trypsin-Ethylenediaminetetraacetic Acid (EDTA) to obtain a cell suspension for seeding onto the scaffolds. The PEGDMA scaffolds were soaked in sterile 1× PBS (pH 7.4) for 12 h, followed by 30 min UV sterilization, and finally soaked in a culture medium for another 12 h under sterile conditions. Cells were seeded at a density of 40 000 cells per well onto the scaffold samples directly using a 24-well plate.
Cell viability and attachment analysis were done according to the manufacturer's LIVE/DEAD assay kit protocol. Three time points on days 3, 7, and 14 were chosen for the in vitro experiment to evaluate the biocompatibility of the prepared scaffolds. The fluorescence microscopy images depicted the cell attachment scenario using the green fluorescence of calcein for live cells obtained by the enzymatic changeover of the cell-permeable acetomethoxy derivative of Calcein. Concurrently, dead cells were imaged using the red fluorescence emitted due to the nucleic acidbinding of the ethidium homodimer-1. The fluorescence microscope and cell imaging system used was Cytation 5 (by BioTek Inc., USA).

Statistical analysis
DOE analysis was done using Origin Pro 2020 software. All parametric testing was conducted after performing the Shapiro-Wilk normality test at a 95% confidence interval. Data transformation was carried out using Minitab version 21.
1. An alpha value of 0.05 corresponding to a 95% confidence interval was selected throughout this study, with a p-value <0.05 signifying statistical significance unless explicitly mentioned as a different p-value criterion.
Among the various plots generated during the analysis of the DOE model, the main effects plot displays the mean response values at each level of a design variable or parameter. This plot compares the relative strength of the effects of various factors on the response variable [56]. As such, a horizontal line in the main effects plot signifies no prominent effect on the mean response variable. In contrast, a deflection from the horizontal in the main effect plot displays an effect on the response by the factor or parameter. Furthermore, a steeper slope in this line highlights the increasing effect of the factor on the response [57]. Thus, when the mean effects plot has an increasing slope in the deflection from the horizontal, it infers that the mean response is higher at the factor-level setting that corresponds to the increasing slope of the main effects plot. Hence, in contrast, the downward slope of the main effects plot would signify a higher value of the mean response at the factor-level setting that corresponds to the decreasing slope of the plot.
Subsequently, the fitted plot displays the scatter plot of the experimental measurement along with the predicted value as per the regression model used in progression with the experimental observation order, also known as row order. One of the fitted plot's critical information is identifying outliers in the statistical model being studied. Along with the fitted plot, residuals are an essential metric for implying the strength of the predictive model. The residual plots are generated using the error terms, i.e. the difference between the measured value and the predicted value in progression with the row order [56]. A random distribution of the residuals vs. row order scatter plot implies a random distribution of the error term, thereby inferring independent error terms with no significant correlation. The residuals vs. fitted values plot detects unequal error variances and nonlinearity. A robust regression model would display a random distribution of the residuals with progression in the fitted values on the x-axis. Any specific trends in the data may infer a weak or biased statistical regression model, requiring further data transformation [58].
Moreover, the histogram of the residuals infers the normality of the data from the statistical model to test the assumption of data normality for validating the statistical analysis generated from the DOE model. Further, the normal probability plot verifies the normality of the residuals where the plotted data would lie close along with the fitted probability straight line for normally distributed residuals with minimal deviations. Moreover, two-way interaction plots are presented to infer the intensity of interaction between factors for a response variable. When the plotted lines for the factors cross each other, it infers a strong interaction between the two factors in that plot, while a parallel line signifies weak or no interaction between the factors.
Further, a full factorial RSM model was generated on Design Expert 13 using the independent variables used in the previous DOE model to generate other three-dimensional surface plots and regression equations for the compressive modulus and strength of the fabricated scaffold.

Variation in compressive modulus
The effect of the three factors on the mean of the uniaxial compressive modulus is shown in figure 1. It was highlighted that adding nHA to the PEGDMA blend has a strong uphill effect on the compressive modulus. As seen in figure 1(A), the enhancement in the mean compressive modulus is strongly affected by adding nHA. It can be explained as the additional reinforcement attained by the hydrogel network walls [59]. Also, during the preliminary study, as shown in supplementary A1, adding 10% (w/v) nHA showed a significantly reduced compressive modulus. This result may be highlighted by the interference of the nHA to the path of the blue light during photocuring. Hence, 2% (w/v) was chosen for the DOE study. nHA incorporation also adds to the well-known benefits during osteogenesis [60].
The effect of change in pH on the rigidity of the material exhibits a downward trend, as shown in figure 1(B), where more acidic pH environments transform the polymer network to render more stiffness. This may be attributed to the mineral deposition triggered in the nHA at low pH [61], thus, increasing the matrix stiffness of the polymer [62,63]. The trend in immersion time as an effect on the modulus, figure 1(C), has been reported as expected as the polymer swelling increases with time, reducing the material stiffness.
The scatter plot of the regression model's predicted values and the compressive modulus's experimental values, as shown in figure 2, shows a good fit with certain outliers. Such outliers may have been due to the experimental variations contributed by scaffold swelling and microstructural variations due to curing.
The behavior of the residuals for the regression model is shown in figure 3. As shown in figure 3(A), the residuals plotted as per the sequence they have been generated show a random distribution of the residuals with no recognizable pattern. This scenario confirms the assumption of independence of the residuals from one another for the compressive modulus testing regime. Moreover, figure 3(B) shows that the residuals have a constant variation through the fitted values except towards the higher end of the predicted value scale. This converging pattern of the residuals may be attributed to the reducing variation between observed vs. expected values at larger values of the predicted response. Also, the normality of the data is depicted in figure 3(C) with a considerably smooth bell curve for the data distribution in the histogram, further confirmed by figure 3(D) with a relatively straight and uniform distribution of the residuals along the fit line.
The two-way interaction plot, as shown in figure 4, depicts the significance of the interaction between the factors. The factor interaction subplot

Variation in compressive strength
The addition of nHA had a strong positive effect on the compressive strength, as shown in figure 5(A), like the trend followed by compressive modulus vs. biopolymer blend. The addition of the filler enhances the reinforcement of the polymer network, which ultimately enhances the strength. However, the change in pH from acidic to more alkaline had a positive effect, highlighted in figure 5(B), especially a steep rise from pH 7.4 to 9. This may be attributed to the enhancement of ductility at a high pH that may have occurred due to the reduction in hydrogen bonds within the polymer network. As can be seen in figure 5(C), the effect of the immersion time on the mechanical strength has a comparatively mild effect. This scenario may be explained due to the significantly minimal effect of swelling of the polymer network on the change in ductility of the polymer matrix.
The scatter plot of the regression model's predicted values and the compressive strength's experimental values, as shown in figure 6, is a random distribution with certain outliers. Such outliers can again be attributed to the microstructural variations during the curing of the scaffolds, along with swelling differences. The residual plots, as shown in figure 7, show the random distribution of the data points without any recognizable pattern. Based on the original data collected, the regression model points to strong twoway interactions between several factors, as shown in figure 8. The interaction between biopolymer blend, pH, and immersion time has significant interactions, especially for the levels pH 5 and 7.4.
However, the normality test of the data was reported as skewed, as shown in figure 7(D). The nonnormal distribution of the data can cause accuracy errors in the p-value and the overall regression model. Therefore, a Johnson transformation operation was performed on the original data to obtain normality, as shown in figure 9. The p-value for the transformation was around 0.3, which rejected the null hypothesis at a confidence interval of 95%, thereby confirming the normality of the transformed data. These transformed data were further analyzed with the regression model to obtain the confirmation of randomness in residual distribution, as shown in figures 10(A) and (B). The histogram skewness has also been improved with the data transformation, as shown in figure 11. The main effects plot for the transformed data, shown in figure 12, replicates the trend of interactions of the factors with   the average compressive strength recorded. The twoway factor interaction plot of the transformed data, shown in supplementary data figure A2, was identical in the inference of significant interactions between the factors-biopolymer blend, pH, and immersion time, as compared to figure 8.

Variation in weight loss
The loss in weight is affected strongly by adding the nHA, as shown in figure 13(A), with more average weight loss reported with nHA incorporation. This may be attributed to the nHA adhered to the polymer network walls being leached off to the surrounding   environment. As shown in figure 13(B), the change in pH had a milder effect on weight loss. However, the weight loss of the scaffold was seen to be affected strongly by the increase in immersion time. This suggests that the increase in swelling over time adds to enhanced loss in the weight of the polymer matrix.
The scatter plot of the predicted values for the weight loss and the experimental values, as shown in figure 14, substantially shows a good fit with one outlier. An outlier may be attributed to measurement variations along with scaffold fabrication and curing variations arising from a random assignment.
The behavior of the residuals for the regression model is shown in figure 15. As shown in figure 15(A), the residuals show a random distribution with no recognizable patterns. Thus, it highlights the assumption of independence of the residuals from one another to be valid. Moreover, figure 15(B) shows that the residuals have a constant variation through the fitted values except for some outliers. Also, the normality of the data is depicted in figures 15(C) and (D) reports the normal distribution of the weight loss data as recorded.
The two way interaction plot, as shown in figure 16, shows the significance of the interaction between the factors. The loss in weight is strongly affected by the significant interaction of pH with the biopolymer blend and immersion time, especially for the pH 7.4 and 9, along with the immersion time. The addition of nHA has a strong effect on weight loss, while the pH of the immersion medium changes from 7.4 to 9. The immersion time has a positive linear trend with increasing weight loss with time. This may be noted as an expected scenario in weight loss due to increased swelling with time at all levels of the factor pH.
Moreover, keeping in mind the timeline for bone defect healing, the mean weight loss was measured for 45 days, n = 3, and extrapolation was done for 120 and 300 days. The PEGDMA-only scaffold degradation experiment showed 8.1% ± 0.7% degradation on day 45. Upon linear extrapolation considering day 0 and day 45 degradation data, at day 120 and day 300, the predicted degradation is 21.5 ± 2.1% and 53.9 ± 4.6%, respectively. The PEGDMA + nHA scaffold degradation experiment showed 7.30% ± 0.89% degradation on day 45. Upon linear extrapolation     considering day 0 and day 45 degradation data, at day 120 and day 300, the predicted degradation is 18.9 ± 2.5% and 48.7 ± 5.7%, respectively.

Regression equation for compressive modulus and strength using RSM
The 3D surface plot generated from the RSM model for optimization of the compressive modulus of the fabricated scaffolds is shown in figure 17.
The variation in the compressive modulus of the PEGDMA-based scaffolds in response to the independent variables can be mathematically represented by equations (1) and (2) The 3D surface plot obtained from the RSM model for the compressive strength is shown in figure 18.
The variation in the compressive strength of the PEGDMA-based scaffolds in response to the independent variables can be mathematically represented by equations (3) and (4) Compressive Strength = 3.83 + (0.004 * A)

Biocompatibility of scaffold
The biocompatibility study of the PEGDMA-based scaffolds was done using fluorescence imaging on OB-6 seeded scaffolds, as shown in figure 19. The day 3 results showed comparatively more cell attachment for the nHA incorporated scaffolds as compared to the PEGDMA-only scaffold group. This may be attributed to the nanofiller, nHA, that creates comparatively more surface roughness for better cell attachment. Also, nHA acts as a biological and chemical cue for the pre-osteoblasts already triggered further down the osteogenesis pathway. Similarly, the day 7 results also followed a similar trend as day 3; however, on day 14, the scenario of cell attachment was similar among both groups. This may be due to the subsequent proliferation of attached cells. Moreover, the biocompatibility of the PEGDMA-based scaffolds reports a promising scope for cell activities with implanted scaffolds.

Limitations and prospects
This study highlighted the appropriate information for evaluating the intensity of interaction of the factors in the DOE model with the variation of uniaxial compressive mechanical properties. However, the maximum mechanical modulus and strength may be a tunable criterion based on the requirement of the application. The prospects in optimizing the compressive mechanical properties of the PEGDMAbased injectable scaffold may further see studies in the incorporation of viscosity-enhancing biopolymers, such as gelatin, which also has a reputation as a bioactivity enhancer for composite scaffold materials. Also, the combination of biopolymer and nanofillers may be further optimized to tune the loss in weight from the scaffolds over time. Further, a broad scope opens up in conducting similar DOE-based porosity and bioactivity studies with both chemical and biological cues and as factors affecting the mechanical behavior of the photocured scaffolds.

Conclusion
The DOE-based approach for evaluation of the significantly affecting factors for PEGDMA mechanical property variation and weight loss yielded crucial information, especially in terms of the modulus vs. strength. The increasing crystalline nature of the polymer network may explain the increase in matrix stiffness at low pH. This creates more rigidity in the polymer network that ultimately enhances the modulus. On the contrary, at a higher pH, the ductility and hence strength has been reported to have been enhanced. The inclusion of the nanofiller, nHA, has a positive effect on the polymer matrix's mechanical properties, including weight loss. This loss in weight of the scaffolds with nanofiller incorporation may be explained as the leaching out of the weakly adhered nanofillers on the walls of the matrix. Weight loss as a response variable of the scaffold also increases significantly with the increase in swelling degree when the immersion time increases. Although the compressive strength reported a milder aftereffect of the increase in immersion time, the material's rigidity was distorted with enhancement in polymer swelling. This study creates a conductive environment for zoning in on the factors that strongly affect the mechanical property of photocurable PEGDMA-based hydrogel matrix. Moreover, the increase in rigidity at an acidic pH, in a way, mimics the probable scenario of defect healing with a comparatively lower pH environment surrounding it. Although more alkaline pH is rare around a tissue defect site, in extraordinary circumstances of high pH environments, the scaffold modulus may drop; however, it would primarily depend on the exposure time to the alkaline environment. The RSM study generated the regression equations for predicting the compressive modulus and strength of the PEGDMA-based scaffolds at different conditions of the independent variables. Future studies may include adding more nanofillers components to balance the mechanical property variation across the pH range and obtain a viable polymer formulation that can resist drastic mechanical property variation in any environment. The biocompatibility of the minimally invasive and photocurable scaffolds was promising, with enhancement in cell activity profile over 14 days. The bioactivity study confirmed the biocompatibility of the active biopolymer blend components, especially for the photoinitiator LAP. In addition, the photocuring method and time were confirmed to have produced an osteoconductive substrate for cell activity.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).