An improved trabecular bone model based on Voronoi tessellation

Background and objective: Accurate numerical and physical models of trabecular bone, correctly representing its complexity and variability, could be highly advantageous in the development of e.g. new bone-anchored implants due to the limited availability of real bone. Several Voronoi tessellation-based porous models have been reported in the literature, attempting to mimic the trabecular bone. However, these models have been limited to lattice rod-like structures, which are only structurally representative of very high-porosity trabecular bone. The objective of this study was to provide an improved model, more representative of trabecular bone of different porosity. Methods: Boolean operations were utilized to merge scaled Voronoi cells, thereby introducing different structural patterns, controlling porosity and to some extent anisotropy. The mechanical properties of the structures were evaluated using analytical estimations, numerical simulations, and experimental compression tests of 3D-printed versions of the structures. The capacity of the developed models to represent trabecular bone was assessed by comparing some key geometric features with trabecular bone characterized in previous studies. Results: The models gave the possibility to provide pore interconnectivity at relatively low porosities as well as both plate-and rod-like structures. The mechanical properties of the generated models were predictable with numerical simulations as well as an analytical approach. The permeability was found to be better than Sawbones at the same porosity. The models also showed the capability of matching e.g. some vertebral structures for key geometric features. Conclusions: An improved numerical model for mimicking trabecular bone structures was successfully developed using Voronoi tessellation and Boolean operations. This is expected to benefit both computational and experimental studies by providing a more diverse and representative structure of trabecular bone.


Introduction
Low-density, porous structures are present broadly in both natural and engineered materials, e.g.trabecular bone, wood, implants, solid foams, or honeycomb structures, to name a few (Gibson and Ashby, 1999).However, the mechanical behavior of these structures is commonly difficult to investigate analytically due to their complex and typically irregular geometries, as is the case of trabecular bone (Fig. 1).Experimental methods are therefore typically used to test the mechanical behavior of these structures, but in the case of trabecular bone the highly variable structure makes it difficult to cover a full populationstructural variations are found between anatomical sites, individuals, as a function of pathology etc. (Rincón-Kohli and Zysset, 2009), and covering the whole spectra through experimental testing is virtually impossible.This can lead to issues when developing new implant designs for example.One specific case is the mechanical evaluation of new screw designs, which is commonlyand according to the ASTM (American Society for Testing and Materials International) standard (ASTM F543-07, 2007) -done using randomly structured polyurethane foams (e.g.Sawbones (Poukalova et al., 2010)).However, this choice of material may lead to different results than testing in real bonee.g., insertion angle has been found to have an effect on screw pullout strength in real bone, implying that the anisotropy found therein is of non-negligible importance (Wilmes et al., 2008).Therefore, lately, 3D-printing has been investigated as a tool to reproduce the trabecular bone structure for implant testing (Wu et al., 2019;Grzeszczak et al., 2021).However, the resolution is not yet satisfactory to this end (Grzeszczak et al., 2021).Also, experimental evaluation would still require testing many different structures, and the tested samples would still consist in examples only, i.e. there is no certainty that a representative selection of structures are tested.Numerical simulations could provide more control as well as a higher time-and cost-efficiency.
Many porous structure designs have been reported as potential representations of trabecular bone (Kou and Tan, 2010;Fantini et al., 2016;Gómez et al., 2016;Cheah et al., 2003;Belda et al., 2023;Bucklen et al., 2008;Lee et al., 2020;Köll and Hallström, 2014;Li et al., 2012;Chen et al., 2021).Several are based on a periodical appearance of unit cells, which can be analysed with both analytical and numerical methods (Cheah et al., 2003;Belda et al., 2023;Bucklen et al., 2008).However, as mentioned above, the mechanical properties of the bone and its mechanical behavior in combination with implants, e.g.resistance to screw pull-out, are sensitive to the testing region chosen due to local variations in the structure (Bennani Kamane, 2012;Roshan-Ghias et al., 2011).It would therefore be valuable to achieve a representation of the porous structure which is closer to reality.
Recently, several stochasticity-based approaches have been reported for the generation of porous structures that could be utilized to represent a trabecular bone structure, many of them based on Voronoi tessellation techniques (Kou and Tan, 2010;Fantini et al., 2016;Gómez et al., 2016;Köll and Hallström, 2014;Li et al., 2012;Chen et al., 2021;Frayssinet et al., 2022).However, the generated structures contain either entirely isolated (closed) pores (Köll and Hallström, 2014) or consist of open-cell structures of truss frameworks (Fantini et al., 2016;Gómez et al., 2016;Li et al., 2012;Chen et al., 2021;Frayssinet et al., 2022).The truss frameworks used are limited to high porosity models in order to adequately simulate the rod-like trabeculae, which is also stated by (Li et al., 2022).Indeed, the rod thickness is in these models changed to control the porosity, i.e. lower porosity structures require rods much thicker than real trabeculae (Fantini et al., 2016;Gómez et al., 2016;Wang et al., 2018).However, at a lower porosity, more plate-like structures are typically present in trabecular structures (Wang et al., 2012), structures which present a smooth transformation between rod-like and plate-like trabeculae, and this cannot be found in the currently available Voronoi tessellation-based algorithms in the literature (Wang et al., 2018).Some studies have simplified the edges/walls to structural rods, beams, or plates (Wang et al., 2012;Liu et al., 2009).However, none of these assumptions correspond to natural trabecular bone structures, which consist of a mixture of rod-and plate-like struts with smooth transfers and open pores (Hildebrand and Rüegsegger, 1997).
In this study, a porous structure-generating algorithm based on Voronoi tessellation was introduced, with the aim of being able to provide more realistic representations of trabecular bone structures for a larger range of porosities, and for different anatomical sites.We hypothesized that by using the surfaces of Voronoi tessellation rather than the edges, a smooth transition from plate-like at lower porosity, to rodlike at higher porosity could be achieved.In this way, the structure could be morphologically more similar to bone for a larger porosity range.The morphometrical parameters of the generated structures were compared to human trabecular bone data reported in the literature.The performance in terms of representing the trabecular bone from a structural integrity perspective was evaluated through numerical and experimental uniaxial compression tests on 3D-printed structures.

