Prerpints.org logo
Preprint
Article

Material Design of Porous Hydroxyapatite Ceramics via Inverse Analysis of an Estimation Model for Bone-Forming Ability Based on Machine Learning and Experimental Validation of Biological Hard Tissue Responses

Altmetrics

Downloads

66

Views

29

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

26 December 2023

Posted:

27 December 2023

You are already at the latest version

Alerts
Abstract
Hydroxyapatite and b-tricalcium phosphate have been clinically applied as artificial bone materials due to their high biocompatibility. The development of artificial bones requires the verification of safety and efficacy through animal experiments; however, from the viewpoint of animal welfare, there is necessary to reduce the number of animal experiments. In this study, we utilized machine learning to construct a model that estimates the bone-forming ability of bioceramics from material fabrication conditions, material properties, and in vivo experimental conditions. We succeeded in constructing the two models as follows: ‘Model 1’, which predicts material properties from their fabrication conditions, and ‘Model 2’, which predicts the bone-formation rate from material properties and in vivo experimental conditions. By an inverse analysis of the two models, we proposed candidates for material fabrication conditions to achieve target values of bone-formation rate. Under the proposed conditions, the material properties of the fabricated material were consistent with the estimated material properties. Furthermore, a comparison between bone-formation rates after 12 weeks of implantation in the porcine tibia and the estimated bone-formation rate. This result showed that the actual bone-formation rates existed within the error range of the estimated bone-formation rates, indicating that machine learning consistently predicts the results of animal experiments using material fabrication conditions. We believe that these findings will lead to the establishment of alternative animal experiments to replace animal experiments in the development of artificial bones.
Keywords: 
Subject: Chemistry and Materials Science  -   Biomaterials

1. Introduction

As the global population is rapidly aging, an increase in healthy life expectancy is required [1]. In order to prolong healthy life expectancy, the maintenance and recovery of locomotor functions are crucial, for which orthopedic bone regeneration is effective [2]. Bone grafting is necessary for the treatment of bone defects caused by tumors and various diseases. Autologous bone grafting is the golden standard [3]; however, numerous complications and secondary invasions occur after extraction from the donor site [4,5]. Therefore, artificial bone grafts have become more common, and artificial bones based on bioceramics have been studied.
Hydroxyapatite (Ca10(PO4)6(OH)2; HAp) is one of bioactive ceramics and b-tricalcium phosphate (b-Ca3(PO4)2; b-TCP) is one of bioresorbable ceramics, which have been clinically applied as an artificial bone material [6,7,8,9]. These porous bioceramics allow biological bone to penetrate the pores of implanted material and integrate with biological bone tissue. Controlling microstructure, such as pore size, porosity, and pore volume, can affect the bone-forming ability [10,11,12].
We have successfully synthesized and succeeded in the single-crystal apatite fiber (AF) and calcium phosphate fiber (CPF) with various Ca/P molar ratios by homogeneous precipitation using urea [13,14,15,16,17]. Furthermore, apatite-fiber scaffolds (hereafter, AFS) and b-tricalcium phosphate-fiber scaffolds (TFS) with macro- and microporous structures have been successfully developed by adding carbon beads (CB) as a pore-forming agent and sintering [17,18,19]. We have also succeeded in developing porous apatite ceramics that include bone minerals (BoneHAp) by adding six types of ions to HAp [20,21]. During the development of these artificial bone materials, it is necessary to verify the efficacy and safety of the materials through animal experiments. However, these processes are costly, time-consuming and require animal sacrifice [22].
The use of machine learning to predict the properties and activity (objective variable Y) of compounds and materials is common. Next, using the experimental and simulation results, a regression model (Y = f(X)) is constructed between Y and the explanatory variables X, which represent the experimental conditions. This model makes it possible to predict the Y values under new experimental conditions without experimentation, thereby contributing to the development of more efficient and cost-effective materials. Research is being conducted to predict the physical properties of materials such as ceramics [23], polymers [24], and metals [25], and to propose experimental conditions for materials with properties that reach the target properties, thus facilitating the development of materials. Although there have been studies using machine learning to predict biomaterial structure and properties [26,27] and predict bone formation in animal experiments with materials [28], it has been impossible to predict the in vivo responses in animal experiments from bioceramic material synthesis conditions using machine learning. Although we have previously constructed a model that predicts bone-formation rates using the results of animal experiments on developed bioceramics and proposed material fabrication conditions based on desired bone-formation rates [29], we have not performed experiments based on the proposed experimental conditions, not validate the validity and practicality of the models.
In this study, we constructed a model to predict the bone-forming ability of bioceramics and designed the material fabrication conditions by inverse analysis of the constructed model, in which various candidate for the material fabrication conditions and in vivo experimental conditions were input into the model, and the material properties and bone-formation rate were predicted. Additionally, porous HAp ceramics designed by inverse analysis were fabricated. We verified whether the properties of the material and the bone-forming ability of the material implanted in the porcine tibia were consistent with the estimated values of the models.

2. Materials and Methods

2.1. Dataset

