Personalized 3D Printing of Artificial Vertebrae: A Predictive Bone Density Modeling Approach for Robotic Cutting Applications

Abstract

Robotic vertebral plate cutting poses significant challenges due to the complex bone structures of the lumbar spine, which consist of varying densities in cortical and cancellous regions. This study addresses these challenges by developing a predictive model for robotic vertebral plate cutting force and bone quality recognition through the fabrication of artificial vertebrae with controlled, consistent bone density. To address the variability in bone density between cortical and cancellous regions, CT data are utilized to predict target bone density, serving as a foundation for determining the optimal 3D printing process parameters. The proposed methodology integrates a Response Surface Methodology (RSM), Back Propagation (BP) neural network, and genetic algorithm (GA) to systematically evaluate the effects of key process parameters, such as the filling density, material flow rate, and layer thickness, on the printed vertebrae’s density. A one-factor experimental approach and RSM-based central composite design are applied to build an initial bone density prediction model, followed by Sobol’s sensitivity analysis to quantify the influence of each parameter. The GA-BP neural network model is then employed to rapidly and accurately identify optimal printing parameters for different bone layer densities. The resulting optimized models are used to fabricate personalized artificial lumbar vertebrae, which are subsequently validated through robotic cutting experiments. This research not only contributes to the advancement in personalized 3D printing technology but also provides a reliable framework for developing patient-specific surgical planning models in robot-assisted orthopedic surgery.

1. Introduction

Degenerative diseases of the lumbar spine are increasingly prevalent among middle-aged and elderly individuals, often necessitating laminectomy as a treatment. This surgical procedure poses significant challenges to surgeons, as extended operating times can lead to wrist fatigue or tremors, increasing the risk of inadvertent damage to adjacent blood vessels or nerves. The rapid advancement in robot-assisted vertebral plate cutting technology offers a promising solution to this issue [1,2]. During robotic cutting of the vertebral plate, the bone’s composition—comprising cortical and cancellous sections with varying densities—directly influences the cutting force. To emulate the discernment of an experienced surgeon in real time, it is essential to develop a model that mirrors the patient-specific bone density of the vertebral plate. This model serves as a cornerstone for the predictive modeling of cutting forces and bone quality recognition.

With the swift progress in medical imaging and 3D printing technologies, fused deposition modeling (FDM) stands out due to its high precision, simplified process, and ability to fabricate complex internal porous structures, making it particularly suitable for irregular bones such as the spine [3,4]. In orthopedic healthcare, 3D printing enables the customization of implants and prosthetics tailored to a unique patient anatomy [5,6], assists surgeons in the pre-surgical planning and simulation of intricate orthopedic procedures [7,8], and holds promise for the direct printing of living tissues and bone structures, offering novel prospects for future orthopedic therapies [9,10]. However, despite these advancements, the current 3D printing methods face significant limitations in achieving the necessary precision and personalization required for patient-specific vertebrae, particularly when attempting to replicate the complex microstructure and varying bone densities found in human vertebrae. These limitations stem from the challenges in precisely controlling the mechanical properties and microstructural details, which are critical for robotic-assisted surgeries.

For example, Yang et al. [11] leveraged CT data to design artificial vertebrae of specified dimensions, followed by the fabrication of a titanium alloy porous prosthesis using 3D printing technology, resulting in favorable postoperative recovery for the patient. Hasan et al. [12] coated 3D-printed prostheses with an elastin-like recombinant (ELR) to promote the ordered deposition of apatite on their surfaces, facilitating a non-cellular mineralization process critical for hard tissue regeneration. Hu et al. [13] successfully implanted 3D-printed, freestanding artificial vertebrae into patient spines, achieving robust fusion in every case with excellent fusion rates and positive short-term outcomes. Kim et al. [14] harnessed the capabilities of FDM printing to manufacture polycaprolactone (PCL)/hydroxyapatite (HA) biocompatible filaments, creating a series of filaments with variable HA content and assessing their mechanical properties. Elsen et al. [15] integrated artificial intelligence into 3D printing, enabling the system to predict, adapt, and autonomously control parameters, facilitating the development of patient-specific bone scaffolds to meet diverse requirements and minimizing error risks. Han et al. [16] created a bioskeletal scaffold with a trabecular structure to mimic the cellular microenvironment during growth. Jiao et al. [17] incorporated zirconia into hydroxyapatite as a matrix material, designing and fabricating 15 sets of bone scaffolds with varying degrees of irregularity and porosity by 3D printing, and exploring the impacts of structural parameters and materials on mechanical properties and biocompatibility. Kantaros et al. [18,19] explored the potential of biomimetic smart materials in 4D printing technology, highlighting the transformative role of these materials in a number of fields, including regenerative medicine, robotics, architecture, and aerospace, which will drive breakthroughs in areas such as personalized medicine, adaptive robotics, and sustainable architecture.

Despite these promising developments, current studies focus primarily on improving biomechanical properties and biocompatibility, rather than addressing the critical need for the accurate replication of patient-specific bone densities and microstructures in the context of robotic-assisted surgeries. The lack of precise bone density matching between artificial models and real vertebrae significantly hinders the application of these models in robotic osteotomy and surgical planning, as bone density directly influences the cutting force and tactile control in robotic surgeries. It is evident that while 3D printing has been widely applied to produce medical models, surgical guides, artificial bone implants, and bioskeletal scaffolds, the primary focus has been on achieving adequate biomechanical properties and biocompatibility, rather than precisely replicating patient-specific bone densities and microstructures for robotic surgical applications. In robotic osteotomy, the cutting force is influenced by physical attributes such as bone density, underscoring the necessity for a bone surrogate material with equivalent bone density for robotic surgical planning and tactile control research. Thus, a novel approach is required to address these challenges by developing artificial vertebrae that precisely mimic both the geometric structure and bone density of natural vertebrae for use in robotic-assisted surgeries.

Given that the bone density of artificial vertebrae is significantly influenced by the parameters of the 3D printing process, achieving artificial vertebrae with a targeted bone density proves challenging. Consequently, there is a pressing need to investigate a method for personalized 3D printing of artificial vertebrae based on bone density prediction. Through a bone density prediction model, one can derive the 3D printing process parameters tailored to individual patients, ensuring that the density and geometric accuracy of artificial vertebrae closely match those of actual human bone. This forms a solid foundation for predictive modeling of robotic vertebral plate cutting force and bone identification. Utilizing the response surface method to construct a predictive model between bone density and printing process parameters can offer extensive data support. By incorporating new sample data, this method facilitates deep learning and training, gradually enhancing the accuracy of model predictions. As a result, this approach overcomes the limitations of conventional prediction methods, substantially reduces the number of required experiments, improves operational efficiency, and presents a novel solution for bone density prediction [20].

