Inverse design of stochastic, voxelated thermo-viscoelastic digital materials

Domain: physics
Published: 2025-11-28
Editorial outcome: minor_revision

Abstract
Polymermaterialjettingenablesthefabricationofvoxelated,multi-materialstructureswithmaterialcontrolatthemicroscale.However,currentworkoftenneglectsviscoelasticeffectsanddesigningvoxelateddigitalmaterialsremainschallengingduetothecomplexityofthevastdesignfreedomandintractabilityofefficientlymodelingmacroscalevoxelstructures.Wepresentanefficientrepresentationofstochasticallymixed,voxelateddigitalmaterialsanddevelopageneralizedviscoelastictemperature-dependentmaterialmodeltodesignandsimulatedigitalmaterialsmixedfromtwoconstituentpolymers.Thematerialmodelisbasedonanextendedpercolationtheoryconsideringfrequencyandtemperature.Anartificialneuralnetworkistrainedonthematerialmodeltodirectlyestimatetargetmaterialbehaviorgivenarbitrarynon-linear,userrequirements.Theapproachisvalidatedusingtwocasestudiesrequiringtailored,non-linearmaterialbehavior:apersonalizedwristorthosisandamachinedamper.Theseshowthenewlyunlockedpossibilitiesforthedesignandfabricationoftuned,stochasticdigitalmaterials.

Introduction
Architected materials have transformed material design due to the unique interplay between the architected mesoscale geometry and the resulting macroscopic mechanical properties1,2. The resulting design flexibility enables advances in different applications such as high stiffness lattices3, optical waveguides4, origami5, and multi-state structures6,7. The realization of these designs is often supported by Additive Manufacturing (AM) due to the process escapability of producing free-form structures8–11.

Recently,polymermaterialjettinghaspushedarchitectedmaterialsfurtherbyenablingmulti-materialAMbymaterialassignmentpervoxelonamicroscale10,12,13.Priorworkhasextensivelyexploredthisprocesstodesignmetamaterialswithfeaturessignificantlylargerthanthevoxelscale. The resulting meso-scale features are designed as deterministic elements, often optimized for shape morphing6,14, stiffness tuning3, and topology control15. These design approaches often assume homogenous or deterministic voxel layouts16, which enable accurate property estimation by simulation at the meso-scale. In contrast, material design on the voxel-scale enables smaller features but is governed by voxel-scale interaction effects of the base materials, causing the meso scale material behavior to be dominated by micro-scale material interfaces rather than only by the base material properties17. This effect is often leveraged to create digital materials exhibiting emergent behavior through the stochastic mixing of base materials18. These stochastic digital materials enable, for example, the direct tuning of material stiffness1,8,14, the polymer shape memory effect19,20, and the design of continuously graded material interfaces15,21. While it has been shown that modeling and simulation of stochastic digital materials can be tackled by coarse-graining and homogenization1,8, most work relies on an assumption of only elastic behavior8,9,15,19,22–26 or neglects the dependence on a wide temperature range18,21 and loading rate17–20,27. This impedes the use of voxelated polymer material jetting for use cases with more complex mechanical behavior, such as damping structures and applications requiring temperature-dependent stiffness modulation.

Moreover, the combinatorial scaling behavior of voxel-based material allocation makes inverse material design, in this context, intractable when using conventional computational methods. While various data-driven and optimization-based inverse design frameworks have been proposed for mechanical metamaterials, they often target unit cell design or fixed lattices, are often constrained to linear behavior and are limited to shape or stiffness matching of a single state of a single structure23,28–31. Additionally, the currently proposed data-driven approaches rely heavily on understanding the microscale structure for the automated generation of large datasets29,30,32–36. For polymer voxel AM, the lack of voxel-level microscale simulation prevents the systematic generation of artificial datasets and calls for empirical-based modeling3,14,17. To date, there is no generalizable and scalable inverse material mapping between desired viscoelastic, time-dependent, and temperature-dependent material response and voxel-scale material composition in stochastically mixed, multi-material systems.

In this work, we present a unified modeling and inverse material design framework for stochastically mixed, voxelated digital materials, combining a generalized viscoelastic material model with a data-driven inverse material design framework (see Fig. 1). The material model homogenizes stochastic digital materials of two base materials using an extension of the physics-based de Gennes percolation theory37 to account for the time-dependent and temperature-dependent stiffness as well as voxel-scale interaction effects. The model is empirically tuned and validated on experimental data. Based on this model, an Artificial Neural Network (ANN) is trained to directly infer voxelated digital materials from user-defined mechanical design requirements. The proposed inverse material design framework is geometry-independent, does not require problem-specific re-training, and enables the direct design of viscoelastic, voxelated digital materials across a range of time scales and temperatures. The approach is validated through two representative case studies. The first case study considers the modeling of a monolithic, voxelated machine damper under cyclic loading, thus showing the simulation of hysteresis and smoothly graded material interfaces that are used to prevent stress concentrations. The second case study investigates the inverse material design of a personalized, multi-material wrist orthosis that can limit wrist deflection and compensate for swelling and unexpected impacts, thus demonstrating consideration of material requirements at multiple timescales.

Results
Efficient representation of voxelated digital materials
The representation of voxelated digital materials quickly exceeds computational limits. While a conventional dogbone type specimen can contain distinct regions of hard and soft materials38, a voxelated specimen can include graded boundaries to prevent stress concentrations and yields a complex material distribution of 25 million individual voxels (Fig. 2a). Further, a close-up inspection of a voxelated, AM material mixed from two polymers reveals a high amount of material smearing on a scale significantly larger than the theoretical AM resolution of 84µm (Fig. 2a). Since computational limits and complex intertwining of material interfaces both prevent the simulation of idealized, voxelated digital materials on the microscale, a representation on an intermediate mesoscale using a coarse grid is proposed, effectively separating the actual voxel layout from the mesoscopic behavior of the digital material. More specifically, macroscopic geometry is approximated by grid points on the mesoscale, where each grid point is assigned a designated material (Fig. 2b). For AM, voxels are assigned a base material stochastically based on the defined digital material mixture. The presented representation approach separates the scales of material assignment from the voxel scale by tuning the resolution of the coarse grid to the finest feature size and enables direct material assignment to Finite Element (FE) meshes for simulation when using a homogenized material model on the mesoscale. The coarse grid resolution is thus adapted to the size of the geometric features in the structure, the desired material assignment resolution, and the FE mesh resolution required for simulation (see Section A.2.4 in the supplementary information).

Material model: stochastic digital materials and de gennes percolation theory
Previous work shows that the homogenized properties of stochastic mixtures of two base materials can be modeled using de Gennes percolation theory37. This theory describes how the homogenized digital material behaves based on the volume fractions of the individual base materials. As the volume fraction of material one ρ increases, the mixture morphology transitions from dispersed droplets to an interconnected network18, thus changing the ratio of the effective properties non-linearly according to Eq. (1).

r = (ρ1/ρc)μ (1)