The dataset consists of data from material preparation to in vivo experiments on biomaterials developed in our laboratory, including AFS, TFS, and porous bone HAp ceramics. Figure 1 shows a schematic of the model that predicts material properties and bone-formation rate from material fabrication and in vivo experimental conditions. In order to predict bone-formation rate, models were constructed using material fabrication conditions, material properties, and in vivo experimental conditions. First, Model 1 (Y1 = f(X1)) was constructed to predict the material properties (Y1) from the material fabrication conditions (X1). Next, Model 2 (Y2 = f(Y1, X2)) was constructed to predict the bone-formation rate based on material properties (Y1) and in vivo experimental conditions (X2). Material fabrication conditions, material properties, and in vivo experimental conditions are listed in Table 1. The bone-formation rate (Y2) was calculated from histological images of bone tissue removed after implantation of the porcine tibia or rat calvaria using the formula in Figure 2. The data were stained with toluidine blue (TB) or villanueva bone (VB) staining. TB staining is a staining method that stains new bone blue, while VB staining is a method that differentiates the calcification stage of the bone. From stained histological images, the area of ceramics remaining after transplantation and the area of bone formation were determined, and the bone-formation rate was calculated by subtracting the area of ceramics remaining from the area of bone loss created when the specimen was transplanted.
Autoscaling is performed as preprocessing of the dataset for each variable to transform the mean to 0 and the SD to 1.

2.2. Construction of Models and Verification of Predictability