The density of artificial vertebrae is contingent upon numerous 3D printing process parameters, and acquiring the specific densities of artificial vertebrae models demands precise control over these parameters. Traditional prediction methods, such as univariate and orthogonal testing, typically require substantial experimental resources and exhibit low prediction accuracy. In this regard, the integration of intelligent computing techniques like neural networks and genetic algorithms holds significant potential for enhancing precision. Methods that integrate neural networks and genetic algorithms require only a small dataset, effectively constructing models linking process parameters to the densities of finished products. These methods boast significant advantages in parameter optimization and performance prediction. Deswal et al. [21] employed Response Surface Methodology–Genetic Algorithm (RSM-GA), Artificial Neural Network (ANN), and ANN–Genetic Algorithm (ANN-GA) hybrid tools to develop mathematical models. Their aim was to establish indirect correlations between various FDM process parameters and part accuracy for process parameter optimization, thereby enhancing the precision of printed products. Dwivedi et al. [22] explored the relationship between process parameters and the mechanical properties of the lattice structure in 3D-printed Polylactic Acid (PLA+) by employing Particle Swarm Optimization (PSO) to refine the settings of Random Forest and XGBoost models, thereby augmenting their predictive capabilities.

In this study, we aim to overcome the limitations of current artificial vertebra fabrication methods by integrating modern 3D printing technology, image processing, and intelligent prediction techniques. Our research focuses on developing a personalized 3D printing method based on bone density prediction to prepare artificial vertebrae that closely replicate patient-specific bone density. By applying methods such as the Response Surface Methodology, Back Propagation neural networks, and genetic algorithms, we can predict 3D printing parameters for artificial vertebrae more accurately, leading to significant improvements in the preparation and performance of artificial vertebrae for robotic-assisted surgeries. This approach addresses the current gap in the preparation of artificial vertebrae that accurately replicate both the geometric and bone density characteristics of natural vertebrae, which is critical for robotic-assisted surgeries. By doing so, we aim to provide a reliable method for producing artificial vertebra models that support both surgical planning and robotic tactile feedback, ultimately contributing to safer and more effective spinal surgeries.

This study holds broad implications for addressing the challenge of obtaining isolated artificial vertebra models with matching bone density through a 3D printing preparation method grounded in bone density prediction.

2. Three-Dimensional Printing Program of Artificial Vertebrae and Analysis of Influencing Factors

2.1. Three-Dimensional Printing Program of Artificial Vertebrae

Based on the anatomical structure of the human lumbar vertebrae and adhering to the principles of the 3D printing process, a comprehensive strategy has been developed to fabricate an artificial vertebral model with targeted bone density. As illustrated in Figure 1, this approach is divided into two key components: the 3D modeling of the lumbar vertebrae, which involves digital reconstruction and density mapping, and the prediction of critical printing process parameters, which ensures precise control over bone density distribution during the fabrication process.

(1)

Lumbar spine 3D modeling program

The image acquisition phase represents the initial step toward crafting a lumbar spine model. Initially, CT data of the human lumbar spine are acquired through imaging, with particular attention given to optimizing image resolution and selecting appropriate slice spacing during the scanning procedure. Following image acquisition, the process of image segmentation ensues, with the primary objective of isolating the region of interest amidst the intricate array of anatomical structures. Utilizing Materialise Mimics 21.0 medical modeling software, the lumbar spine’s 3D model is generated by judiciously setting the threshold range to delineate the desired area. Through image segmentation, the target region is extracted, culminating in the production of a preliminary model. Ultimately, to refine the surface quality and geometric precision of the model, Geomagic Studio 2012 reverse engineering software is employed. This software facilitates wrapping, smoothing, and reducing the number of triangular facets on the model’s surface, thereby enhancing its overall fidelity.

(2)

Three-Dimensional printing process parameter prediction program

Initially, a single-factor experiment is conducted to identify the pivotal printing process parameters that significantly impact the specimen’s density. Subsequently, Response Surface Methodology is applied, opting for a central composite design coupled with a quadratic polynomial response surface model. Here, principal printing process parameters serve as independent variables, while the density of the printed specimen acts as the dependent variable or response. Distinct prediction models are then formulated for both cortical and cancellous bone densities. These models provide a theoretical foundation for forecasting the printing process parameters required to achieve a specified bone density, thus enabling the more accurate and efficient production of bone-mimetic structures.

Based on the established bone density prediction model, data are randomly divided into training and test samples. A Back Propagation (BP) neural network model is then constructed using these samples. This model serves as the foundation for deriving the fitness function for the genetic algorithm (GA). To establish a bone density prediction and optimization model, we combined several algorithms. First, we collected data on how different printing parameters (such as filling density, layer height, and material flow rate) affect the density of printed samples. Next, we used Response Surface Methodology (RSM) to analyze the relationship between these parameters and bone density, helping us find the best settings. Then, we employed a Back Propagation (BP) neural network to learn from these data and identify patterns between the input parameters and bone density. Finally, we applied a genetic algorithm (GA) to further optimize the printing parameters, simulating the process of natural selection to find the most effective combinations. By integrating these methods, we can accurately predict and optimize the bone density of 3D-printed artificial vertebrae, ensuring their quality. Ultimately, the integration of the GA with the BP neural network allows for the effective prediction and optimization of the 3D printing parameters.

2.2. Effect of 3D Printing Process Parameters on Bone Density

The process parameters significantly impact the quality and density of printed parts. Among these, factors influencing density during the printing process include the filling density, material flow rate, layer thickness, print speed, and nozzle temperature, among others [23,24].

In this paper, specimens are prepared under various process parameters for comparative density experiments. Initially, the sample model is created using 3D design software SolidWorks 2021. The sample is designed as a small cube measuring 30 × 30 × 10 mm. Subsequently, slicing software UltiMaker Cura 5.4 is employed to process the layering information, enabling the samples to be printed.

2.2.1. Effect of Filling Density on Specimen Density

To investigate the impact of filling density on the density of the printed specimens, filling density is selected as a variable parameter while all other parameters are held constant at their default values, which are listed in

Table 1. Default values for other process parameters when selecting fill density as a variable parameter.

Parameters	Layer Thickness (mm)	Material Flow (%)	Print Speed (mm/s)	Molding Temperature (°C)
Numerical value	0.2	100	60	210

. Specimens are then printed with varying filling densities of 50%, 60%, 70%, 80%, 90%, and 100%.

To visualize the relationship between filling density and specimen density, the filling density vs. specimen density curve is depicted in Figure 2. As observed from the figure, there is a significant increase in specimen density as the filling density rises. Specifically, when the filling density increases from 50% to 100%, the specimen density escalates from 0.625 g/cm³ to 1.246 g/cm³.

2.2.2. Effect of Layer Thickness on Specimen Density

To investigate the influence of layer thickness on the density of the specimen, six distinct layer thicknesses are chosen: 0.10 mm, 0.15 mm, 0.20 mm, 0.25 mm, 0.30 mm, and 0.35 mm. During the experiments, all other parameters are kept constant, with the specific settings detailed in

