Inverse-designed spinodoid metamaterials

After a decade of periodic truss-, plate-, and shell-based architectures having dominated the design of metamaterials, we introduce the non-periodic class of spinodoid topologies. Inspired by natural self-assembly processes, spinodoid metamaterials are a close approximation of microstructures observed during spinodal phase separation. Their theoretical parametrization is so intriguingly simple that one can bypass costly phase-field simulations and obtain a rich and seamlessly tunable property space. Counter-intuitively, breaking with the periodicity of classical metamaterials is the enabling factor to the large property space and the ability to introduce seamless functional grading. We introduce an efficient and robust machine learning technique for the inverse design of (meta-)materials which, when applied to spinodoid topologies, enables us to generate uniform and functionally graded cellular mechanical metamaterials with tailored direction-dependent (anisotropic) stiffness and density. We specifically present biomimetic artificial bone architectures that not only reproduce the properties of trabecular bone accurately but also even geometrically resemble natural bone.


INTRODUCTION
Tailoring the architecture of cellular materialsincluding random foams as well as deterministic truss-, plate-, and shell-based latticeshas produced a wide variety of metamaterials (also referred to as architected materials) whose macroscale physical and mechanical properties can be controlled by a careful design at the microstructural level.Supported by optimization techniques [1][2][3][4][5][6][7] and advances in additive manufacturing 8 , truss-and plate-based mechanical metamaterials [9][10][11][12][13][14][15] have gained prominence as, e.g., lightweight cellular solids with engineered stiffness, strength, or energy absorption.Unfortunately, they also suffer from detrimental stress concentrations found at all intersections of structural members, which generally leads to low strength and poor recoverability [16][17][18] .Metamaterials based on smooth topologies such as triply periodic minimal surfaces (TPMS)address this shortcoming 19 by avoiding points of stress concentration while also showing excellent scaling of stiffness and strength with respect to density [20][21][22] (leveraging the benefits of doubly curved surfaces 23 that engage slender structures primarily in stretching rather than bending 24,25 ).Most recently, spinodal topologiesextensively studied in the 1960s-1990shave found a revival in metamaterials [26][27][28][29] .They are naturally created by self-assembly processes such as spinodal decomposition during phase separation 30,31 in, e.g., nanoporous metal foams [32][33][34] , microemulsions 35,36 , and polymer blends 37,38 .Similar to TPMS, spinodal topologies consist of smooth, non-intersecting surfaces having nearly zero mean curvature 29,39 with unique topological properties 26,28,29,[39][40][41][42] .Spinodal topologies can be realized as solid networks (obtained from removing one of the two phases after phase separation, such as in nanoporous foams 26,29 ) and as shell networks (retaining only the interfaces between the two phases 27,28 ).Overall, spinodal topologies cover an extreme and seamless range of anisotropic (direction-dependent) mechanical properties, as verified, e.g., by their tailorable anisotropic elastic moduli 26 creating materials that are stiff in some directions and soft in others.By contrast to unit-cell-based truss-, plate-, or TPMS-type architectures, spinodal designs are non-periodic, which greatly enhances the design space and achievable directionality.Moreover, non-periodic architectures do not suffer from the fabrication-based introduction of symmetry-breaking defectsby contrast to periodic metamaterials whose sensitivity to defects deteriorates their mechanical properties 22,27,28,43,44 .As we further demonstrate here, the non-periodicity of spinodal networks (without geometric unit cell limitations or tessellation requirements)counter-intuitivelyenables the effortless creation of functionally graded materials with spatially varying properties.Finally, the emergence of spinodal topologies during selfassembly (such as by arresting spinodal decomposition in polymeric emulsions 27 ) promises scalable manufacturing of centimeter-scale samples with nanoscale features, in contrast to conventional additive manufacturing-based designs.
While spinodal topologies hence promise an intriguingly large design space as well as advantages over most competing architectures, key challenges persist.It is practically impossible to theoretically explore the full anisotropic design space, since topologies are obtained by simulating the time-dependent Cahn-Hilliard phase separation process 26,28 (see Fig. 1a), which can take hours per simulation on a modern computer.The only reported alternative 29 is limited to isotropy, which strongly limits the achievable property space.To overcome this limitation, we here deploy a cellular design strategy that bypasses expensive simulations and leverages the statistics of spinodal topologies, resulting in a class of metamaterial architectures, which we term spinodoid (spinodal-like) topologies, with enormous anisotropic design and property spaces.Exploring those, however, requires another approach for the inverse design to be reported here.
Structure-property relations of all existing metamaterials have primarily been explored in a forward fashion: given a microstructure, one extracts the effective properties by methods of homogenization.The inverse challengeidentifying a microstructural topology that meets the mechanical property requirementshas often been addressed by inefficient trial and error, requiring a designer's intuitive understanding of the structure-property relations.Systematic approaches such as topology optimization and genetic algorithms 2,45,46 are not only beneficial but also computationally expensive (relying on repeated sampling and/or computation of the effective properties) and highly dependent on initial guesses.In addition, while the forward problem, i.e., mapping from topological parameters to property space, is well defined, the inverse problem is ill-posed (multiple topologies can have the exact same or similar effective properties).Recent advances in machine learning and deep learning have shown success in overcoming this challenge.Machine learning-based surrogate modelstrained in a data-driven 47,48 (supervised) or physics-driven [49][50][51][52] (unsupervised) setting, which can bypass expensive simulations and/or experimentshave been developed for applications, including metamaterials 3,53 , composites 54,55 , two-dimensional materials, and nanotubes 56,57 .More recently, surrogate models based on autoencoders and variational autoencoders 58,59 have attracted interest due to their ability to obtain a low-dimensional latent space parametrization of the design space.Once trained, surrogate models accelerate the search for possible candidates that fit the design requirements from a large design space.However, such approaches require solving an optimization problem based on the (forward-only) surrogate model, which prevents an on-demand inverse design framework.
We here deploy a machine learning strategy for the inverse design of the newly introduced spinodoid topologies, which efficiently predicts an optimal topology for a given set of sought properties.Our approach, in contrast to surrogate optimization methods, provides a computationally inexpensive two-way relationship between topological parameters and mechanical properties.Although we focus on the anisotropic elasticity of spinodoid architectures (and demonstrate their potential for, e.g., spatially varying patient-specific bone replacements), our inverse design approachinspired by similar problems in the nanophotonic and plasmonic community [60][61][62][63][64][65] is sufficiently general to apply, in principle, to any physical material properties and to any finite set of design parameters.
Fig. 1 Design of spinodoid metamaterials.a Conventional methodology of generating spinodal topologies by computationally expensive phase field simulations in combination with a level set.The anisotropic lamellar topology shown here was obtained by such simulations with anisotropic surface energies (see ref. 26 for details).b Schematic view of the design parameters (cone angles) θ 1 , θ 2 , and θ 3 in our anisotropic GRF approach.Black dots indicate the N wave vectors n i sampled from within the cones.c Spinodoid topology generation strategy using the GRF approach, followed by computational homogenization for effective property extraction.Here the example of an isotropic GRF (θ 1 = θ 2 = 0°, θ 3 = 90°) and its corresponding topology (for ρ = 0.5) and elastic surface (obtained by finite element simulations) are shown.The elastic surface illustrates the effective (normalized) directional Young's modulus E(d) for all directions d ∈ S 2 in the Cartesian basis fê 1 ; ê2 ; ê3 g.Indicated in light and dark gray are the elastic Voigt bound (ρE s ) and the Hashin-Shtrikman upper bound, respectively, for comparison.
In the following, we first introduce the class of spinodoid topologies, followed by a discussion of the tunable anisotropic material behavior with opportunities for functionally graded spinodoid metamaterials.Next, we present the data-driven inverse design framework, whose application to spinodoid metamaterials results in topologies with as-designed anisotropic stiffness.We close by demonstrating the (inverse) design of synthetic bones based on spinodoids.