Theory and algorithm
The porous structure generation algorithm is summarized in the flowchart in Fig. 2. First, a number of seeds are used as input to generate a corresponding Voronoi tessellation.The Voronoi cells are then prescaled to separate the design space and subdivided to acquire smooth individual polyhedrons.Then, the subdivided polyhedrons are scaled, or tuned, to acquire a target porosity and intersections between different polyhedrons are assessed.If intersections are found, a Boolean operation is utilized to create new surfaces and to merge the intersected polyhedrons.Due to the intersection of polyhedrons, the porosity is not simply computable by doing a summation of the individual Voronoi cell volumes without volumetric meshing.To acquire the actual porosity, the structure is discretized using voxels and the porosity is evaluated by computing the voxel volumes.If a target porosity is not achieved, the Voronoi cells are tuned and merged again in an iterative process until the target porosity is achieved.The procedure is described in more detail in the following paragraphs.
A Voronoi diagram, or a so-called tessellation, is a partitioning of a field based on N points in Euclidean space.The first step in generating a Voronoi diagram is to create initiation points (seeds) in a Cartesian reference frame (x, y, z).For example, the seeds can be positioned by uniformly distributed random numbers (Fig. 3a).If a homogenous structure is required, the points are redistributed, by minimizing the areal differences between the cells (Fig. 3b).This is achieved by iteratively recalculating the seeds as the central point of the cells in the last  where S x,y,z,m,n+1 is the Cartesian coordinate of seed m in iteration n + 1, and the C x,y,z,m,n is the centroid of polyhedron m in iteration n.The centroid is calculated by discretizing the convex hull to several tetrahedrons.The seeding and tessellation strategy can also be applied to construct anisotropic structures to mimic natural trabecular structures according to e.g.Wolff's law (Frost, 1994), as shown in Fig. 3c.When generating anisotropic structures, the probability density function (PDF) is set differently in different directions.
After this, N spheres around the predefined seeds with vanishing small radii are created.The radii are then increased simultaneously until neighbour spheres come into contact (Kou and Tan, 2010).Thereby, the boundaries create surfaces between different seeds.Mathematically it can be defined by a set of points with (Okabe et al., 2000): where, V(p m ) is the defined volume for cell m, d(p, p m ) is the Euclidean distance between point p and the seed p m of cell m.After traversing all the seeds, an isolated-cell structure is created.Some three-dimensional prototypes with different seed distributions are shown in Fig. 4. The generated polyhedrons were initially pre-scaled (Fig. 5b).The pre-scaling was only utilized to facilitate the identification of the elements for each Voronoi cell and give a clear definition of the porous structure before the subdivision.The reference points of the pre-scaling are the centroids of each cell, and all the Euclidean distances between all vertices in a specific cell and the centroid are scaled with the same scale factor f. A cube was selected as the design domain for convenience.
Every surface that exceeded the domain was cut by the boundary plane, i.e. the boundary surfaces were plane.
Until here a porous structure based on Voronoi tessellation is created.Then, the polyhedrons from the Voronoi diagram were subdivided to create smoother surfaces as shown in Fig. 5c.There are different subdivision techniques to generate smooth surfaces from a polyhedron, e.g. the Catmull-Clark scheme, Doo-Sabin, Loop (Hakenberg et al., 2014).Here, a B-spline curve continuity-based algorithm was utilized (Muir, 2020).
After the pre-scaling and subdivision processes, a porous structure was created with smooth cells.In the next step, the porosity was controlled by scaling all individual cells with the same algorithm used for pre-scaling.With the scale factor f, the porosity φ can be controlled with before the first contact between two arbitrary spheres, where V is the volume of the whole design domain, and V m is the volume of the cell m.
The scale factor f is computed as, where φ is the predefined porosity and φ is the present porosity of the structure, which is calculated by summing the volume of every cell.The upper limit of the porosity for a non-intersection condition between voids is around 55%-65% for different seed distributions -the sphere packing theory gives an average 64% for random packing (Conway and Sloane, 2013).Above that limit, the boundaries of the polyhedrons will start to intersect each other.The intersection of polyhedrons is achieved by a Boolean operation (Schneider and Eberly, 2002;Jacobson, 2018).
The intersection can provide permeability to the structure, and the permeability can be increased as the porosity increases.This is an additional parameter that is very important to the correct representation of a trabecular structure, e.g. for the evaluation of implant augmentation materials such as bone cements.
The overlapping of the polyhedrons also results in the inconsistency of Eq. ( 3).Specifically, the summed value of each V m does not account for the entire space occupied by the polyhedrons.To determine the actual porosity of the porous model, the structure is discretized (using voxels) (Amanatides and Woo, 1987;Min, 2021).A voxelized structure is shown in Fig. 6b, corresponding to the discretized structure in Fig. 6a.After that, the porosity can be simply calculated by counting the volume occupied by the voxels.If the calculated volume does not match the targeted porosity, a modification of the scale factor f is applied as introduced in Eq. ( 4), to the original structure, i.e. the triangular meshed structure.This procedure is repeated until the desired porosity of the structure is obtained.