Table 2. Default values for other process parameters when selecting layer height as a variable parameter.

Parameters	Filling Density (%)	Material Flow (%)	Print Speed (mm/s)	Printing Temperature (°C)
Numerical value	100	100	60	210

.

In order to more intuitively see the relationship between fill density and specimen density, Figure 3 shows the layer thickness–specimen density curve; it can be seen that the specimen density with the increase in layer thickness changed in a small manner but did have an impact, the layer thickness increased from 0.10 mm to 0.35 mm, and the specimen density increased from 1.240 to 1.249.

2.2.3. Effect of Material Flow on Specimen Density

In investigating the influence of the material flow rate on specimen density, six distinct flow rate values are selected: 80%, 85%, 90%, 95%, 100%, and 105%. All other process parameters are maintained at constant levels to ensure the integrity of the experimental design. The baseline settings for the remaining molding parameters are detailed in

Table 3. Default values for other process parameters when selecting material flow as a variable parameter.

Parameters	Filling Density (%)	Layer Thickness (mm)	Print Speed (mm/s)	Printing Temperature (°C)
Numerical value	100	0.2	60	210

. The density of the specimen increases with an increasing material flow rate, as visualized from the experimental results in Figure 4.

2.2.4. Effect of Printing Speed on Specimen Density

To investigate the impact of printing speed on specimen density, five distinct printing speeds (in mm/s) are chosen: 30, 50, 70, 90, and 110. Throughout the experiment, all other process parameters are held constant to isolate the effect of printing speed. The baseline settings for the remaining molding parameters are provided in

Table 4. Default values for process parameters when selecting printing speed as a variable parameter.

Parameters	Filling Density (%)	Layer Thickness (mm)	Material Flow (%)	Printing Temperature (°C)
Numerical value	100	0.2	100	210

.

Based on the measurement outcomes presented in

Table 5. Measurements of specimen density and printing time at different extrusion speeds.

Printing Speed (mm/s)	$\mathbf{Specimen} \mathbf{Density} \mathbf{\left(}{\mathbf{g}/\mathbf{c}\mathbf{m}}^3\mathbf{\right)}$	Printing Time (min)
30	1.246	66
50	1.245	41
70	1.246	30
90	1.244	24
110	1.242	21

, it can be concluded that the printing speed exerts a negligible impact on the density of the printed specimens. Any slight variations in density observed are attributable to inherent machine errors. Moreover, the effect of printing speed on printing duration is notably substantial.

The impact of various process parameters, including the filling density, material flow rate, layer thickness, and printing speed, on bone density has been analyzed through a series of printing experiments. As a result, the primary influential parameters have been identified: the filling density, material flow rate, and layer thickness. This analysis provides a solid foundation for defining the critical factors that will be used in the predictive modeling of bone density.

3. Establishment of Bone Density Prediction Models and Analysis of Parameter Effects

3.1. Bone Density Prediction Modeling

The response surface method (RSM) finds broad application in engineering disciplines, encompassing areas such as structural optimal design, reliability evaluation, dynamic analyses, and the development of predictive models [25]. When conducting optimization designs, it is essential to select representative discrete sampling points. Among various experimental design methods, the central composite design (CCD) is the most commonly employed technique. This design incorporates three categories of test points. The specific composition of these test points is as follows:

$$N=m_c+m_r+m_o,$$

(1)

In this equation, $m_0$ denotes the center point, $m_c$ denotes the cubic point, and $m_r$ denotes the axial point, where the distances of the axial points from the center point are all $\alpha$. This $\alpha$ is called the asterisk arm.

In this paper, a bone density prediction model is developed, with bone density serving as the response variable and printing parameters, specifically the filling density (%), material flow rate (%), and layer thickness (mm), acting as the influencing factors. The model considers three influencing factors (k = 3). An experimental design is adopted with eight levels for each factor ($m_{\mathrm{c}}$ = 8) and six center points ($m_r=6$). Since the center point in this modeling serves primarily to estimate the experimental error, a total of six center points ($m_0$ = 6) are selected, and the star arm length is set to unity ($\alpha$ = 1), facilitating the estimation of the experimental error and the robustness of the model.

After gathering the data samples, selecting an appropriate response surface model is essential to derive the regression model that relates inputs to outputs. Due to its straightforward structure and superior capability in handling both linear and non-linear relationships, the quadratic polynomial method is favored in this context. Accordingly, this method has been adopted as the form of the bone density prediction model in this study. The mathematical representation of the model is as follows:

$$y={\beta }_0+{\displaystyle \sum _{i=1}^k{\beta }_i{x}_i}+{\displaystyle \sum _{i=1}^k{\beta }_{ii}{x}_i^2}+{\displaystyle \sum _i{\displaystyle \sum _j{\beta }_{ij}{x}_i{x}_j}}+\epsilon ,$$

(2)

In this equation, $y$ represents the response value of the model; $\beta$ denotes the matrix of regression coefficients, which includes terms for linear, squared, and interaction effects. Specifically, ${\beta }_i$, ${\beta }_{ii}$, and ${\beta }_{ij}$ are the coefficients corresponding to linear, squared, and interaction factor effects, respectively. $x_i$ and $x_j$ stand for the values of the design parameters, and $\epsilon$ represents the error term, capturing any unexplained variability in the response.

According to the RSM, selecting the central composite design requires determining the range and levels of printing parameters first. The data for the filling density, material flow rate, and layer thickness are input into the experimental design software Design-Expert12.0. Based on the chosen experimental design method and the form of the response surface, the software automatically generates 20 sets of parameter combinations for the printing process. The software then analyzes these combinations, ultimately yielding the experimental data results along with predictive models for the parameters’ effects on the output.

3.2. Results and Discussion

3.2.1. Cortical Bone Density Prediction Modeling and Analysis

The Design-Expert software is used to fit the experimental data to obtain a predictive model of cortical bone density:

$$\begin{array}{l}{Y}_1=0.4445+0.0026A-0.0118B+0.0417C+1.23\times {10}^{-4}AB+1.25\times {10}^{-4}AC\\ +7.5\times {10}^{-4}BC-1.5\times {10}^{-5}{A}^2+6.2\times {10}^{-5}{B}^2-0.2455C^2\end{array},$$

(3)

In this equation, $Y_1$ denotes the quadratic regression equation for cortical bone density, A represents filling density, B denotes material flow, and C represents layer thickness.

The ANOVA results presented in

Table 6. Analysis of variance for $Y_1$ regression equation.

