Uncovering a high-performance bio-mimetic cellular structure from trabecular bone

The complex cellular structure of trabecular bone possesses lightweight and superior energy absorption capabilities. By mimicking this novel high-performance structure, engineered cellular structures can be advanced into a new generation of protective systems. The goal of this research is to develop an analytical framework for predicting the critical buckling load, Young’s modulus and energy absorption of a 3D printed bone-like cellular structure. This is achieved by conducting extensive analytical simulations of the bone-inspired unit cell in parallel to traverse every possible combination of its key design parameters. The analytical framework is validated using experimental data and used to evolve the most optimal cellular structure, with the maximum energy absorption as the key performance criterion. The design charts developed in this work can be used to guide the development of a futuristic engineered cellular structure with superior performance and protective capabilities against extreme loads.

structure, an open-cell strut-like structure or a hybrid of the two. Plate-like structures (Fig. 1a,b) can be found in denser bones that can support high loading demands, such as the human femur. Hybrid plate-like and strut-like structures typically exist at the interface between the compact (cortical) shell and the cellular (trabecular) core of bone. Trabecular bone typically has a relative density less than 70% 17 . More specifically, the relative density of open cell rod-like and closed cell plate-like trabecular structures is typically less than 0.13 and greater than 0.2, respectively. At intermediate relative densities, trabecular bone is a hybrid of rod-and plate-like elements 17 . Concave (CCV) and convex (CVX) cell geometries can be observed in the hybrid (HYB) plate-like cellular structure of trabecular bone (Fig. 1b). Mimicking these novel characteristics of trabecular bone can lead to the development of engineered cellular structures with superior energy absorption capabilities, by controlling their buckling and collapse mechanisms. An overall snapshot of the analytical and experimental framework developed in this work is shown in Fig. 1, from the biomimetic conceptual level to the most optimal engineered cellular structure. The bone-like cellular structure (Fig. 1b) is mimicked using a periodic Voronoi diagram, which is an organised set of points (sites) that control the shape of its polygons (Fig. 1c) 40 . A unit cell (Fig. 1d) is extracted from the Voronoi diagram and manipulated to yield a bone-like CAD model with defined design parameters (Fig. 1e). The 3D printed bioinspired prototype (Fig. 1f) is generated from this CAD model and tested under compression. The analytical model is developed to predict the collapse load and Young's modulus of the bone-like structure, and validated using the experimental force-displacement curve (Fig. 1g,h). Extensive simulations are executed in parallel to  Table 1. Element lengths and angles of the bone-like half unit cell (see Fig. 2 for a visual). * a = H www.nature.com/scientificreports/ exhaustively manipulate the design parameters and identify the best bone-like unit cell that exhibits the highest energy absorption capacity (Fig. 1h,i).
Design of bone-like unit cell. Conventional cellular structures, including re-entrant and honeycomb, are typically built from a single periodic unit cell 30,41 , which limits its design parameters. In contrast, the bioinspired unit cell (Fig. 2) is constructed with an upper and lower sub-cell, which provides more freedom in design. Comprehensive details on generating the 3D printed bioinspired unit cell using Voronoi diagrams and its demonstrated advantages over traditional hexagonal and re-entrant structures are provided in 15 . The design parameters of the bio-inspired unit cell are summarised as follows: H and W are the height and width of the unit cell, respectively; α and β are the angles of the upper and lower sub-cell, respectively; l ut and l lt are the lengths of the upper and lower tie, respectively; t w is the thickness of the unit cell wall; and t is the out-of-plane thickness of the unit cell. Note that the width of the unit cell w s can be determined from the above parameters, as demonstrated hereafter.
The lengths and angles of each element of the unit cell are tabulated below, noting that only half the cell will be analysed due to symmetry. The angles are determined with respect to the global horizontal axis.
Constraints are placed on the key design parameters to obtain a valid unit cell within the bounds (0, W) : 0 < w s < W ; 0 < b < W ; 0 < l lt < W ; 0 < g < W ; and 0 < l ut < W.
The relative density can be inferred from the lengths of the elements in Table 1 as follows: where ρ is the density of the cellular structure; ρ S is the density of the solid material; t w is the thickness of the cell walls; L e is the length of element e; and W and H are the width and height of the unit cell, respectively.
Analytical model development. The deformation of a cellular structure can be predicted by applying the direct stiffness method of structural analysis to the bio-inspired unit cell design shown in Fig. 2 42 . The forces and moments in each element are related to the displacements and rotations in the global coordinate system as follows: where T e is the transformation matrix from the local to global coordinate system; K e is the stiffness matrix of (1) (2) T e K e T T e � u Ge = � f Ge Figure 3. Compressive force-displacement curves for the 3D printed bone-like cellular structure 16 . The deformations in the linear-elastic, plastic and densification regions are also illustrated. www.nature.com/scientificreports/ The boundary conditions are determined based on the symmetry of the unit cell: horizontal sway-rigid support at node A (v Ga = θ a = 0) ; and vertical sway-rigid support at nodes F to J (u G = θ = 0) . The stress, strain and Young's modulus of the unit cell can be written as follows: where the force P is arbitrarily chosen in the linear-elastic region (Fig. 3), n is the number of unit cells across the width of the structure, W and H are the width and height of a single unit cell, t is the out-of-plane thickness of the unit cell, and v Gf is the global vertical displacement at node F.
Buckling of the unit cell is dictated by the critical column, where each column is supported by elastic restraints. The critical buckling load of an elastically supported column can be inferred by solving the Euler buckling ordinary differential equation of a pinned-pinned column and applying elastic boundary conditions to the pinned joints 43,44 . This well-known equation and its solution are represented as follows: where v is the lateral displacement, 2 = P/E S I and E S I is the flexural rigidity of the column; P is the compressive load; and A and B are constants of integration. By differentiating Eq. (4) twice and imposing elastic boundary conditions at the joints of the column (at x = 0 and x = L ), the following expressions are obtained: where M(0) and M(L) are the moments at the joints of the column, and A and B are constants. By combining Eqs. (5a-f), the following expression can be obtained, which relates the rotations at the joints of the column with the bending moment at the joint at x = L:  Table 3. Euler buckling load for each column of the analysed half unit cell (Fig. 2).