Since there are multiple objective variables in Model 1, it is desirable to construct a model using a machine learning method that consider the relationship between the objective variables, y, and the explanatory variables, x. Therefore, Models 1 and 2 were constructed with Gaussian mixture regression (GMR) [30], which is a nonlinear machine learning model. GMR is a method to perform regression analysis based on the Gaussian mixture model (GMM) [31], which is a model that represents a dataset with a mixture of Gaussian distributions. z connects x and y. In the GMM, the probability distribution function of the multivariate Gaussian distribution p(z) is given as follows:
p z = ∑ i = 1 n π i N   ( z ❘ μ i ,   ∑ i ) ,                   ∑ i n π i = 1
where n is the number of Gaussians, mi and Si are the mean vector and variance-covariance matrix in the ith Gaussian. By explicitly separating x and y in Equtaion 1 in the GMM, the joint probability distribution of x and y can be expressed as:
p ( x ,   y ) = ∑ i = 1 n π i N   x       y ❘ [ μ x ,     i       ∑ y , i ] , ∑ x x , i ∑ y x , i ∑ x y , i ∑ y y , i
where mx,i, my,i, Sxx,i, Syy,i, Sxy,i, and Syx,i, are the mean vectors of x and y, variance–covariance matrices of x and y and the covariance matrices of x and y in the ith Gaussian, respectively. Estimating y from x corresponds to calculating the posterior distribution of y given x (p(y|x). p(y|x) can be transformed by the multiplication theorem of probability and Bayes’ theorem as follows:
p y | x = ∑ i = 1 n p   ( y | x , μ x , i ,   ∑ x x , i ) p ( μ x , i ,   ∑ x x , i   | x ) = ∑ i = 1 n p   ( y | x , μ x , i ,   ∑ x x , i ) p ( x | μ x , i ,       ∑ x x , i ) p ( μ x , i ,           ∑ x x , i ) ∑ j = 1 n p ( x | μ x , j ,       ∑ x x , j ) p ( μ x , j ,           ∑ x x , j )      = ∑ i = 1 n p   ( y | x , μ x , i ,   ∑ x x , i ) π i p ( x | μ x , i ,       ∑ x x , i ) ∑ j = 1 n π j p ( x | μ x , j ,       ∑ x x , j )
where p(y|x, mx,i, Sxx,i) is the multivariate Gaussian distribution of the predicted y values in the ith Gaussian. wx,i, which is the weight of the Gaussian, is given by
w = π i p ( x | μ x , i ,       ∑ x x , i ) ∑ j = 1 n π j p ( x | μ x , j ,       ∑ x x , j )
Substituting Equtaion 4 into Equtaion 3 gives
p y | x = ∑ i = 1 n w x , i p   ( y | x , μ x , i ,   ∑ x x , i )
In p(y|x, mx,i, Sxx,i), the mean vector mi(x) is given by
m i ( x ) = m   + ( x   -   m )   ∑ x x , i − 1 ∑ x y , i  
In this study, the predicted y values are those with the highest weights. GMR can perform regression analysis and inverse analysis without combining other machine learning methods, even with multiple y, while taking their relationships into account. Additionally, since the GMM is constructed for the entire dataset, by changing the ranges of x and y, regression and inverse analyses that consider the constraints can be performed. Therefore, GMR is suitable for this study and was used to construct models.
Due to the small sample size of 50 in the dataset, double cross-validation (DCV) [32] was used to evaluate the predictive ability of models. DCV is a method to evaluate the predictability of a model with a nested structure of cross-validation (CV). The dataset is divided into two groups. One group is used for verification, while the other is used for hyperparameter optimization by CV and model construction with the obtained hyperparameters. Then, the model is then used to estimate the data for verification, and the estimated y values are obtained. This is conducted until all the groups are used for verification, and the estimated y values are compared with the measured y values to evaluate the predictability.
This study used a dataset consisting of 50 samples, and 49 samples to construct the model and predict the other sample. Furthermore, this step was repeated 50 times so that all samples were test data not involved in model construction. Therefore, the model can be tested using all 50 samples.
The coefficient of determination (r2DCV) and the root mean squared error (RMSEDCV) are used as evaluation indicators after DCV. r2DCV and RMSEDCV are given by:
r 2 D C V = 1 − ∑ i = 1 n ( y ( i ) − y D C V i   ) 2 ∑ i = 1 n ( y ( i ) − y m e a n   ) 2
R M S E D C V = ∑ i = 1 n ( y ( i ) − y D C V i   ) 2 n
where n is the number of samples, y(i) is the actual value of ith y, y(i)DCV is the estimated value of ith y, and ymean is the average value of y. Higher r2DCV and lower RMSEDCV indicate higher predictive accuracy.

2.3. Feature Importance

We calculated the permutation feature importance (PFI), which represents how much each explanatory variable (material properties, in vivo experimental conditions) of Model 2, which predicts bone-formation rate, affected the accuracy of the model. However, PFI calculation requires training and validation data, and a small number of samples makes the calculation of feature importance unstable. Additionally, the PFI of strongly correlated features is estimated to be lower than that of other independent features. Therefore, to solve these problems and properly calculate the importance of features, we used cross-validated permutation feature importance (CVPFI) [33].
Assuming that the number of iterations is J, the algorithm to calculate for the CVPFI can be described as follows.
  • Calculate the correlation coefficients for all the features.
  • Calculate the absolute correlation coefficient for all the feature and set the coefficient to zero when there is no correlation.
  • During CV, for n = 1, 2, ⋯, N, where N is the number of CV folds, the following procedures were performed using the training and validation data (VD) at each fold.
3-1
Construct a model using the nth training data.
3-2
Estimate y for the nth VD, that is VDn, using the model.
3-3
For each feature I, which is the ith column of VDn, and for each repetition j in 1, 2,⋯, J, a randomly sampled column, I, of the original data set without duplication, and for each feature, m (the mth column of VDn), for which the absolute correlation coefficient between features I and m, that is, ri,m, after statistical testing is higher than 0, use a randomly sampled column, m, of the original data set without duplication with a probability of ri,m to generate a corrupted version of the data set CVDn,I,j and estimate the y values of the model.
  • Integrate the y values estimated during the CV for VD1, VD2, ⋯, and VDN and calculate the reference score, rscv, with the integrated y. This score represents the determination coefficient, r2, for a regressor.
  • Integrate the y estimated during the CV for CVD1,I,j, CVD2,I,j, ⋯, and CVDN,I,j and calculate the score, scvi,j, with the actualy.
  • Calculate the importance, CVPFIi, for the ith feature as follows:
    CVPFIi,j = rscv – scvi,j
    C V P F I i = 1 J ∑ i = 1 J C V P F I i , j
As CVPFI integrates the estimated y values and since all samples can be used to calculate r2 for a regressor, the importance can be calculated stably using CVPFI.
Therefore, we calculate the CVPFI using the GMR used in the construction of Model 2 to interpret the precision of each explanatory variable in predicting bone-formation rates.

2.4. Materials Designs for New Bioceramics by Inverse Analysis of the Model

New bioceramics were designed using the constructed models. Virtual samples were generated under the following conditions: The AFS and TFS fabrication conditions (e.g., heating conditions and firing conditions) were set to known conditions. 924 samples were generated by combining 3 types of molding pressure (30–50 MPa in a increments of 10 MPa), 11 types of 150 mm CB addition ratios (0–100% in a increments of amounts of 10%), and 14 types of CB addition (50–700 mass% in a increments of amounts of 50 mass%). The virtual samples were input into Model 1 and the material properties were calculated. The calculated material properties and three in vivo experimental conditions (porcine tibia 12 weeks; rat cranial crown 2, 4 weeks) were then combined to generate 2,772 samples. The virtual samples generated were input into Model 2, and the bone-formation rate was calculated. Five of the virtual samples were selected and fabricated under these conditions.

2.5. Fabrication of Materials

AFs were synthesized according to an earlier report [13,15]. The following aqueous solutions were prepared: 0.167 mol･dm-3 Ca(NO3)2･4H2O, 0.100 mol･dm-3 (NH4)2HPO4, 0.500 mol･dm-3 (NH2)2CO, and 0.100 mol･dm-3 HNO3. The Ca/P ratio was adjusted to 1.67. The AFs were mixed with a spherical CB of ~150 mm and ~20 mm diameter in a mixed solvent (ethanol/water = 1:1 [v/v]). Based on the results of the inverse analysis, the amount of CB added and the ratio of CB with different particle sizes were determined, and uniaxial pressure molding was performed at 30 MPa. AFS were prepared by baking the compacts in a steam atmosphere at 1,300°C, with a temperature increase rate of 5°C･min-1, for 5 h. AFSs with dimensions of 3.8 mm (diameter) and 7 mm (height) were used in this study. The samples were abbreviated as AFS X(Y-Z), where X is the amount of CB added (mass%), Y is the ratio of 150 mm CB added (%), and Z is the ratio of 20 mm CB added (%). Five types of AFS were fabricated.

2.6. Material Characterization and Validation

The AFSs were characterized as follows. Phase composition was determined by X-ray diffractometry (XRD; MiniFlex, Rigaku Co., Japan). The X-ray generator with CuKa was operated at 30 kV and 15 mA. Data were collected in the range of 2θ 4.0-50.0° with a step size of 0.04° and a count speed of 4 ° ･min-1. We achieved phase identification by comparing the diffraction patterns with the International Centre for Diffraction Data standard. Fourier transform infrared (FT-IR; IR Prestige21, Shimadzu Co., Japan) spectra were measured using the KBr tablet method. Spectra were obtained at a 4 cm-1 resolution over the range of 400-4,000 cm-1 with a scan speed of 4 mm·s-1. The microstructures were observed by scanning electron microscopy (SEM; JSM6390LA, JEOL Ltd., Japan). Compressive strength was measured using a universal test machine (Autograph AGS-J, Shimadzu Co., Japan) at a crosshead speed of 0.5 mm･min-1. Porosity was determined as follows. Firstly, the bulk density was determined from the dimensions and mass of the fabricated ceramics specimens. This value was divided by the theoretical density (for intense, 3.16 g/cm3 in the case of HAp) and multiplied by 100 to calculate the relative density (%). Porosity (%) was determined by subtracting the relative density (%) from 100%. The actual porosity and compressive strength were compared with the estimated values of the model.

2.7. In Vivo Evaluation of Materials Using a Pig Tibia Model and Validation

We used a miniature pig weighing 83.8 kg to evaluate the bone-forming ability of AFSs in vivo. Subsequently, the right tibia of the pig was exposed and a 3.8 mm diameter cylindrical defect was opened in the tibial epiphysis. A dry, heat-sterilized AFS was placed into the defect. The pig had access to a standard laboratory diet and water throughout the study. Surgical procedures were performed according to the Animal Care and Use Committee of Meiji University.
After implantation for 12 weeks, the pig was sacrificed, and the tibiae were removed. Noncalcified polished sections were prepared and observed under an optical microscope after VB staining. Bone-formation rates were measured using image analysis software (WinROOF, Mitani Co., Japan). Subsequently, we verified whether the measured bone-formation rates were in agreement with the estimated values of the model.

3. Results and Discussion

3.1. Construction of Model 1

The Porosity, compressive strength, and four peaks at full width at half maximum (FWHM) of the XRD pattern were predicted from the material fabrication conditions. FWHM was calculated from the XRD data as the (2 1 1) and (3 0 0) peaks for HAp and the (2 2 0) and (0 2 10) peaks for b-TCP. The plots of the actual material property versus the predicted material property in DCV are shown in Figure 3. The calculated values of the r2dcv and RMSEdcv are shown in Table 2. DCV results show that approximately 90% of porosity, 70% of the compressive strength, and 90%–99% of FWHM were predicted for the new data. Outliers were identified in the scatter plots for compressive strength. Furthermore, in compressive strength tests, there tends to be a large variation in the measured values even under the same experimental conditions. In this study, the model was constructed using average values without considering the variation. This suggests that some samples were off the diagonal due to variations caused by individual differences in materials. These results show that GMR can predict multiple material properties simultaneously with high predictability.

3.2. Construction of Model 2

The bone-formation rates were predicted on material properties and in vivo experimental conditions. We compared the predictability of Model 2 with and without FWHM. The plots of the actual bone-formation rate versus the predicted bone-formation rate in DCV are shown in Figure 4. The calculated values of the r2dcv and RMSEdcv are also shown in Table 3. Without FWHM, the predictability was approximately 50%. With the addition of the FWHM, we found that the model could predict approximately 70% of the new data. This result suggests that FWHM is an important factor in predicting bone-formation rates. When FWHM is an explanatory variable, and information on the crystal structure of the material can be included in the model. Since the crystalline structure of the material affects its bone-forming ability, the addition of FWHM may have improved accuracy. As a result, Model 2 was to predict complex in vivo reactions with 70% predictability.

3.3. Feature Importance

In order to determine how much each feature of Model 2, which predicts bone-formation rate, affects the accuracy, the feature importance was calculated. The feature importance was calculated as the CVPFI using GMR. Model 2 was constructed using material properties, in vivo experimental conditions, vascular endothelial growth factor (VEGF), and FWHM 4 peaks for bone formation rate, and CVPFI values were calculated. The calculated feature importance is shown in Figure 5. The FWHM 4 peak was found to be the most important. The four FWHMs were equally important. FWHM’s high importance suggests that the identification of the crystalline phase and the slight differences in crystallinity affect the model. Next, using the half width, the crystallite size can be calculated using Scherrer’s equation [34], where a large crystallite size indicates a small FWHM. However, when the FWHM is small, Scherrer’s equation cannot be applied to determine the crystalline size due to increase in the error. The FWHM of the material included in these data is also small and Scherrer’s equation could not be applied. Therefore, by using half-width instead of crystallite size, the crystallinity of the material could be expressed as the explanatory variable X. The ratio of crystallite sizes of (2 1 1) to (3 0 0) for HAp or (0 2 10) to (2 2 0) for b-TCP can represent the anisotropy of the particles. Furthermore, the equal importance of the four FWHMs suggests that the relationship between these two HAp and b-TCP peaks contributed to the model. It is also suggested that the model could include the crystalline morphology of the material. When the material properties and the in vivo experimental conditions, the material properties were found to be more important. This shows that the material properties contribute significantly to the model. In other words, material properties have a greater influence on bone-forming ability than in vivo experimental conditions (the week of implantation and implantation animal). Porosity is a property that describes the microstructure, and has been reported in the past to affect bone-forming ability, which is consistent with its importance [11,18,19]. Compressive strength is a property that depend on porosity. It can be inferred that the importance of compressive strength increases with the importance of porosity.

3.4. Materials Designs for New Bioceramics by Inverse Analysis of Model

A total of 2,772 virtual samples were generated under various constraints. Material properties were calculated by inputting the generated samples into Model 1. Here, FWHM was not used and the crystalline phase was added to the variables as a dummy variable (0: HAp, 1: b-TCP) to classify the crystalline phase. Specifically, AFS was set to 0 (HAp) and TFS to 1 (b-TCP) and calculated as the crystalline phase, and the calculated material properties and in vivo experimental conditions were input into Model 2 to calculate bone-formation rate. The applicability domain (AD) [35] was set to investigate the reliability of the virtual sample. The distance between the existing and designed material fabrication conditions was calculated using the k-nearest neighbor (k-NN) method [36], a method of averaging the distance to new data. k-NN calculations used Euclidean distance, and k was set to 1 due to the small number of training samples. The distance between the 2,772 generated samples and the closest known distance was calculated, and the 99.7% value was used as the threshold for AD. Of the 2,772 samples generated, 2,733 were identified as samples within AD. Five AFSs were selected among the samples within AD. The experimental conditions, calculated material properties, and calculated bone-formation rates for the selected samples are presented in Table 4.

3.5. Experimental Validation of Material Properties

In this section, we describe the results of AFSs prepared of the experimental conditions presented in the inverse analysis and compare them with the properties of their actual and predicted materials. Material XRD patterns showed that all AFSs were single-phase HAp. The FT-IR spectra of AFS showed absorptions of PO43- groups at 1300-900, 600, and 570 cm-1 and absorptions of OH- groups at 3,570 cm-1. SEM images of AFSs that contain many macropores because of the burning of CBs and many interconnected micropores derived from entangled AFs.
Figure 6(a) shows the comparison of the actual and estimated porosities of AFS. The error bar of the estimated porosity is the SD of the error between the actual and predicted porosities for each sample in DCV (SD; ±3.09). If the actual value exists within the error range of the estimated value, the actual value is defined as consistent with the estimated value. The porosity increased with increasing amount of CB added. The porosity also increased with the addition of 150 mm CB. This trend was similar to that of the estimated porosities of the model. The actual porosities of all AFSs were within the error range of the estimated porosities and the values were generally in agreement. In the DCV, the predictability for approximately 90% of the unknown data suggests that the actual and estimated values were mostly in agreement.
Figure 6(b) shows the comparison of the actual compressive strength and the estimated compressive strength of the AFSs. The error bar of the estimated compressive strength is the SD of the error between the actual and the estimated compressive strengths for each sample in DCV (SD; ±1.32). The compressive strength decreased as the amount of added CB increased. Additionally, the compressive strength also increased with the ratio of CB with a smaller particle size. This trend was the same as the estimated compressive strength of the model, and the actual compressive strength of all AFSs was within the error range of the estimated compressive strength. However, for the two AFSs with a high compressive strength, there was a large error between the actual and estimated values. The actual values of the material with a high compressive strength exhibited significant variation due to individual differences in the materials and measurement errors. Additionally, assuming the same trend in the dataset, individual differences in materials and measurement errors could be the cause of the prediction errors of Model 1.
A comparison of porosity and compressive strength showed relatively good agreement between the actual and estimated values. The dataset used in this study contains a large amount of data on porous bioceramics with different conditions for the addition of pore-forming agents. In the inverse analysis, virtual samples were generated by varying the pore-forming agent addition conditions and the molding pressure. This indicates that the virtual samples were generated close to the samples included in the dataset. In fact, 98% of the virtual samples generated were within AD. It restricted inverse analysis based on the dataset result in good agreement between the estimated and actual values. These results suggest that porous ceramics with the desired strength and porosity can be designed and fabricated.

3.6. In Vivo Evaluation of Materials and Validation

Among the five AFSs fabricated, three kinds of AFSs were used in animal experiments: 1) AFS100(100-0), 2) AFS100(70-30), and 3) AFS350(50-50). Histological images of VB staining of pig tibia in visible light and fluorescence are shown in Figure 7. It was observed that all AFSs were directly bound to the newly-formed bone of the host. Additionally, good bone formation was observed around and inside the material, which confirmed its excellent osteogenic potential. Fluorescence imagery also showed that two AFSs with a porosity of approximately 70% formed a large amount of calcified bone.
Figure 8 shows a comparison of the bone-formation rates calculated by image analysis with estimated bone-formation rates. The error bar of the estimated bone-formation rate is the SD of the error between the actual and the estimated bone-formation rate for each sample in DCV (SD; ±12.93). The AFS350(50-50) bone-formation rate was the highest among the porous AFS specimens. The trend between actual and estimated bone-formation rates was mostly consistent. In particular, the actual bone-formation rate of AFS100(100-0) was almost in agreement with the estimated bone-formation rate. Actual bone-formation rates in all AFSs were within the error range of estimated bone-formation rates. This error could be attributed to errors in predictability and individual differences among experimental animals i.e., the results showed that the predictions were correct within the error range due to individual animal differences.
In contrast, there were no differences in the estimated bone-formation rate between AFS100(100-0) and AFS100(70-30). Only porosity was included as an explanatory variable for the microstructure. It has been reported that not only porosity, but also pore diameter and pore volume affect bone-forming ability [10,12,37], and it is possible that porosity alone does not explain the bone-formation rate. In the future, we aim to further improve the predictability by including more detailed data on the microstructure, such as pore size distribution, and by extracting features through image analysis using SEM images and adding them as explanatory variables.