Typology	Sum of Squared Deviations	Degrees of Freedom	Mean Square Error	F-Value	p-Value	Significance
Model	0.1680	9	0.0187	11,958.13	<0.0001	significance
A—Filling density	0.1364	1	0.1364	87,412.05	<0.0001	significance
B—Material flow	0.0311	1	0.0311	19,950.51	<0.0001	significance
C—Layer thickness	0.0001	1	0.0001	43.31	<0.0001	significance
AB	0.0003	1	0.0003	192.30	<0.0001	significance
AC	1.250 × 10−7	1	1.250 × 10−7	0.0801	0.7829
BC	1.125 × 10−6	1	1.125 × 10−6	0.7208	0.4157
A2	5.818 × 10−6	1	5.818 × 10−6	3.73	0.0823
B2	6.568 × 10−6	1	6.568 × 10−6	4.21	0.0673
C2	0	1	0	10.62	0.0086
Residual	0	10	1.561 × 10−6
Misfit term	0	5	2.055 × 10−6	0.9306	0.2446	insignificance
Pure error	5.333 × 10−6	5	1.067 × 10−6

indicate that the p-values for all three parameters (filling density, material flow rate, and layer thickness) are below 0.0001, indicating a highly significant impact of these parameters on the regression equation for cortical bone density. Additionally, the p-value for the interaction between the filling density and material flow rate is also less than 0.0001, highlighting a substantial effect on the predictive model for cortical bone density.

The lack-of-fit F-statistic of 0.9306 for this regression suggests that the cortical bone prediction model fits the observed data well. Moreover, the coefficient of determination, $R^2$ = 0.9967, signifies excellent predictive power, indicating that the model explains 99.67% of the variability in cortical bone density. The adjusted $Adj{R}^2$ = 0.9953 further confirms the model’s robustness, suggesting a reliable predictive capacity for cortical bone density.

As depicted in Figure 5a, the normal probability plot of the residuals for the $Y_1$ prediction model reveals that the residuals are roughly aligned along a straight line. This alignment suggests that the residuals adhere to a normal distribution, thereby affirming the validity and applicability of the $Y_1$ model. Moreover, Figure 5b illustrates the scatter plot comparing the predicted and actual values of the $Y_1$ model. It is evident that the data points closely align with the reference line, indicating a high degree of correlation between the model’s predictions and the observed values. This observation underscores the high predictive accuracy of the cortical bone density prediction model.

Figure 6a displays the surface plot illustrating the interaction effect of the filling density and flow rate on density. From this, we ascertain the impact of four boundary conditions on the specimen’s density: a high filling density and high flow rate (100%, 100%) yield the highest density, followed by a high filling density and low flow rate (100%, 90%), then a low filling density and high flow rate (80%, 100%), and lastly, a low filling density and low flow rate (80%, 90%).

Figure 6b presents the surface plot depicting how filling density interacts with layer thickness to affect the density. Here, we observe the effect of four extreme scenarios on the print block’s density: a high filling density and thick layer (100%, 0.3 mm) produce the densest outcome, succeeded by a high filling density and thin layer (100%, 0.1 mm); subsequently, there is a low filling density and thick layer (80%, 0.3 mm), and finally, a low filling density and thin layer (80%, 0.1 mm).

Figure 6c is the surface plot revealing the interaction between the layer thickness and material flow rate on density. It becomes evident that the interplay between the layer thickness and material flow rate exerts a negligible effect on the density of cortical bone, implying that this interaction is of minor significance.

3.2.2. Predictive Modeling and Analysis of Cancellous Bone Density

A predictive model of cancellous bone density can be obtained by Design-Expert software:

$$\begin{array}{l}{Y}_2=0.0749-0.0004A-0.0017B+0.0573C+1.24\times {10}^{-3}AB+2.50\times {10}^{-3}AC\\ +1.25\times {10}^{-3}BC+6.4647\times {10}^{-6}{A}^2+9.5455\times {10}^{-6}{B}^2-0.1545C^2\end{array},$$

(4)

In this equation, $Y_2$ denotes the quadratic regression equation for cancellous bone density, A represents filling density, B denotes material flow, and C represents layer thickness.

Similarly, the analysis of variance (ANOVA) for the regression equation of cancellous bone density is presented in

Table 7. Analysis of variance for $Y_2$ regression equation.

Typology	Sum of Squared Deviations	Degrees of Freedom	Mean Square Error	F-Value	p-Value	Significance
Model	0.2944	9	0.0327	14,799.1	<0.0001	significance
A—Filling density	0.2819	1	0.2819	12,7500	<0.0001	significance
B—Material flow	0.0097	1	0.0097	4403.35	<0.0001	significance
C—Layer thickness	0.0001	1	0.0001	7.64	0.0200
AB	0.0028	1	0.0028	1255.33	<0.0001	significance
AC	1.125 × 10−6	1	1.125 × 10−6	0.5089	0.4919
BC	1.250 × 10−7	1	1.250 × 10−7	0.0565	0.8168
A2	5.818 × 10−6	1	5.818 × 10−6	2.63	0.1358
B2	2.506 × 10−6	1	2.506 × 10−6	1.13	0.1321
C2	6.568 × 10−6	1	6.568 × 10−6	2.97	0.1155
Residual	0	10	2.211 × 10−6
Misfit term	9.273 × 10−6	5	1.855 × 10−6	0.7226	0.6349	insignificance
Pure error	0	5	2.567 × 10−6

. The ANOVA table shows that the p-values for the filling density and material flow rate are both less than 0.0001, indicating a highly significant influence of these two parameters on the regression equation for cancellous bone density. The p-value for layer thickness, which is less than 0.05, also signifies that its effect on the regression equation is statistically significant.

Among the interaction terms involving these parameters, the p-value for the interaction between the filling density and material flow rate is less than 0.0001, emphasizing that this interaction term significantly affects the predictive model of cancellous bone density. However, the other interaction terms do not demonstrate statistically significant effects.

The lack-of-fit F-statistic for this regression, F = 0.7226, suggests that the cancellous bone prediction model fits the observed data well. The coefficient of determination $R^2$ = 0.9991 indicates exceptional predictive power, signifying that the model explains 99.91% of the variability in cancellous bone density. Moreover, the adjusted $Adj{R}^2$ = 0.9976 confirms that the cancellous bone prediction model effectively fits the data and can predict changes in the dependent variables.

As illustrated in Figure 7a, the normal probability plot of the residuals from the $Y_2$ prediction model reveals a near-linear distribution. This indicates that the residuals are normally distributed, supporting the applicability and validity of the $Y_2$ model. Additionally, as depicted in Figure 7b, the scatter plot of predicted versus actual values for the $Y_2$ model shows that the data points closely align with the reference line. This close alignment signifies that the model’s predicted values are in close proximity to the observed values, confirming the high predictive accuracy of the cancellous bone density model.

Figure 8a presents a surface plot illustrating the interaction effect of the filling density and flow rate on the density of the print block. The analysis reveals that among the four boundary conditions, the combination of a high filling density and high flow rate (40%, 100%) results in the highest density. This is followed by a high filling density with a lower flow rate (40%, 80%), then a low filling density with a high flow rate (10%, 100%), and finally, a low filling density with a low flow rate (10%, 80%).

