Inner design of artificial test bones for biomechanical investigations using topology optimization

Artificial or human test bones are used for the biomechanical testing of implants. Human test bones are rare and not always available. These must, therefore, be substituted with artificial test bones. However, current artificial test bones are only available with specific characteristics (e.g., age groups or disease characteristics). Additionally, their mechanical properties are only comparable to a limited extent to those of a human bone. This paper presents a methodology for designing additively manufactured artificial test bones for biomechanical testing that replicate the mechanical behavior of a human bone. Topology optimization methods are used to generate the artificial test bone's internal structure. The geometric model is based on a computed tomography dataset of a human bone. The input data can be manipulated in advance to reproduce defects or disease patterns. The bone was fixed at the distal diaphysis and loaded with different biomechanical forces for topology optimization. Boundary conditions due to possible additive manufacturing processes were incorporated into the optimization to ensure manufacturability. The optimization result is compared with experimental data from a human bone. A bone-like internal structure and increased compliance of the topology-optimized test bone model compared to the commercial model were observed.


Introduction
The examination of bone plates, also called osteosynthesis plates, is performed with biomechanical tests using test bones. The availability of suitable human bones is limited, especially for bones from young-to medium-aged donors or particular disease types. In addition, the use of human test bones is associated with high costs and complex preservations along with high fluctuations of their mechanical properties. Due to the disadvantages of natural test bones, artificial test bones are increasingly being used in the orthopedic biomechanics research and surgical training. These artificial bones, mostly composite bones, are undergoing steady improvement. They are less expensive than cadaveric femora, are sufficiently available, and are ethically justifiable [1,2]. However, their availability is currently limited to a few sizes and geometries.

Conventional approach
The market leader in supplying artificial test bones is the company Sawbones, a division of Pacific Research Laboratories, Inc. (Pacific Research Laboratories, Inc., Vashon, USA) [3]. The so-called Sawbones, named as the company itself, consist of a polyurethane foam, intended to represent a cancellous bone, which is encased in epoxy resin to replicate the cortical bone. Short-glass-fiber-reinforced (SGFR) epoxy, used since the third generation of Sawbones, is intended to have similar mechanical properties as the cortical bone. The epoxy component has been optimized for the current fourth generation of Sawbones [1]. Former studies concluded that artificial bones represent the femora of healthy adult men below 80 years of age [4]. They are primarily available in three sizes: large, medium and small [3].
Heiner [6] investigated the mechanical properties of the fourth-generation artificial femora and tibiae. An overview of her findings, including the Young's modulus, the tensile strength, and the compressive strength, as well as similar values for the natural femur, is given in Table 1. However, the comparison was not investigated in a biomechanical test. Gardner et al. [4] provided comparable measurement data for fourth-generation Sawbones of size large. The results of the mechanical tests of the artificial bones are in the range of their natural counterparts, which themselves show a wide range. Basso et al. [9] observed a high stability in synthetic bones of the fourth-generation, which was not in accordance with human bones. They also observed differences in the fracture behavior.
Despite the availability of synthetic bones with mechanical properties comparable to those of natural bones, these bone replicas cannot reproduce the anisotropy of the cancellous bone of their natural counterparts. This is due to the polyurethane foam with isotropic properties [3,10]. As a result, the polyurethane foam can only reproduce the mechanical properties of the natural cancellous bone to a limited extent. Such structures, which replicate the anisotropic behavior, can be created by additive manufacturing (AM) combined with topology optimization.

Topology optimization of bone-like structures
According to 'Wolff's law', the bone is able to adapt to new boundary conditions and remodel itself to an optimal internal structure relative to its mechanical requirements [11].
Using topology optimization, optimal structures can also be found with respect to a defined objective function. When combined with AM, complex components and even individualized geometries can be produced cost-effectively.
An approach based on the strain energy density function is one method of topology optimization. Huiskes et al. [11] used this approach to simulate bone-remodeling. According to their theory, the bone density is increased in regions where a certain threshold of the strain energy is reached and is decreased when it is below a certain threshold. Jang et al. [12] compared this approach with the compliance-based topology optimization. In two case studies, they obtained similar results for both algorithms. They concluded, that a bone is capable of self-optimization and that topology optimization also has the potential for the generation of artificial bone structures. Park et al. [13] used a topology approach with a perimeter control to generate the internal bone structures. The perimeter can be used to constrain the design space and to enable the generation of complex structures. One-leg stance, abduction, and adduction were chosen as load cases on the femur. By randomly selecting suitable initial conditions, trabecular structures could be generated.
Another algorithm for structural optimization is the design space optimization. Here, the material is iteratively distributed into areas of high stresses while simultaneously varying the design space. Structures are thereby given minimum compliance and thus maximum stiffness. This algorithm was used by Boyle and Kim [14] to verify Wolff's law using the femoral head as an example. For loading, they used the data from Bergmann et al. [15] and Heller et al. [16] of stair climbing and walking. Their results show trabecular structures similar to those of the human femur. Wu et al. [17] introduced a local density upper bound to avoid larger material accumulations and to achieve more evenly distributed structures. The density around each voxel must not exceed an upper limit. This prevents regions where a large amount of material accumulates and thus favors the formation of fine structures, such as those found in bones.