Column
Critical buckling load ( P CR N) Effective length ( L e mm) k factor ( L e /L) www.nature.com/scientificreports/ The rotation at the joint at x = L can be written in terms of the bending moment M(L) and the stiffness of the restraining member: Table 4. Optimal unit cells obtained by modifying all five design parameters simultaneously, namely the cell wall angles and tie lengths. The re-entrant auxetic topology is also presented as a benchmark case.
Maximum specific properties (MPa/g) Other specific properties (MPa/g) Angles Tie lengths (mm) Unit cell www.nature.com/scientificreports/ where K L is the stiffness of the joint at x = L . The negative sign signifies that the bending moment in the bracing element opposes the bending moment in the column. Substituting Eq. (7) into Eq. (6) gives the following relationship: Figure 5. Influence of the sub-cell angles (a) and tie lengths (b) on the Euler buckling stress. The unit cell corresponding to the peak Euler buckling stress is also illustrated. Table 5. Maximum buckling stress obtained for each design parameter, and its corresponding critical column and L t /L ratio. The reader may refer to Fig. 2 to identify the critical columns.  where α is the stiffness coefficient of the joint; L t is the length of the restraining tie; and L is the length of the column element being analysed for buckling. The smallest root for Eq. (9) can be obtained graphically for each combination of L t /L and recorded in a design table.