4. Conclusion

We used machine learning to construct Model 1 (Y1 = f(X1)) to predict material properties (Y1) from material fabrication conditions (X1) and Model 2 (Y2 = f(Y1, X2)) to predict bone-formation rate (Y2) from material properties (Y1) and in vivo experimental conditions (X2). The predictability of Model 2 was improved by adding FWHM to the explanatory variables. Based on the feature importance, the FWHM was the most important in predicting the bone-formation rate, as far as we were able to determine in this work. This suggests that the crystalline phase and crystallinity affect the bone-forming ability. Through an inverse analysis of the constructed models, we proposed the fabrication conditions for new porous hydroxyapatite ceramics and then fabricated the materials under the conditions. The actual material properties were relatively close to the estimated material properties. The created materials were implanted into the porcine tibia and the bone-formation rate was calculated. The actual bone-formation rates existed within the error range of the estimated bone-formation rates. We believe that this research will enable us to propose artificial bone materials that match the individual patients in some experiments. Furthermore, we propose a novel process of artificial bone material development that consistently predicts the in vivo response of material fabrication to animal testing with the help of machine learning, leading to the establishment of an alternative method of animal testing to replace animal testing.

Author Contributions

Conceptualization, M.A.; methodology, M.A, H.K. and S.H.; software, S.H., K.M. and H.K.; –formal analysis, S.H.; investigation, S.H., K.N., M.N. and K.S.; resources, H.N. and M.A.; writing original draft preparation, S.H.; writing, reviewing and editing, K.S., H.K. and M.A.; supervision, H.K. and M.A.; funding acquisition, M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by a Grant-in-Aid for Scientific Research (KAKENHI) (Grant Number 20H04538) from the Japan Society for the Promotion of Science.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data supporting the findings of this study are available from the corresponding author (M.A.) upon reasonable request.