Approach of this paper
This study aims to develop a methodology for creating artificial test bones. It utilizes the approach of topology optimization for the design of additively manufacturable artificial test bones, which are intended to represent the anisotropic behavior of human bones. The outer geometry of a bone is used as the input data. The internal structure of the bone is generated in a next step with the use of the topology optimization. By defining boundary conditions, the manufacturability by additive manufacturing is ensured. With this approach, it is possible to manipulate bone geometries in advance and use the results of that as initial geometries for the topology optimization. This allows the creation of a wide variety of individual bone geometries for the respective biomechanical test. Moreover, for the first time, the developed methodology allows the comparison between results from different topology optimization approaches and experimental data of human bone.

Methodology
The methodology used in this study for generating the biomechanical test bone through topology optimization is illustrated in Fig. 1. The geometry of the bone was derived from the computed tomography (CT) data of Kluess et al. [18] using segmentation and reverse engineering. The advantage of this data is that a round-robin study was already conducted to compare simulated and experimental results [18]. The model was divided into the cortical and the cancellous bone. Based on this geometry, a simulation model was created for topology optimization. A variation of the optimization boundary conditions was performed to investigate the internal structures (Sect. 2.4). All topology optimization results were simulated in a reanalysis step. The reanalysis step must be performed because the topology optimization solution is a density distribution. However, for a realistic approach, only values of 0 or 1 can be assumed. Finally, the resulting deformations of the femoral head were compared to the experimental data of Kluess et al. [18].

CT data and geometry reconstruction
The CT data from the study by Kluess et al. [18] were used as the basis for generating the geometry of the femur. The data were processed with ImFusion Suite 2020 (ImFusion GmbH, Munich, Germany). The segmentation was performed semi-automatically using the interactive watershed algorithm. Initial regions of the background, the cortical bone, and the cancellous bone were marked. The algorithm is used to perform automatic segmentation of areas with the same gray value within a certain tolerance. The outer surfaces can be exported as a stl file. The software Geomagic Design X 2019 (3D Systems Corporation, Rock Hill, USA) was used to prepare the surface data for simulation. Since a stl file is not suitable for meshing, the surface was approximated by patches using non-uniform rational basic spline surfaces. Once this step was completed, the volume bodies of the cortical bone and cancellous bone could be generated.

Material selection
Because of the good availability and its high design freedom, laser-based powder bed fusion of polymers was chosen as a suitable AM technology. Material jetting or vat photopolymerization are two techniques that could also be used but are not considered further in this paper due to the high production costs compared to laser-based powder bed fusion of polymers. Due to high deformations with the use of polyamide 12 in preliminary tests, Alumide® (from EOS GmbH, Krailling, Germany) material-a mixture of polyamide 12 and aluminum granules [19]-was used to create the internal structures. An overview of the used materials is shown in Table 2. After the printing process, the internal structure is casted with an outer shell of epoxy resin. The Young's modulus of the epoxy resin representing the artificial cortical bone can be adjusted by adding short glass fibers. The optimizations were performed with a Young's modulus of 15 GPa. All materials are considered to be isotropic for the optimization step. Nonlinearities of the manufacturing process can be implemented in the nonlinear reanalysis. Within the scope  of this work, no predefined building direction was known. Therefore, the anisotropy of the building process cannot be considered.