Evaluation of the generated structures 2.2.1. Numerical simulation
To evaluate the influence of different model parameters on the mechanical properties, several models were numerically analysed in the commercial finite element code Abaqus (ABAQUS, 2016; Vélizy-Villacoublay, France) (Hibbitt and Sorensen, 1992).The structural response under two simple loading cases was evaluated (Fig. 7).
To evaluate the influence of porosity, nine models were generated with the same number of seeds, i.e. 100 seeds, but with a varying porosity, from 10% to 90%.To test the influence of the number of seeds, five models were generated with the same porosity, but with a varying number of seeds from 5 to 150.To evaluate the anisotropy, an anisotropic model was generated with a density ratio of seeds 0.5:1:2 in the x-, y-, z-directions.The anisotropic model had a porosity of 80%, with 100 seeds.The approximate edge size of tetrahedral elements in these models was 1% of the edge length of the cubic design domain (10 mm).
The number of elements to be used was determined by a sensitivity analysis of the effect of element size on the structure stiffness.The element type used was standard quadratic tetrahedron elements (C3D10).The material properties used were those of the material to be printed, i.e. polylactic acid (PLA), with a Young's modulus E = 2 GPa and a Poisson's ratio of 0.3, as assessed by uniaxial compression tests on homogeneous specimens in a testing machine (giving E = 1.98 ± 0.14 Gpa).The material was assumed to be isotropic and elastic, and small deformations were assumed to prevail.
A compression test was simulated using a displacement of 0.1 mm in the vertical z-direction on the top surface, while the bottom surface was locked in the z-direction, as indicated in Fig. 7a.A shear test was simulated using a 0.1 mm displacement on the top surface in a horizontal (x-or y-) direction with the bottom surface fixed, Fig. 7b.