Design parameter Critical column
The energy absorbed by the half unit cell can then be calculated from Eqs. (3) and (9) as follows, assuming elastic buckling: where σ CR is the critical buckling stress (P CR /A) and A is the cross-sectional area of the column element. The mass of the unit cell can be calculated as follows: where ρ S is the density of Markforged Nylon, namely 1.1 kg/m 3 , t w is the thickness of the cell walls, t is the outof-plane thickness of the unit cell and L i is the length of each element.
Results and discussion experimental validation. The analytical model is validated using the results from our compressive tests on the bone-like cellular structure, which are presented in Fig. 3 16 . The goal is to verify that buckling is the primary collapse mechanism of the structure as opposed to plastic deformation in the material. The following parameters are used for validating the analytical model, which are identical to the tested cellular structure (see Fig. 3): W = H = 20 mm ; α = 63.43 • , β = 25 • , γ = 41.14 • ; l ut = 10 mm , l lt = 5 mm ; t w = 2 mm ; t = 30 mm ; E S = 320 MPa (Markforged Tough Nylon, 3D printing material). A force of 1,000 N is chosen arbitrarily within the linear elastic region of the force-displacement curve (Fig. 3), and divided over the number of half unit cells  www.nature.com/scientificreports/ (six); the force per half unit cell is thereby 166.67 N . A vertical displacement of v Gf = 0.84 mm at node F v Gf was predicted from the analytical model by solving Eq. (2). The vertical displacement at node F v Gf is used to calculate the strain and Young's modulus from Eq. (3). The critical buckling load was then calculated using Eq. (9) for each column and presented in Table 3. It can be inferred from Table 2 that the predicted Young's modulus of the half unit cell is in agreement with the experiment within 7%. The experimental critical buckling load is also within the range of the predicted results. Although column AB provided the closest buckling load (Table 3) with the experiment, it is difficult to predict that this column will first buckle, as the condition of instability depends on imperfections that are inherent in the 3D printing process, as well as the load transfer through the structure. Regardless, the analytically predicted Young's modulus and critical buckling load are conservative for design purposes. Importantly, the predicted minimum critical buckling stress (in column BC) is well within the linear elastic region of the material curve of the Markforged Tough Nylon 3D printing material 16 , which verifies that buckling is the primary collapse mechanism of the structure as opposed to plastic collapse. The effective lengths and k factors are calculated from the Euler buckling load and presented in Table 3, which are slightly larger than 1 as expected, noting that k = 1 represents a sway rigid-to-rigid column with one point of inflexion or a braced pinned-to-pinned column. The critical buckling loads listed in Table 3 signify that buckling can occur in either column AB, BC, CD or DE. Column EF is ignored because the critical buckling stress exceeds the linear elastic limit of the material (see 16 ). Hence, for the optimisation studies conducted hereafter, the critical buckling load is calculated for columns AB to EF for each unit cell configuration and the minimum buckling load is used to compute the strain energy density. optimal design parameters. A computer algorithm was developed to execute many parallel simulations to obtain the design parameters ( α, β, γ , l ut , l lt ) that produce the biomimetic unit cell with the optimal energy absorption. To this effect, the sub-cell angles and tie lengths were modified as follows: α, β, γ : [5, 175 • ] in increments of 1 • ; and l ut , l lt : [1, 19 mm] in increments of 1 mm. This equates to 1.7 billion cases, which were run on a supercomputer using 170 threads (total runtime of 2.3 h). The results are reported in Table 4 and analysed herein. Case 1 ( α = β = γ = 90 • , l ut = l lt = 1 mm ) results in vertical columns and produces the stiffest structure with the lowest energy absorption. This is expected because the columns are more significantly rigid under axial deformation as opposed to flexure. Case 2 ( α = 79 • , β = 30 • , γ = 34 • , l ut = l lt = 9 mm ) results in a structure with a significantly low stiffness and a high buckling load compared to Case 1. This effectively increases the energy absorption relative to Case 1 according to Eq. (10). The design parameters are analysed individually in the following sub-section to gain further insight into this behaviour. Case 3 ( α = 122 • , β = 20 • , γ = 29 • , l ut = 14 mm, l lt = 18 mm ) produces the most optimal structure with the lowest stiffness and highest energy absorption. Although the Euler buckling load is nearly identical to Case 1, the stiffness is significantly lower. According to Eq. (9), the buckling load depends on the ratio of the length of the tie restraint L t to the length of the restrained column. The mechanical properties of the re-entrant auxetic structure, which is extensively investigated in the literature for protective applications, are also presented in Table 4 for comparison. It is evident that the bioinspired structure shows a significantly higher energy absorption capacity with a slight increase in mass (around 4%) under quasi-static loading. This promising result can motivate further numerical investigations to obtain the optimal bioinspired structure under extreme loadings in future work. The design parameters are also analysed individually in the following sub-section to gain further insight into this behaviour.