Modeling
The finite element analyses and the topology optimizations were performed using the software package from Altair Engineering 2021 (Altair Engineering Inc., Troy, USA). This includes the preprocessing in Hypermesh, the utilization of Optistruct as a solver, and the post-processing in Hyperview. The finite element mesh was created with a target element size of 1 mm. A mesh independence study was performed, but the average element size was chosen lower to allow the generation of fine structures in the topology optimization. It includes only tetra-elements of the first order, because the solution does not change significantly with the use of second order elements. The minimum element size was set to 0.8 mm and the growth rate to 1.0 to ensure the creation of uniform elements throughout the entire bodies. To validate the model, a test run was created using the material properties of the natural bone described by Wirtz et al. [10]. The validation of the simulation model was performed by the evaluation of the mean deviations of the deformations compared to the results of the round-robin study by Kluess et al. [18]. The comparison of the mean deviations to the experimental data is shown in Fig. 2. After the validation, the materials were given properties suitable for AM. In Fig. 3, the model setup of the topology optimizations is shown in a simplified way with the points of force application (green and orange) and the fixing boundary condition (red). Five load cases were defined to represent the daily loading of the bone in the human body. For this purpose, a combination of the most common loads on the bone throughout the day was collected from current literature. Two linear static load cases for the topology optimizations were implemented according to Heller et al. [16] (see Tables 3 and  4) with a body weight (BW) of 847 N. For further analysis, the load cases one-legged stance, abduction, and adduction, that were used in former studies, were implemented as well (see Tables 5, 6 and 7) [20][21][22]. In the reanalysis step, an additional geometric nonlinear quasi-static load case was used, which was defined according to the test conditions of Kluess et al. [18]. A geometric nonlinear load case was calculated as large deformations occur. All optimization results are reanalyzed with a density threshold of 0.5 and compared to the experimental data of Kluess et al. [18]. This means, that relative density values below 0.5 are truncated, and densities above 0.5 are considered solid material.
The forces were applied to the nodes by means of rigid body elements (RBE) of the type RBE3. These types of RBEs only distribute acting forces. The nodes of the attachment surface thus remain movable relative to each other. Consequently, the application of the RBE of this type does not induce artificial stiffness [23]. The femur was constrained with a fixed boundary condition up to 70 mm above its distal end analogue to the study of Kluess et al. [18], which is comparable to a distal diaphyseal fixation in the current literature [24][25][26][27].

Topology optimization
The topology optimization was performed with the objective of compliance minimization for the two load cases and afterward with the additional three load cases. This  approach has the advantage of not requiring previously determined experimental data to create internal structures of the artificial bone. Optistruct uses the Solid Isotropic Material with Penalization (SIMP) method to formulate the optimization problem. As a constraint for the optimization, the parameters controlling the geometry of the structures to be generated were defined in a design variable in Hypermesh. These were incorporated to ensure the manufacturability of the solution. The parameters mindim, maxdim, and mingap control and penalize the minimum and maximum diameters of the structures as well as the minimum gap between them [28]. A mindim value of 3 mm was chosen, which is three times the average mesh size. The maxdim and the mingap parameters were set to 6 mm. It was found in previous studies, that a discretization factor of six gives suitable results and was, therefore, used in this study. To prevent the formation of checkerboard patterns, the checkerboard control option was enabled [29].  Table 4 Applied forces of the load case stair climbing in percent of body weight from Heller et al. [16] and points of application (see Fig. 3)  Table 5 Applied forces of the load case one-leg stance from Beaupré et al. [22] and points of application (see Fig. 3)  Table 6 Applied forces of the load case abduction from Beaupré et al. [22] and points of application (see Fig. 3)  Table 7 Applied forces of the load case adduction from Beaupré et al. [22] and points of application (see Fig. 3) To explore further possibilities for improving the topology optimizations, variations of the boundary conditions were conducted. An increase in the penalty parameter for the SIMP method was performed with values of 8 and 10. A further optimization was modeled with a limitation of the maximum displacement of the force application point in the corresponding load case according to the experimental study of Kluess et al. [18]. This constraint was also coupled with a minimization of the volume as the objective function.

Printable geometry
The geometries of the reanalysis can be used for AM. The resulting bodies can be exported as stl files and must subsequently go through the steps of building job preparation.