Validation -experimental testing
2.2.2.1.Additive manufacturing.Additively manufactured, or 3Dprinted, structures were used to verify the validity of these types of models in terms of mechanical properties.The STL files were used as   input into the software Cura (Braam, 2022) and printed with Fused Filament Fabrication (FFF) (Ultimaker S5, Ultimaker BV, Utrecht, Netherlands).The structure was scaled to a 2 × 2 × 2 cm 3 cube to improve the printing quality and printed with 0.25 mm nozzle size and 0.1 mm layer height.A print rate of 10 mm/s was utilized to print PLA filament (3D4Makers B.V, Haarlem, Netherlands) at a printing temperature of 170 • C. The infill density was 100% and the infill pattern was "Grid".No support structure was used.An example of a printed structure is shown in Fig. 8.

Mechanical testing.
The printed structures were tested using a standard laboratory materials testing machine (AGS-X, Shimadzu, Japan) at a loading rate of 1 mm/min.A lubricant was applied to the rigid plates between the structure and the testing machine to reduce any influence of friction.A load cell with a capacity of 5 kN was used.The tests were conducted either until a displacement of 5 mm, i.e. 25% of the specimen height, was reached, or until 4.8 kN.The stiffness was calculated from the initial linear part of the measured stress-strain curve, where the stress was calculated by averaging the reaction force on the structural top surface.In each group, 6 specimens were printed and tested.The mean value for each porosity level was used for comparison with analytical and numerical results.Two force-displacement curves are shown to illustrate the typical experimental processes (Fig. 9), where the elastic regions were picked to calculate the stiffness.

Verification -analytical analysis
To be able to compare the numerical, experimental, and theoretical stiffness of the same porous structure, the global Young's modulus of the porous material was estimated considering the porosity and the anisotropy, as given by the simplified expression where E s is the Young's modulus of the solid material, φ is the average porosity in the whole domain.The ρ i is the density ratio of seeds in direction i for anisotropic structure, and E i is the Young's modulus of the porous material in direction i = x, y, z in the Cartesian system.The density ratio of seeds of x-, y-, z-direction is normalized so that ρ x ρ y ρ z = 1.(Rueden et al., 2017) were used to determine the ratio between bone volume and total volume (BV/TV), degree of anisotropy (DoA), trabecular thickness in μm (Tb.Th), trabecular separation in μm (Tb.Sp), trabecular number in 1/μm (Tb.N), and ellipsoid factor (EF) (Frayssinet et al., 2022) to enable comparison with other porous structures from literature (Bouxsein et al., 2010).The estimation of Tb.Th and Tb.Sp is shown in Fig. 1, i.e. the Tb.Th and Trabecular separation (Tb.Sp = 1/ Tb.N − Tb.Th).The DoA was estimated with the mean intercept length (MIL) method (Harrigan and Mann, 1984).The closed cell porosity was measured to estimate the permeability of the models from this data.

Permeability
A blood mimicking fluid (BMF) was utilized to investigate the longitudinal permeability, similar to a previous study (Diez-Escudero et al., 2020).The BMF was prepared with water (39%), glycerol (58%), and dextran (3.3%), resulting in a density of 1168 kg/m 3 and dynamic viscosity of 4.69 mPa⋅s (as measured on a Discovery Hybrid Rheometer 2 (TA Instruments, Sollentuna, Sweden)).Beakers were utilized to contain the BMF and maintain the liquid height during the tests.A BMF height was set at 26.6 mm to provide the same pressure as in the previous study, i.e. 305 Pa (Diez-Escudero et al., 2020).A hole was drilled in the bottom of the beaker to fit the 3D printed scaffolds.The models used for the permeability test were the same as the ones utilized for mechanical testing.BMF was flowed through the bone models (porosity 60%, 70%, 80%, 90%) and Sawbone structures (porosity 70%, 85%, 89%, 92%).To ensure only longitudinal flow of the BMF the scaffolds were laterally sealed with parafilm.The gaps between the beaker hole and the scaffolds were also filled with parafilm.The bottom of the scaffold was covered with parafilm until the liquid level reached 26.6 mm, then the bottom parafilm was removed and the BMF was continuously added to maintain the liquid level, and the time was recorded.The time was recorded until 300 mL BMF had passed through to calculate the fluid flow rate.The longitudinal permeability was calculated by K = μqL z /[AΔp], where K is the longitudinal permeability, μ is the BMF viscosity (4.69 mPa⋅s), q is the fluid flow rate, L z is the height of the model, A is the cross-section area of the models (4 cm 2 ), and Δp is the pressure resulting from the BMF (305 Pa).