RESULTS
From spinodal to spinodoid topologies with tunable anisotropy Spinodal topologies form naturally when a homogeneous solution decomposes into two spatially separated uniform phases in a diffusion-driven fashion.During the early stages of phase separation, small fluctuations in the phase distribution are described accurately by the linear Cahn-Hilliard model 66 for an isotropic evolution of the underlying phase field φ(x) representing the concentration fluctuation of one phase at position x ∈ Ω in a domain Ω & R 3 (Fig. 1a).Cahn 30 showed that the phase field may be described by a superposition of a large number N ≫ 1 of standing waves of constant wavenumber β > 0, mathematically known as a Gaussian random field (GRF) 67 : (1) where S 2 ¼ fk 2 R 3 : k k k ¼ 1g denotes the unit sphere in threedimensional (3D).n i and γ i denote, respectively, the directions and phase angles of the i th wave vector, which are random variables independently sampled from uniform probability distributions.While the original Cahn-Hilliard equation applies only to isotropic systems, it has been exploited for simulating anisotropic phase separation [68][69][70][71][72][73][74][75] by introducing directional-dependent interface energy or phase mobility/diffusivity, both of which result in cellular architectures with anisotropic properties.Toward the same objective but without the need to simulate the phase separation process, we here propose an anisotropic extension of Cahn's GRFbased solution (1) as a means to efficiently generate smooth spinodal-like topologies.
Our point of departure is the observation that the wave vectors in (1) are isotropically sampled from the unit sphere.To introduce anisotropy, we restrict those n i directions by biasing their statistical sampling with a non-uniform orientation distribution function (ODF), favoring some directions and neglecting others.Since the resulting topologies approximate the products of spinodal decomposition, we speak of spinodoid topologies.As a particular example, inspired by crystallographic texture poles, we introduce the distribution (2) where fê 1 ; ê2 ; ê3 g denotes the Cartesian basis and angles θ 1 ; θ 2 ; θ 3 2 f0g ∪ ½θ min ; π=2 represent the allowed spread of wave vectors about each of the three orthogonal base directions (Fig. 1b).Choosing, e.g., θ 1 = θ 2 = 0 restricts wave vectors to the cone, which forms an angle less than θ 3 > 0 with the x 3 -axis (and the interfaces of the topologies align preferentially perpendicular to those n i vectors).Controlling the value of θ 3 allows a seamless transition from a phase field with an isotropic structure (θ 3 = π/2, reducing to the unit sphere) to that with a lamellar structure (θ 3 < π/2).As representative examples, Figs 1c and 2 show specific choices of θ 1 ; θ 2 ; θ 3 ð Þ and the resulting isotropic, lamellar, columnar, and cubic symmetries.
From the above phase field, a bicontinuous topology is constructed by computing level sets or isosurfaces of the phase field 29,30,40,41 .To this end, we define a binary indicator function χ(x) that denotes the presence of material vs. void at x.This applies equally to the generation of solid-and shell-type architectures 27,28 : For solid networks, threshold φ 0 is computed as the quantile (of normally distributed random variable φ) evaluated at the average relative density ρ ¼ E½χ of the solid phase: φ 0 ¼ ffiffi ffi 2 p erf À1 ð2ρ À 1Þ.The relative density ρ provides a measure of the porosity of a given topology.Since small relative densities ρ ≪ 1 may contain disjoint solid domains, we restrict the design space to ρ ∈ [0.3, 1] in the following (for the same reason, we choose the minimum cone angle as θ min ¼ 15 ); for isotropic topologies, ρ > 0.159 was recently shown to ensure bi-continuity 29 .For shell topologies, the zero isosurface (φ 0 = 0) is chosen.We note that the area per unit volume of the isosurface of stochastic fields of type ( 1) is constant regardless of the anisotropy 76 , and the area per unit volume of the zero isosurface of an isotropic GRF bounded by a unit cell is 40,41 so that the relative density of all slender shell topologies is obtained as ρ = sh with shell thickness h (note that the relative density of shell architectures can be orders of magnitude below that of their solid counterparts).The complete set of design parameters to characterize our anisotropic spinodoid topologies is hence Θ = (ρ, θ 1 , θ 2 , θ 3 ).We stress that the thus generated topologies are not mere mathematical constructs but have physical relevance: they may be interpreted as an approximate solution to the modified Cahn-Hilliard equation modeling spinodal decomposition in case of an anisotropic diffusive mobility (see Supplementary Information, Section 1).It is for this reason that the spinodoid architectures presented here may indeed be expected from phase separation with direction-dependent interface mobility.While the self-assembly based fabrication has been demonstrated only for isotropic structures 27 , phase separation in block copolymers and polymeric micro-emulsions is a promising avenue for the selfassembly of anisotropic spinodal topologies [36][37][38]77,78 .
The above architectures are generally non-periodic, which is expected from foams but is a departure from conventional metamaterial designs and comes with advantages.First, nonperiodic topologies are more resilient to fabrication-related, symmetry-breaking imperfections 27,28 .The lack of symmetry also affects the buckling behavior and suppresses bifurcation modes of intermediate wavelengths (spanning single to multiple unit cells) often exploited in periodic metamaterials [79][80][81][82][83] .Second, despite the complex architecture (applying to an, in principle, infinite design space without periodicity or symmetries), the mathematical description is simple and can be encoded in an ODF and parametrized by a finite set of design parameters (e.g., Θ).Third, the topology can be varied smoothly and seamlessly across a given body for the purpose of functional grading without issues arising from unit cell tessellation and discontinuous topologies (which, in contrast, is an active challenge for periodic topologies, including trusses, plates, and shells).Although functional grading with smooth transitions was demonstrated for TPMS topologies [84][85][86] , the triple-periodicity constraint results in a limited and discrete design space (including the well-known Schwarz Primitive 87 and Schoen Gyroid 88 ) and included limited examples of anisotropy in TPMS topologies [89][90][91] .A unified continuous design space, particularly in terms of a tunable anisotropy of the resultant metamaterial, has remained a challenge.Here a spatially variant spinodoid topology can be generated by the spatial superposition of two or more random fields (Fig. 3 and Supplementary Information, Section 4), which offers unprecedented opportunities for spatially variant architectures.In addition to locally varying the topology, this approach can also create architectures with spatially variant microstructural length scales (e.g., to vary locally the characteristic pore size); an example is shown in Supplementary Information, Section 4. Finally, the non-periodic field-based formulation (1) is space-filling and can be leveraged to create spinodoid solids of arbitrary macroscopic shapes (without the need for adaptive tessellation and conforming algorithms near the body's boundary, as is the case for most periodic unit cell-based designs).
Tunable anisotropic material behavior As a representative example of effective properties, we characterize the 3D elastic stiffness of solid topologies (shell topologies can be treated analogously, which have great potential due to their optimal scaling of stiffness and strength with density 28 and their extreme resilience verified experimentally in ceramic thin-shell architectures 27 ).We apply computational homogenization by the finite element method (FEM), using a representative volume element (RVE) for each chosen topology made of a homogeneous, isotropic, linear elastic base material (Young's modulus E s , Poisson's ratio ν s = 0.3).The resulting effective (homogenized) constitutive behavior of the metamaterial is linear elastic with an effective fourth-order elastic modulus tensor C extracted from the RVE via applying average strains 〈ϵ〉, computing volume-averaged stresses 〈σ〉, and solving hσ ij i ¼ C ijkl hϵ kl i.Since the topologies lack periodicity, we apply affine displacement boundary conditions to the RVE, so that the homogenized stiffness is generally an overestimate of the true stiffness (for sufficiently large RVEs relative to microstructural features, assuring statistical homogeneity and a separation of scales, this estimate presents a good measure of the effective stiffness of the metamaterial).Numerical details are contained in Supplementary Information, Section 2. The microstructural length scale is inversely proportional to the wavenumber β in (1).For a cubic RVE of size l × l × l, we found β = 10π/l to be sufficient for an effective separation of scales; similar observations were reported previously for spinodal topologies 28,29 .
To visualize the elastic anisotropy, the elastic surface is computed, showing the effective Young's modulus along a direction d ∈ S 2 as EðdÞ ¼ Figure 2 illustrates how design parameters (θ 1 , θ 2 , θ 3 ) control the elastic surface.In addition, the relative density ρ controls the absolute stiffness scaling (see Supplementary Fig. 2).Note that the principal stiffness directions here are by choice always aligned with coordinates fê 1 ; ê2 ; ê3 g.However, a coordinate rotation (i.e., Ĉpqrs ¼ P i;j;k;l R pi R qj R rk R sl C ijkl for a rotation tensor R ∈ SO(3)) can further expand the design space to obtain stiffness tensors with arbitrary principal directions.
All information related to the structure of the elastic surface of a metamaterial is now encoded in the effective stiffness matrix C, whose componentsby exploiting symmetries and the orthotropic nature of all architectures defined by Eq. ( 2)reduces to the following vector of independent elastic moduli:

Data-driven inverse design
Toward the creation of cellular solids with as-designed properties, we need to address the inverse design question: how can we systematically and efficiently find a topology from within the design space that has the target anisotropic elastic moduli?In contrast to existing methods, we use offline training of a machine learning model that provides a computationally inexpensive twoway relationship between topological parameters and mechanical properties.We propose a machine learning technique based on deep neural networks (NNs) 92 , which require the a-priori creation of a sufficiently large and representative training dataset D ¼ fΘ i ; S i g; i ¼ 1; ; n f g consisting of n pairs of design parameters Θ and corresponding stiffness S, which are computed by FEM homogenization.Note that a large dataset is made possible only because of the computationally inexpensive GRF-based formulation of the spinodoid topology.Creating a dataset by simulating the time-dependent Cahn-Hilliard phase separation process is impractical as each simulation can take several hours on a modern computer (in contrast to few seconds for the GRF-based formulation; see Supplementary Information, Section 5).
Let F ω be a NN that maps design parameters Θ onto stiffness S in a forward fashion and is hence referred to as f-NN.We choose a multi-layer perceptron (MLP) architecture (Fig. 4a) whose parameter set ω contains the weights and biases of all hidden layers (see Supplementary Information, Section 5).Training the f-NN requires minimizing the (mean squared error) loss between true values and predictions with respect to the f-NN-parameters ω, i.e., This problem is well posed, and we use a back-propagation algorithm 93 to perform the optimization in Eq. ( 6).We leverage automatic differentiation 94 (which is also the implementation basis of back-propagation in most modern machine learning packages) to compute the gradients F 0 ½Θ i ¼ ∂F ω ½Θ i =∂Θ i (the latter will be important for the inverse problem).
To address the inverse problem, we introduce another MLP NN G τ (the i-NN), which maps stiffness S onto design parameters Θ.
Here the challenge lies in defining a measure (a loss function) that identifies the correctness of an answer while ignoring the multiplicity of correct answers (since various design parameters may yield the same or similar effective stiffness).As a trivial example, the spinodoid topologies corresponding to Θ = (ρ, 90°, 0, 0), (ρ, 90°, 90°, 0), and (ρ, 90°, 90°, 90°) all have the same (isotropic) stiffness.This prevents the straightforward application of machine learning tools, which rely on direct computation of errors.For example, the naive approach to minimize Fig. 3 Functional grading.A spatially variant architecture interpolating between columnar (θ 1 = 0°, θ 2 = θ 3 = 30°) and lamellar (θ 1 = 30°, θ 2 = θ 3 = 0°) topologies (indicated by orange and blue colors, respectively).The cellular solid is generated by the superposition of the corresponding anisotropic GRFs φ 1 and φ 2 in a linearly graded fashion: with λ(x 1 ) ∈ [0, 1] (see Supplementary Fig. 4); the elastic surfaces corresponding to the two GRFs are shown with their design parameters.The level set φ 0 is set to zero for a homogeneous relative density of ρ = 0.5.As a consequence of the heterogeneous stiffness distribution, loads are carried and deformation is observed inhomogeneously across the sample (as demonstrated in Supplementary Information, Section 4 for axial compression).
ill-posed and may not converge to a correct solution (see Supplementary Information, Section 7).To overcome this challenge, we leverage the f-NN and propose to train the i-NN using the loss function The reconstruction loss computes the error between the stiffness of a predicted topology and the true stiffness that has been queried, the rationale being that the predicted topology must have the correct stiffness irrespective of the design parameters used to create that topology.Ideally, the reconstruction should be computed using FEM.For efficiency, we instead leverage the above f-NN as a computationally inexpensive approximator of the FEM homogenization scheme.The advantage is two-fold.First, the numerous evaluations of the reconstruction loss via the pre-trained f-NN are several orders of magnitude less expensive than if FEM was used for reconstruction (see Supplementary Information, Section 5) during the i-NN training.Second, the backpropagation algorithm 93 requires computing the derivatives of the loss function with respect to the NN parameters (so-called sensitivities).Numerical differentiation by FEM (perturbing the design parameters and re-calculating the effective stiffness for each perturbation) would incur prohibitive expensesthe f-NN as an approximator instead provides analytical gradients via automatic differentiation 94 (see Supplementary Information, Section 5).Note that during the i-NN training, the pre-trained f-NN remains unchanged; i.e., no further do we train the parameters ω of the f-NN.The prediction loss in Eq. ( 7) acts as a soft regularization by measuring the error in the predicted design parameters.Although this regularization is not necessary (equivalent to choosing λ = 0), we observe that it accelerates the training of i-NN and the convergence of the reconstruction loss as long as λ is sufficiently small (e.g., λ ~0.01-1.0).For example, in the initial epochs of training when the parameters of the i-NN are untrained, the prediction loss ensures that the design parameter predictions from the i-NN are not nonsensical.For this reason, we activate the prediction loss only during the first few epochs of the training stage before deactivating it again (details in Supplementary Information, Section 5).
For a quantitative assessment of our machine learning technique, we generated a set of training and test data containing, respectively, 19,170 and 2130 pairs of topologies and their corresponding effective elastic stiffnesses (details about datasets and training protocols are in Supplementary Information, Section 5).For each queried stiffness S i from the test dataset, the i-NN predicts a set of design parameters Θ Ã i expected to yield the sought effective stiffness S i .For verification, we reconstruct the stiffness S Ã i of the identified design Θ Ã i (i) exactly via FEM simulations and (ii) approximately using the f-NN (recall that we cannot compare Θ Ã i to Θ i from the test data, since the mapping between design parameters and stiffnesses is non-unique).Ideally, in both cases the reconstructed stiffness of the identified topology should agree with the target stiffness, i.e., S Ã i % S i .In a plot of predicted (or reconstructed) vs. queried stiffness components, we hence expect each data point to ideally lie on a line with zero intercept and unit slope, and we may use the coefficient of determination (R 2 ) with respect to the aforementioned line as a measure of accuracy.
Representative example results of reconstruction vs. true (queried) stiffness component C 1111 are shown in Fig. 4b, c.Since the i-NN has been trained against the f-NN (which only serves as an approximation of FEM-based homogenization), the FEM reconstruction of stiffness shows a lower accuracy (R 2 = 0.997) than when using f-NN (R 2 = 0.999), as expected.Nevertheless, we generally observed excellent agreement of queried (true) vs. achieved (reconstructed) stiffness values across all stiffness components tested (we verified R 2 ≥ 0.995 and 0.998, respectively, across all FEM-and f-NN-based reconstructions).
As expected, the predicted topological parameters Θ * for a queried stiffness S may vary significantly from those in the dataset due to non-uniqueness (Fig. 4e).When queried, e.g., with a stiffness S originally obtained from a sample with Θ = (0.35, 70°, 70°, 0°), the i-NN predicted an architecture with Θ * = (0.34, 47.99°, 55.88°, 31.79°) with only negligible differences in stiffness.In addition, for close-to-isotropic stiffness queries (obtained from a sample with e.g., at least one of θ 1 , θ 2 , or θ 3 >80°), the i-NN predicts all three θ 1 , θ 2 , and θ 3 in the range of 45°-65°( contributing to the decrease in prediction accuracy for high values of θ 1 in Fig. 4e).This highlights the advantages of our approach in overcoming the ill-posedness of the inverse problem.Note that the stiffness dependence on relative density ρ is so dominant that the i-NN generally recovers the relative density of architectures from the training set for a given stiffness (Fig. 4d).
While the inverse model is successful on test queries similar to the training dataset, a natural question arises: is the model applicable to arbitrary stiffness queries?There are two requirements: (i) a queried stiffness must be thermodynamically admissible, i.e., it must satisfy major and minor symmetries (C ijkl ¼ C klij ¼ C ijlk ; i; j; k; l ¼ 1; 2; 3) and strong ellipticity for reasons of stability (using Einstein's summation convention: C ijkl u i v j u k v l > 0 for all non-zero u; v 2 R 3 ).Furthermore, it must not violate the Voigt upper bound that limits the maximally achievable stiffness for a given relative density (e.g., one cannot ask for an architecture that is stiffer than the base material).(ii) The query S must lie within the design space (i.e., there must be at least one Θ * within the bounds of the design parameters whose stiffness matches the query S).If these requirements are not satisfied, the i-NN will likely predict a Θ * , which does not lie within the design space (an indication that the S does not satisfy the aforementioned conditions).In the unlikely event Θ * does lie within the design space, a quick reconstruction of S * via f-NN (and comparison with the query S) will invalidate the i-NN prediction.
Once trained, the f-NN and i-NN together provide a computationally inexpensive two-way structure-property map as a design tool for spinodoid metamaterials (Supplementary Information, Section 5 estimates computational costs).Moreover, the machine learning framework presented here can be readily integrated into the functional grading approach of Fig. 3 to inverse-design spatially variant solids.The i-NN can be used to generate multiple independent GRFs that locally satisfy the anisotropic stiffness requirements (queries to the i-NN), followed by the straightforward GRF interpolation approach outlined in Supplementary Information, Section 4.
Inverse design of artificial bone Synthetic bones or bone-mimetic scaffolds and implants are prime examples that have benefited from advances in additive manufacturing [95][96][97] .Matching the topological and mechanical properties of bones in 3D-printed architectures (important for successful long-term compatibility 98 ) has remained a challenge though.Although bone properties have been mimicked by truss and TPMS architectures 95,98,99 , the thus available design space is highly limited.For example, those designs did not cover the required high level of elastic anisotropy and heterogeneity as well as porosity found 100 across different patients or even across bones within the same patient.Leveraging the approach introduced here can overcome those challenges and create bone-mimetic structures with properties closer to those of natural bone than any of the aforementioned designs.
As an example, we consider trabecular bone from bovine femur samples, whose relative density and anisotropic stiffness components were measured experimentally 101 (see Supplementary Information, Section 5).Using the measured directional modulus variations as a query (and assuming the base material of bone tissue is isotropic), the i-NN predicts the spinodoid architectures and their stiffnesses and relative densities summarized in Fig. 5. Remarkably, the inverse model accurately predicts the target relative density and matches the anisotropic elastic stiffness of the bone specimens despite no prior information about trabecular bone during the learning stage.The small differences between the reconstructed and experimentally measured stiffnesses are attributed to small but non-zero normal-shear coupling terms in the measured stiffness tensor (deviations diminish when the corresponding stiffness components are set to zero).Moreover, unlike previous approaches 95,101 , the resulting columnar spinodoid architectures bear topological resemblance to natural bone specimens (micro-computed tomographic images of the original bones are included in Fig. 5).The columnar topology makes structural features align with the load-bearing directions, akin to the trabecular alignment in natural femoral bones.

DISCUSSION
Bone datasets are scarce, particularly from a data-driven inversedesign perspective.Our approach bypasses this challenge and applies an inverse model trained on spinodoid metamaterials to the in silico generation of artificial bone, preserving the anisotropic stiffness and level of porosity (while the characteristic length scale is controlled by the wavelength).If a small but representative dataset for trabecular bones becomes available, the accuracy of the inverse model may be improved even further by first training on the spinodoid design space followed by finetuning on bone data.Finally, spatially varying architectures matching patient-specific mechanical properties across an individual bone (if known) can be realized readily by the GRF interpolation approach of Fig. 3.This capability of creating cellular solids with as-designed spatially varying stiffness and density is key to, e.g., acoustic cloaking 102 or tissue engineering 103 .We close by noting that, with possible integration into multiscale topology optimization or as a standalone framework to explore the design space (e.g., via genetic algorithms), our combination of forward and inverse maps of the structure-property relation accelerates the design process of metamaterials with a wide range of tunable mechanical response (anisotropic stiffness only being the tip of the iceberg).

METHODS
Details of the anisotropic spinodal decomposition theory (Section 1), the simulation procedures and post-processing of data (Sections 2 and 3), the spatially variant architecture simulations (Section 4), the machine learning Fig. 5 Inverse design of bone-mimetic topologies.Measured anisotropic stiffnesses of bovine femoral specimens #2 and #3 from Colabella et al. 101 (see Supplementary Table 1) were used as queries; shown are micro-computer tomographic (CT) scans of both bone samples as well as the predicted topologies (with β = 15π/l for visual comparison) and the corresponding design parameters Θ. Elastic surfaces of the bone specimens (dashed lines) and the predicted topologies (solid lines) are illustrated via projections onto the ê1 -ê 3 -, ê2 -ê 3 -, and ê1 -ê 2 -planes.Also shown are adjusted elastic surfaces of both the bone samples (obtained from setting all off-diagonal normal-shear and shear-shear coupling terms to zero in the measured stiffness tensors).Micro-CT images are adapted from ref. 101 with permission from Springer Nature.

Fig. 2
Fig. 2 Anisotropic spinodoid topologies and their properties.a-f Example solid topologies and their corresponding elastic surfaces, demonstrating the wide range of anisotropy controlled by the design parameters (θ 1 , θ 2 , θ 3 ) for ρ = 0.5.Elastic surfaces illustrate the effective (normalized) directional Young's modulus E(d) for all directions d ∈ S 2 in the Cartesian basis fê 1 ; ê2 ; ê3 g.Indicated in gray is the elastic Voigt bound (ρE s ) as an upper bound.

Fig. 4
Fig. 4 Data-driven inverse design.a Schematic of the inverse design process.The inverse model (i-NN) takes a queried stiffness S as input and outputs design parameters Θ defining a topology according to Eq. (1).The predicted design parameters are fed to the forward model (f-NN) to reconstruct the stiffness and verify the predictions from the i-NN.Both NN models consist of six hidden layers with the indicated numbers of nodes.b, c Reconstructed vs. true stiffness component C 1111 in the test dataset.The reconstructed stiffness (i.e., stiffness of the topologies predicted by the i-NN) is computed using f-NN and FEM.d, e Predicted vs. true design parameters ρ and θ 1 in the test dataset.All dashed lines represent the ideal line with zero intercept and unit slope; the corresponding R 2 deviations are indicated.Further details are provided in Supplementary Information, Section 5.