Figure 8b shows a surface plot demonstrating the interaction effect of filling density and layer thickness on the density of the print block. The findings indicate that a high filling density combined with a thick layer (40%, 0.3 mm) yields the highest density, followed by a high filling density with a thin layer (40%, 0.1 mm). This is succeeded by a low filling density with a thick layer (10%, 0.3 mm), and lastly, a low filling density with a thin layer (10%, 0.1 mm).

Figure 8c depicts a surface plot illustrating the interaction effect of the layer thickness and material flow rate on the density of cancellous bone. The data indicate that the interaction between the layer thickness and material flow rate has a negligible effect on the density, suggesting that this factor does not significantly influence the overall density.

3.3. Parameter Sensitivity Analysis of Bone Density Prediction Models

The Sobol method is a variance-based sensitivity analysis method that enables a sensitivity analysis of strongly non-linear, nonmonotonic models [26] and relies on the estimation of the following values:

$$S={\displaystyle \frac{V_x[E(y|x)]}{V(y)}},$$

(5)

In this equation, $y$ denotes the output variable,$x$ denotes the input variable, $E(y|x)$ represents the expected value of $x$ conditional on $y$ taking a constant value, and $V(y)$ is the total output variance of the output $y$.

Setting the model as $y=f(x),(x=x_1,x_2,\dots ,x_k)$, $x_i$ obeys a uniform distribution, and $f^2(x)$ is integrable in the domain of definition; it can be decomposed as

$$f(x_1,\dots ,x_k)=f_0+{\displaystyle \sum _{i=1}^k{f}_i(x_i)}+{\displaystyle \sum _{1\le i<j\le k}{f}_{ij}(x_i,x_j)}+\dots +f_{1,2,\dots ,k}(x_1,x_2,\dots ,x_k),$$

(6)

Employing the Sobol method, the total variance of model $y$ is decomposed into components attributable to individual parameters and interactions among them. Consequently, the total variance $V$ of model $y$ can be dissected as follows:

$$V={\displaystyle \sum _{i=1}^n{V}_i}+{\displaystyle \sum _{1\le i<j\le n}{V}_{ij}}+\dots +V_{1,2\dots ,n},$$

(7)

In this equation, $V_i$ denotes bias variance produced by parameter $x_i$, $V_{ij}$ presents variance produced by the interaction of parameters $x_i$ and $x_j$, and $V_k$ denotes variance generated by the joint interaction of the k parameters.

The specific steps for conducting a global sensitivity analysis using the Sobol method are as follows: initially, the parameters to be analyzed are selected. In this particular study, the printing parameters—namely, the filling density, material flow rate, and layer thickness—are designated as the factors under investigation. Subsequently, the variation intervals for each of these parameters are established. For cortical bone, the filling density is varied from 80% to 100%, the material flow rate ranges from 90% to 100%, and the layer thickness is adjustable between 0.1 mm and 0.3 mm. Similarly, for cancellous bone, the filling density spans from 10% to 40%, the material flow rate is adjustable within the 80% to 100% interval, and the layer thickness remains consistent, varying between 0.1 mm and 0.3 mm.

Each input print parameter follows a uniform distribution. The formulas for calculating the mean S and total variance $V$ of each parameter using the Monte Carlo method are as follows:

$$f_0={\displaystyle \int f(x)dx\approx \frac{1}{N}}{\displaystyle \sum _{k=1}^{N}f(x_k)},$$

(8)

$$V={\displaystyle \int f^2(x)dx-f_0^2\approx \frac{1}{N}}{\displaystyle \sum _{k=1}^N{f}^2(x_k)}-f_0^2,$$

(9)

In this equation, the mean value $f_0$ is a constant and N is the size of the Monte Carlo sampling, which in this study is 4000, and by fixing a parameter while adjusting the values of the other parameters, the bias or first-order effect $D_i$ of each parameter is computed; that is, the effect of a single parameter on the model output is

$$V_i={\displaystyle \int f_{i}d{x}_i},$$

(10)

The extent to which the interaction between parameters influences the model outcomes is quantified by the subsequent equation:

$$V_{i_1,i_2,\dots ,i_k}=\int {{\displaystyle f}}_{i_1,i_2,\dots ,i_k}^{2}d{x}_{i_1}d{x}_{i_2}\cdots d{x}_{i_k},$$

(11)

The sensitivity index S of the objective function model for a single parameter and its interaction is denoted as

$$S_{i1,\dots ,i_s}={\displaystyle \frac{V_{i1,\dots ,i_s}}{V}},(1\le i_1<\dots <i_s\le k),$$

(12)

Therefore, for parameter $x_i$, the first-order sensitivity coefficient is obtained as $S_i=V_i/V$, which represents the main effect of parameter $x_i$ on the output model result y; i.e., parameter $x_i$ partially contributes to variance $V(y)$. Similarly, for $i\ne j$,$V(y)$ the second-order sensitivity coefficient is $S_{ij}=V_{ij}/V$, which calculates the effect of the cross-talk between parameter $x_i$ and parameter $x_j$ on the model result y.

Normalizing the above equation yields the following equation:

$${\displaystyle \sum _{i=1}^k{S}_i}+{\displaystyle \sum _{1\le i<j\le k}{S}_{ij}}+\dots +S_{1,2,\dots ,k}=1,$$

(13)

The global sensitivity of each input parameter is the sensitivity of all included parameters $x_i$. The global sensitivity definition index can be expressed by the following equation:

$$S_{Ti}=S_i+{\displaystyle \sum _{i\ne j}{S}_{ij}}+\cdots ,$$

(14)

The global sensitivities of the individual printing parameters have been computed, and the outcomes of these sensitivity analyses are summarized in

Table 8. First-order sensitivity and global sensitivity of predictive models.

Predictive Model	Parameters	First-Order Sensitivity	Global Sensitivity
Cortical bone density	Filling density (%)	0.8700	0.8100
	Material flow (%)	0.1878	0.1892
	Layer thickness (mm)	0.0004	0.0005
Cancellous bone density	Filling density (%)	0.9738	0.9658
	Material flow (%)	0.0340	0.0371
	Layer thickness (mm)	0.0000	0.0001

. To provide a more intuitive visualization of the relative impact of various printing parameters on the printed density, the first-order and total-order sensitivities of the prediction model are depicted in bar charts, as illustrated in Figure 9.

As shown in Figure 9, the disparity between the global sensitivity and the first-order sensitivity of the predictive model for bone density across distinct bone layers is negligible, suggesting a minimal role of parameter interactions. The quantitative outcomes indicate that filling density exerts the predominant influence, followed by the material flow rate, while the impact of layer thickness appears to be inconsequential.

4. Prediction of 3D Printing Process Parameters Based on GA-BP Neural Network

4.1. Bone Density BP Neural Network Prediction Modeling

A robust regression model of bone density has been established using RSM. The neural network training and test sample data were generated through quadratic regression modeling. A total of 600 different sets of data were selected, with 450 sets used as training samples and 150 sets as test samples.