Porous models
To achieve different porosities, the volumes of individual cells were scaled as described above.Four structures with the same seeds but different porosities are shown in Fig. 10 to illustrate the influence of porosity on the topology of the structure.As the porosity increases, more intersections occur.When the porosity is low, the smooth polyhedrons are small and isolated.At higher porosities, the individual cells are intersected, and as the porosity reaches 80%-90%, the individual cells are hard to discern.The shape of the struts turns from plate-like to rodlike.
Another important factor controlled in the structures is the number of seeds.Four structures with different number of seeds are plotted in Fig. 11.The number of seeds varies from 20 to 1000, and the porosity remains constant at 80%.As the number of seeds increases, the volume of one individual cell decreases and more struts are generated.When the number of seeds is low, the influence of every cell on the whole structure is high.On the other hand, a structure with 1000 cells is more homogenous, and the influence of one individual cell on a global mechanical property is naturally lower.For human vertebrae (Arlot et al., 2008), we can estimate the number of trabecular struts in a 1 cm 3 cube.If we aim to build a trabecular bone model for a 1 cm 3 cube, and assume there are 6-12 trabeculae to construct a cell, around 250 to 500 cells are required to construct the 1 cm 3 trabecular structure cube.
The distribution of the seeds is another important factor that influences the properties of the model.Fig. 12 shows three structures generated with different distributions of seeds.The model (b) with uniformly distributed cells is more isotropic and homogenous compared to the model with a random distribution of seeds with uniform PDF (a).Different distributions can be used for trabecular bone from different anatomical sites, e.g. the trabecular structure in human vertebrae is more homogenous and isotropic than the femoral part, and a model with a uniformly distributed PDF could be used to mimic the vertebral structure.

Morphology
The morphological measurements from the micro-CT analysis of two of the bone models are shown in Table 1.

Stiffness evaluation and comparisonverification and validation
The volume differences between the 3D-printed structure and the original model were less than 5% for all models, as evaluated by the micro-CT analysis described above.The global stiffnesses of different structures with different porosities are shown in Fig. 14a, contrasted