Influence of individual design parameters.
In this set of studies, the influence of each individual design parameter of the biomimetic unit cell (see Fig. 2) on its energy absorption, stiffness and Euler buckling stress is analysed. To this effect, the parameters of the baseline unit cell are modified individually from the reference baseline configuration ( α = 63.43 • , β = 25 • , γ = 41.14 • , l ut = 10 mm and l lt = 5 mm ). It is evident in Fig. 4 that the sub-cell angle β has the most significant influence on the stiffness of the unit cell, followed by the length of the lower tie l lt . It is expected that β = 90 • results in the stiffest unit cell because the lower sub-cell columns become parallel such that they transmit the load via axial deformation predominantly, as opposed to flexure.
It can be observed in Fig. 5 that both sub-cell angles β and γ , and the tie lengths l ut and l lt produce the unit cell with the highest Euler buckling stress. According to Eq. (9), the Euler buckling load is governed by the ratio of the length of the tie restraint ( L t ) to the length of the restrained column ( L ). The ratio L t /L and the critical column are listed in Table 5 for each individual design parameter that produced the maximum buckling stress. It can be deduced that the plateau region e.g. for 45 • < α < 103 • gives the same L t /L ratio and thereby produces an identical Euler buckling stress.
It is evident in Fig. 6 that the sub-cell angle γ = 31 • and upper cell wall angle of l ut = 7.8 mm result in unit cells with the highest strain energy density. According to Eq. (10), the optimal energy absorption can be obtained by minimising the Young's modulus E whilst maximising the critical buckling stress σ CR . According to the fiveparameter optimisation in Table 4, modifying each design parameter individually will significantly underestimate the unit cell with the maximum strain energy density ( U max = 18.8 mJ/mm 3 ). It can be deduced from Figs. 4 and 5 that minimising the sub-cell angle γ produces a relatively constant low stiffness ∼ 13 MPa and a maximum buckling stress of ∼ 19 MPa , respectively. The sub-cell angle β and lower tie length l lt can overcome this restriction to reduce the stiffness further and thereby obtain the unit cell with the optimal energy absorption (see Fig. 4).