Conflicts of Interest

The funders had no role in the design of this study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Alajlouni, A.D.; Bliuc, D.; Tran, S.T.; Blank, D.R.; Center, R.J. Muscle strength and physical performance contribute to and improve fracture risk prediction in older people: A narrative review. Bone 2023, 172, 116755. [Google Scholar] [CrossRef] [PubMed]
  2. Gibon, E.; Lu, Y.L.; Nathan, K.; Goodman, B.S. Inflammation, aging, and bone regeneration. J. Orthop. Translat. 2017, 10, 28–35. [Google Scholar] [CrossRef]
  3. St John, A.T.; Vaccaro, R.A.; Sah, P.A.; Schaefer, M.; Berta, C.S.; Albert, T.; Hilibrand, A. Physical and monetary costs associated with autogenous bone graft harvesting. Am. J. Orthop. 2003, 32, 18–23. [Google Scholar]
  4. Gristina, G.A. Biomaterial-centered infection: Microbial adhesion versus tissue integration. Science 1987, 237, 1588–1595. [Google Scholar] [CrossRef]
  5. Banwart, C.J.; Asher, A.M.; Hassanein, S.R. Iliac crest bone graft harvest donor site morbidity. A statical evaluation. Spine 1995, 20, 1055–1060. [Google Scholar] [CrossRef] [PubMed]
  6. Jiao, C.; Xie, D.; He, Z.; Liang, H.; Shen, L.; Yang, Y.; Tian, Z.; Wu, G.; Wang, C. Additive manufacturing of bio-inspired ceramic bone scaffolds: Structural design, mechanical properties and biocompatibility. Mater. Des. 2022, 217, 110610. [Google Scholar] [CrossRef]
  7. Wang, W.; Yeung, K.W.K. Bone grafts and biomaterials substitutes for bone defect repair: A review. Bioact. Mater. 2017, 2, 224–247. [Google Scholar] [CrossRef]
  8. Yamasaki, N.; Hirao, M.; Nanno, K.; Sugiyasu, K.; Tamai, N.; Hashimoto, N.; Yoshikawa, H.; Myoui, A. A comparative assessment of synthetic ceramic bone substitutes with different composition and microstructure in rabbit femoral condyle model. J. Biomed. Mater. 2009, 91B, 788–798. [Google Scholar] [CrossRef]
  9. Bohner, M.; Döbelin, N.; Santoni, B.L.G. β-tricalcium phosphate for bone substitution: Synthesis and properties. Acta Biomater. 2020, 113, 23–41. [Google Scholar] [CrossRef]
  10. von Doernberg, M.C.; von Rechenberg, B.; Bohner, M.; Grünenfelder, S.; van Lenthe, G.H.; Müller, R.; Gasser, B.; Mathys, R.; Baroud, G.; Auer, J. In vivo behavior of calcium phosphate scaffolds with four different pore sizes. Biomaterials 2006, 27, 5186–5198. [Google Scholar] [CrossRef]
  11. Dias, M.; Fernandes, P.; Guedes, J.; Hollister, S. Permeability analysis of scaffolds for bone tissue engineering. J. Biomech. 2012, 45, 938–944. [Google Scholar] [CrossRef]
  12. Shibahara, K.; Hayashi, K.; Nakashima, Y.; Ishikawa, K. Effects of channels and micropores in honeycomb scaffolds on the reconstruction of segmental bone defects. Bioact. Mater. 2022, 6, 490–502. [Google Scholar] [CrossRef]
  13. Aizawa, M.; Porter, E.A.; Best, M.S.; Bonefield, W. Ultrastructural observation of single-crystal apatite fibres. Biomaterials 2005, 26, 3427–3433. [Google Scholar] [CrossRef]
  14. Aizawa, M.; Ueno, H.; Itatani, K.; Okada, I. Synthesis of calcium-deficient apatite fibers by a homogeneous precipitation method and their characterizations. J. Eur. Ceram. Soc., 2006, 26, 501–507. [Google Scholar] [CrossRef]
  15. Aizawa, M.; Howell, S.F.; Itatani, K.; Yokogawa, Y.; Nishizawa, K.; Toriyama, M.; Kameyama, T. Fabrication of porous ceramics with well-controlled open pores by sintering of fibrous hydroxyapatite particles. J. Ceram. Soc. Jpn. 2000, 108, 249–253. [Google Scholar] [CrossRef]
  16. Aizawa, M.; Ito, M.; Takeoka, Y.; Rikukawa, M.; Okada, I. Fabrication of porous tricalcium phosphate ceramics from calcium-phosphate fibers for a matrix of biodegradable ceramics/polymer hybrids. Phosphorus Res. Bull. 2004, 17, 209–210. [Google Scholar] [CrossRef] [PubMed]
  17. Kawata, M.; Uchida, H.; Itatani, K.; Okada, I.; Koda, S.; Aizawa, M. Development of porous ceramics with well-controlled porosities and pore sizes from apatite fibers and their evaluations. J. Mater. Sci: Mater. Med., 2004, 15, 817–823. [Google Scholar] [CrossRef] [PubMed]
  18. Yamada, Y.; Inui, T.; Kinoshita, Y.; Shigemitsu, Y.; Honda, M.; Nakano, K.; Matsunari, H.; Nagaya, M.; Nagashima, H.; Aizawa, M. Silicon-containing apatite fiber scaffolds with enhanced mechanical property express osteoinductivity and high osteoconductivity. J. Asian Ceram. Soc. 2019, 7, 101–108. [Google Scholar] [CrossRef]
  19. Shigemitsu, Y.; Nagashima, H.; Matsunari, H.; Aizawa, M. In vivo evaluation of calcium-phosphate ceramics with highly-interconnected pores using porcine tibia defect model. Solid State Phenomena 2022, 340, 113–117. [Google Scholar] [CrossRef]
  20. Yokota, T.; Miki, T.; Honda, M.; Ikeda, T.; Ishii, K.; Matsumoto, M.; Aizawa, M. Fabrication and biological evaluation of hydroxyapatite ceramics including bone minerals. J. Ceram. Soc. Jpn. 2018, 126, 99–108. [Google Scholar] [CrossRef]
  21. Okitsu, S.; Yokota, T.; Aizawa, M. Effect of ball-milling treatment on sinterability of hydroxyapatite ceramics including bone minerals. Phosphorus Res. Bull. 2022, 38, 25–31. [Google Scholar] [CrossRef]
  22. Fontana, F.; Figueiredo, P.; Martins, J.P.; Santos, H.A. Requirements for animal experiments: Problems and challenges. Small. 2020, 17, e2004182. [Google Scholar] [CrossRef] [PubMed]
  23. Kaufmann, K.; Maryanovsky, D.; Mellor, M.W.; Zhu, C.; Rosengarten, S.A.; Oses, C.; Toher, C.; Curtarolo, S.; Vecchio, S.K. Discovery of high-entropy ceramics via machine learning. NPJ Comput. Mater. 2020, 6, 42. [Google Scholar] [CrossRef]
  24. Huang, D.; Li, Z.; Wang, K.; Zhou, H.; Zhao, X.; Zhang, R.; Wu, J.; Liang, J.; Zhao, L. Probing the efffect of photovoltaic material on Voc in ternary polymer solar cells with non-fullerene acceptors by machine learning. Polymers 2023, 15, 2954. [Google Scholar] [CrossRef] [PubMed]
  25. Guo, Y.; Rui, S.; Xu, W.; Sun, C. Machine learning method for fatigue strength prediction of nickel-based superalloy with various influencing factors. Materials 2023, 16, 46. [Google Scholar] [CrossRef] [PubMed]
  26. Yu, J.; Wang, Y.; Dai, Z.; Yang, F.; Fallahpour, A.; Nasiri-Tabrizi, B. Structual features modeling of substituted hydroxyapatite nanopowders as bone fillers via machine learning. Ceram. Int. 2021, 47, 9034–9047. [Google Scholar] [CrossRef]
  27. Kwaria, J.R.; Mondarte, Q.A.E.; Tahara, H.; Chang, R.; Hayashi, T. Data-driven prediction of protein adsorption on self-assembled monolayers toward material screening and design. ACS Biomater. Sci. Eng. 2020, 6, 4949–4956. [Google Scholar] [CrossRef] [PubMed]
  28. Wu, C.; Entezari, A.; Zheng, K.; Fang, J.; Zreiqat, H.; Steven, G.P.; Swain, M.V.; Li, Q. A machine learning-based multiscale model to predict bone formation in scaffolds. Nat. Comput. Sci. 2021, 1, 532–541. [Google Scholar] [CrossRef]
  29. Motojima, K.; Shiratsuchi, R.; Suzuki, K.; Aizawa, M.; Kaneko, H. Machine learning model for predicting the material properties and bone formation rate and direct inverse analysis of the model for new synthesis conditions of bioceramics. Ind. Eng. Chem. Res. 2023, 62, 5898–5906. [Google Scholar] [CrossRef]
  30. Shimizu, N.; Kaneko, H. Direct inverse analysis based on gaussian mixture regression for multiple objective variables in material design. Mater. Des. 2020, 196, 109168. [Google Scholar] [CrossRef]
  31. Bishop, C.M. Pattern Recognition and Machine Learning; Springer: New York, 2006. [Google Scholar]
  32. Filzmoser, P.; Liebmann, B.; Varmuza, K. Repeated Double cross validation. J. Chemom. 2009, 23, 160–171. [Google Scholar] [CrossRef]
  33. Kaneko, H. Cross-validated permutation feature importance considering correlation between features. Anal. Sci. Adv. 2022, 3, 278–287. [Google Scholar] [CrossRef]
  34. Rabiei, M.; Palevicius, A.; Monshi, A.; Nasiri, S.; Vilkauskas, A.; Janusas, G. Comparing methods for calculating nano crystal size of natural hydroxyapatite using X-ray diffraction. Nanomaterials 2020, 10, 1627. [Google Scholar] [CrossRef]
  35. Tetko, V.I.; Sushko, I.; Pandey, K.A.; Zhu, H.; Tropsha, A.; Papa, E.; Oberg, T.; Todeschini, R.; Fourches, D.; Varnek, A. Critical assessment of QSAR models of environmental toxicity against tetrahymena pyriformis: Focusing on applicability domain and overftting by variable selection. J. Chem. Inf. Model. 2008, 48, 1733–1746. [Google Scholar] [CrossRef]
  36. Cover, T.; Hart, P. Nearest neighbor pattern classification. IEEE Trans. Inf. Theory 1967, 13, 21–27. [Google Scholar] [CrossRef]
  37. Kasten, P.; Beyen, I.; Niemeyer, P.; Luginbühl, R.; Bohner, M.; Richter, W. Porosity and pore size of b-tricalcium phosphate scaffold can influence protein production and osteogenic differentitation of human mesenchymal stem cells: An in vitro and in vivo study. Act. Biomater. 2008, 4, 1904–1915. [Google Scholar] [CrossRef] [PubMed]