In this study, Matlab was utilized to execute the model design. The built-in normalization function ‘mapminmax’ was employed to normalize both the input and output data, transforming the learning and testing sample data into the [−1, 1] interval. After the model training and testing were completed, the inverse normalization procedure was applied to retrieve the actual output values. The normalization formula used is as follows:

$$x\prime ={\displaystyle \frac{x-x_{\min }}{x_{\max }-x_{\min }}},$$

(15)

In this equation, $x$ and $x\prime$ represent the data before and after normalization, respectively; $x_{\max }$ and $x_{\min }$ represent the maximum and minimum values of the data for a sample within a given dimension.

(1)

Design of the input and output layers

The objective of this paper is to construct a network model that correlates printing bone density with processing parameters. As determined in the preceding section, the processing parameters include the filling density (%), flow rate (%), and layer thickness (mm). Consequently, the prediction model formulated in this study features three nodes in the input layer, corresponding to the aforementioned parameters, and a single node in the output layer, representing the bone density.

(2)

Design of hidden layers

To streamline the model architecture, this paper employs a single-hidden-layer BP neural network. Through systematic network construction, we sought to determine the optimal number of nodes within the hidden layer. Initially, the node count was bracketed between 3 and 20. By training the network and evaluating the Euclidean distance errors of the training set, the optimal configuration was identified as the one yielding the minimum error.

As depicted in Figure 10a, setting the hidden layer to 14 nodes minimizes the Euclidean distance error to 0.020709. Conversely, Figure 10b reveals that a five-node hidden layer reduces the error to 0.017854. Consequently, the structure of the BP neural network prediction model for cortical bone density is determined to be 3-14-1, while that for cancellous bone density is 3-5-1, as illustrated in the model structures shown in Figure 11.

(3)

Selection of activation function

The activation function plays a crucial role as the transfer mechanism between the hidden layer and the output layer within a neural network, with its selection directly impacting the network’s output. Given that the input and output data underwent normalization in the preceding section, their values now fall within the range of [−1, +1]. This ensures that the sample data entering the hidden layer via the input layer align with the value domain of the chosen activation function, tansig.

For the neural network prediction model developed herein, the tansig function has been selected as the activation function for the hidden layer, facilitating the processing of data with a bipolar sigmoidal distribution. Meanwhile, the output layer utilizes the purelin function, allowing for the generation of outputs over a continuous, unbounded range.

(4)

Selection of training algorithm

Upon defining the neural network’s structure and activation functions, the selection of an appropriate training algorithm becomes crucial. To identify the most suitable training approach for our model, we assessed the primary training methodologies utilized in Back Propagation neural networks. Our evaluations revealed that the trainlm algorithm offers superior performance, characterized by the fewest iterations and the lowest mean square error on the validation set. Consequently, we elected to employ the trainlm algorithm as the training protocol for our neural network model.

(5)

BP neural network prediction results

The parameter settings utilized in the BP neural network are detailed in

Table 9. BP neural network parameter settings.

Number of Training Sessions	Learning Rate	Training Goal	Momentum Factor	Training Algorithms	Activation Function
1000	0.01	0.00001	0.01	trainlm	Tansig, purelin

.

The comparison between the predicted and actual values of the BP neural network test set, obtained through training, is presented in Figure 12. The relative error of the BP neural network’s prediction for the test set is depicted in Figure 13. Upon examining these figures, it becomes evident that the predicted values closely align with the actual test sample values, with a maximum relative error of 0.15% and a minimum relative error of 0.00005% observed in the prediction of cortical bone density.

To further validate the model’s reliability, three key performance indicators are employed for the analysis: the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE), and the Mean Absolute Percentage Error (MAPE). For the cortical bone density predictions, the metrics reveal MAE = 0.00035, RMSE = 0.00045, and MAPE = 0.033%. In the case of cancellous bone density predictions, the maximum relative error reaches 0.96% with a minimum of 0.00166%, and the associated error metrics are MAE = 0.00046, RMSE = 0.00063, and MAPE = 0.17255%.

Based on these comprehensive data points, it is demonstrable that the established neural network prediction model effectively captures the intricate non-linear relationship between the molding process parameters and the resulting densities.

4.2. GA-BP Neural Network Prediction Method for 3D Printing Process Parameters

The essence of applying the GA-BP neural network algorithm to predict parameters lies in ensuring that the constructed neural network accurately captures the underlying relationship between the input and output variables. This critical aspect is visually depicted in the GA-BP parameter prediction flowchart, as illustrated in Figure 14. The flowchart illustrates the iterative optimization process that refines the network’s architecture and parameters to achieve high prediction accuracy.

(1)

Initial population

The BP neural network prediction model developed in this study for correlating molding process parameters with cortical bone density features an input layer with three nodes, corresponding to the filling density (%), material flow rate (%), and layer thickness (mm). The hidden layer comprises 14 nodes, and the output layer consists of a single node representing the density of the printed component. For the prediction model concerning cancellous bone density, the hidden layer contains five nodes. The length of the chromosome, which is essential for the genetic algorithm to manipulate during the optimization process, is determined as follows:

$$S=n\times l+l\times m+l+m,$$

(16)

In this formula, $n$ represents the number of nodes in the input layer, $m$ represents the number of nodes in the output layer, and $l$ represents the number of nodes in the implicit layer, and the length of the chromosome can be found to be 71 and 26 by bringing the corresponding values into the formula.

Determining the population size, a key parameter, significantly impacts the algorithm’s performance. In this paper, we set the initial population size, M = 30.

The initial iteration number is usually set to a lower value, and then adjusted based on repeated trials to determine the ideal number. In this paper, the maximum iteration number, G = 100, is used after multiple rounds of debugging to find the most suitable figure.

(2)

Adaptation function

The focal parameter in the process parameter prediction addressed in this paper is bone density. When designing the parameter prediction framework, the objective is to engineer a suitable ensemble of parameter combinations that enable the printed bone density to closely approximate, or ideally coincide with, the pre-established target value. In line with this rationale, the BP neural network prediction model detailed in the preceding section serves as the fitness function, guiding the optimization process toward achieving the desired bone density outcomes.

(3)

Genetic operator

In this paper, the roulette method is chosen for the selection operation, in which the probability $P_i$ that individual $i$ is selected is

$$P_i=F_i/{\displaystyle \sum _{i=1}^N{F}_i}=F_i/F_{sum},$$

(17)

In this equation, N is the population size, $F_i$ is the fitness of individual $i$, ${\displaystyle \sum _{i=1}^N{F}_i}$ is the total fitness value of the population, and $P_i$ is the selection probability for individual $i$.

The general crossover probability takes the value of [0.4, 0.9]. In this paper, the crossover probability $P_c$ = 0.75. In this paper, a simulated binary crossover method is used to perform crossover operation, which is commonly used in genetic algorithms for real number coding, and the final expression is