The percolation threshold, ρc, is set to 0.163 based on previous literature19, validation by visual inspection of different stochastic digital materials and an optimization procedure based on sample II from Fig. 3a (see supplementary information). The percolation coefficient, µ, for two materials with different stiffnesses is evaluated by testing the two base materials, namely a soft material and a rigid material, and five stochastic digital materials (Fig. 3a) using a dynamic mechanical analysis (DMA). The Time-Temperature Superposition Principle39,40 (TTSP) is assumed and the materials are subsequently modeled using a generalized Maxwell model represented using the Prony series41,42 in combination with a Williams-Landel-Ferry40 (WLF) and Arrhenius function43. The resulting master curves show a non-linear dependency between the volumetric mixture ratios and resulting stiffness in the frequency range of 10−10 Hz and 1010 Hz. The resulting stiffnesses E’ at different frequencies are biased toward the stiffer material at higher frequencies and toward the softer material at low frequencies. This indicates that the percolation coefficient µ is dependent on frequency and due to the TTSP also on temperature. The percolation coefficient µ is thus fit per frequency and temperature using four of the seven experimentally determined Prony series curves (Fig. 3b, sample I, III, V, and VII). The percolation threshold, ρc, is fit using sample II and the final percolation model is validated using samples IV and VI. Additionally, all other mixtures are also compared to the percolation model to assess the quality of the fitting procedure (Fig. 3d).

Material behavior estimation for inverse design
The percolation-based interpolation enables the stiffness estimation of voxelated digital materials at each timescale and temperature. This is used to create a mapping M between mixture ratio, timescale, temperature and stiffness (Fig. 4a) as well as to the material viscoelasticity, expressed as a stiffness decrease over a characteristic time period (Fig. 4b). While the stiffness is maximal for a higher percentage of the rigid material, the viscoelasticity ΔEτ shows an interesting switching behavior across temperature.

When designing, it is desirable to directly determine the material mixture ratio based on given requirements. In this case, each requirement can be defined as a tuple of target stiffness at a given timescale and temperature. For multiple requirement tuples, the determination of the target mixture becomes ambiguous. Here, the ambiguity is bypassed by training an ANN with a multi-layer perceptron architecture on the mapping M for both estimating mixture ratios (ANNmix) and stiffness values (ANNE) from the remaining three dimensions, respectively and chaining the ANNs into an inverse design framework. Detailed information on the ANN architecture and training can be found in the subsequent methods section and in the supplementary information. In this framework, the user-provided requirement tuples are extended by linearly spaced requirement tuples and the network ANNmix estimates target material mixtures from the extended set of requirement tuples. These estimations are averaged, providing the target material mixture ρmix (Fig. 4c). The target material mixture determines the glass transition temperature Tg,mix by percolation from the experimentally determined Tg,test of the tested samples (Fig. 4d). The target mixture ratio is provided as input to the network ANNE, which estimates the material stiffness over the full time range from 10−15 s to 1015 s for the determined mixture at room temperature. The resulting stiffness over time is the master curve for the target mixture. A generalized Maxwell model with hyperelastic stiffness elements is fitted to the master curve using the previously determined Tg,mix and an interpolated hyperelastic curve based on experimental data (Fig. 4e). Depending on the diversity of requirement tuples, the degree to which the estimated stiffness of both ANN and the fit model fulfill the required stiffness can vary. The agreement between the fit model and the network estimation, however, is a matter of the fitting performance and varies depending on the temperature and timescale (Fig. 4f). The observed differences can be attributed to the Arrhenius shift factor of the experimental curves not being equal, which cannot be depicted using a single Prony series with a single shift factor in the fit model.

Since the inverse material design framework is trained on the theoretical material model, the influence of real-world variability on the estimated properties is evaluated experimentally and analytically (detailed description in supplementary information). Experimental mechanical characterization of identical samples revealed a manufacturing variability of 3% in force response. The analytical analysis of the impact of stochastic base material assignment is conducted using a cube of a given feature size, the number of voxels within the cube determines the effective sample size of the binomial mixing process. Assuming voxel-wise material assignment with a certain probability, the standard deviation of the number of voxels allocated can be quantified. This analysis reveals that the stochastic variability is larger than the manufacturing variability of 3% for features below 1–2 mm, depending on the target material mixture. Similar values of around 5 mm can be found in literature18. This suggests that the deterministic ANN model behaves accurately for homogenized mesoscale structures.

Case study: machine damper
A promising example to demonstrate the validity of the proposed material model is the inverse design of a standard machine damper, which is used in a range of applications from soft robotics to aerospace44,45. There is a need to tune machine dampers for a particular application to improve performance46,47, but current simulation approaches fail to depict the typical hysteresis. Further, multi-material structures are plagued by stress concentrations during loading due to large stiffness gradients at the material boundaries38. Voxelated digital materials offer a potential solution since the damping material can be tuned for a target application and graded boundaries can be introduced between single, solid materials in the structure. The presented damper design is based on a Vibrostop AAT 20/N damper (Vibrostop SRL, Milan, Italy) with the interfaces between the compliant dampening material and the rigid housings showing a linear gradient to reduce stress concentrations. The mesoscopic material formulation enables the representation of the full structure on a 405,769-point grid instead of using 273.3 × 106 individual voxels. The discussed mesoscopic representation, the tetrahedral finite element mesh and the AM damper are shown in Fig. 5a. The viscoelastic material formulation in this work enables the prediction of the voxelated damper’s hysteresis, which was not feasible in previous work on voxelated digital materials due to the typical adherence to linear elasticity and a quasi-static time scale19. The hysteresis is simulated and tested at three different frequencies (0.1, 1, and 10 Hz), each at two temperatures (15 °C and 25 °C). The simulation is performed in Abaqus CAE, with runtimes ranging from 100 to 145 min for a 253,777-element mesh assigned to 28 individual percolated digital material sections. The experimental and simulated force-displacement hysteresis curves are shown in Fig. 5b. The best agreement is observed at 0.1 Hz, with a mean reaction force (RF) error of just 5.6% across both temperatures. At 1.0 Hz, the mean error increases to 20.0% and at 10 Hz, it reaches 26.6%. The overall accuracy of the predicted hysteresis is significant, given the challenges of simulating a complex 3D structure with graded material transitions subject to cyclic loading across multiple frequencies and temperatures.

Case study: personalized wrist orthosis
The direct inverse design of the material distribution of a personalized, multi-material wrist orthosis is now shown. Personalization in healthcare often requires versatile design methods that can provide a design estimate for multiple local requirements, such as long-term compliance and short-term stability at the same location. While patient-specific design requirements vary significantly across patients and pathologies, and are not available in this study, the orthosis design is based on literature-derived design requirements to demonstrate the method’s potential. The design requirements are formulated for three design regions all having a material thickness of 3mm: the thumb region, the stabilization splint, and the bulk shell (Fig. 6a). The stabilization splint provides resistance against wrist flexion and extension allowing a maximum displacement of 2mm under a force of 15N48. The thumb region must permit free thumb motion, and the bulk shell should provide impact resistance while still being compliant to long-term changes of the underlying anthropometry, such as swelling (Fig. 6a). The requirements for the stabilization splint and the thumb region are thus formulated as the maximum and minimal material stiffness available. The bulk shell is subject to more than one requirement, impacts are approximated as a concentrated force of 120N over a timescale of 0.01s for impact resistance and a stiffness decrease from 900MPa to 150MPa within 86400s (see supplementary information for further details).