Influence of combined design parameters.
In this set of studies, the combined effects of two design parameters of the biomimetic unit cell (see Fig. 2 www.nature.com/scientificreports/ values are also illustrated near each surface plot. It is evident from Fig. 7a,c that the sub-cell angle β (see Fig. 2) is a key design parameter that dictates stiffness, whereby the lower columns of the unit cell are relatively vertical. This trend is consistent with that observed in Fig. 4, which showed that modifying β alone produces the stiffest unit cell by far. The unit cell with the highest stiffness ( E max = 50.3 MPa ) is obtained when β = 99 • and γ = 69 • (Fig. 7c). This is an expected result, given that the columns of both sub-cells are relatively vertical, which promotes axial deformation over flexure. A relatively stiff unit cell ( E = 44.7 MPa ) can also be obtained by other combinations of sub-cell angles e.g. α = 76 • and β = 88 • (Fig. 7a), which further reinforces that the stiffness is controlled by β as a governing parameter. In contrast, Fig. 7b,d show that the unit cells with the highest stiffness have more oblique members such that flexural deformation becomes more prominent, which consequently reduces the stiffness of the structure. The highest critical buckling stress ( σ CR = 19.7 MPa ) is obtained from α = 80 • and β = 33 • (Fig. 8a). This is followed by σ CR = 19.3 MPa , which results from the tie lengths l ut = l lt = 9 mm (Fig. 8d). From the discussion preceding Fig. 5 and Table 5, and referring to Eq. (9), the critical buckling load depends on the ratio of the restraining tie length to the length of the restrained column L t /L . It is evident that the unit cell in Fig. 8a has the shortest tie (attached to longer columns CD and DE), which effectively yields a small L t /L ratio and thereby attracts a higher buckling stress. Longer plateau regions can be observed in the plots of Fig. 8b,c because the unit cells that produce the highest buckling stress have a more uniform structure, whereby the relative lengths between the struts and ties are similar.
The unit cell with the highest strain energy density ( U max = 14.4 MPa ) is obtained when α = 88 • and γ = 32 • (Fig. 9b). Recalling from Eq. (10) and the discussion preceding Fig. 6 that the peak strain energy density is obtained when the critical buckling stress is maximised and the stiffness is minimised, it is evident that the unit cell illustrated Fig. 9b is designed to balance axial and flexural deformation, whereby the oblique orientations of its columns reduce the stiffness. The unit cell also has relatively similar column and tie lengths, which thereby produces a relatively high critical buckling load ( σ CR = 18.9 MPa ). Note that the critical buckling load does not vary within a large interval compared to the stiffness, given that the L t /L is constrained by the size of the unit cell. Importantly, other design parameters can be modified to produce unit cells with high energy absorption e.g. U = 14.2 MPa is obtained from β = γ = 153 • (Fig. 9c). This reinforces that γ is a governing parameter that has a significant influence on the strain energy density of the bone-like unit cell. The unit cells in Fig. 9a,d that produce the peak strain energy density have relatively short tie lengths compared to the struts, which has the effect of increasing the critical buckling stress. The oblique struts of these unit cells produce more prominent flexural deformations, which consequently reduces the Young's modulus and thereby increases the strain energy density according to Eq. (10).

conclusions
An analytical framework was developed to predict the critical buckling load, stiffness and energy absorption of a cellular structure based on trabecular bone. The framework was validated using experimental data obtained from compressive tests on a 3D printed bone-like cellular structure. The validated framework was used to run extensive simulations in parallel to obtain the most optimal unit cell configuration with high computational efficiency, using the energy absorption as the key performance criterion. It was observed that all design parameters, namely the sub-cell angles and tie lengths of the unit cell played a significant role in obtaining the highest energy absorption. Effectively, oblique sub-cell angles control the stiffness of the unit cell by balancing axial and flexural deformations. Specifically, unit cells with vertical column elements produced the highest stiffness whereas oblique unit cells exhibited a lower stiffness and higher energy absorption. The relative lengths of the columns and ties were found to control the buckling stress of the unit cell by modifying the joint stiffness of the column elements. It was observed that the optimal bioinspired structure showed re-entrant geometries in both sub-cells, which resulted in a low stiffness and high energy absorption. The optimal bioinspired structure was benchmarked against a re-entrant auxetic topology, which is extensively investigated in the literature, and showed an increase in energy absorption by a factor of four. The next step is to manufacture the optimal unit cells and assess their energy absorption capacity under static and dynamic loads. The analytical model will also be extended further to account for dynamic loading conditions. Machine learning techniques, including artificial neural networks and genetic algorithms, will also be employed in future work to reduce the number of simulations required under dynamic loading, which can be computationally expensive. The size effect of the unit cell can also be important under dynamic loading conditions and will also be investigated in future work.