$$X_A^{t+1}=\alpha X_B^t+(1-\alpha )X_A^t,$$

(18)

$$X_B^{t+1}=\alpha X_A^t+(1-\alpha )X_B^t,$$

(19)

In this equation, $X_A^t$ and $X_B^t$ denote that the two chromosomes A and B crossover at position t, and $\alpha$ is a random number taken between [0, 1] for the crossover parameter.

In this section, the basic bitwise mutation method is employed, which involves executing a random mutation operation based on a predefined mutation probability. Typically, the mutation probability falls within the interval [0.001, 0.2]. For the purposes of this study, the mutation probability is set to A = 0.2.

Prediction process: Initially, individuals with high fitness values from the initial population are selected to advance to the subsequent generation. Subsequently, new individuals are generated through crossover and mutation operations, designed to prevent the algorithm from converging on a local optimum. This evolutionary process is iterated until the optimal combination of process parameters is attained.

Assuming the target parameter to be the cortical bone density, denoted as ρ = 0.963 g/cm³, the program is executed in accordance with the established GA-BP neural network model. The variation in individual fitness values throughout the prediction process is graphically represented in Figure 15.

The fitness of the optimal individual exhibits a progressive decline, trending towards stabilization, suggesting that the algorithmic optimization is increasingly converging towards the optimal solution. Upon reaching the 12th generation, the fitness of the optimal individual plateaus at a low level and remains consistent until the algorithm concludes at the 100th generation. At this juncture, the optimal combination of process parameters is identified: a filling density of 81.346%, a material flow rate of 94.436%, and a layer thickness of 0.209 mm.

To further validate the model, ten sets of cortical and cancellous bone density values are randomly selected and fed into the program to derive the corresponding process parameter values. These parameter combinations are then utilized to 3D print specimens, enabling us to measure the actual densities. The relative errors between the predicted and actual densities are calculated, and to visualize these discrepancies more vividly, the relative error rates are plotted on a line graph, as depicted in Figure 16. The relative errors are found to be within 4%, indicating that the deviations fall within an acceptable margin. This validates the suitability of the BP neural network for serving as the fitness function in a genetic algorithm, thereby facilitating accurate predictions of process parameters.

5. Personalized 3D Printing of Artificial Vertebrae

5.1. Data Acquisition and Modeling

The lumbar spine data utilized in this study were sourced from a 40-year-old male patient at a partner hospital using a Siemens SOMATOM Definition Edge 128-slice spiral CT scanner (Munich, Germany). The data acquisition was conducted with the patient’s informed consent and in full compliance with medical ethics standards. A total of 530 CT slice images were obtained and exported in DICOM format. Initially, the DICOM data were imported into Mimics software for 3D model reconstruction. However, the resulting models were typically coarse and lacked the desired quality. To enhance the model’s surface smoothness and structural integrity, the reverse engineering software Geomagic was subsequently employed for optimization and refinement. Finally, the cortical and cancellous bone components were assembled into a unified model using SolidWorks. The workflow for lumbar spine data acquisition and model development is depicted in Figure 17. The lumbar spine modeled according to the above process is shown in Figure 18.

5.2. Bone Density Measurement and 3D Printing

Bone density includes several types such as apparent wet density, apparent dry density, gray density, and bone mineral density [27]. In this paper, the Quantitative Computed Tomography (QCT) method is utilized to quantify the CT value (expressed in Hounsfield Units, Hu) of bones, upon which the bone density value is calculated. With the apparent density serving as the principal investigative parameter, previous studies have conducted linear regression analyses to explore the relationship between CT values and apparent densities. This has led to the derivation of a linear regression equation linking the two variables, as referenced in [28].

$${\rho }_{app}=0.896CT+134,$$

(20)

In this equation, ${\rho }_{app}$ is the apparent density unit (mg/cm³) and CT is the CT value in Hounsfield Units (Hu). The D′ensity in Ellipse′ tool within the Mimics 21.0 software is employed to designate regions of interest (ROIs) within each image slice. Upon selection, the system autonomously computes and displays the CT values, as exemplified in Figure 19. In this illustration, A‘rea’ denotes the extent of the delineated ROI, whereas M′ean (HU)′ signifies the arithmetic mean of the density readings in Hounsfield Units (HU). Specifically, the mean grayscale value of the cancellous bone within the current slice is quantified as 215.51 HU.

In the realm of medical imaging, HU serves as a standardized metric to depict the relative density of tissues in CT scans. The computation of M′ean (HU)′ involves aggregating a series of density measurements and dividing the sum by the total number of data points, yielding the average density. Conversely, S′td.Dev (HU)′ represents the standard deviation of the density data, expressed in HU.

Should the image quality be subpar, marked by the presence of streaks or other imaging artifacts, the affected layer must be excluded from the analysis. In such instances, the conversion of CT values should adhere to the methodology outlined by Bradley [29]. CT values for each slice are meticulously recorded in correspondence with their respective layer numbers. By averaging these values, a representative CT value for each bone segment is derived, facilitating a comprehensive assessment of bone density across the entire structure.

The CT values of the cancellous and cortical bone for each layer were recorded following the aforementioned procedure, resulting in average CT values of 231 HU for cancellous bone and 1156 HU for cortical bone. Upon substituting these average CT values into the previously derived equation, the computed densities were 340.98 mg/cm³ for cancellous bone and 1169.78 mg/cm³ for cortical bone.

After determining the densities of both cortical and cancellous bone, the anticipated printing process parameters are estimated through a hybrid approach combining a BP neural network analysis with genetic algorithm prediction. The resultant printing process parameters are tabulated in

Table 10. Model Printing Parameters.

Bone Type	Bone Density (mg/cm3)	Filling Density (%)	Material Flow (%)	Layer Thickness (mm)
Cancellous bone	340.98	32.83	84.28	0.14
Cortical bone	1169.78	100	93.85	0.19

.

For the fabrication of the vertebrae models, slicing software is employed to define and optimize the 3D printing parameters, ensuring accurate configuration and proper alignment within the build platform. The artificial vertebral bodies are subsequently constructed using 3D printing technology, as shown in Figure 20. Specifically, Figure 20a displays the cancellous bone model, and Figure 20b depicts the cortical bone model, while Figure 20c,d present the composite models that integrate both the cancellous and cortical bone structures.

6. Experiments on Robotic Cutting of Artificial Vertebrae

To validate the efficacy of our method for fabricating artificial vertebral bodies based on bone density predictive modeling, a robotic vertebral plate cutting experiment was devised. The experimental setup comprises a 6-DOF EC66 collaborative robot (Elite Corporation, Shanghai, China), a type of KWR75B six-axis force sensor (Kunwei, Changzhou, China), a ball end milling cutter (China Dezhou Lansuo Medical Equipment Co., Dezhou, China), and a computer for control and data acquisition (HP ProBook, Beijing, China). During the experiment, the 3D-printed artificial vertebrae are securely mounted on a fixture. An overview of the experimental platform is illustrated in Figure 21.