The design requirements for the three orthosis regions are provided to the ANN and the resulting digital materials are used to generate an orthosis with distinct material regions and linear material gradients in-between. The coarse grid representation reduces a total of 15.25 billion voxels to 28,175 grid points. The proposed material model enables the simulation of orthosis flexion and extension, long-term compliance, and impact resistance in Abaqus CAE using a 5587-element mesh assigned to 19 individual percolated digital material sections. The wrist flexion and extension are simulated as a radially exerted force of 15N and tested experimentally on a custom rig using an additively manufactured wrist and 15N of weights attached to the palm. The more critical case of flexion shows an experimental and simulated displacement of 1mm and 0.923mm, respectively (Fig. 6c), indicating high accuracy considering the complex load case. The long-term compliance is experimentally evaluated by measuring the normalized RF of a half pipe of the shell material on an indenter initially loaded with 6N over a period of 24h. The simulation is performed by applying a radial displacement on the shell for 24h. The normalized RF is compared to the normalized simulated reaction pressure and shows good correspondence (Fig. 6d). To further validate long-term stability, the shell material is tested under cyclic loading for 10,000 cycles, and no significant change in stiffness is observed (see supplementary information). The impact resistance is evaluated in simulation only with a ramped concentrated force to the shell (Fig. 6e) and shows very low displacement results, indicating that the orthosis can absorb impacts while still exhibiting long-term compliance to account for swelling. Additional analysis of these four load cases at 0°C (cold weather) and 37°C (body temperature) is included in the supplementary information. The aggregated results across all three temperatures indicate that the multifaceted requirements placed on the orthosis have been fulfilled, both in the simulated and experimental behavior. The ability to predict complex time-dependent and temperature-dependent behavior in voxelated digital materials enabled by the presented modeling approaches is a crucial step towards highly personalized wearables.

Methods
Materials
All samples and designs are printed on a Stratasys J850 Prime (Stratasys ltd., Eden Prairie, MN, USA) using the two build materials Agilus 30 Black and VeroUltraWhite, SUP706 is used for support structures. The prints are printed using the Voxel Print™ utility included in the GrabCAD Print Pro™ Package. The designs are generated in MATLAB R2024a51 and imported in GrabCAD as a set of PNG images, where each print material is encoded by a separate color. The J850 machine has a layer height of 27µm for voxelated manufacturing and a resolution of 600 dpi and 300 dpi in x- and y-direction, respectively. To align the planar resolutions, the pixels in the imported images are doubled in x-direction to form quasi-voxels with a size of 84.66×84.66×27µm. The prints are post processed by removing the support material mechanically by hand, followed by cleaning the part at a water jet cleaning station.

Mechanical testing
The two base materials and five stochastic digital materials are tested with a DMA and a tensile test. The sample geometry is a type IV dogbone geometry according to ASTM D638 with the dogbone ends printed in VeroUltraWhite, the core printed in the digital material under test and the transition between the two being a 2mm long linear material gradient. The seven material mixtures are AB/VW mixtures with an AB content of 0, 16.5, 33.3, 50, 66.6, 83.5, and 100%. The DMA testing is performed on an Instron ElectroPuls E3000 (Instron, Norwood, MA, USA) using an Instron Dynacell 2527 load cell with a range of ±5kN, Instron fatigue-rated mechanical wedge action tensile grips and an Instron 3119-600 series temperature-controlled chamber. Instron WaveMatrix is used to record all relevant data and to calculate the storage modulus E’, loss modulus E’’ and dissipation factor tan(δ) for a sinusoidal displacement input for frequencies between 1Hz and 18Hz with 1Hz increments. The measurement temperatures are set for each sample to force a change from the glassy to the rubbery state. The sample is pre-strained to always be under tension and the pre-strain is increased for an increasing temperature to compensate for thermal expansion. An additional tensile test is performed at temperatures where the samples are in the rubbery state with a speed of 3mm/min to ensure quasi-static conditions. The resulting force is measured, yielding the hyperelastic stress-strain data for all samples.

Material model
The material is modeled using a generalized Maxwell model initially consisting of a single, linear Hookean element k0 and n Maxwell elements consisting of a linear Hookean element of stiffness kn and a Newtonian dashpot, or viscous element, with a relaxation time of τn. As shown in Eq. (2), the relaxation time is modulated using a scaling factor α based on the WLF and Arrhenius equation.

