Boosted by additive manufacturing, architected materials have opened new opportunities to extend the perfor- mance of engineering materials. Yet, their development is held back by the intense efforts required to understand their com- plex property -structure -process -performance relationship. Therefore, data -driven biomimetic approac hes are becoming in- creasingly popular to unveil such relationships . Here we mimic the functionally graded structures found in Coscinodiscus sp. diatom to understand the role of their shapes and define new guidelines for the design of novel architected honeycombs with tunable mechanical properties. Finite element simulations, validated on the outcome of a testing campaign performed on 3D- printed elastomeric samples, are used to build a dataset for machine learning algorithm training. Different machine learning techniques are used to link the geometric features of the designed biomimetic structures to their energy absorption proper- ties and, in particular , to the specific absorbed energy divided by the peak force , here used as the performance index . The proposed approach leads to a novel design, which features a performance increase of 250% w .r.t. conventional honeycomb s. Architected materials are a class of materials in which chemical and geometrical features are carefully controlled on different length scales to ensure the optimal accomplish- ment of specific tasks. Typical examples are honeycombs, strut -based lattices , and triply periodic minimal surface lat- tices1,2. Their power lies in the ability to tune their perfor- mance by acting on the hierarchical and mutual interplay between material and structure. Thus , they find application in a multitude of fields, such as lightweight structures, wave propagation control, bio -scaffolds, and sensing2–4. Advances related to architected materials have been mainly driven by the mimicry of natural materials like bone, wood, nacre, and diatoms5–11, whose remarkable capabili- ties and multifunctionality are mainly attributed to the so- phisticated multiscale hierarchical structure12,13, despite the limited palette of building blocks they are made of. In Nature, a few universal building blocks, combined into com- mon structural arrangements, functional gradients, and in- terfaces with recurring patterns over multiple length scales5,14 –20 lead to a wide diversity of structures and materials, able to meet different functions and needs. To- day, thanks to imaging techniques, numerical simulations, and additive manufacturing (AM), the application of the “universality -diversity” paradigm found in Nature21–23 to the design of architected materials has turned out to be a very fruitful guide23–27. Yet, unraveling the property -struc- ture -process -performance relationship28 and transferring it to a man -made material is generally a non -straightforward highly resource -intensive process. In recent years, data - driven approaches, especially machine learning (ML), have proven to be valuable methods to overcome this issue29–34. ML algorithms have been successfully used for the data - driven design of new materials by: (i) replacing constitutive equations for the analysis of new material properties, (ii) identifying new material architectures, (iii) guiding the choice of process p arameters, and (iv) improving perfor- mance such as fatigue life, buckling resistance, and fracture toughness29–34. Latest cutting -edge research has demon- strated that Materials Informatics tools at the forefront , such as attention -based diffusion models35 and multi -modal 2 language models36, are capable of both providing forward predictions ( i.e., obtaining nonlinear structural behaviors as a function of material hierarchical structure) and effectively solving inverse design problems ( i.e., discover new architec- tures that have specified response curves), giving increas- ing room for customization of material representations and offering ever -expanding generalization capabilities ( i.e., the same models can be applied to different materials, dimen- sional scales, and design tasks). While biomimicry has been successfully used to provide design solutions for materials with superior mechanical properties, ML -based methods have further accelerated the design process and the discov- ery of new materials37. Inspired by the fascinating diatom morphology38–42, in this research we propose a biomimetic data -driven ap- proach for the design of novel architected honeycombs with tunable energy absorption performances. What caught our attention on these aquatic algae is the organization of the amorphous biosilica that forms the exoskeleton, aka frus- tule. The frustule of Coscinodiscus sp. diatoms exhibits a hi- erarchical architected design, consisting of three main lay- ers: (i) cribrum, (ii) areolae, and (iii) foramen ( Figure 1). From a structural point of view, this protective shell is de- signed to withstand predatory attacks that result in compressive loads. It provides both lightness and mechani- cal protection at the same time, and the central honeycomb - like layer is the one that affects such performance the most41,43 –49. Previous AI diatom -based designs focused on a publicly available diatom dataset to translate their intri- guing structural features into 3D printable material models, showing the great potential offered by the intersection of human and artificial intelligenc e to generate innovative bi- ologically inspired design solutions50. Here, we start from solutions offered by Nature and investigate the relationship between their shapes and mechanical properties. Specifi- cally, an intriguing feature of Coscinodiscus diatoms is the cell wall profile of the areolae , which shows varying thick- ness gradients ( Figure 1 (b)), not observed in other n atural and synthetic honeycombs. This research work arises from the need to answer the following questions: what role do es it play in determining the overall properties of the frustule? Can we tune the shape of the cell walls to improve the me- chanical performance of bioinspired architected materials? The deformation of honeycombs under out -of-plane com- pression is mainly governed by buckling -induced collapse. Thus, this failure mode, typically detrimental, could be har- nessed to adapt the material’s energy absorption properties according to the specific application of interest. Figure 1. Hierarchical architecture of the studied diatom frustule, schematization, geometric, and numerical model. (a) Hierarchical architecture of the Coscinodiscus sp. diatom frustule with its four different levels: “Petri dish” con- figuration shared by all diat oms (modified from39, CC BY 4.0); sandwich -like structure distinctive for this family of diatoms (modified from51); intricate porous network of amorphous silica that covers the entire outer surface of the 3 frustule (modified from51); building blocks that make up the frustule (modified from52, CC BY 4.0). (b) Cross -section of the Coscinodiscus sp. diatom multilayer exoskeleton: view of the thickness gradient present in the cell walls of its honeycomb substructure (modified from39, CC BY 4.0). (c) Simplified geometric model of the sandwich -like structure of the Coscinodiscus sp. frustule with its three layers: cribrum, areolae, and foramen. (d) Schematic representation of the thickness gradient in the analyzed honeycomb -like struct ures (e) Example of the stress -strain response curve of elastomeric honeycombs under out -of-plane compression. Before the advent of AM , it was impossible to fabricate structures with such a level of complexity, so the effect of functional gradients in honeycomb walls has not been fully investigated yet. Previous studies have highlighted that the spatial tailoring of honeycomb geometry can substantially enhance their performance and be exploited for fine -tuning their properties. Regarding the in -plane mechanics of hon- eycombs, in particular, the use of deep learning models has enabled control over the material distr ibution in th eir cellu- lar walls through encoded patterns to generate innovative designs that provide designer -specified nonlinear constitu- tive relationships35,36. The out-of-plane behavior has been partially analyzed, focusing on simple material distribution laws 53,54. In this study, we leverage ML techniques, finite element analysis (FEA), and AM to deeply delve into the effect of thickness gradients on the out-of-plane honeycomb me- chanical performance. Given the diversity of shapes and gradients occurring in the Coscinodiscus sp. diatom frus- tule38,40,41 , we focus the study on a single hexagonal pris- matic cell. A large dataset of bioinspired digital models of diatom geometry is created. A subset is adopted to generate a numerical dataset describing the relationship between the cell shape and its mechanical performance ( i.e., the out -of- plane compression properties calculated by FEA). Based on simulation outcomes , supervised learning algorithms are trained for regression and the design space is extensively explored by using the best metamodel. A comprehensi ve parameter, based on critical load, absorbed energy, and density of each bioinspired honeycomb cell, is used for per- formance analysis and topology ranking. Finally, AM and testing are exploited for the validation of the proposed data - driven designs. This bioinspired ML -driven approach al- lows one to i dentify an ideal structure for protective gear applications, featuring a 250% boost in energy absorption capability compared to a classic honeycomb unit cell with constant wall thickness. Moreover, the proposed frame- work, based on ML and AM, represents a v ersatile, powerful, and sustainable approach, which accelerates the design phase and expands the design space by orders of magnitude compared to conventional methods, considerably reducing the resources, time, and costs involved in simulations, man- ufacturi ng, and testing. Numerical model validation . Driven by the growing in- terest in “multi -hit” devices for mechanical energy absorption55, we focus on the geometry of the honeycomb - like architectures found in Coscinodiscus sp. diatoms and use a thermoplastic polyurethane (TPU), instead of the nat- ural biosilica, as parent material. This choice is driven by the fact that (i) elastomeric mater ials are typically used in en- ergy absorption applications, (ii) TPU is easily available for AM technology. To have a reliable dataset for the ML algo- rithm training, the finite element (FE) model, used to gener- ate the synthetic data, is validated analytical ly and experi- mentally. The analytical formulae based on honeycomb me- chanics g ive us a first approximation, and then the numeri- cal model is validated on a single honeycomb unit cell with constant wall thickness taken as a benchmark. According to experimental test results of 3D -printed samples, the Blatz - Ko hyperelastic model is calibrated to replicate the non -lin- ear behavior of the material. Figure 2 (a) shows the comparison among data obtained from the three approaches used. The unit cell is analyzed without considering the periodicity of the architected mate- rial. Figure 2 (b-d) depicts the theoretical, simulated, and physically observed deformation modes. In the elastic re- gion there is an optimal match between the numerical and experimental curves, with a difference of 11.3% in terms of compression modulus. The experimental test results are very similar to the one predicted by the FEA and the analyt- ical formulae, showing an error of 2.1% and 0.8% on the maximum stress, respectively. In general, there is a good comparison in terms of elastic modulus of compression, peak stress, an d absorbed energy. A more quantitative com- parison is provided in Table 1. It should be noted that the calculation of elastic moduli is very sensitive to the selected strain range, so the errors re- lated to it are less meaningful than those associated with the other mechanical properties evaluated. Beyond the maxi- mum stress, the fini te element model becomes stiffer and deviates from what is measured experimentally. This differ- ence is primarily due to the anisotropies and imperfections introduced by the AM process in the testing specimens, and it also affects the structure deformation mode. The height of the benchmark cell is twice its side. From the theory of elas- tic stability we expect, in each of its walls, two whole defor- mational half -waves in the vertical direction and only one in the direction perpendicular to it56. The 3D printed cell (Figure 2 (d)) shows a shift of such deformation waves and shapes that are not perfectly harmonic w .r.t. the theoretical and numerical counterparts . 4 Figure 2. Finite element model validation. (a) Analytical, numerical, and experimental comparison of the response curve of the benchmark honeycomb cell with constant -thickness walls in out -of-plane compression. Cell deformation modes according to (b) honeycomb mec hanics, (c) numerical simulations based on FEA, and (d) experimental com- pression test. Table 1. Mechanical properties of the honeycomb cell with a constant -thickness wall calculated analytically, numer- ically, and experimentally. Elastic modulus is calculated between 0.01 and 0.03 strain. Percent errors are calculated with respect to experimental results (“ground truth”). Parameter Unit Experimental Numerical err % Analytical err % Elastic Modulus MPa 33.27 29.50 11.3% 24.74 25.7% Maximum stress MPa 2.42 2.37 2.1% 2.40 0.8% Specific absorbed energy J/mm3 0.432 0.45 3.6% 0.49 12.2% Metamodel Selection . After the FE model validation, a numerical dataset is built. This dataset, which includes about one tenth of all the generated structures populating the design space, is solved numerically via FEA to determine the mechanical behavior of such structures. This dataset is used for training, tuning, and testing three different super- vised machine learning algorithms to determine an optimal metamodel capable of predicting the maximum force F max (i.e., the critical load) and absorbed energy , EA, of each cell mod el within the design space, based on knowledge of the material distribution along their walls. Specifically, a Linear Ridge Regression (LRR), a Kernel Ridge Regression (KRR), and a Deep Neural Network (DNN) are used. Figure 3 shows the comparison between the distributions of the target me- chanical properties calculated numerically and predicted using ML for all the structures present in the test set. LRR is too simple to capture the nonlinear phenomena involved in the out -of-plane compression of elastomeric honeycomb cells. Better results are achieved with KRR and DNN with R2 scores of 0.964 and 0.981, respectively (see Table 2). Fur- thermore, it is worth noting that both KRR and DNN predict the same best -performing structure on the validation and test sets. Although the prediction accuracy for KRR and DNN are similar, with KRR being 250 times faster than DNN in terms of training ti me (0.3 seconds compared to 75 sec- onds), the DNN is selected to explore the remaining portion of the design space ( i.e., predicting the performance of structures not included in the dataset). This decision is mo- tivated by the DNN’s stronger correlation in distribution tails, meaning its ability to confidently predict the outliers. Furthermore, once trained, it is ~105 faster than FEA in cal- culating the honeycomb cell mechanical properties, thus it significantly facilitates the design space exploration pro- cess. Details regarding its tuning, performance, and the hardware resources employed for its utilization are de- scribed in the Supporting Information ( SI). Effect of thickness gradients on energy absorption properties. In energy absorption applications the choice of material is usually guided by (i) the amount of energy to be dissipated, (ii) the maximum stress associated with it, and (iii) the weight limits imposed by the product design57. To rank the analyzed structures, a performance parameter, based on the energy absorbed divided by the maximum force of each cell and by the relative density ρrel, is used. Fig- ure 4 (a) shows the results obtained for our design space. A significant variability in mechanical properties is observed from the graph. From the theory of constant -wall honey- combs, we expect an increase of the critical load according to the mass of the cell walls , thus their stiffness. Intuitively one might assume a monotonic trend between EA/F max and ρrel. However, such trend cannot be found in the results be- cause the shape variations, featured in the bioinspired cells, significantly affect both the buckling load and the energy 5 absorption mode following the collapse of the cell due to in- stability ( i.e., the post -buckling behavior). To understand the effect of thickness gradients, three structures, representative of this variability, are selected and compared with the benchmark model. Specifically, the honeycomb cell models with geometric tailoring are charac- terized by : (i) the best mechanical performance, (ii) the worst mechanical performance, and (iii) an average struc- ture representing the 50th percentile of the overall dataset. Figure 4 (b) depicts the wall shape of these models. Table 3 shows a quantitative comparison of their performance. For a detailed understanding of the energy absorption mode re- lated to the selected geometries, FEA and experiments are performed on them . Figure 4 (c-d) shows a comparison be- tween the stress -strain response curves of the best -per- forming and worst -performing cell and the stress -strain curve of the benchmark model. The stresses are normalized to the absolute maximum value. Comparing the numerical results with the experimental ones, their similarity first emerges, confirming the effective- ness of the approach used. On the other hand, analyzing the differences between the mechanical response of the indi- vidual structures, it can be seen that the honeycomb cell with the best performance is characterized by a low peak stress and a post -buckling deformation, which is almost perfectly plastic, ensuring energy absorption at a constant force. The worst honeycomb structure maximizes ene rgy storage predominantly in the linear elastic regime and ex- hibits a noticeable drop after buckling. These results are in line with our expectation . In fact, to efficiently absorb a given amount of energy, the maximum force must be mini- mized, and the behavior of the honeycomb cell must be as close as possible to an elastic -perfectly plastic behavior57. The outcome of the proposed framework shows how this can be effectively achieved through the honeycomb cell wall shape tailoring. Remarks . In this study, a data -driven biomimetic ap- proach is proposed to customize the mechanical properties of architected materials by controlling buckling and post - critical behavior . Inspired by honeycomb -like structures found in diatom exoskeletons, which exhibit variable mate- rial distributions along the out -of-plane direction, we ana- lyze the effect of their functional gradients on the energy ab- sorption properties of millimeter -sized elastomeric struc- tures that can be 3D printed. Despite the limited ran ge of thickness distributions examined, by tailoring the shape of honeycomb cell walls, a performance increase of 250%, compared to the benchmark model with a constant wall thickness, is achieved. In contrast to most numerical data - driven approaches, which typically require significant com- putational power for both dataset creation and metamodel- ing training, this result is obtained using a standard implicit numerical solver and a commercial laptop. Therefore, the adopted framework shows great potential for t he rapid and cost -effective design of innovative structures and materials for energy -absorbing applications such as helmets, sports impact protection, and packaging. 6 Figure 3. Absorbed energy and peak force distributions obtained from FEA analysis and predicted by three different machine learning algorithms on the test set. (a) Comparison between FEA and Linear Ridge Regression (LRR). (b) Comparison between FEA and Kernel Ridg e Regression (KRR). (c) Comparison between FEA and Deep Neural Net- work (DNN). Table 2. R2 score of the Kernel Ridge Regression (KRR) and the Deep Neural Network (DNN). The table includes the overall R2 calculated considering both the absorbed energy EA and the peak load F max, the R2 on the single metric on the 90th and 95th percentile to show the better tail reconstruction of the DNN in distribution tails. ML algorithm R2 score Overall 90th percentile 95th percentile EA Fmax EA Fmax KRR 0.964 0.770 0.572 0.711 0.519 DNN 0.981 0.842 0.762 0.700 0.653 7 Figure 4. Predictions of the Deep Neural Network (DNN) and comparison with the benchmark solution. (a) Mechani- cal performance of bioinspired structures contained in the design space predicted by the DNN. Highlighted in the diagram: the benchmark model and three st ructures representative of the property variability induced by thickness gradients. (b) Representation of the geometry of the analyzed models with a focus on their cell wall thickness pro- files. (c) Numerical and (d) experimental stress -strain curves of out -of-plane compression loading for the best, worst, and benchmark honeycomb cells. Table 3. DNN prediction of absorbed energy, critical load, and performance parameter EA/(ρ relFmax) of the bench- mark, best, average, and worst honeycomb cell; performance gain with respect to the benchmark structure. Honeycomb ρrel (-) EA (J) Fmax (N) EA/(ρ relFmax) (J/N) Performance gain Benchmark 0.33 7.88 1099.3 0.022 --- Best 0.19 1.56 149.5 0.055 254% Average 0.29 2.37 303.4 0.027 122% Worst 0.45 9.99 1697 0.013 60% Although there is no doubt that geometric gradients in di- atom exoskeletons are responsible for mechanical protec- tion, this approach d oes not identify a match between the profile of the optimized cell and that of the algae. This mis- match can be mainly associated with the following reasons: (i) we focus on ly on the thickness gradient effec t, neglecting other morphological aspects and using TPU as base material while diatom frustule is composed of amorphous biosilica; (ii) the frustule shapes are extremely complex and require a more sophisticated approach to be represented; (iii) the diatom exoskeleton is optimized to perform multiple func- tions at the same time and not only to absorb mechanical energy. To generalize this study, the base material can be as- sumed as a design variable and tailored according to the specific application requirements ( e.g., biocompatible met- als for biomedicine, biodegradable polymers for sustainable design, nanocomposites for smart capabilities) or to further understand the biological role of the frustule ( i.e., assuming silica -based ceramics). The interplay between different ma- terial failure modes and thickness gradients could provide additional insights to enhance honeycomb st ructural per- formance and broaden the range of functions they can fulfill. Also, alternative ways to describe the natural thickness gra- dients should be investigated. Since the problem of honey- comb out-of-plane elastic instability can be traced back to harmonic -type solutions56,58, a promising idea is to assume 8 a space of parametric functions ( e.g., using Fourier series) to represent the functional gradients and build an AI model that can fit them to a set of diatom images to : (i) discover the most recurrent material distribution patterns, (ii) build a thickness profile generator that faithfully represents bio- logical designs, and (iii) investigate the shape -property re- lationship in more detail. However, the success of such ap- proach is highly dependent on the availability and quality of diatom data, which is very l imited to date. A more versatile solution, on the other hand, is offered by language -based representations, as they apply to a variety of high -dimen- sional data and are independent of a priori knowledge of material design logic36. Nevertheless, we can still conclude that the thickness gra- dients featured in diatom cell walls allow one to improve the structural performance of the biological material. In- deed, one can notice the similarity of the natural cell shape with the shape of the average cell predicted by our frame- work, which already shows a 122% performance gain w.r.t. the benchmark honeycomb. It should also be noticed that periodic boundary conditions (PBCs) are not applied to the vertical edges of the honeycomb cells. Their abs ence does not allow us to recreate the mechanical behavior of the pe- riodic material from a single cell analysis. Yet, it allows the study to be generalized in that: (i) by optimizing the perfor- mance of the single cell, a global improvement in the perfor- mance of the corresponding architected material can also be expected; (ii) the results obtained, besides being exploit- able for the class of materials mentioned, can also be ex- tended for modeling columns and rods designed against buckling. There are several examples in the engineering field where a pointwise material design would improve dramatically the performance. A common example is given by prosthe- ses, where a local material tuning is crucial to provide a more homogenous stress and stra in distribution, reducing local concentrations. Another example is provided by sports protective equipment, such as helmets, pads, and other in- jury prevention devices ( e.g., body armor), where protec- tion and lightness must coexist. These generally mutually exclusive characteristics can be combined by a targeted op- timization of the structure, to have a mechanical response that at the same time limits the peak load and absorbs as much energy as possible to avoid severe injuries. Similar needs are encountered in the field of soft robotics. Espe- cially when dealing with fragile objects it is fundamental that grippers do not exceed certain load limits. This can be achieved by implementing optimized architected materials that can significantly deform after reaching the critical load but, thanks to the use of elastomeric polymers, recover their shape when the load is removed. To fully understand the shape -performance relationship in diatoms, a multidisciplinary approach is needed. As fu- ture developments we plan to expand the study by evaluat- ing the multifunctionality of diatom -inspired architected materials, including fluid dyn amical and optical perfor- mances in a multi -objective optimization framework. EXPERIMENTAL SECTION Biomimetic Design. We focus on the structural arrange- ment of the diatom exoskeleton frustule characteristic of the micrometer scale ( Figure 1) and reminiscent of sandwich panels. To characterize the effect of shape heter- ogeneities occurring within the areolae on the out -of-plane compression loading scenario, different geometries of unit cells are studied : (i) a geometric benchmark model with constant -thickness walls and (ii) ~18000 models with var- ying wall thickness gradients. The former provides a bench- mark with known mechanical behavior58,59 for comparison with the geometrically tailored structures. The dimensions are defined on the millimeter scale, while maintaining the same proportions as the natural structure. Specifically, the edge length side of the mean perimeter defining the section of the hexagonal prismatic cell, the height of the latter , and the thickness of its walls are 14 mm, 28 mm, and 2 mm, re- spectively. This scale -up is necessary to cope with the limi- tations imposed by the AM technology used. From this model, we derive the oth er bioinspired structures through an algorithm that randomly defines the material distribu- tions in the honeycomb walls, keeping the other geometric features unchanged ( Figure 1 (d)). To limit the number of possible design solutions, we divide the benchmark unit cell into seven equal partitions along its height, then assign each one a thickness value belonging to the following vector t = [1.0 1.5 2.0 2.5 3.0] mm . This range of variation has been chosen so that the benchmark model represents an average configuration. We also impose a constraint on the gradient to avoid unphysical configuration. This is achieved by limit- ing to two steps the possible increments between two adja- cent partitions. More details are available in the Supporting Information (SI ) file. Analytical Modeling. Honeycombs mechanical response under out -of-plane compression is characterized by three different stages (see Figure 1 (e)): (i) a first linear elastic re- gime; (ii) a region characterized by constant stresses known as the “plateau”, triggered by buckling cell failure; and (iii) a densification phase in which the material becomes in- creasingly stiff due to the filling of voids contained within it associated with cell failure58,59. In the case of an elastomeric material, this mechanical behavior is also influenced by vis- coelastic effects, thus on the loading rate. This study is lim- ited to a quasi -static loading scenario . The honeycomb unit cell is assumed to have an elastic -perfectly plastic behavior; thus, it is described by the elastic modulus and the maxi- mum force F max estimated from the linearized buckling the- ory58,59. The honeycomb out-of-plane elastic properties are functions of its relative density ρrel58,59: ρrel=Ahoney Aapp≅2 √3𝑇 𝐿=Ehoney Ebulk (1) where T and L represent the values of the wall thickness and the edge length, respectively, Ehoney and Ebulk are the Young’s moduli of the unit cell and of the bulk material, Ahoney and Aapp are the effective and the apparent cross -sec- tional area of the honeycomb cell. F max can be calculated as: 𝐹𝑚𝑎𝑥 =3K Ebulk (1-νbulk2) 𝑇2 𝐿 (2) where K is a constraint factor dependent on the boundary conditions applied to the plates constituting the honeycomb cell walls56 and νbulk is the Poisson’s coefficient of the parent material. The value of the maximum stress is derived from Fmax using the Aapp section. 9 Finite Element Analysis (FEA). All the modeled cells have the same average surface area. To minimize the com- putational effort, this area is discretized through four -node shell elements with an average size of 0.5 mm, then the thickness associated with each element is locally adapted through a custom macro. The Blatz -Ko60 hyperelastic mate- rial model, calibrated on the outcome of uniaxial tensile and compressive tests performed on 3D -printed samples, is adopted to describe the behavior of the TPU. The calibration phase is carried out using the commercial software PolymerFEM . Simulations are performed using an implicit approach, based on the Newton -Raphson method for calcu- lating the nonlinear solution. The bottom and top bases of each cell are placed in contact with two rigid surfaces (rep- resenting the compression plates of t he testing machine). The top surface is moved axially at a slow speed to induce quasi -static compression, while the bottom surface is kept fixed. Contact interactions are modeled using a pair -based definition using target segment and line -to-surface ele- ments in ANSYS. A friction coefficient of 1.2 is used to simu- late sliding between the rigid surfaces and the FE mesh. For each cell, a single load step is applied, and the vertical dis- placement of the upper base is ca lculated at each substep, obtaining the fo rce-displacement curve for large displace- ments. From the latter, the stress -strain relationship is de- rived, and the energy absorption properties are calculated. The number of substeps is adjusted by the program accord- ing to the solution gradient. Supervised Learning. The numerical dataset used for the ML algorithm tuning is divided into three parts: (i) train- ing set (1000 data), (ii) validation set (500 data), and (iii) test set (500 data). The training set is used to make the ML algorithms learn the parameters that d efine the regression models ( i.e., the link between the input thickness distribu- tion, x i, and the output mechanical properties, y i). The vali- dation and the test set are used to tune their hyperparame- ters and evaluate their goodness of approximation, respec- tively. The optimal metamodel is defined and then used to calculate the mechanical properties of the non -simulated geometric model s (about 16000). Figure 5 depicts the adopted workflow. To enhance the performance of ML algo- rithms, both the labels and features are normalized and cen- tered. These transformations ensure that all data points have a mean value of zero and are normalized to one. Also, the 2000 structures are selected randomly to avoid any bias. Three different supervised learning approaches are used and compared: • Linear Ridge Regression (LRR). Thanks to its simplicity and computational efficiency, a linear regression is adopted as initial approximation tool. For its imple- mentation, we use the related Scikit -learn package choosing a regularization parameter α=0.1 to prevent underfitting. • Kernel Ridge Regression (KRR). Nonlinear kernel algo- rithms allow capturing complex input -output relation- ships while maintaining high computational efficiency. Here we define a kernel regression model based on ra- dial basis functions by exploiting Scikit -learn. In this case, the regularization -related hyperparameter is op- timized by cross -validation, varying it on a logarithmic scale in the interval [10−6, 101] and choosing the opti- mal one (Figure S3). • Deep Neural Network (DNN). DNNs have no limits in terms of capability to approximate functions61, how- ever, compared to previous approaches they require significantly more computational effort for training and hyperparameters selection, due to the large number of parameters and hyperparameters to be tuned . We cre- ate our own DNN using Keras, a user -friendly front -end API for TensorFlow. About fifty different DNN architec- tures are analyzed, with different combinations of hy- perparameters ( i.e., number of hidden layers, number of neurons in each hidden layer, dropout fraction for each layer) to maximize the fitting. Detailed infor- mation is reported in the SI file. The selected DNN ar- chitecture consists of an input layer with seven n eu- rons (one for each thickness partition), an output layer with two neurons (one for each label: EA and F max), and five fully connected hidden layers. The first and last hidden layers have 32 and 128 neurons, respectively, while the three central layers have 256 neurons. We use the ReLU activation function for all the layers ex- cept for the output one, where we imp lement the ‘lin- ear’ activation function commonly employed for re- gression tasks. Finally, to minimize the loss function, we utilize the Adam optimi zation algorithm. Training the DNN and its parameters involves running a total of 50 epochs with a batch size of 200. For all three ML models, we evaluate the goodness of the fit using the R2 score, which represents the fraction of infor- mation reconstructed by the model. The algorithms are run using Google Colab . 3D Printing and Testing. The Fused Filament Fabrica- tion (FFF) technology 3D printer Original Prusa i3 MK2s is used to manufacture both the specimens needed for TPU mechanical characterization and the honeycomb cell sam- ples with geometric and non -geometric tailoring. A mini- mum num ber of three samples per test and geometry is con- sidered, to ensure repeatability of results. The process pa- rameters are provided in SI. For all the samples, the com- mercial filament FILOALFA® ALFA TPU hard is used. To cal- ibrate the hyperelastic model of TPU, both tensile and com- pression tests are carried out, according to ASTM D412 and ASTM D575 standards62,63, respectively. To experimentally assess the behavior of honeycomb cells, we carry out quasi - static compression tests in displacement control (with a crosshead speed v=2mm/min), using a ZwickRoell ProLine universal testing machine endowed with a 10 kN load cell. 10 Figure 5. Schematic of the workflow used to build the numerical dataset, train the machine learning algorithms, and predict the behavior of geometrically tailored honeycomb cells. The design space, consisting of approximately eight- een thousand different structures, is defined by creating seven -component vectors, x i, describing the wall thickness distribution. Two thousand structures are randomly selected and analyzed by FEA to obtain the 2D vector, y i, char- acterizing the energy absorption capacity of the structure: energy absorbed (EA) and peak load (F max). The meta- model built by supervised learning algorithms is trained, tuned, and tested on the dataset. The mechanical proper- ties of the geometric mod els, not included in the dataset, are predicted, and those with the optimal performance identified. ASSOCIATED CONTENT Supporting Information The Supporting Information (SI) file contains g eometric details of the bioinspired model used as a benchmark for numerical model validation ; 3D printing process parame- ters; additional information on m achine learning algo- rithms . This material is available free of charge via the In- ternet at http://pubs.acs.org AUTHOR INFORMATION Corresponding Author *Flavia Libonati - University of Genoa, Polytechnic School, Department of Mechanical, Energy, Management and Trans- portation Engineering, Via all’Opera Pia 15/A, 16145, Ge- nova, Italy; orcid.org/0000 -0001 -6490 -1922; Email: flavia.libonati@unige.it ACKNOWLEDGMENTS The authors acknowledge support from the University of Genoa, under the Curiosity Driven Starting Grant initiative. ABBREVIATIONS AM, Additive Manufacturing ; ML, Machine Learning ; FEA , Finite Element Analysis ; TPU , thermoplastic polyurethane ; LRR, Lin- ear Ridge Regression; KRR, Kernel Ridge Regression; DNN, Deep Neural Network; PBCs, Periodic Boundary Conditions.