The 6-DOF collaborative robot, equipped with a ball end milling cutter, executes transverse feed cutting along the Y-axis. The transverse feed rate and the longitudinal depth of the cut are controlled by the collaborative robot’s controller and the host computer. The rotation speed of the end mill, set at 5000 rpm, is achieved by commanding the integrated servo motor through the host computer interface. Cutting trials are conducted on both the printed cortical and cancellous bone sections. Throughout the process, a six-axis force sensor is employed to monitor and record the cutting forces exerted by the end mill. For this specific experiment, parameters are set as follows: a rotation speed of 5000 rpm, a feed rate of 3 mm/s, and a depth of the cut of 0.2 mm.

A porcine spine, exhibiting a density akin to that of the artificial vertebra, is chosen for comparative cutting experiments. Utilizing CT scanning, the bone densities of various layers within the porcine spine are determined. The comparative density analysis, summarized in

Table 11. Comparison of artificial vertebral density with porcine spine density.

Bone Type	Artificial Vertebra Density (mg/cm3)	Porcine Spine Density (mg/cm3)	Relative Error (%)
Cortical bone	340.98	353.26	3.60
Cancellous bone	1169.78	1152.36	1.15

, reveals relative errors all confined within a 4% margin. This substantiates the suitability of the porcine spine for implementation in the cutting comparison tests.

A porcine spine, exhibiting a density akin to that of the artificial vertebra, is chosen for comparative cutting experiments. Utilizing CT scanning, the bone densities of various layers within the porcine spine are determined. The comparative density analysis, summarized in Table 11, reveals relative errors all confined within a 4% margin. This substantiates the suitability of the porcine spine for implementation in the cutting comparison tests. To corroborate the hypothesis that distinct bone qualities yield differential force signals during the cutting process, the robot performs cuts on differing bone layers of both the artificial vertebral plate and the porcine spine. The outcomes demonstrate that the force signals indeed vary in response to the varying bone properties. Illustrations of the robotic cutting procedure on the porcine spine are depicted in Figure 22, while the analogous procedure on the artificial vertebral plate is showcased in Figure 23.

Considering the y-direction milling force for comparison, Figure 24 displays the plots of the cutting force signal data for the cortical bone and cancellous bone models, respectively. Similarly, Figure 25 illustrates the cutting force signal data for the cortical and cancellous bone segments of the actual porcine spine.

The data analysis reveals that the average cutting force for the cortical bone model hovers around 1.3 N, whereas for the cancellous bone model, it is approximately 0.5 N. In contrast, the average cutting force for the cortical bone of the porcine spine is around 1.5 N, and for the cancellous bone, it is roughly 0.8 N. These findings validate the cutting force predictions derived from the bone density models, confirming the effectiveness of the predictive model in simulating real bone cutting characteristics.

Despite the overall success of the cutting experiments, several unexpected challenges arose. During some trials, minor discrepancies in force measurements were observed, particularly when cutting regions with complex microstructures. These variations could be attributed to the inherent differences in the microstructures of artificial vertebrae compared to biological bone. Moreover, the surface roughness and material composition of the 3D-printed models occasionally resulted in deviations from expected cutting forces.

The comparison of cutting forces between the artificial models and the porcine spine provides crucial insights into the accuracy of the predictive model. The relatively small differences in force magnitudes between the two suggest that the model can reliably predict cutting behaviors. However, the observed discrepancies, though minimal, highlight areas for potential refinement in material selection and model parameterization. Further investigation into how material anisotropy and microstructural properties affect cutting force signals is necessary to enhance the reliability of the robotic cutting process.

From the experimental findings, it is evident that a significant disparity exists between the cutting forces of the cortical and cancellous bone models, mirroring the substantial differences observed in real bone. This similarity in cutting force characteristics between the models and real bone validates the practicality of the method for preparing artificial vertebrae, which is grounded in bone density predictive modeling. Overall, while the cutting force data validate the accuracy of the bone density-based predictive model, further refinement is needed to address the minor discrepancies and ensure greater consistency in the cutting outcomes for future applications.

7. Conclusions

In this study, we addressed the challenges associated with personalized 3D printing for the preparation of artificial vertebrae, focusing on predictive bone density modeling. This comprehensive approach highlights the feasibility and effectiveness of utilizing predictive modeling techniques, particularly GA-BP neural networks, to enhance the precision and efficiency of personalized artificial vertebrae fabrication based on bone density predictions. The investigation encompassed several key areas:

(1)

The development of programs for 3D printing artificial vertebral bodies, including an analysis of influencing factors: A 3D modeling technique for lumbar vertebrae was established, alongside methods for predicting printing process parameters. Through one-way experiments, primary parameters affecting bone density—such as the filling density, material flow rate, and layer thickness—were identified.

(2)

The systematic development of bone density prediction models using RSM: Experiments were designed to establish predictive models under varying 3D printing conditions (filling density, material flow rate, layer thickness). Sobol’s sensitivity analysis quantitatively assessed each parameter’s influence on bone density, offering critical insights into process optimization.

(3)

The implementation of a GA-BP neural network for predicting 3D printing process parameters: This involved constructing BP neural network models to correlate printing parameters with bone density, optimizing model architecture, activation functions, and training algorithms using genetic algorithms. The resulting GA-BP framework provided robust predictions of optimal process parameters.

(4)

The construction and preparation of lumbar spine models, including the L4 segment and composite vertebral models suitable for 3D printing: Bone density values derived from CT scans of the L4 segmental lumbar spine were integrated into the GA-BP neural network to determine precise 3D printing parameters, facilitating accurate model preparation.

(5)

The execution of robotic cutting experiments on artificial vertebrae to validate the proposed preparation technique: A specialized setup enabled the collection of cutting force data, demonstrating the practicality of the method based on bone density predictions.

Although the current research confirms the effectiveness of using predictive bone density models for personalized vertebral fabrication, several areas warrant further investigation. Future studies could focus on refining the GA-BP neural network by incorporating additional physiological parameters such as anisotropic mechanical properties or incorporating stochastic elements to account for individual variations in bone quality. These enhancements would lead to a more comprehensive and robust predictive model for diverse patient conditions.

Additionally, conducting preclinical and clinical trials is crucial to evaluate the safety, reliability, and long-term stability of the fabricated vertebral models under realistic physiological conditions. Such trials will offer deeper insights into the clinical utility of these artificial vertebrae, providing data to guide potential modifications to the predictive modeling approach.

Beyond vertebral applications, exploring the applicability of this method to other anatomical structures, such as the pelvis or cranial bones, would expand its utility. Further, integrating smart materials into the 3D printing process could enhance the adaptability of the printed models, allowing them to respond dynamically to changes in their biomechanical environment, thereby paving the way for more advanced and patient-specific orthopedic and surgical solutions.