Table 1
Morphological parameters of our porous model.DoA is the degree of anisotropy; Tb.Th is the trabecular thickness; Tb.Sp is the trabecular separation (Bouxsein et al., 2010 with different observation methods, i.e. numerical simulation, experiments, and the analytical expression, Eqn. ( 5).Here, the porosities are equal in all three material directions, i.e. the structures are isotropic on a global scale.The experimental data is naturally variable, due to e.g.print quality, particularly notable at lower porosities.As one can observe, the three methods give very similar results.The maximum relative difference between the simplified analytical estimation and the numerical simulation was approximately 10% (Fig. 13a).In Fig. 14b, the influence of the number of seeds on the compressive and shear stiffness is illustrated.When the number of seeds is low, there is a stiffness oscillation due to the stochastic positions of the seeds caused by the limitation of the random number generator.As the number of seeds increases, both compressive and shear stiffness are more stable and showed a convergence after approximately 75 seeds.In Fig. 14c, the influence of structure anisotropy on the compressive stiffness in different directions is shown.The stiffness in the y-direction did not change compared with Fig. 14b.The stiffness in the x-direction increased to around three times compared to the y-direction.The stiffness in the z-direction decreased to around 10% compared to the y-direction.The differences between numerical and experimental results were 5% in the x-direction and around 15% in the other two directions.
The differences between the analytical estimations and the results from the experiments and simulations in x-and y-directions were around 3% and 15%, however much higher in the z-direction, about five times.

Permeability
The longitudinal permeability measurement results are shown in Fig. 15.The permeability of the 3D printed bone models was typically higher than the Sawbones scaffold for the same porosity, with the exception of the 90% porosity.The bone model showed higher permeability at 60% porosity than the Sawbones model at 70% porosity.

Discussion
An algorithm for generating porous structures based on Voronoi tessellation was introduced, with an aim to better represent trabecular bone structures than currently available models.Compared to many other Voronoi diagram-based algorithms (Kou and Tan, 2010;Fantini et al., 2016;Gómez et al., 2016;Köll and Hallström, 2014;Li et al., 2012;Chen et al., 2021), this algorithm generates porous structures with partially intersected smooth polyhedrons rather than scaffolds consisting of rods only.The substantial difference of these previously reported, truss-like structures to human trabecular bone has been highlighted repeatedly in literature (Li et al., 2022).The algorithm in this study generates structures with mainly completely interconnected pores, which is also the case of trabecular bone (Gibson and Ashby, 1999).It subdivides the surfaces to a predefined smoothness, rather than introducing straight edges and planes (Köll and Hallström, 2014).The final structure is determined by the positions of the original seeds and the porosity.Furthermore, we can achieve a global target porosity of the structure, rather than using additional approximate parameters (Fantini et al., 2016).
To evaluate the degree of similarity between our model and real trabecular structures, three key morphology parameters were compared with rabbit trabecular bone and human trabecular bone as a function of porosity, namely the Tb.Th, the Tb.N, and the DoA (Fig. 16).Porosity is a design parameter in our algorithm and is hence not an indicator.The human trabecular bone porosity has been reported to range from 60% to more than 90%, therefore, we also created models starting from 60% (Ding et al., 1999).The DoA can represent the structural anisotropy within a volume by comparing the mean intercept length in different directions (Odgaard, 1997).Trabecular bone data collected from five previous studies (Rincón-Kohli and Zysset, 2009;Arlot et al., 2008;Mueller et al., 2009;Fields et al., 2009;Chen et al., 2010;Turunen et al., 2013) for different parts of human trabecular bone, i.e., radius, vertebra, and the femoral neck were used for comparisons (Fig. 16).One group of rabbit femur trabecular bone (38 samples), and 8 models generated with the algorithm introduced in this study were plotted in the same figure (Fig. 16).When the porosity was in the range of that of human trabecular bone, the Tb.Th, Tb.N, and DoA of our models were also in the same range.It can be noted that the DoA in the human femoral bone was not much different to the human vertebral bone values, which could be due to sampling very locally in a specific region.However, the DoA can be further influenced by changing the PDF, as described above.The comparison suggests that our algorithm could mimic human trabecular bone   to a certain extent for a large range of porosities.Furthermore, it is flexible in terms of ability to represent different porosities, anisotropy, and stochasticity requirements of the structure.However, it is the combination of the above parameters that constitutes a trabecular bone, and the relative importance and occurrence of different parameters is not entirely mapped.
To further compare human bone (Arlot et al., 2008), a Sawbones (Poukalova et al., 2010) specimen, rod-like Voronoi models from a previous study (Frayssinet et al., 2022), and the new bone models, two models were generated to have similar morphological parameters to two human vertebral bone structures.One of them had a high porosity, i.e. 90%, to represent human osteoporotic bone (Fields et al., 2009).Another one was generated with a lower porosity, i.e. 65%, to represent healthy bone.The parameters are shown in Table 2.As we can scale our model with any number to suit different bone situations, the Tb.Th of all models were normalized to 17 pixels.Therefore, we can compare other parameters under the same scale, i.e. same Tb.Th.When we controlled the porosity and DoA with our input parameters and seeds, similar Tb.Sp parameters were found between our models and human bone (Table 2).The differences in DoA can be controlled to below 10% for both porosities, and the differences of Tb.Sp were 20% for the model of 65% porosity and 10% for the 90% porosity model.
Comparing to the rod-like Voronoi model, our model also showed improved fitting for DoA and Tb.Sp.Another advantage of our model is the larger range of EF.As reported in previous rod-like Voronoi tessellation-based models (Frayssinet et al., 2022), the EF is almost a constant number (approximately 0.23) despite varying the porosity (0.6-0.8), which indicates a constantly more rod-like structure.However, as reported in (Doube, 2015), the EF can be varied for different trabecular bone structures from − 1, highly oblate (discus-shaped or plate-like), to +1, i.e. highly prolate (rod-like) structure.In our models, we found that the EF factor could increase from approximately − 0.1 to +0.3 for models with increasing porosity.
Comparing to Sawbones, the algorithm is more flexible in changing different parameters to suit different situations, e.g. the different porosities of osteoporotic or healthy bone, the anisotropy and distribution of seeds for different anatomical sites and bodies.In addition, the strut size can also be controlled and the bone model can be manufactured with 3D printing to generate repeatable experimental results.The morphological parameters of Sawbones are predefined by the manufacturing process, and not adapted to the bone site to be tested.In contrast, the generated structure can be optimized to better fit a specific trabecular bone structure, i.e. the positions and number of seeds are the design parameters and the differences in morphological parameters between the generated structure and the destinated real trabecular bone structure are the objective functions to be minimized.
Furthermore, the algorithm can be used to develop different unit cells and be combined to build bone structures with varying properties in their domains.If the material properties change within the bone, e.g.due to local pathological differences, unit cells with different parameters can be used to accommodate for the corresponding properties (Zheng et al., 2020).
To further visualize the similarity between a generated structure and real trabecular bone structure, an example is shown in Fig. 17.As the porosity is relatively low, more plate-like structures are found.Between the plate-like structures individual or connected smooth voids can be found.Anisotropic irregular voids can also be observed in the three structures, which leads to anisotropic struts, and structure.
There were some variations in the experimental results, which could be due to several factors.First, there is a layering effect of the 3D-printed material, i.e. the material properties are not isotropic.Second, the printing resolution is insufficient, i.e. the slicing algorithms in the 3D printing software must omit some details of the structure in order to create a continuous printing path.These could be contributing reasons to the slightly lower experimental stiffness in comparison to that found in the numerical analysis.However, the average values were very similar (Fig. 14a,c), demonstrating the predictability of the global mechanical properties of the models.The numerical simulation can hence be utilized to predict the stiffness of these types of structures.The stiffness differences between the analytical, numerical, and experimental results of the anisotropic structure in the z-direction were larger than other directions, as in the sparse direction there is possibly high amounts of bending of thin slender structural parts which influence the global stiffness in a substantial way, which is not considered in the theoretical expression, Eqn. ( 5).With that in mind, the overall agreement between numerical simulations, experiments and the simplified theoretical expression is surprisingly good.A tuning of parameters in Eqn. ( 5) may provide more accurate prediction results.However, the optimized exponent may only suit one type of seed distribution.(Rincón-Kohli and Zysset, 2009;Arlot et al., 2008;Mueller et al., 2009;Fields et al., 2009;Chen et al., 2010;Turunen et al., 2013), and our porous bone models (blue).

Table 2
Morphological parameters of different porous structures.DoA is the degree of anisotropy; Tb.Th is the trabecular thickness; Tb.Sp is the trabecular separation (Bouxsein et al., 2010) Therefore, it was not further investigated in the study.
The pore interconnectivity was estimated from the closed cell porosity measurements and longitudinal permeability.For the bone model 1 shown in Table 2, the closed cell porosity was 0%, i.e. no closed cells were present in the model of 66% porosity.The pore interconnectivity and permeability will determine how well blood and nutrients spread within the structure in the living material, but also how well an augmentation material such as a bone cement spreads, and are hence important parameters for accurate modelling of bone structures.Compared with Sawbones, which contains closed cells at 30 PCF (75% porosity) (Hoffmeister et al., 2018), our model can generate open-cell structures in a lower porosity range (≤66%) as well.Both Sawbones and our model showed increased longitudinal permeability when porosity increases, as expected (Fig. 15).However, the permeability value of our model was always higher than the Sawbones for the same porosity, except at 90% porosity, where it was similar.In particular, our model showed permeability at porosity 60%.We can hence conclude that our models provide a better interconnection, and that we have interconnection for a large porosity range.
If this model is utilized for bone substitution, the high porosity models (≥70%) are suggested rather than the low porosity models, as the lower porosity models can result in excessive stiffness and low permeability, both of which may limit bone tissue ingrowth.However, the porosity of human trabecular bone is indeed usually higher than 60% -the reported range is from 60% to more than 90% (Ding et al., 1999).When the porosity is higher than 97%, small, isolated trabecular structures may be found.On the other hand, these can be eliminated by identifying the largest volume and delete the smaller ones.
The main weakness of the present algorithm is its relative complexity.The intersection between different polyhedrons must be checked by the algorithm ergodically, which is time consuming.Furthermore, the Boolean operation could create small surface triangles with poor quality, especially when the number of seeds is high, which complicates FEM meshing.
Another limitation lies within the possibilities of manufacturing these models at the same scale as human bone tissue, i.e. similar exact values of Tb.Th, Tb.Sp, etc. Indeed, in the present study the printed 3Dmodels were scaled up by a factor of 2-4.Previous studies have worked on the relationship between the printing quality and the scale factor when using FFF and SLA techniques (Grzeszczak et al., 2021;Wu et al., 2020).We believe that this problem can be improved by the use of other additive manufacturing techniques, such as laser beam powder bed fusion and/or further development of existing 3D printing technologies.
The influence of the distribution of the seeds on the morphometrical parameters of the structures could be further studied (Zheng et al., 2020).As mentioned before, the distribution of the seeds will not only influence the distribution of cells, but also the geometry of the polyhedrons and the mechanical properties of the structure.The influence of the distribution can be further evaluated to generate anisotropic structures, for example, such as in the femoral bone (Li et al., 2012).
Biocompatibility as a function of morphology could be investigated in a future study, if the model is aimed to be used as a bone scaffold in vivo.Indeed, both pore size and shape have been found to influence cell response (Diez-Escudero et al., 2020, 2021).Therefore, the scaffolds could potentially be optimized to the biological response by changing the seed distribution.

Conclusions
In this study, a numerical algorithm was developed to design heterogeneous partial open-cell porous geometries in order to provide an improved representation of trabecular bone.The algorithm was verified with numerical analyses and validated by experimental results, which showed differences of approximately 5% for the structures' global stiffness.Several morphological parameters were compared between the models generated with the algorithm and actual trabecular bones and showed reasonable agreements for certain key morphological parameters.
The developed models are not limited to rod-like, high-porosity truss frameworks as in classical Voronoi tessellation-based porous structures, i.e. the algorithm can provide a more bone-like structure, for a larger porosity range.The structure can indeed transfer smoothly from platelike at lower porosity, to rod-like at higher porosity, and display both at medium porosity.Several control parameters were evaluated to create structures to suit different scenarios, such as different anatomical sites.
The algorithm provides an improved way to represent the trabecular bone structure compared to earlier models.the work reported in this paper.

Fig. 3 .Fig. 4 .
Fig. 3. 2D Voronoi diagrams with uniformly distributed random seeds (a), averaged cells (b), and a structure with the main axis in the horizontal direction (c).

Fig. 6 .
Fig. 6.A discretization (a) and the voxelization of the same structure (b).
Y.Zhou et al.

Fig. 9 .
Fig. 9. Two typical experimental force-displacement curves (the test was interrupted close to the load cell capacity since the elastic region only was of interest).

Fig. 12 .
Fig. 12.Three different models with different distributions of seeds, i.e. a) randomly distributed with uniform PDF; b) uniform distribution; and c) normal distribution with the mean at the center of the structure.The number of seeds in each model is 100 and they have the same porosity (80%).
Y.Zhou et al.
Trabecular bone preparation.To compare the model with real trabecular bone structures, 38 rabbit distal femurs were acquired from a local butcher.The soft tissue was removed with scalpel, and the bone was stored in the freezer at − 20 ∘ C until use.Two artificial bone models were generated for scanning, to mimic trabecular bone structures of different conditions, i.e. low porosity (~65%) and high porosity (~90%).Skyscan 1172, Bruker, Belgium) to verify the printing quality and the geometrical parameters of the structure.The voltage and current used were 60 kV and 167 μA with an Al filter for the printed structures, 100 kV and 100 μA with an Al + Cu filter for the lapine bones.The voxel size was 11.95 μm, and the exposure time 2350 ms.The rotation step was 0.4 • with 180 • rotation.The micro-CT images were reconstructed with NRecon (Bruker, Kontich, Belgium) with automatic threshold.CTAn (Bruker, Kontich, Belgium) and ImageJ2 ).
Y.Zhou et al.