Results and discussion
In all optimizations, the medullary canal is mapped in the design space. The more material the optimization algorithm can distribute in the design space controlled by the volume fraction boundary condition, the more material accumulates at the caput femoris, in the medullary canal (especially distally), and at the greater trochanter (see Fig. 4 marked in red). Moreover, no large spaces without material form in the proximal region of the femur above a volume fraction of 0.35. Figure 5 shows the resulting relative density in the design space after the optimization for the changed boundary conditions. Increasing the SIMP parameter to a value of 8 shows approximately the same material distribution as the initial value of 6. A further increase to a value of 10 leads to a deterioration of the convergence to discrete 0 and 1 values. This is evident from higher density values in the medullary canal. According to the load case of Kluess et al. [18], the introduced limitation of the displacement does not lead to a more prominent differentiation of the density distribution. If this boundary condition is coupled with a minimization of the volume, a strongly deviating density distribution is shown (see Fig. 5d). Here, the generated structures were distributed over the entire design space. In all topology optimizations, a small amount of material accumulated proximal to the caput femoris at low values of the volume fraction boundary condition. However, this region has a high density in the CT data used in this work. According to the analogy of human bone self-optimization and topology optimization, the material should also form in this area by the objective function of compliance minimization.
For further analysis of the internal structures, cross-section cuts of the different geometries are shown in Fig. 6. It can be seen that the internal structures show only minor adjustments when optimized with a volume boundary Fig. 4 Results of the topology optimization with a varying volume fraction with the proximal region marked in red Fig. 5 Results of the topology optimization with varying boundary conditions condition and the minimum compliance objective function (see Fig. 6a). The optimizations with the modified objective function or load cases show considerable bone-like orientations of the internal structures (see Fig. 6b, c). Furthermore, the combination of all load cases leads to a less pronounced alignment of the internal structures (see Fig. 6d).
The displacement values of the reanalysis at the force application point of the femoral head under load for different volume fraction constraints can be seen in Fig. 7. Even when solid material is used in the cancellous bone (volume fraction = 1), the horizontal displacement is larger than in the measurements performed by Kluess et al. [18]. In turn, the vertical displacement is too small even with a completely hollow bone compared to the measured values. For classification purposes, it should be noted that the deviations of the laboratories participating in the round-robin procedure are significantly higher in some cases. Three of the seven laboratories obtained deviations of more than 120% on average [18]. The results of this study can be assigned to the more accurate simulations with average deviations of less than 40% [18]. With increasing volume fractions, the displacements decrease in both directions. The resulting displacement values for the study of different optimization boundary conditions are shown in Fig. 8. For a better comparability, the result of the volume fraction condition of 0.35 is also displayed. The lowest displacement values were achieved with a combination of the displacement constraint, the objective of volume minimization, and the adjustment of the load case. A higher penalization factor increases the displacements, leading to higher deviations from the experimental   [18] data. Implementing the displacement constraint without the change of the objective shows no significant advantages compared to the standard formulation.
However, this boundary condition does not consider the natural formation approach of a bone. Therefore, the optimizations represent the self-optimization mechanism of the bone less accurately than those modeled solely based on the load cases [16,[20][21][22]. The displacement boundary condition was additionally non-active in every case and thus cannot be used for an improved displacement behavior. This is due to the segmented geometry showing larger displacements in the reanalysis, even with full material [18]. Such a boundary condition also requires data that can only be obtained in an experimental setup. If a bone from a living patient needs to be modeled, these data cannot be collected. Either approximated data from comparable experiments must be used or a different optimization modeling must be performed. As a result, the generated structure has a lower stiffness.
All topology optimizations were performed with distal fixation according to the mechanical tests of the study by Kluess et al. [18]. This boundary condition still does not represent the real mechanical environment of the femur in the human body. Rossmann et al. [30] and Speirs et al. [31] showed that the fixing boundary conditions significantly influence the resulting mechanical behavior. It can be assumed, that these boundary conditions also affect the topology optimizations. With a better biomechanical representation of the model, the topology optimization could reproduce the self-optimization of the human bone even better. The use of alternative load cases or the integration of patients' gait analysis data also offers room for improving the methodology.

Conclusions and outlook
To enable the development of individualized implants, they must be subjected to mechanical testing. An approach to design individualized artificial test bones with bone-like structures was pursued in this work. First, topology optimization was performed to replicate the trabecular structure of a cancellous bone. A geometry was derived from CT data and different optimization models were set up. These contain load cases that reproduce the actual loads on the bone occurring in the human body and boundary conditions to generate the fine structures. All results of the topology optimizations were subsequently verified in the reanalysis and compared to the experimental data of Kluess et al. [18]. The best results in terms of the alignment of the internal structures and the resulting displacements were obtained using the volume boundary condition approaches in conjunction with the displacement boundary condition and the modified load cases.
In the future, experimental verification of the relationship between the generated structures and the Young's modulus is of great interest. The anisotropy of the manufacturing process will be included in the reanalysis for further examination. If the processes are sufficiently refined and validated, the findings can be applied to the design of artificial bones. Besides the femora, there are many other relevant bones in the human body. The transferability of the topology optimization approach needs to be confirmed. In particular, the transfer to bones with similar properties, such as the tibia, is of high interest.
Funding Open Access funding enabled and organized by Projekt DEAL. The presented study was part of the research project "FOMIPU-Implant and minimally invasive surgical technique for repositioning osteotomy in children "(13GW0293F), which was funded by the Federal Ministry of Education and Research (BMBF). We thank the BMBF for its excellent support.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.