log10 αT = { C3 / (T - Tg) if T < Tg; -C1(T - Tg) / (C2 + T - Tg) if T ≥ Tg (2)

The glass transition temperature Tg is derived from the observed peak in tan(δ), the WLF parameters C1 and C2 are fixed to 17.44 and 51.6, respectively according to Ferry40. The Arrhenius parameter C3 is determined based on the TTSP. The shifting procedure is conducted using a numerical optimization scheme based on a Nonuniform Pattern Search algorithm. The objective is to shift the individual frequency sweeps to create a master curve that is as smooth as possible by minimizing the second derivative along its entire length. This is achieved by fitting a spline through the individual data points gathered during the mechanical testing procedure and numerically determining the second derivative. Thermal expansion is also considered using a thermal element with an expansion coefficient of α. The experimental data is fit with a model consisting of n=31 Maxwell elements, the estimated material in the design framework is also fit with a model with n=31 Maxwell elements. The fitting procedure is conducted using a numerical optimization scheme based on a nonlinear least-squares algorithm. The fitted models are expressed using a normalized Prony series, the stiffness of the Maxwell elements is thus defined as a fraction of the long-term stiffness k0. In the final simulation model, the Hookean element k0 that models the long-term stiffness is replaced with a Mooney–Rivlin hyperelastic element. In combination with the normalized Prony series, this enables the simulation of larger strain ranges. The fitting of the Mooney–Rivlin hyperelastic model is conducted at runtime in the FE simulation software.

Neural network setup
The neural networks ANNmix and ANNE are both trained on a dataset created using the temperature and frequency-dependent percolation. The stiffness E is calculated for each temperature between −30°C and 130°C with 1°C increments, each AB/VW mixture with 1% increments and each time between 10−15 s and 1015 s with 300 regular increments in the logarithmic space, yielding a dataset of 4.878 million data points. The training of ANNE is performed on the full dataset and the training of ANNmix is performed on a reduced dataset, where all data for mixture ratios above the percolation threshold is culled to prevent one-to-many mapping. The full and reduced dataset are individually normalized to values between 0 and 1. The networks are trained on a training set with 80% of the data points and tested on a test set with the remaining 20%. Both networks are multi-layer perceptrons with a constant width and depth, an activation function globally defined for all neurons, an input size of three, and an output size of one. The network ANNmix has a width of 180, depth of seven (five hidden layers), GELU activation functions, and is trained using the ADAM optimization algorithm52 with a learning rate of 0.001 and a batch size of 256 over 150 epochs. The loss function is the mean squared error, and the test loss of the trained network is around 10−5. The network ANNE has a width of 150, depth of seven (five hidden layers), GELU activation functions, and is trained using the ADAM optimization algorithm52 with a learning rate of 0.0001 and a batch size of 256 over 150 epochs. The loss function is a weighted mean squared error with increased weights for lower stiffness labels to counteract the decreasing relative accuracy for labels around zero, a penalty term for stiffness estimations below zero is also added. The test loss after training is around 10−4. All hyperparameters are evaluated using a full factorial hyperparameter evaluation, including the ANNs width, depth, the activation function, learning rate, batch size, and number of training epochs (see supplementary information).

Finite element simulation
All finite element simulations are performed in Abaqus CAE 2022 using a viscoelastic model with hyperelastic elements. The temperature shifting behavior is modeled using a custom UTRS user subroutine. The simulation of the damper and the orthosis is performed using a standard implicit solver. The details on the mesh generation and the boundary conditions for both case studies can be found in the supplementary information.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents a computational framework for the inverse design of architected polymer structures with voxel-based microstructures, focusing on thermo-viscoelastic behavior. It combines percolation-inspired mixture modeling, Prony-series-based viscoelasticity, and a ANN model for inverse design. Applications to damping components and orthotic devices are included. While the manuscript is well written and the workflow is well implemented, the core methodology lacks novelty, and the proposed contributions seem incremental. Many similar frameworks already exist in the literature, and this work does not provide new theoretical insights, modeling techniques, or significantly more powerful design capabilities. Thus, the proposed work is not appropriate for publication in the journal Nature Communications, and a computation-focused journal would be more suitable.\n1. The paper reuses well-established techniques—homogenized mesostructural modeling, time–temperature superposition, and ANN-based inverse design, without introducing significant innovation or theoretical advancement. In addition, the authors appear to miss the majority of the vast literature on the inverse design of architected structures and soft materials. There are a large number of studies that this manuscript ignores.\n2. There is little to no quantitative comparison with existing frameworks for inverse design of soft architected materials, especially those using voxel representations and ML-based and gradient-based design tools. The authors should first conduct a comprehensive literature review and potentially revise and improve their approach.\n3. There are two example applications (a damper and a wrist orthosis) that are rather simple and don’t demonstrate the need for computational inverse design, and it is unclear how the method works for more technically challenging case (such as involving geometric constraints, multiple conflicting objectives, or large deformations).\n4. The neural network is used as a black-box inverse predictor, but its architecture, training dynamics, generalization behavior, and failure cases are not well analyzed.\n5. The ANN and mixture models assume deterministic relationships, but the actual AM processes and viscoelastic properties exhibit high variability.\n6. The percolation-based mixture model is used, but not deeply analyzed or validated. It is unclear whether this offers any benefit over more traditional homogenization models in terms of prediction accuracy or computational cost.", "reviewer_2": "It is a great honor to receive an invitation from [Redacted] to review this manuscript from Nature Communications.\nThe overall recommendation is to accept it with minor revisions.\nThe detailed recommendation is as follows:\nIn this paper, an innovative reverse design framework, combining osmosis theory and neural networks, is proposed to provide a valuable tool for the design and simulation of multi-material voxel structures. It is very interesting. However, the academic rigor and engineering utility need to be further enhanced by additional experimental validation, quantifying the effect of fabrication errors, and refining the biomechanical assumptions. I recommend this manuscript to be accepted after a slight revisions.\na. The current model is based on a blend of two polymers (Agilus30Black and VeroUltraWhite) and the permeation threshold (ρc=0.15 ρc=0.15), but independent experimental validation against the material combinations of the present study is not provided. Furthermore, whether the frequency dependence of the permeation index μμ is applicable to other material systems (e.g., thermoplastic elastomers versus rigid plastics) has not been discussed.\nRecommendation:\n-Add experimental validation of permeation thresholds (e.g., microstructural characterization) for the selected material combinations.\n-Discuss the applicability of the model to other material systems, or clarify the scope of application in the “Limitations” section.\nb. The training data is generated by the theoretical model (1.6M data points) and not directly based on experimental data. Although the theoretical model has been partially validated, it may result in a network that is not robust to actual manufacturing errors (e.g., material “fuzziness”).\nRecommendation:\n- Introduce test data from actual print samples as part of the training set to enhance the generalization of the model.\n- Provide a sensitivity analysis of the network under noisy data (e.g., prediction error after randomly perturbing the mixing ratio).\nc. The design requirements for the wrist orthosis (e.g., maximum displacement of 2mm at 15N force) are based on literature 2929 but do not incorporate actual biomechanical data from the patient (e.g., joint stiffness for different pathologies). Fatigue properties for long-term use (e.g., material relaxation after cyclic loading) were also not validated.\nRecommendation:\n- Supplement patient-specific biomechanical data or cite broader clinical studies to support design criteria.\n- Add long-term fatigue testing (e.g., change in stiffness after 10,000 cycles of loading).\nd. The paper points out the “blurring” phenomenon during material injection (Figure 2a), but does not analyze its impact on simulation accuracy. For example, the actual printed gradient transition region may deviate from the linear mixing assumption of the design.\nRecommendation:\n- Quantify manufacturing errors (e.g., analyze actual material distribution by micro-CT) and introduce tolerance parameters in the model.\n- Discuss the extent to which manufacturing errors affect case study performance (e.g. damper hysteresis curves).\ne. Simulation parameters (e.g., mesh sensitivity analysis, time step selection) are not specified, and no mention is made of the convergence of multi-threaded calculations (the supplementary material mentions a single-threaded Abaqus run). In addition, the parameter fitting procedure for the hyperelastic model (Mooney-Rivlin) is not disclosed.\nRecommendation:\n- Add grid-independent validation (e.g., comparison of displacement errors at different grid densities).\n- Disclosure of the fitting method and experimental data support for the parameters of the hyperelastic model.\nf. The paper emphasizes the “high efficiency” of the method (7 orders of magnitude reduction in data points) but does not discuss the practical limitations in its engineering applications, e.g:\n- The demand on computational resources for high-resolution voxel design (e.g., preprocessing time for 280,000 grid points).\n- Cost and scalability of multi-material printing (e.g., difficulty in removing the support material SUP706).\nRecommendation:\n- Add an analysis of the trade-off between computational efficiency and accuracy (e.g., mesh resolution vs. simulation error).\n- Discuss economic challenges in industrial-scale applications (e.g., material waste rate, print time).\nAll the best to [Redacted] and everyone! Have a good day! Yuze Nian May 1, 2025", "reviewer_3": "This manuscript presents a novel inverse design framework for voxelated thermo-viscoelastic materials. This approach effectively reduces computational complexity via mesoscale homogenization and demonstrates promising predictive accuracy in two functional case studies. While this work offers a valuable contribution to data-driven design of architected soft materials, its limited validation is a significant weakness. Therefore, I would like to recommend it for publication after addressing some concerns. Detailed comments are provided below.\n• The manuscript relies heavily on simulation results without providing adequate experimental validation, which may raise questions about the practical applicability of the findings. The Prony series fitting is based on only four experimental datasets for training and a single dataset for validation, which may be insufficient to robustly capture the full range of material behavior, as acknowledged by the authors. Moreover, the inverse design framework relies solely on simulations. If possible, incorporating such experimental validation would significantly enhance the credibility of the approach.\n• Although the authors emphasize the importance of temperature-sensitive material responses, most of the work are conducted at room temperature. Extending the analyses to include elevated or variable temperatures, especially those relevant to the case studies, would better demonstrate the robustness and applicability of the model.\n• I would recommend replacing the term “graduated” with “graded” to reduce potential confusion among interdisciplinary readers, as it is the more commonly accepted term in the related field.\n• Page 7, line 186: please check the error.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to the Reviewers’ Comments for Manuscript:\n“Inverse Design of Stochastic, Voxelated Thermo-Viscoelastic Digital\nMaterials”\nDear Reviewers,\nThank you very much for the helpful comments to further improve the quality of our\nmanuscript. We have carefully considered every comment and address them here, in\naddition to the changes made in the revised manuscript that are marked using a yellow\nhighlight and the corresponding identifier [R*C*]. Following your comments, we conducted\nextensive additional experimental validation and literature review. This additional validation\nled to refitting the material model and complete retraining of the ANN to include the newly\ngenerated experimental data. The main paper has been updated with the new results, and\nthe supplementary information has been significantly expanded to include additional analysis\nand results. To better emphasize the scope of the work, the title has been changed to\n“Inverse Design of Stochastic, Voxelated Thermo-Viscoelastic Digital Materials”.\nWe further made additional minor changes to the manuscript to improve the readability and\nclarity.\nThe Authors\nPage 1 of 12\nReviewer #1: This manuscript presents a computational framework for the inverse design of\narchitected polymer structures with voxel-based microstructures, focusing on thermo-\nviscoelastic behavior. It combines percolation-inspired mixture modeling, Prony-series-based\nviscoelasticity, and a ANN model for inverse design. Applications to damping components\nand orthotic devices are included. While the manuscript is well written and the workflow is\nwell implemented, the core methodology lacks novelty, and the proposed contributions seem\nincremental. Many similar frameworks already exist in the literature, and this work does not\nprovide new theoretical insights, modeling techniques, or significantly more powerful design\ncapabilities. Thus, the proposed work is not appropriate for publication in the journal Nature\nCommunications, and a computation-focused journal would be more suitable.\nR1C1: The paper reuses well-established techniques—homogenized mesostructural\nmodeling, time–temperature superposition, and ANN-based inverse design, without\nintroducing significant innovation or theoretical advancement. In addition, the authors appear\nto miss the majority of the vast literature on the inverse design of architected structures and\nsoft materials. There are a large number of studies that this manuscript ignores.\nAn additional literature review has been conducted, and the introduction has been\nsignificantly rewritten. The key novelty of the presented work is the material modeling of\nhigh-resolution voxelated digital materials that are geometry-independent and use physics-\nbased percolation material modeling. Specifically, the novelty lies in the definition of a time-\ndependant and temperature-dependent adaption of the de Gennes percolation theory. This\nmodeling approach enables the simulation of complex behavior such as hysteresis, shape\nmemory effect, and stiffness modulation across a wide range of time scales and\ntemperatures. Additionally, the combination of the efficient mesoscale representation and the\nability to predict the properties of any digital material composed of the two base materials\nenables the simulation of complex graded transitions in a computationally efficient manner\nwhile undergoing complex loading. The experimentally validated ANN-assisted inverse\ndigital material design method enables designers to determine a digital material that satisfies\ngiven requirements at various temperatures and time scales. In summary, these findings\nrepresent a significant step towards a generalized modeling approach to stochastic\nvoxelated digital materials and enable the computational design of high resolution voxelated\nstructures with tuned behavior across multiple time and temperature scales for a variety of\napplications in engineering and science.\nR1C2: There is little to no quantitative comparison with existing frameworks for inverse\ndesign of soft architected materials, especially those using voxel representations and ML-\nbased and gradient-based design tools. The authors should first conduct a comprehensive\nliterature review and potentially revise and improve their approach.\nMany existing inverse design frameworks in the field of metamaterials and soft materials\nfocus on inverse geometric design, such as optimizing unit cell topologies and lattice layouts.\nIn contrast, our framework determines digital material allocation for a given geometry. Thus,\nwe changed the terminology in the paper to “Inverse Design of Stochastic, Voxelated\nThermo-Viscoelastic Digital Materials”. The challenge addressed in this paper lies in creating\nstochastic voxelated digital materials to meet viscoelastic, time-dependent and temperature-\ndependent requirements. The goal of our work is not to generate local shapes or structures,\nbut rather to determine a fitting digital material that satisfies the user requirements. This\nPage 2 of 12\nmakes existing geometric optimization frameworks not directly comparable to the problem\nwe address, but rather complementary. To emphasize this distinction and to highlight the\ncontributions of this work, we extended the introduction and related this work more tightly to\nthe existing body of research. We hope this clarification resolves the concern and helps to\nposition our unique contribution.\nR1C3: There are two example applications (a damper and a wrist orthosis) that are rather\nsimple and don’t demonstrate the need for computational inverse design, and it is unclear\nhow the method works for more technically challenging case (such as involving geometric\nconstraints, multiple conflicting objectives, or large deformations).\nThe damper and orthosis case study have been adapted to better show the unique\ncapabilities of the model. The damper is now validated at two different temperatures with\nloading at three different frequencies for each temperature. The orthosis is validated at two\nadditional temperatures. These additional results have been added to Fig. 5 and Fig. 6 in the\nmain text and are further discussed in the supplementary information.\nThe context of the inverse digital material design approach has been further clarified as\ndescribed in [R1C2]. Since the goal of our work is not to generate local shapes or structures,\nbut rather to determine the digital material composition that satisfies the user requirements,\ngeometric constraints cannot be addressed directly. In the case of conflicting material\nrequirements, the current method determines a digital material that represents a tradeoff\nbetween the different objectives. The inverse digital material design method uses the dataset\nof all percolation-based material master curves to determine the mechanical response of the\nstochastic digital material. When there is large deformation or complex, nonlinear load\ncases, the global mechanical response of the macroscale structure is defined by a strong\ninterplay between the digital material composition, the geometry of the structure and the load\nhistory. To address this interplay, the proposed modeling approach could be integrated into\na numerical optimization scheme that controls coarse grid digital material assignment to find\nan optimal digital material distribution for the given global load case. This is outside the\nscope of the presented work.\nR1C4: The neural network is used as a black-box inverse predictor, but its architecture,\ntraining dynamics, generalization behavior, and failure cases are not well analyzed.\nWe have added a detailed analysis of the ANNs hyperparameter selection, training behavior,\nand ANN accuracy to the supplementary information (section C, Figure S14-S17). The\nhyperparameter selection is performed by a full factorial analysis and results in ANNs with\nlower training loss and test loss. For the ANNs with the lowest loss, the convergence of the\ntests show that the data is not overfit and the ANN acts as a good surrogate for the\nunderlying dataset. Further, we provide accuracy maps to visualize the spatial accuracy of\nthe ANN to the test data. The accuracy is interpreted and set into the context of this work.\nThe ANN estimations show a high accuracy over the full range of all parameters, this\ndemonstrates that the ANNs operate as transparent, well-posed interpolators, rather than a\nblack-box predictor. We further show that the failure cases, i.e. where the ANN matches the\ndataset the least, occur in regions with limited mechanical sensitivity and have negligible\nimpact on the mechanical properties of the resulting digital material.\nPage 3 of 12\nR1C5: The ANN and mixture models assume deterministic relationships, but the actual AM\nprocesses and viscoelastic properties exhibit high variability.\nWe fully agree that the experimental viscoelastic behavior of AM materials can be affected\nby both the variability of the AM process and variability of the stochastic material\nassignment. We added a quantification of both variabilities by a combined experimental and\nanalytical analysis, where details are given in the supplementary information (Section A.4)\nand summarized in the main text. Manufacturing variability is characterized by uniaxially\ntesting 20 gradated dogbone samples with the same stochastic seed. We show that the\nmanufacturing variability is 3% at peak reaction force. Variability due to stochastic\ndistribution is experimentally characterized by uniaxially testing 20 gradated dogbone\nsamples with unique stochastic seeds but having the same mesoscale material allocation.\nThese show a variability of 2.5%, which is below the manufacturing variability. We further\nextend the analysis of the design variability by an analytical approach, where we show that\nstochastic material allocation variability only becomes significant for features sizes below 1-\n2mm. This supports the validity of the deterministic modeling approach for the mesoscale\nstructures targeted in this work.\nR1C6: The percolation-based mixture model is used, but not deeply analyzed or validated. It\nis unclear whether this offers any benefit over more traditional homogenization models in\nterms of prediction accuracy or computational cost.\nThe analysis of the percolation-based mixture model has been significantly expanded. Two\nadditional digital materials have been tested experimentally. The model is now fit in a two-\nstage approach based on five digital materials and is subsequently validated with two digital\nmaterials instead of one (Fig. 3d). The percolation exponent that results from the fitting\nprocedure is also analyzed in more detail in Fig. 3c and the corresponding text.\nThe key benefit of the presented approach over more traditional homogenization models is\nthe relatively low number of samples needed to conduct the material fitting procedure due to\nthe use of a physics-based percolation model. While previous methods use seven1 to\neleven2,3 different digital materials to fit a model valid for a specific temperature and\nfrequency, we use only five to fit a model over the full temperature and frequency range. The\nlow sample number is especially significant since the data gathering techniques that allow\nhigh-value data extraction, such as DMA, are associated with high experimental effort.\nPage 4 of 12\nReviewer #2: In this paper, an innovative reverse design framework, combining osmosis\ntheory and neural networks, is proposed to provide a valuable tool for the design and\nsimulation of multi-material voxel structures. It is very interesting. However, the academic\nrigor and engineering utility need to be further enhanced by additional experimental\nvalidation, quantifying the effect of fabrication errors, and refining the biomechanical\nassumptions. I recommend this manuscript to be accepted after a slight revisions.\nR2C1: The current model is based on a blend of two polymers (Agilus30Black and\nVeroUltraWhite) and the permeation threshold (ρc=0.15 ρc=0.15), but independent\nexperimental validation against the material combinations of the present study is not\nprovided. Furthermore, whether the frequency dependence of the permeation index μμ is\napplicable to other material systems (e.g., thermoplastic elastomers versus rigid plastics)\nhas not been discussed.\nRecommendation:\n-Add experimental validation of permeation thresholds (e.g., microstructural characterization)\nfor the selected material combinations.\n-Discuss the applicability of the model to other material systems, or clarify the scope of\napplication in the “Limitations” section.\nThis comment addresses an important point regarding both experimental validation and\nmodel generalizability. To address the experimental validation, we have added an additional\nvalidation of the percolation threshold for the Agilus30Black (AB) and VeroUltraWhite (VW)\nsystem used in this study (Supplementary Information A.3). In this qualitative validation, a\nset of digital material samples with varying base material ratio is fabricated and visually\ninspected, confirming that the transition from connected to dispersed material regions occurs\nbetween 15 and 20% AB content. Further, we have characterized two further digital material\nsamples using DMA testing to validate the percolation-based material model. The two new\ndigital material samples have base material mixture ratios of 16.5% AB and 83.5% VW and\nvice versa. These two digital material samples and the previous validation sample with a\n50/50 AB/VW base material mixture are placed between the four initial samples to set up the\nmodel. The percolation threshold 𝜌 is now determined based on the 16.5% VW and 83.5%\n𝑐\nAB sample using an additional fitting step on top of the fitting of the percolation coefficient 𝜇 .\nThe accuracy of the model prediction to the validation sample data of two samples, as well\nas the accuracy of the model prediction to the samples used for fitting, shows that the\nproposed model is valid across the whole design domain.\nTo clarify the generalizability of the frequency-dependent percolation coefficient across\ndifferent material systems, we have expanded the discussion to explicitly state that the\ntemperature and frequency dependence may vary with different material systems. This\nvariation is likely influenced by interfacial interactions, polymer chain entanglement, and\nmaterial smearing during fabrication. We thus argue that applying the model to other base\nmaterial combinations would require a new fitting procedure that requires only five samples\n(2 base materials and 3 digital materials) to be tested.\nR2C2: The training data is generated by the theoretical model (1.6M data points) and not\ndirectly based on experimental data. Although the theoretical model has been partially\nvalidated, it may result in a network that is not robust to actual manufacturing errors (e.g.,\nmaterial “fuzziness”).\nRecommendation:\nPage 5 of 12\n- Introduce test data from actual print samples as part of the training set to enhance the\ngeneralization of the model.\n- Provide a sensitivity analysis of the network under noisy data (e.g., prediction error after\nrandomly perturbing the mixing ratio).\nWhile the ANN is trained on the theoretical model data, we have taken steps to assess its\nrobustness in realistic settings. We analytically quantified the standard deviation in effective\nmixture ratio of the base materials due to the stochastic base material assignment and\nevaluated the expected difference in estimated stiffness in Section A.4 of the supplementary\ninformation. This quantifies the stiffness misestimation compared to experimental conditions\nand also with respect to feature size. The ANN is trained on model data since the\nexperimental data from five samples leaves significant data gaps between the measured\nmaterial mixes. The model data thus fills these gaps by physically modeling the behavior\nwith the de Gennes percolation theory. Inclusion of raw experimental data is avoided since\nmixing model and experimental data leads to a noisier dataset and duplicate entries. This\nwill effectively lead to worse ANN training.\nR2C3: The design requirements for the wrist orthosis (e.g., maximum displacement of 2mm\nat 15N force) are based on literature 2929 but do not incorporate actual biomechanical data\nfrom the patient (e.g., joint stiffness for different pathologies). Fatigue properties for long-\nterm use (e.g., material relaxation after cyclic loading) were also not validated.\nRecommendation:\n- Supplement patient-specific biomechanical data or cite broader clinical studies to support\ndesign criteria.\n- Add long-term fatigue testing (e.g., change in stiffness after 10,000 cycles of loading).\nA current state-of-the-art wrist orthosis design is either manually manufactured as a plaster\ncast or additively manufacturing as a rigid orthosis. The stiffness of the orthosis is typically\ntuned heuristically on a patient-by-patient basis. To our knowledge, patient-specific\nbiomechanical data does not currently exist in the public domain and has not been quantified\nin a generalized approach in literature. The design requirements in this paper are derived\nfrom the body of available literature. This assumption and the choice of design variables is\nclarified in the revised manuscript. We hope that this clarification highlights our focus on\nvalidating the design and modeling frameworks under generalizable conditions, while future\nwork may incorporate patient-specific data when available.\nTo address the concern regarding long-term material durability, additional experimental\nvalidation of fatigue performance under cycling loading is added to the supplementary\ninformation (section E.5). A standardized dogbone sample made from orthosis shell material\nis tested and showed no significant reduction in stiffness after 10’000 cycles of loading. This\ndemonstrates the materials' suitability for long-term cyclic applications, such as orthoses.\nR2C4: The paper points out the “blurring” phenomenon during material injection (Figure 2a),\nbut does not analyze its impact on simulation accuracy. For example, the actual printed\ngradient transition region may deviate from the linear mixing assumption of the design.\nRecommendation:\n- Quantify manufacturing errors (e.g., analyze actual material distribution by micro-CT) and\nintroduce tolerance parameters in the model.\nPage 6 of 12\n- Discuss the extent to which manufacturing errors affect case study performance (e.g.\ndamper hysteresis curves).\nWe have added both experimental and analytical analyses to quantify the effect of blurring,\nstochasticity, and manufacturing variability. Detailed analysis is shown in the supplementary\ninformation (A.3 and A.4). The findings are also summarized in the main text.\nFirst, we analyze the additively manufactured morphology of stochastic digital material\nacross varying base material ratios. A visual inspection confirms a clear transition from\nconnected to dispersed regions between 15% and 20% VW content. This partially validates\nthe literature-derived assumption of the percolation threshold. A new dogbone sample at\n16.5% VW content is added to the percolation fitting procedure to more accurately model the\ncritical percolation threshold (Fig. 3). Further, it provides an empirical understanding of how\nbase material transitions evolve during the additive manufacturing process, including blurred\nboundaries between regions.\nSecond, we conduct a repeatability study, where we distinguish between two sources of\nvariability: manufacturing process variation and stochastic voxel layout. Mechanical testing\nof 20 samples with identical voxel layout shows a manufacturing-induced variability of 3% in\nreaction force. Characterization of 20 samples with randomized voxel layouts, but with the\nsame base material mix ratio, shows variability of 2.5%. This shows that blurring and\nstochasticity have limited effects for feature sizes of the tested magnitude relative to\nmanufacturing accuracy.\nTo further assess the influence of stochasticity at smaller scales, we conduct an analytical\nstudy. We model the stochastic assignment as a binomial process and quantify the standard\ndeviation in the effective mixture ratio relative to the target base material mix ratio as a\nfunction of feature size. We propagate this error through the ANN to evaluate the mismatch\nin estimated material stiffness. We deduce that the stiffness error increases significantly for\nsub-millimeter features and shadows the manufacturing variability of 3% for features below\n1-2mm, depending on the base material mix ratio.\nWhile quantification of actual base material distributions within a digital material with 3D\nvisualization methods is an excellent next step, our current combination of morphological\ninspection, experimental testing, and analytical modeling provides a validated feature size\nthreshold for digital material assignment. Importantly, all features in the damper and orthosis\ncase studies are above the critical threshold of 1-2mm, ensuring that base material\nallocation variability does not critically impact the case study results.\nR2C5: Simulation parameters (e.g., mesh sensitivity analysis, time step selection) are not\nspecified, and no mention is made of the convergence of multi-threaded calculations (the\nsupplementary material mentions a single-threaded Abaqus run). In addition, the parameter\nfitting procedure for the hyperelastic model (Mooney-Rivlin) is not disclosed.\nRecommendation:\n- Add grid-independent validation (e.g., comparison of displacement errors at different grid\ndensities).\n- Disclosure of the fitting method and experimental data support for the parameters of the\nhyperelastic model.\nPage 7 of 12\nAdditional simulation parameters are added to the Supplementary Information for both case\nstudies. As mentioned in the Supplementary Information, multi-threaded runs using the\nAbaqus viscoelastic step show significant convergence issues. We observe that models\nexhibiting stable convergence in a single-threaded environment often failed to converge\nunder multi-threaded execution even when subject to significant time step reduction or\nrelaxed viscoelastic strain error tolerances. This could also be caused by the custom\nsubroutine used to model Arrhenius and WLF behavior. However, to the best of our\nknowledge, the current version of the subroutine is thread safe.\nA specific mesh convergence study is not conducted as there is an interaction between\nfeature size, coarse grid spacing, and element sizing in the FE simulation. The coarse grid\nneeds to be chosen such that it can accurately assign material at the minimum feature size\nof the current structure. Material is assigned to the FE mesh based on the materials\nassigned to coarse grid points. The sizing of the FE mesh is thus directly related to the sizing\nof the coarse grid. A full mesh sensitivity analysis exceeds the scope of this work as the\nfocus of this work is mainly on the material model and less on the specifics of the FE model.\nThe Mooney-Rivlin hyperelastic model is fit at runtime by the Abaqus solver. The\nhyperelastic data is passed to Abaqus in the form shown in Table S2 with the setting\n“uniaxial test data”. The hyperelastic curve of the material mixtures is determined based on\nthe hyperelastic curves of the two base materials (AB/VW in Table S2) and a mixture ratio\nbased on the long-term response of the percolated master curve relative to the long-term\nresponse of the base materials. The online methods section has been adapted to better\nreflect this step.\nR2C6: The paper emphasizes the “high efficiency” of the method (7 orders of magnitude\nreduction in data points) but does not discuss the practical limitations in its engineering\napplications, e.g:\n- The demand on computational resources for high-resolution voxel design (e.g.,\npreprocessing time for 280,000 grid points).\n- Cost and scalability of multi-material printing (e.g., difficulty in removing the support\nmaterial SUP706).\nRecommendation:\n- Add an analysis of the trade-off between computational efficiency and accuracy (e.g., mesh\nresolution vs. simulation error).\n- Discuss economic challenges in industrial-scale applications (e.g., material waste rate,\nprint time).\nWe added a paragraph to the discussion containing an assessment of the trade-offs\ninvolved. We now clarify the computational bottleneck and how to address it. Together with\nthe response to [R2C4], where we analyzed feature size influencing manufacturing\naccuracy, we provide some practical guidelines for voxel-based design. On the\nmanufacturing side, we have placed the proposed framework in the current additive\nmanufacturing environment and provided recommendations on when to use it.\nFurther, we added a section (B.11, including Fig. S13) to the supplementary information that\nshows the runtime for generating the designs of both case studies. The digital material\nPage 8 of 12\nassignment on the coarse grid is very fast. This shows the framework's capability for fast\ndesign iterations.\nPage 9 of 12\nReviewer #3: This manuscript presents a novel inverse design framework for voxelated\nthermo-viscoelastic materials. This approach effectively reduces computational complexity\nvia mesoscale homogenization and demonstrates promising predictive accuracy in two\nfunctional case studies. While this work offers a valuable contribution to data-driven design\nof architected soft materials, its limited validation is a significant weakness. Therefore, I\nwould like to recommend it for publication after addressing some concerns. Detailed\ncomments are provided below.\nR3C1: The manuscript relies heavily on simulation results without providing adequate\nexperimental validation, which may raise questions about the practical applicability of the\nfindings. The Prony series fitting is based on only four experimental datasets for training and\na single dataset for validation, which may be insufficient to robustly capture the full range of\nmaterial behavior, as acknowledged by the authors. Moreover, the inverse design framework\nrelies solely on simulations. If possible, incorporating such experimental validation would\nsignificantly enhance the credibility of the approach.\nAdditional experimental validation samples have been measured and incorporated into both\nthe fitting and validation procedures. Fig. 3 has been reworked to include these new\nsamples. The percolated material model is now compared to five digital materials (Fig. 3d),\ntwo of these digital materials are used to fit the percolation coefficient 𝜇 , one is used only for\nthe fitting of the critical percolation threshold 𝜌 , and two are used only for validation\n𝑐\npurposes. The entire material response is now split into six equally spaced bands by the two\nbase materials and five stochastic digital materials. The dataset used to train both ANNs is\ngenerated from this fitted percolation model and is based on experimental data and validated\nby the samples not included in the fitting procedure.\nAdditionally, the damper case study in Fig. 5 has been expanded to include testing at\nmultiple temperatures and frequencies, thus further validating the model and simulation\napproach based on the ANNs. The orthosis case study has also had additional simulations\nintroduced at temperatures relevant to its application.\nR3C2: Although the authors emphasize the importance of temperature-sensitive material\nresponses, most of the work are conducted at room temperature. Extending the analyses to\ninclude elevated or variable temperatures, especially those relevant to the case studies,\nwould better demonstrate the robustness and applicability of the model.\nThe damper case study and orthosis study have been extended to better show the\ncapabilities of the model. The damper is now validated at two different temperatures (15°C\nand 45°C) with loading at three different frequencies (0.1Hz, 1Hz, 10Hz) for each\ntemperature. This shows the predictive capabilities of the modeling approach for the\ncomplex load case of hysteresis. The orthosis is also validated at expected ambient\ntemperatures and skin temperature. A cold environment simulation at 0°C is conducted, and\na skin temperature validation at 37°C is also performed.\nR3C3: I would recommend replacing the term “graduated” with “graded” to reduce potential\nconfusion among interdisciplinary readers, as it is the more commonly accepted term in the\nrelated field.\nPage 10 of 12\nWe agree on using “graded” to enhance clarity for a broader audience. All terms have been\nadjusted accordingly in the manuscript and the supplementary information.\nR3C4: Page 7, line 186: please check the error.\nThe error is a reference with the wrong hyperlink. We corrected it accordingly.\nPage 11 of 12\nReferences\n1. Saldívar, M. C., Doubrovski, E. L., Mirzaali, M. J. & Zadpoor, A. A. Nonlinear coarse-\ngraining models for 3D printed multi-material biomimetic composites. Addit Manuf 58,\n103062 (2022).\n2. Dong, L. & Wang, D. Optimal Design of Three-Dimensional Voxel Printed\nMultimaterial Lattice Metamaterials via Machine Learning and Evolutionary Algorithm.\nPhys Rev Appl 18, 054050 (2022).\n3. Yuan, C., Wang, F., Rosen, D. W. & Ge, Q. Voxel design of additively manufactured\ndigital material with customized thermomechanical properties. Mater Des 197, 109205\n(2021).\nPage 12 of 12"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have revised according to the reviewer's comments and the revised manuscript is improved. Thus I have no other comments for the paper.", "reviewer_2": "Cheers! All my comments have been responded to. The paper contains useful results and is ready for publication. Congratulations and looking forward to further research results! It is a great honor to receive an invitation from [Redacted] to review this manuscript from Nature Communications.\nBased on my current capabilities, the code in this paper is highly feasible, reproducible, and accurate, and can support the publication of this paper. Best, Yuze Nian. August 12, 2025.", "reviewer_3": "The authors have revised the manuscript in response to the initial review, addressing several of the previously raised concerns. The integration of experimental data to validate the simulation results of the storage modulus of voxelated materials under frequency- and temperature-dependent conditions enhances the credibility of the proposed approach. The inclusion of multiple temperature conditions in the case study analyses is appropriate, as it highlights the method’s ability to account for different loading scenarios. Changing the term graduated to graded has also minimized potential confusion for readers, aligning the terminology with common usage in the field. Overall, the manuscript has improved meaningfully, and the authors have been responsive to reviewers feedback.\nThat said, I have a few remaining or new points for consideration before recommending acceptance:\n• After the mixture ratio is determined, it would be helpful to clearly describe how the mesoscale unit size (coarse grid size) for applying the stochastic distribution is selected in practice. Providing explicit criteria or a decision process for choosing the size and shape of these units, preferably in a way that allows readers to visualize the workflow, would help future researchers replicate the proposed method.\n• In the damper experiment, the error increased to 26.6% under the 10 Hz condition. It would strengthen the manuscript to include an analysis of the cause of this reduced accuracy in the high-frequency range and to discuss potential improvement strategies (e.g., augmenting ANN training data or adjusting model parameters).\n• In Figure 2a, the term graduated is still used. Please replace it with graded for consistency with the rest of the manuscript.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to the Reviewers’ Comments for Manuscript:\n“Inverse Design of Stochastic, Voxelated Thermo-Viscoelastic Digital\nMaterials”\nDear Reviewers,\nThank you very much for the helpful comments to further improve the quality of our\nmanuscript. We have carefully considered every comment and address them here.\nWe further made additional minor changes to the manuscript to improve the readability and\nclarity based on editorial input.\nThe Authors\nPage 1 of 3\nReviewer #1: The authors have revised according to the reviewer's comments and the\nrevised manuscript is improved. Thus I have no other comments for the paper.\nReviewer #2: Cheers! All my comments have been responded to. The paper contains useful results and is ready for publication. Congratulations and looking forward to further research\nresults! It is a great honor to receive an invitation from [Redacted] to review this manuscript\nfrom Nature Communications.\nBased on my current capabilities, the code in this paper is highly feasible, reproducible, and\naccurate, and can support the publication of this paper.\nReviewer #3: The authors have revised the manuscript in response to the initial review,\naddressing several of the previously raised concerns. The integration of experimental data to\nvalidate the simulation results of the storage modulus of voxelated materials under\nfrequency- and temperature-dependent conditions enhances the credibility of the proposed\napproach. The inclusion of multiple temperature conditions in the case study analyses is\nappropriate, as it highlights the method’s ability to account for different loading scenarios.\nChanging the term graduated to graded has also minimized potential confusion for readers,\naligning the terminology with common usage in the field. Overall, the manuscript has\nimproved meaningfully, and the authors have been responsive to reviewers feedback.\nThat said, I have a few remaining or new points for consideration before recommending\nacceptance:\nR3C1: After the mixture ratio is determined, it would be helpful to clearly describe how the\nmesoscale unit size (coarse grid size) for applying the stochastic distribution is selected in\npractice. Providing explicit criteria or a decision process for choosing the size and shape of\nthese units, preferably in a way that allows readers to visualize the workflow, would help\nfuture researchers replicate the proposed method.\nThe corresponding description in the main manuscript has been rewritten to better highlight\nhow the coarse grid resolution is chosen. In practice, the coarse grid resolution is adapted to\nthe size of the geometric features in the structure, the desired material assignment resolution\nand the FE mesh resolution required for simulation.\nAdditionally, Section B.4 in the supplementary information has been expanded to discuss\nthese three key considerations in more detail.\nR3C2: In the damper experiment, the error increased to 26.6% under the 10 Hz condition. It\nwould strengthen the manuscript to include an analysis of the cause of this reduced\naccuracy in the high-frequency range and to discuss potential improvement strategies (e.g.,\naugmenting ANN training data or adjusting model parameters).\nThe discussion has been expanded to include additional information on what drives this error\nand what measure could be used to address it.\nAt higher frequencies the materials are generally in a region of larger stiffness gradients.\nThis means that small inaccuracies in the material model fitting procedure and specifically\ntemperature shifting can cause a more significant mismatch between the experiment and the\nPage 2 of 3\nsimulation. The results at 10 Hz and 25C are also less accurate than the results at 15C.\nThis is driven by the limitations of the Maxwell model and shifting functions. The WLF and\nArrhenius function are limited in their modeling of frequency-dependent shifting biases. The\nANN response can model more complex behavior in terms of temperature shifting than the\nfinal Maxwell model as shown in Fig. 4f. A possible solution that could be used to address\nthis issue is the use of a custom temperature shifting function for the Maxwell model based\non the percolated material data.\nR3C3: In Figure 2a, the term graduated is still used. Please replace it with graded for\nconsistency with the rest of the manuscript.\nFigure 2a has been updated to the new terminology.\nReviewer #4: I co-reviewed this manuscript with one of the reviewers who provided the\nlisted reports. This is part of the Nature Communications initiative to facilitate training in peer\nreview and to provide appropriate recognition for Early Career Researchers who co-review\nmanuscripts.\nThe repository provides clear and sufficient instructions, and I was able to run the code\nwithout major issues. The results presented in the paper are reproducible using the provided\nimplementation.\nPage 3 of 3"}]

Decision letter