Architected implant designs for long bones: Advantages of minimal surface-based topologies

(cid:1) Beam-and surface-based lattice topologies are all suitable for bone implant design. (cid:1) In surface-based topologies, speciﬁc area and implant stiffness can be decoupled. (cid:1) Surface-based topologies are more efﬁcient at promoting bone growth than beam-based topologies. (cid:1) Stochastic spinodal shell topologies require the smallest safety factor.

the bone defects.More specifically, the implants must be porous to allow efficient transport of nutrients and removal of waste and promote bone ingrowth (which is controlled by the surface area and permeability of the implant) [12][13][14][15][16]; they must have sufficient strength to support daily physical loadings and minimize the occurrence of fatigue failure (which is governed by the peak stress experienced in the implant) [13,17]; and they should have similar stiffness (characterized by the effective Young's modulus of the implant) to the host bones, to prevent bone loss due to post-surgery stress shielding [13,[18][19][20].
Architected cellular materials, i.e. porous solids consisting of a (generally) periodic repetition of an optimally designed unit cell, have been extensively investigated over the past couple of decades.Topology optimization of their unit cell provides a unique opportunity to tailor the effective properties of the material to maximize multifunctional performance for a wide range of applications.In the case of implants, the amount of porosity u (defined as the volume fraction of the open space, and equal to u ¼ 1 À q À , with q À the relative density of the cellular material) and its topological arrangement (i.e., the unit cell design) can be designed to facilitate bone growth.Optimal implant designs match the stiffness of the surrounding bone while providing constant mechanical stimulation to evolving tissue in the pores, maximize surface area for bone/implant adhesion [21,22], optimize fluid transport [14,23] through the pores, and minimize stress intensification for optimal fatigue life [24,25].With the advancement of metal additive manufacturing technologies, in particular selective laser melting (SLM) and electron beam melting (EBM), it is now possible to accurately and cost-effectively fabricate porous metallic implants made of unit cells with nearly arbitrary topological complexity [26][27][28][29][30][31].
Strut-based architected materials have garnered significant attention in the development of candidate bone implants [13,[31][32][33][34][35]; they are relatively well understood (at the unit cell level) and provide a clear pathway to fill volumes through repetition of unit cells.These advantages are considerable since implant performance can be anticipated from well-established design principles [36][37][38].However, the intersections of struts lead to local stress concentration that can severely limit fatigue life [39][40][41].Furthermore, these intersections may play a role in inhibiting homogeneous bone growth [42][43][44], although additional in vivo studies are required to assess the magnitude of this problem.
Surface-based topologies, characterized by very uniform nearzero mean curvature and the lack of strut intersections (or ''nodes"), offer an attractive alternative to strut-based geometries, and may circumvent these shortcomings.The most notable examples of these surface-based topologies are triply periodic minimal surfaces (TPMSs), such as Schwarz P (primitive) and Schwarz G (gyroid) surfaces [45].Over the past few years, TPMSs have received significant attention and have been investigated for their multifunctional performance in multiple fields [46][47][48].For bone implant applications, TPMSs are rapidly becoming topologies of choice, thanks to their uniform curvature, high specific surface area, high strength, and topological similarity to human trabecular bones [49,50].As thoroughly demonstrated in [51,52], a wide range of TPMS-based topologies (primitive, I-WP, Gyroid, and F-RD) possess ideal combinations of Young's modulus, yield strain, fatigue behavior, permeability and morphological characteristics, rendering them excellent candidates for the fabrication of implant scaffolds that best mimic the properties of human trabecular bone.
A recently introduced minimal surface-based topology is the spinodal shell, a stochastic surface derived from the interface of two spinodally decomposing phases [53].While naturally stochastic, these spinodal shells exhibit fully isotropic properties, exceptional specific stiffness and strength (almost on par with those of topologically optimized TPMSs), combined with high imperfection insensitivity (unique for surface-based systems) [53].Spinodal surface-based architected materials can be computationally designed and additively manufactured with commercially available approaches (e.g., laser powder bed fusion or electron beam melting), with pore sizes comparable to those of existing implants.Unlike all regular implant topologies, though, this novel class of shell metamaterials is uniquely amenable to self-assembly via spinodal decomposition of suitable precursors, followed by material conversion techniques [54,55]; this emerging manufacturing approach allows cost-effective fabrication of cm-scale implants with micro-scale porosity, yielding combinations of specific surface area and pore topology optimization that are unavailable in any regular metamaterial topology.While the mechanical performance (stiffness, strength, toughness and failure mechanisms) [53,[56][57][58][59] and the transport properties [60,61] of these spinodal surface-based architected materials have been recently investigated for structural and energy applications, to the best of our knowledge their effectiveness as metamaterials for bone implants (and in particular, their potential to promote bone growth) has not been explored.
In this work, we numerically investigate and compare the performance of three distinct implant topologies for long bones, at three relative densities (q À r ¼ 0:1; 0:2 and 0.3): the octet trussbased lattice (one of the most mechanically efficient periodic truss-based lattices), the Schwartz P surface-based lattice (a periodic and minimal surface-based design) and the spinodal surface-based lattice (a stochastic near-minimal surface-based design).Each implant topology (reinforcement) and its surrounding soft tissue (matrix) are modeled as a composite system via finite element analysis, for calculation and analysis of stress fields in the metallic implant and strain fields in the soft evolving tissue.Based on the stress and strain fields, three performance metrics are defined to assess the relative performance of different implant topologies: (i) the effective elastic modulus of the implant (relevant to stiffness matching to connecting bone), (ii) the load safety factor (relevant to fatigue life in service conditions), and (iii) a new metric (strongly related to the well-established osteogenic index [62][63][64]), which captures the complex distribution of strains in the tissue into a scalar parameter and provides a heuristic but rigorous estimate of the implant's ability to promote bone growth.Naturally, there are a host of factors that impact bone growth, and the metric should be regarded as a foundational attempt to assess the relative performance of implant designs.
Inspired by previously-published work on minimal surfacebased implants [51,52], this study introduces two novel concepts: (i) the definition of a simple scalar parameter, which can be readily extracted from finite elements calculations of the strain fields in the evolving soft tissue and provides a basis for the assessment of implant topology effect on bone growth, and (ii) a class of stochastic surface-based architected materials, spinodal shells, which are amenable to scalable fabrication via self-assembly, allowing controlled introduction of microscale porosity and unprecedented surface area per unit volume.
For q À r ¼ 0:1 À 0:3, we show that all topologies fabricated from common constituent materials for metallic implants (Ti-6Al-4 V, 316L SS, and Co-Cr-Mo) can be stiffness-matched to long bones.However, surface-based topologies can be realized with larger and nearly shell thickness-independent surface area per unit volume, and may promote significantly more bone growth than conventional truss-based designs; this makes them optimally suited for long bone implants.Spinodal shell designs share the same advantages as all TPMS-based topologies, while exhibiting less sensitivity to fatigue failure and being compatible with fabrication processes that enable orders of magnitude increase in surface area.

Synopsys of approach
We use four measures of biomechanical efficiency to assess implant performance: the Young's modulus of the implant, the safety factor of the implant under service load, the specific surface area of the implant and an index to assess the relative bone growth rate between implants (to be defined subsequently).Whereas the first three measures could be extracted by finite elements analyses of the porous metallic implant alone, the fourth one requires modeling of the tissue occupying the pores of the implant material.Hence, we model each implant as a composite material consisting of a metallic phase (the actual implant structure) with the desired topology, surrounded by a softer phase simulating the boneforming tissue.Finite element meshes for each implant design (varying implant topology and relative density) are constructed as detailed in Section 2.2, and interfacial surface areas per unit volume are extracted as explained in Section 2.3.The models are then loaded in uniaxial compression, and the Young's modulus, the safety factor and the bone growth rate index are calculated as described in Section 2.4 and Appendices A and B.

Generation of composite topologies
Each of the three distinct implant topologies is generated numerically to fill a cubic volume of length L s with 3 Â 3 Â 3 unit cells.(The number of unit cells was chosen based on a previously reported convergence study that ensures minimal edge effects [65]; incidentally, a larger number of unit cells in a typical implant would require printing resolution beyond state-of-the-art powderbed fusion capabilities.)The void space in the cubic unit cells is subsequently filled with a soft matrix to build the topologyreinforced composite.The generation process is illustrated in Fig. 1 (a)-(c) and is explained in detail below.
Octet lattices of relative densities q À r = 0.1, 0.2, and 0.3 (defined as the material-filled volume divided by the unit cell volume) are generated by manually merging cylindrical beams of equal radius d and length l, connected through their identical end points at common intersections (as shown in the red box of Fig. 1 (a)), in the Dassault Systèmes CAD software Solidworks.The relative densities can be approximated by [66]: where C 1 is a constant depending on nodal geometry and determined to be 5.75 for the current study by fitting the Solidworkscalculated q À r to d=l.The octet lattices are exported to the Dassault Systèmes finite element software Abaqus where a Boolean operation is performed to fill the open space with a matrix.The resulting octet-reinforced composites are subsequently meshed with approximately 1,200,000 to 1,700,000 tetrahedron elements in Abaqus, following a mesh sensitivity analysis (Appendix C).Schwartz P shells-based lattices of q À r = 0.1, 0.2, and 0.3 are generated using the 3D image processing tool SimpleWare ScanIP.First, the mask generator function (a thresholding process that separates a portion of a given volume from the rest by assigning values to its voxels) of the built-in lattice factory is used to generate two masks, for the solid and its inverse, by infilling the two equally divided sub-domains separated by the level set surface equation [43,67,68]: wherex, y, and z (each bounded in [0, L s ], with L s the domain size) represent the coordinates of points on the surface; n x , n y , and n z are the number of unit cells in the x, y, and z-direction respectively; and l is a level-set constant, here set to zero.Second, the shell active mask function is used to generate a thin shell on the boundary of the two masks obtained above; this step is conceptually equivalent to thickening the minimal surface approximated by Eq. ( 2) to a desired shell thickness.Third, Schwartz P shell-reinforced composites are generated using Boolean operations and meshed in ScanIP with about 3,500,000 mixed tetrahedron and hexahedron elements, following a mesh sensitivity analysis (Appendix C), and subsequently exported to Abaqus for mechanical analyses.Spinodal surface-based lattices are generated by solving the Cahn-Hillard equation (the governing equation to model spinodal decomposition) [69]: where À1 u 1is the concentration difference between the material and the void phase (u ¼ 1 denotes full material and u ¼ À1 denotes void space); h is the interfacial width between the two phases; f is a double-well free-energy function, chosen as f ðuÞ ¼ 1 4 ðu 2 À 1Þ 2 ; and time t describes the coarsening process and time evolution of topologies (with t ¼ 0 representing a nearly homogeneous system, and the final value of t chosen to result in a unit cell size, m).Eq. ( 3) is solved over a cubic domain ofMx M Â M nodes using a finite difference approach with periodic boundary conditions, as detailed in [53,70], to obtain the profile of u.Stochastic spinodal solids with various unit cell sizes and relative densities can then be generated by thresholding the u profile (refer to [53] for numerical details).Spinodal cellular solids with q À r = 0.5 at ''early" and ''late" generations (corresponding to relatively early and late coarsening stages of spinodal decomposition with time, and topologically corresponding to M=m = 6 and M=m = 3, respectively) are obtained and resampled into 2-dimensional image slices for further image processing (schematic shown in Fig. 2); detailed topological characterization of early-and late-generation spinodal solids can be found in [53].Note that the images of early-generation (M/m = 6) spinodal solids are cropped in half along x, y and z-directions, to result in 3 Â 3 Â 3 unit cells (see Fig. 2), consistently with all the other topological samples discussed above.Subsequently, we mask these images in ScanIP, where spinodal shells of q À r = 0.1, 0.2, and 0.3 are subsequently generated by adjusting the shell thickness pixel values via the shell active mask function between the two masks; this is equivalent to assign a finite thickness to the interface of spinodal solids at q À r = 0.5.Finally, Boolean operations and image stacking are performed where spinodalreinforced composites are generated, smoothed, and finally meshed with 3,000,000-3,500,000 mixed tetrahedron and hexahedron elements, following a mesh sensitivity analysis (Appendix C).The meshed spinodal-reinforced composites are imported into Abaqus for mechanical analyses.

Investigation of interfacial morphology
The interfacial surface area ðS A Þ between the two phases (implant reinforcement and tissue matrix) of a composite is proportional to the number of unit cells n and the unit cell volume V, and inversely proportional to the unit cell size m in the composite.We can write the topological relation as follows: where D is a constant that uniquely depends on the interfacial surface topology.For octet-reinforced composites, D can be approximately obtained as: where d and l are the diameter and length of the truss members, respectively; C 2 depends on the beam intersection geometry and is determined to be 115.8 by fitting the Solidworks-calculated S A at various reinforcement phase densities (q À r = 0.1, 0.2, and 0.3) in this study.
For Schwartz P shell-reinforced composites, we first construct the interfacial surface according to equation (2) using MiniSurf, a triply periodic minimal surface STL file generating software [71,72].The generated interfacial surface is composed of many triangular patches, rather than tetrahedra or hexahedra as produced in ScanIP, hence allowing a simple calculation of the interfacial surface area.S A can then be obtained as: where H is the total number of triangular patches, i is the patch index that goes from 1 to H, a and b are the two edge vectors of patch i, k k is the norm operator, and Â stands for the cross product.Finally, D is extracted using Eq. ( 4).
For spinodal-reinforced composites, S A can be calculated by similar procedures (MiniSurf cannot be used as it currently cannot generate spinodal topologies): The two-dimensional slices of spinodal images generated in Section 2.2 are imported into Matlab as binary three-dimensional matrices (entries of 0 and 1 only).(ii) The Matlab built-in function isosurface is used to generate a triangularly discretized three-dimensional interfacial surface between the 0's and 1's of the matrices followed by sufficient smoothing to eliminate zig-zag pixel-like surfaces; the triangular discretization provides us with the information about the nodal connectivity.(iii) Eq. ( 6) is used to obtain S A of spinodal-reinforced composites and finallyD is calculated using Eq. ( 4).Note that the number of triangular patches is ensured to be sufficient to obtain a consistent S A value.

Mechanical simulations
All mechanical simulations are performed using quasi-static solvers in the commercial finite element package Abaqus.All constituent materials are modeled as linearly elastic solids.The denser phase (q À m = 0.9, 0.8, and 0.7; subscript m stands for ''matrix") of each composite is modeled as a soft matrix with Young's modulus, E m ¼ 1 MPa and Poisson's ratio, t m ¼ 0:3 (representing the elastic response of weak cartilage tissue, which generally forms at the surgery site), while the less dense phase (q À r = 0.1, 0.2, and 0.3; subscript r stands for ''reinforcement") is modeled as a stiff implant reinforcement with E r ¼ 123 GPa and t r ¼ 0:3 (representing the elastic response of the commonly used titanium alloy Ti-6Al-4 V).
We simulate the mechanical response of each composite under a displacement-controlled uniaxial compression by applying the following boundary conditions (Fig. 3): (1) A downward displacement d (corresponding to an effective strain e ¼ d Ls ¼ 1%, with L s the cube side length) in the zdirection is applied on the top surface, while the bottom surface is constrained from moving in the same direction.We ensure that the applied displacement is smaller than all the strut or shell thicknesses of all implant topologies; hence nonlinear kinematics can be ignored in all simulations under the applied strain.Both surfaces are free to translate in the x and y-directions.
(2) All the side surfaces are free to translate in all directions.
(3) Two corner nodes are constrained to prevent rigid body translations and rotations: one node is fixed in the ydirection and the other node is fixed in the x-direction.
To obtain the effective Young's modulus, E, we extract the reaction force F in the z-direction such that, The maximum allowable load in a porous implant is limited by the resulting peak stress S peak , subjected to a desired safety factor, SF (SF > SF min ¼ r y =S peak , where r y is the yield strength of the porous implant's constituent material, taken to be 932 MPa in this case).For highly complex implant topologies, S peak could be strongly affected by the underlying irregular mesh; hence we adopt a more robust measure of S peak as the upper 1% Mises stress threshold, as described in Appendix A.
Finally, the relative bone growth (a purely mechanical parameter used in this work as an indicator of the qualitative extent of endochondral ossification, the process by which cartilage is gradually turned into bones) in the tissue matrix phase is estimated by calculating the combined effect of maximum shear strain and compressive volumetric strain, following the procedures summarized in Appendix B. While this metric, based on established concepts in the bone community [62][63][64], should not be interpreted as an exact estimate of bone growth, it nonetheless allows us to rigorously and quantitatively compare the bone growth potential of dif-ferent implant topologies which subject the organic matter in their porosity to different states of complex cyclic strains.

Young's modulus
As we assumed linearly elastic small strain response in all mechanical simulations discussed in Section 2.4 and the reinforcement dominates the mechanical response (the modulus of the reinforcement is 10 5 times larger than that of the matrix), the effective Young's moduli of the Titanium implant can be readily extrapolated to other implant materials (such as PEEK) by simple scaling, as long as the ratio of reinforcement-to-matrix moduli remains in a similar range.The Young's moduli, E (summarized in Table 2), of composites with octet lattice, Schwartz P surface, and spinodal surface implant topologies in three different metallic materials (Ti-6Al-4V, 316L SS, and Co-Cr-Mo) commonly used in long bone implants [73] are plotted against the relative density of the implant in Fig. 4 (a).Fig. 5. Relative Young's modulus E=Er versus interfacial surface area constant D for implant topologies: octet, primitive, early-and late-generation spinodal, with relative densities of 10%, 20%, and 30% (q À r ¼0.1, 0.2, and 0.3).Note that an increase in shell and beam thickness increases q À r , and hence the Young's modulus, E, for all topologies; the same increase in thickness only increases the interfacial surface area for truss-based topologies, though, with the relative density of surface-based designs nearly independent on shell thickness.The implication is that interfacial surface area and Young's modulus can be effectively decoupled in surface-based implants.
For reinforcement/matrix moduli ratio approximately equal to 10 5 and q À r in the range of 0.1-0.3,minimal surface-based topologies (spinodal and Schwartz P) are consistently stiffer than the truss-based topology (octet lattice).The superior mechanical efficiency of the surface-based designs (primitive and spinodal) over the octet lattice derives from the uniform double curvature of the surface, which prevents bending deformation without inplane membrane stretching, resulting in stretching-dominated response and uniform stress distribution across the surface.While the octet lattice is also stretching dominated, practical implementations result in significant stress intensifications near the nodes, which depresses the mechanical performance [53,74,75].The lack of long-range periodicity typical of the spinodal design results in further reduction in stress intensifications, as previously noted in [56]; as a result, spinodal shell topologies are the most efficient.These effects are clearly illustrated in Fig. 4(b), which displays the distribution of equivalent (Mises) stresses (limited by a chosen stress limit r limit = 800 MPa to allow a more visible comparison between topologies) across the material in the topologies under investigation.The area fraction of elements that exceed the limit stress under a prescribed uniaxial strain are illustrated in red and quantified near each topology.The larger this value, the larger the fraction of the cellular material that supports the critical stress and hence the higher its contribution to the Young's modulus of the composite.Notice that spinodal topologies have the largest distribution of critical Mises stress (31% and 29%), followed by the primitive (21%) and the octet topology (13%).On the other hand, the tissue matrix phase provides essentially no structural support due to its high compliance (relevant to the earliest stages of bone growth) compared to the implant reinforcement.(The Mises stress is nearly zero everywhere in the tissue matrix phase providing little contribution to the Young's modulus of the composite as shown in Fig. 4b.)As the result, the elastic properties of the composites are nearly the same as those of the metallic cellular materials.
In terms of mechanical biocompatibility (degree of similarity in mechanical properties between the porous implant and the surrounding tissue), all topologies can generate implants that are stiffness-matched to cortical bone, with a proper choice of relative density of the reinforcement phase: q À r ¼ 0:1 À 0:3 for the spinodal shell,q À r ¼ 0:13 À 0:3 for the primitive shell andq À r ¼ 0:2 À 0:3 for the octet, as shown in Fig. 4 (a) and summarized in Table 1.Since E is one of the key factors affecting the design of bio-mechanically compatible implants, it is cross-plotted with the other properties under considerations in the subsequent sections.

Interfacial surface area
The interfacial surface area between an implant topology and its surrounding tissue is another important factor that affects mechanical biocompatibility, as larger interfacial surface area allows more bone/implant contact and osseointegration [90,91].In Fig. 5, the relative Young's modulus E=E r , with E r the Young's   modulus of the reinforcement base metal, is plotted against the interfacial surface area constant D (summarized in Table 2), a non-dimensional measure of surface area defined in Section 2.2.It can be immediately noticed that D of truss-based topologies (octet) increases fast with increasing relative density, while D of minimal surface-based topologies (both primitive and spinodal) is essentially independent of relative density over the relative density range of interest.(Note that a change in the shell thickness slightly increases the external surface area of the cubic samples, but these sample size effects are not considered in the interfacial area calculations.)The implication is that with minimal surfacebased topologies we can decouple E from surface area: the relative density of the implant q À r can be chosen to match the desired long bone modulus, while the implant can still exhibit large surface area; by contrast, with truss-based topologies, the Young's modulus of an implant would need to be increased beyond that of human long bones in order to obtain the same surface area.Importantly, the early-generation spinodal topology has the largest D, about 1.05-1.5 times that of the octet topology.Its surface area can be potentially further increased by increasing the number of pores while reducing the pore sizes to the microscale, an opportunity only provided by self-assembly manufacturing approaches unique to spinodal topologies.

Safety factor: Peak stress
The effective normalized Young's modulus E=E r is plotted in Fig. 6a against the safety factor SF min ¼ r y =S peak (summarized in Table 2) for the implants, resulting from an applied macroscopic stress of 10 MPa (about one order of magnitude lower than the typical yield strength of cortical long bones, thus allowing a reasonable physiological S peak comparison of different implant topologies while ensuring the S peak in these topologies would not reach the yield strength of the constituent material, ~1 GPa).As metallic implants can redistribute highly localized stress intensifications via localized plastic deformation, with no measurable effect on the overall mechanical response of the implant, we define S peak as the stress level that is exceeded in 1% of the structure, and Fig. 8.The maximum shear strain emaxs and the volumetric strain e vol field experienced in the q À m = 0.7 tissue matrix phase (analogous to the encompassed tissue of an implant) of octet, primitive, early-and late-generation spinodal shell-reinforced composites.The direction of the uniaxial compression is indicated by a red arrow.
extract it with the procedure detailed in Appendix A. Using a higher threshold (e.g., 2% or 3% of the structure) would result in a less conservative measure, but have no effect on the conclusions of this work (see Appendix A for details).At q À r ¼ 0:3, all minimal surface-based topologies have similar and much larger SF min (i.e., smaller S peak ) than truss-based topologies; however, this difference slowly diminishes as q À r is decreased to 0.2, where SF min in both minimal surface-based and truss-based topologies are essentially the same.As q À r is decreased to 0.1, SF min in truss-based topologies is about 18% larger than (S peak is 15% smaller than) in minimal surface-based topologies.Such difference in SF min is attributed to the fact that at q À r ) 10% truss-based topologies often have very large nonsmoothed nodal geometry and carry the majority of the load through their connecting nodes resulting in high nodal stress concentration (higher S peak ); as the relative density becomes lower (generally q À r < 10%), the nonsmoothed nodal geometry becomes small and negligible such that the connecting beams contribute to the majority of the load-carrying capacity resulting in a lower S peak than minimal surface-based topologies [75,92,93].In addition, at q À r > 10%, the nodes, rather than the connecting beams, are often the initial point of failure [94][95][96].By contrast, there are no defined nodes in minimal surface-based topologies [53,97], resulting in much more uniform stress distribution.Furthermore, minimal surface-based topologies also have higher Young's modulus, for a choice of SF and implant constituent material, than the trussbased topologies; notably, E of spinodal topologies outperform all others by almost 1.4-1.9times.The implication is that with minimal surface-based topologies (especially spinodal-shell topologies), we can choose a smaller implant relative density q À r (resulting in less material) than with the truss-based topologies, in order to achieve the same required E while satisfying the desired SF and the choice of the implant material.Note that S peak in octets is experienced in significant fractions of the nodal region and not merely at a few points, regardless of the value of q À r (Fig. 6b and  Appendix A); this confirms that truss-based designs are not unfairly penalized in this comparison.

Relative bone growth index
The relative bone growth index I growth (summarized in Table 2), an indicator of the relative endochondral ossification rate in the tissue surrounding an implant topology subject to external periodic loading, is estimated as detailed in Appendix B and plotted against normalized effective Young's modulus E=E r in Fig. 7.In general, I growth scales with E=E r in a power-law fashion, with topologyspecific exponent: I growth of minimal surface-based topologies decreases linearly with increasing E=E r , while it decreases parabolically in truss-based topologies.Furthermore, minimal surfacebased topologies have much higher I growth ; almost by a factor of two, than truss-based topologies.This can be explained by examining two strain fields in the soft tissue matrix: (1) the maximum shear strain e maxs (accelerating bone growth) and (2) the volumetric strain e vol (delaying bone growth when compressive), as shown in Fig. 8. (While the implant with a matrix phase relative density q À m ¼ 0:7 was chosen in this figure as the pronounced strain variations in this simulation enable a clearer comparison, the conclusions apply to all densities under consideration.)Notice that tissue matrices in truss-based topologies experience moderate e maxs , but highly compressive e vol ; conversely, minimal surfacebased composites experience high e maxs and relatively low compressive e vol .(We emphasize that this finding is profoundly important, clearly explaining the advantage of bone growth in minimal surface-based implants over truss-based implants from a purely mechanistic point of view.With a large number of TPMS topologies available, we believe this finding will motivate investigations into additional minimal surface-based implant designs in the bone community.)Due to their stochastic nature, spinodal topologies induce higher e maxs than other topologies; but in early-generation spinodal topologies, this gain is offset by a decrease in e vol , resulting in overall lower I growth than the primitive topologies.

Conclusions
The relative performance of both truss-based (octet) and minimal surface-based (primitive and spinodal) implant topologies was investigated numerically and presented in terms of four key design factors: the Young's modulus E, the interfacial surface area S A , the peak stress S peak and the relative bone growth index I growth .Four key results emerged: (i) All metallic implant topologies, across a wide range of relative densities, can be designed to have similar Young's modulus as human long bones (E = 5-22 GPa).(ii) Unlike in truss-based topologies, E and S A are decoupled in minimal surface-based topologies, allowing the choice of the desired implant stiffness without compromising the specific surface area.(iii) Bone growth in the tissue surrounding the implant is expected to be much higher in minimal surface-based topologies than in truss-based topologies.(iv) For a desired implant stiffness and a given external design load, spinodal surface-based implants require the smallest safety factor, resulting in minimal amount of metal required to meet all design constraints.The advantage is particularly pronounced over truss-based implants.
In summary, considering all the design factors discussed in this study for long bone implants, minimal surface-based topologies clearly emerge as superior candidates over truss-based designs.Among minimal surface-based topologies, one must consider the trade-offs between bone ingrowth (S A and I growth ), safety (S peak ), and mechanical biocompatibility (E) of the implants, with spinodal-shell designs often offering the best combinations of properties.Although not explicitly investigated in the current work, the potential of scalable fabrication via self-assembly of spinodal surface-based micro-architected materials present a unique opportunity for a dramatic increase in specific surface area S A , without compromising any other attribute.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.(3) The volume-Mises stress pairs are then fitted with the weighted kernel density estimator [99][100][101]:

CRediT authorship contribution statement
where f is a probability distribution as a function of Mises stress variable s and probability weight vector p = ½p 1 ; p 2 ; Á Á Á ; p e T ; p j is the volume of j-th element divided by the volume of all elements, and e is the total number of elements.K is a symmetric kernel smoothing function which takes the shape of a normal distribution in this study.h is the bandwidth that controls the smoothness of the resulting probability distribution.
(4) S peak (1% Mises threshold) is extracted by calculating the inverse cumulative distribution of f such that the probability that a random Mises stress variable s take on a value less than or equal to S peak is 99%, P s S peak À Á ¼ 0:99.Note that for a lower Mises stress threshold (for example a ¼ 2% or 3%), S peak can be simply calculated from step (4) with P s S peak À Á ¼ 1 À a.
The effective normalized Young's modulus E=E r is plotted against the safety factor SF min ¼ r y =S peak in the implants in Fig. A2 to show that the effect of the choice of threshold value, while observable, has no impact on the relative performance of the various topologies.Fig. A3 illustrates the regions of the implant that exceed S peak for the three different threshold values.While increasing the threshold values increases the S peak region (red), the locations of S peak practically remain unchanged; notably, S peak is experienced in major load-supporting paths in octet and primitive topologies (prone to catastrophic failure), while S peak is dis-  tributed more randomly in spinodal topologies (allowing a more gradual failure response), as previously reported [56].

Appendix B. Prediction of long bone growth
Most of the bones in human bodies are long bones such as femurs and tibias.These long bones grow longer through endochondral ossification, the process by which cartilage is gradually turned into bone tissue.It was shown that endochondral ossification is related to the osteogenic index I (defined in [62]) as a linear combination of shear stresses (bone growth facilitator) and com-pressive dilatational stresses (bone growth inhibitor) experienced in cartilage cells [62][63][64].We can write a similar relation for the relative bone growth in the tissues surrounding the implants as a linear combination of average maximum shear strain e À maxs > 0 and average compressive volumetric strain e À vol < 0, experienced in the tissue matrix phase of the composite models, as: where V m is the volume of the matrix phase over the volume of the composite and 0 j Àe  and e À vol are assumed to have equal weight (j ¼ 1) in the current study.Note, I growth is only a relative index for comparative purposes.
In order to obtain the bone growth index I growth for the tissue matrix phase, we perform the following procedures: (1) The three normal strains e 11 , e 22 , e 33 , the maximum principal strain e maxp and the minimum principal strain e minp are extracted at the centroid of each element, together with the volume fraction of the corresponding element -where p is the ratio of the volume of the element to the total volume of the matrix phase.
(2) The quantity e maxs is calculated as 2 and e vol is calculated as (3) We then produce two sets of data for all elements, with one consisting of ½e maxs p and the other consisting of ½e vol p, where p = ½p 1 ; p 2 ; Á Á Á ; p e T , e maxs ¼ ½e maxs;1 e maxs;2 ; Á Á Á ; e maxs;e T , e vol ¼ ½e vol;1 e vol;2 ; Á Á Á ; e vol;e T , and e is the total number of elements in the matrix phase.
(4) In order to reduce the effect of extreme high local strains (outliers) on the measures of central tendency, we eliminate all rows in ½e maxs p whose first component is less than the 0.5th or greater than the 99.5th percentile of e maxs .Similarly, we eliminate all rows in ½e vol p whose first component is less than the 0.5th or greater than the 99.5th percentile of e vol .
(5) Weighted kernel density estimators (Eq.(A.1)) are used to fit ½e maxs p and ½e vol p as f maxs and f vol respectively. (6)e À maxs and e À vol can be calculated as the mean of f maxs and the mean of f vol respectively.I growth ¼ V m Á (e À maxs þ j Á e À vol ) then can be found.

Appendix C. Mesh sensitivity study
A mesh sensitivity study is performed on the composite materials samples with each of the four reinforcement topologies (octet, primitive, early-and late-generation spinodal), with sample size 100Â100Â100 mm and reinforcement relative density q À r ¼ 0:3, to ensure that results are not affected by the mesh quality.Uniaxial compression (see Section 2.4 for detailed boundary conditions) is applied via a downward displacement of the top face along the zdirection (d z ¼ À1 mm), and the resulting reaction force in the zdirection is extracted.The convergence of the reaction force (normalized by the converged value) as a function of mesh density is plotted in Fig. C1.To find a suitable balance between accuracy and computational efficiency, the reasonable mesh density for each topology (indicated by a circle) is chosen when the reaction force is within 2% of the final converged value.All subsequent mechanical simulations performed in this work to assess biocompatibility efficiency have a mesh density above the respective reasonable mesh convergence value.Fig. C1.Mesh sensitivity study of octet (blue line), primitive (green line), early-(brown line) and late-generation (red line) spinodal composites.Reasonable mesh convergence for each topology, denoted by a circle in its respective color, is assumed when the reaction force in a compressive test is within 2% of the converged value.All subsequent mechanical simulations performed in this work to assess biocompatibility efficiency have a mesh density above the respective reasonable mesh convergence value.

Fig. 1 .Fig. 2 .
Fig. 1.The numerical generation process of (a) octet lattice, (b) Schwartz P shell and (c) spinodal shell-reinforced composites.The magnified view of a beam intersection is highlighted in the inset in (a), where the red dot represents identical beam end points.By contrast, all surface-based topologies (b and c) have no defined intersections.For illustration purposes, only 10%-reinforced composites are shown.Insets are not drawn to scale.

Fig. 3 .
Fig. 3. Illustration of the applied boundary conditions in the y-z and x-y plane.The triangles represent roller supports.

Fig. 4 .
Fig. 4. (a) Effective Young's modulus E of composites with four different implant topologies (represented by the markers) -octet, primitive, early-and lategeneration spinodal -with three different implant constituent materials (represented by colors) -Ti-6Al-4 V, 316L SS, and Co-Cr-Mo -at implant relative density q À r = 0.1, 0.2 and 0.3.The typical Young's modulus range for human long bones are bounded by the dashed lines [76-89].(b) The Mises stress field in both the encompassed tissue and the implant phases of the composites under 1% uniaxial compressive strain.The stress legend bar is limited by a stress limit r limit (chosen to allow a more visible comparison between topologies) such that all stresses above it are colored in red.Er and Em are the Young's modulus of the reinforcement and matrix constituent material, respectively.

Fig. 6 .
Fig. 6.(a) Relative Young's modulus E=Er versus minimum safety factor SF min ¼ ry=S peak (S peak is the Mises peak stress and ryis the yield strength of Ti-6Al-4 V) for octet, primitive, early-and late-generation spinodal implants with relative densities of 10%, 20%, and 30% (q À r = 0.1, 0.2, and 0.3).(b) Half-cut and zoomed-in views of Mises stress distribution in octet implants at q À r = 0.1, showing that the peak stress S peak is experienced in a significant fraction of the nodal volume (stress values higher than S peak are shown in red).See Appendix A for more details.

Fig. 7 .
Fig. 7. Bone growth index I growth in the tissue (the green matrix phase in each composite inset) surrounding octet lattice, Schwartz P shell and spinodal surfacebased implants versus effective Young's modulus E=Er.

( 1 )
The load-controlled stress field r l is obtained using the relation r d =r l ¼ r de =10 (all in MPa), where r d and r de are displacement-controlled stress field and displacementcontrolled effective stress that can be found from mechanical simulations in Section 2.3.(2)The elemental volume and the Mises stress at the centroid of each element are extracted as pairs.

Fig. A1 .
Fig. A1.Determination of 1% Mises stress threshold as the peak stress in a porous implant.The percentage of the porous implant volume is defined as the volume of partial implant regions that have a specific Mises stress value divided by the total volume of the implant.

Fig. A3 .
Fig. A3.The half-section view of the Mises stress distribution (limited by the 1%, 2%, and 3% Mises upper threshold r threshold ; the regions with stresses above r threshold are shown in red and the regions with stresses below r threshold are shown in blue) experienced in the q À r ¼ 0:3 porous implants with octet, primitive, early and late-generation spinodal topologies.

Table 1
The typical Young's modulus of human long bones.

Table 2
Summary of the performance metrics for all the implant topologies under investigation.