Preprints 94507 g001
Figure 1. Schematic diagram of a model estimating bone-forming ability.
Figure 1. Schematic diagram of a model estimating bone-forming ability.
Preprints 94507 g001
Preprints 94507 g002
Figure 2. Methods for determination of bone-formation rate.
Figure 2. Methods for determination of bone-formation rate.
Preprints 94507 g002
Preprints 94507 g003
Figure 3. Actual Y1 versus estimated Y1 in DCV: (a) Porosity, (b) compressive strength, (c) FWHM211, (d) FWHM220.
Figure 3. Actual Y1 versus estimated Y1 in DCV: (a) Porosity, (b) compressive strength, (c) FWHM211, (d) FWHM220.
Preprints 94507 g003
Preprints 94507 g004
Figure 4. Actual Y2 versus estimated Y2 in DCV: (a) Without FWHM, (b) with FWHM.
Figure 4. Actual Y2 versus estimated Y2 in DCV: (a) Without FWHM, (b) with FWHM.
Preprints 94507 g004
Preprints 94507 g005
Figure 5. Feature importance of Model 2.
Figure 5. Feature importance of Model 2.
Preprints 94507 g005
Preprints 94507 g006
Figure 6. Comparison of actual porosity and compressive strength with estimated ones. Error bars for actual values represent the standard deviation (n = 5); error bars for estimated values represent the standard deviation of error between estimated values and actual ones in DCV.
Figure 6. Comparison of actual porosity and compressive strength with estimated ones. Error bars for actual values represent the standard deviation (n = 5); error bars for estimated values represent the standard deviation of error between estimated values and actual ones in DCV.
Preprints 94507 g006
Preprints 94507 g007
Figure 7. VB staining of AFSs. High magnification: highly magnified images of dashed square area. Visible light: osteoid purplish red; hypocalcified bone: orange; calcified bone: light brown. Fluorescence: osteoid: red; hypocalcified bone: orange; calcified bone: green.
Figure 7. VB staining of AFSs. High magnification: highly magnified images of dashed square area. Visible light: osteoid purplish red; hypocalcified bone: orange; calcified bone: light brown. Fluorescence: osteoid: red; hypocalcified bone: orange; calcified bone: green.
Preprints 94507 g007
Preprints 94507 g008
Figure 8. Comparison of actual bone-formation rate with estimated ones. Error bars for estimated values represent the standard deviation of error between estimated bone-formation rates and actual ones in DCV.
Figure 8. Comparison of actual bone-formation rate with estimated ones. Error bars for estimated values represent the standard deviation of error between estimated bone-formation rates and actual ones in DCV.
Preprints 94507 g008
Table 1. X1, Y1, X2, Y2 used in construction of models.
Table 1. X1, Y1, X2, Y2 used in construction of models.
Variables Variable names Variables Variable names
X1 Ca(OH)2 / mol･dm-3 X1 20 mm CB rate / %
MgCl2･6H2O / mol・dm-3 Molding pressure / MPa
NaCl / mol・dm-3 Firing tempreture / °C
KCl / mol・dm-3 Firing atmosphere
(NH4)2CO3 / mol・dm-3 (0: air, 1: steam, 2: steam in carbonate gas)
H3PO4 / mol・dm-3 Y1 Porosity / %
NH4F / mol・dm-3 Compressive strength / MPa
Ball mill grinting hour / h Half-width of 4 peaks / degree
(NH4)2HPO4 / mol・dm-3 X2 Implantation periods / weeks
Si(OC2H5)4 / mol・dm-3 Implantation animals
Heating tempreture / °C (0: pig, 1: rat)
Heating time / h Vascular endothelial growth factor adding
Amount of carbon beads (CB) / mass% (0: without, 1: with)
150 mm CB ratio / % Y2 Bone-formation rate / %
Table 2. r2dcv and RMSEdcv for Each Feature Model 1.
Table 2. r2dcv and RMSEdcv for Each Feature Model 1.
Material properties r2dcv RMSEdcv
Porosity 0.933 3.02
Compressive strength 0.745 0.12
FWHM211 0.935 0.02
FWHM220 0.996 0.01
Table 3. r2dcv and RMSEdcv for without or with FWHM Model 2.
Table 3. r2dcv and RMSEdcv for without or with FWHM Model 2.
Feature r2dcv RMSEdcv
Without FWHM 0.497 14.33
With FWHM 0.689 11.30
Table 4. Results of inverse analysis.
Table 4. Results of inverse analysis.
X1 Y1 X2 Y2
Amount of CB
[mass%]		 150 mm CB rate
[%]		 Porosity
[%]		 Compressive
strength
[MPa]		 Crystalline phase Implantation
weeks		 Implantation
animal		 Bone-formation rate
[%]
100 100 70.0 4.7 HAp 12 Pig 37.0
100 70 68.1 5.4 36.8
300 50 84.0 0.9 49.6
350 50 85.4 0.7 51.6
600 60 90.5 0.2 59.4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated