Machine learning guided resolution of mechanical trade-off in polymer composites via stress adaptive interface

Domain: physics
Published: 2026-02-24
Editorial outcome: major_revision

Abstract
Developing polymer composites that simultaneously achieve high strength, toughness, and impact resistance remains a fundamental challenge due to inherent trade-offs and brittle interfacial failure. Here, we propose a universal toughening strategy that integrates a bone-inspired trabecular interlock architecture with a thermodynamically driven, stress-adaptive interface to enable efficient energy dissipation under mechanical loading. To address multi-objective optimization in composites design, we further develop a data-driven framework combining Pareto Set Learning and Active Learning, which systematically explores the composition–performance landscape to identify balanced, high-performance formulations. The optimized composites exhibit synergistic mechanical properties: strength up to 250 MPa, fracture toughness exceeding 14 MPa·m1/2, and impact resistance of nearly 4.8 J, surpassing most bioinspired and engineered polymer counterparts. The strategy is scalable, chemically versatile, and broadly applicable, offering a programmable route to next-generation lightweight composites for aerospace, transportation, and protection.

Introduction
Achieving synergistic enhancement of strength, toughness, and impact resistance remains a fundamental challenge in structural polymer composites, primarily due to intrinsic trade-offs and interfacial failure under complex loading conditions1–3. While conventional approaches4–12, such as incorporating stiff fibers, nanofillers, or layered architectures, can improve individual properties, they often induce brittle zones or interfacial mismatches that degrade overall mechanical performance.

Nature offers an alternative paradigm where hierarchical architectures and adaptive interfaces synergistically dissipate energy across multiple scales. Although bioinspired designs mimicking bone, nacre, and other natural systems have emerged13–18, most rely on static interfacial reinforcements, including micro-architectured interlocking19, surface functionalized fillers20, polymer grafting21, supramolecular/covalent/ionic bonding22,23, and nanostructured coatings24. Critically, such chemically constrained interfaces lack dynamic stress responsiveness and broad materials applicability.

Herein, we propose a universal toughening strategy that integrates a bone-inspired trabecular interlock architecture with a thermodynamically driven, stress-adaptive interface. This multiscale design enables dynamic energy dissipation, effectively overcoming the classic strength–toughness trade-off. To address multi-objective optimization challenges, we further develop a machine learning (ML) framework combining Pareto Set Learning (PSL) and Active Learning (AL). This data-driven approach efficiently explores the composition–performance landscape to identify balanced, high-performance formulations. Demonstrated across diverse commercial matrices (Polyketone, Polyamide 6) and reinforcements (1D fibers to 2D nanosheets), our strategy establishes a programmable pathway for advanced composite design. Beyond exceptional properties, this work advances mechanistic understanding of architecture–interface–property relationships and provides a generalizable platform for accelerated materials discovery.

Results
Biomimetic trabecular architecture with thermodynamically adaptive interfaces
Inspired by the hierarchical porosity and mechanical efficiency of natural trabecular bone1,2, we establish a biomimetic composite design paradigm governed by thermodynamic principles (ΔG=ΔH−TΔS)25,26. Central to this strategy is the reduction of interfacial enthalpy (ΔH) during the thermal blending, which lowers the energetic barrier for polymer chain interdiffusion and promotes the formation of stable physical interlocks at soft–hard interfaces. This concept is realized in the trabecular interlocked composites (TICs), characterized by aligned reinforcing phases embedded within a robust three-dimensional porous skeleton (Fig. 1b). The stress-adaptive interfaces in these composites rely on dynamic diffusion and entanglement of polymer chains with rigid constituents under thermal activation, thus transitioning from traditional static bonding toward a thermodynamically governed entanglement strategy.

The effectiveness of this approach is experimentally validated in a ternary system, where thermoplastic polyurethane (TPU) functions as the entanglement promoter, poly(phenylene ether ketone) (POK) serves as the rigid matrix, and aligned basalt fibers (BF) provide structural reinforcement. Analogous to natural bone structures, TPU forms collagen-like entangled networks, while the POK phase provides structural rigidity, and the shear-aligned BF creates the trabecular scaffold. The hierarchical porous architecture is fabricated via gas-assisted pore formation under controlled thermo-shear conditions. Upon thermal activation, TPU chains become highly mobile and gradually diffuse into the POK matrix and along the BF interfaces. This interdiffusion lowers the interfacial enthalpy (ΔH), thereby reducing the energetic barrier for physical entanglement formation. The resulting dynamic, entangled interfaces enhance compatibility and mechanical integrity. These phenomena are substantiated by thermal analyses, structural characterizations (Figure S5), and molecular dynamics simulations (Fig. S6 and Supplementary Note 2).

Machine learning (ML) driven multi-objective optimization
While the feasibility of the trabecular interlock architecture and thermodynamically driven interfacial strategy has been established, the compositional design space for these polymer composites remains extensive and intrinsically high-dimensional. Traditionally, optimized polymer composite formulations are developed through empirical trial-and-error methods27,28, which involve adjusting one variable at a time (e.g., x in Fig. 2a) while keeping others constant. Performance metrics such as f(x) improve gradually through repeated cycles of formulation, testing, and adjustment. However, conventional approaches become inefficient when addressing multiple mechanical objectives with competing trade-offs.

To address the intrinsic trade-offs among strength, toughness, and impact resistance, we introduce a machine learning (ML)–assisted multi-objective optimization framework (Fig. 2b). This data-driven approach systematically explores the composition–performance landscape, identifies approximate Pareto-optimal solutions, and resolves competing performance objectives. Consequently, the framework significantly accelerates the discovery and refinement of high-performance composite formulations.

Pareto set learning (PSL)
An initial baseline dataset is established using a design-of-experiments (DoE) strategy by selecting 50 representative composite formulations spanning the feasible composition space. Each formulation is experimentally fabricated and mechanically characterized to measure tensile strength, fracture toughness, and impact energy dissipation. These data are used to train three independent Gaussian process regression (GPR) surrogate models.

As the training set size increases from 30 to 50, cross-validated predictive performance improves systematically for all three targets, as evidenced by increasing coefficients of determination (R2), decreasing mean absolute errors (MAE), and decreasing mean squared errors (MSE), together with a concurrent reduction in predictive uncertainty quantified by the standard deviation across validation folds (Supplementary Note 3.1 and Fig. S7). Although toughness-related properties exhibit comparatively higher uncertainty at smaller sample sizes, all surrogate models attain statistically reliable accuracy at 50 samples, thereby establishing a robust foundation for subsequent multi-objective optimization.

Building on the validated GPR surrogates, multi-objective optimization is performed using Pareto set learning (PSL). PSL learns a mapping from preference vectors λ, which encode user-defined trade-offs among competing objectives, to optimized formulations x by scalarizing the multi-objective problem into a tractable single-objective form29. By sampling 1000 preference vectors, the trained PSL model generates a diverse and structured set of Pareto-consistent candidate solutions spanning a broad range of engineering priorities (Supplementary Note 3.2 and Fig. S8).

Active learning (AL)-guided optimization
To further improve optimization efficiency and solution quality, we integrate active learning (AL) with multi-objective Bayesian optimization to iteratively refine the learned Pareto set and guide exploration of the trabecular interlocked composites (TICs) design space30,31. Specifically, we adopt an acquisition strategy that combines a lower confidence bound (LCB) criterion with hypervolume improvement (HVI) to select small batches of informative and non-dominated candidates from the PSL-generated solution pool for experimental validation (Supplementary Note 3.3). This uncertainty-aware strategy balances exploitation of high-performing regions with exploration of poorly sampled areas, enabling targeted expansion of the experimentally validated performance envelope under a fixed experimental budget.

Starting from the baseline dataset (Fig. 3a), the experimentally sampled formulations in the strength–fracture toughness space follow an approximately monotonic trend, with impact energy dissipation co-varying with the high-performance region. However, the initial data are concentrated around intermediate strength (≈220–235 MPa) and fracture toughness (≈12–13.5 MPa·m1/2), leaving the extreme high-performance boundary sparsely sampled. With successive AL iterations (Figs. 3b and S9), newly selected formulations progressively populate the upper-right region of the design space, corresponding to tensile strengths exceeding ~245 MPa and fracture toughness values above ~14.5 MPa m1/2, while maintaining high impact energy dissipation (>4.6 J). Furthermore, the selected candidates are not clustered in regions of high sample density but are preferentially drawn toward areas associated with elevated predictive uncertainty and potential Pareto-front expansion, consistent with the intended uncertainty-aware acquisition strategy.

The efficiency of this frontier expansion is quantified by hypervolume (HV) analysis. More than 90% of the total normalized HV gain is achieved within the first ten additional experiments, driven by a dominant early HV increase (ΔHV≈0.012), followed by rapid saturation with marginal gains below 0.001 per iteration. Across five AL iterations, learning-curve analysis further reveals systematic improvements in cross-validated predictive accuracy and progressive stabilization of the surrogate models, particularly for fracture toughness and impact energy dissipation (Fig. S10). Independent hold-out testing confirms robust out-of-sample generalization for all three properties, with no evidence of systematic prediction bias (Fig. S11).

Finally, benchmarking against a classical response surface method (RSM) under identical experimental budgets (50 initial formulations plus 25 additional trials) demonstrates that the combined PSL+AL framework consistently drives the Pareto front toward higher-performance regimes and achieves over 90% of its terminal normalized hypervolume improvement within the first ten additional experiments. In contrast, RSM yields only marginal performance gains even after exhausting the full experimental budget (Fig. S12). These results establish that uncertainty-aware, iterative learning enables substantially higher sample efficiency and reliability than conventional quadratic design-of-experiments baselines.

Preference-guided solution extraction
After completion of five active learning iterations (25 additional experiments), the PSL model is retrained using the expanded dataset, yielding the optimized Pareto set shown in Fig. 3d–f. The resulting Pareto-consistent solutions exhibit a pronounced and organized structure across both composition and performance spaces. In the TPU–BF–POK composition space (Fig. 3d), the solutions are confined to a relatively narrow compositional window rather than being dispersed across the full admissible domain, revealing a well-defined high-performance basin.

A consistent structural organization is observed in the three-dimensional performance space (Fig. 3e), where the Pareto-consistent solutions form a smooth, manifold-like trajectory instead of a diffuse cloud. This behavior indicates that the achievable multi-objective trade-offs are constrained by a common structural regime, within which improvements in tensile strength are accompanied by predictable and coupled variations in fracture toughness and impact energy dissipation, rather than arbitrary performance fluctuations.

This structured Pareto manifold enables preference-guided engineering decision making (Fig. 3f and Table S2). By specifying target preference vectors, representative solutions corresponding to balanced, strength-prioritized, fracture-toughness–prioritized, and impact-energy–prioritized design modes can be directly extracted along the learned manifold. The close agreement between the prescribed preference vectors and their matched counterparts confirms that PSL preserves preference information during solution generation, allowing explicit design priorities to be translated into concrete compositions and predicted performance outcomes.

These results demonstrate that the combined AL–PSL framework not only accelerates the discovery of high-performance Pareto fronts but also constructs a mechanistically admissible and queryable design manifold, enabling practical and transparent composite formulation design under competing performance objectives.

Experimental validation of preference-guided Pareto solutions
To validate the machine-learning–selected Pareto solutions experimentally, we fabricated four representative formulations corresponding to distinct design preferences (balanced, strength-first, fracture-first, and impact-first). Across tensile, notched tensile, and impact loading, the macroscopic response curves remain highly consistent (Fig. 4a–c), with comparable yielding, strain hardening, and terminal failure behavior. This similarity indicates that preference-guided selection operates within a narrowly defined high-performance basin, rather than shifting the material into qualitatively different deformation regimes.

Quantitatively, the absolute properties cluster tightly. Tensile strength spans approximately 248–250 MPa (maximum deviation ~1.2 MPa relative to the balanced formulation, <0.5%, Fig. 4d, g), and fracture toughness remains near 14.8–14.9 MPa m1/2 (differences ~0.02 MPa m1/2, Fig. 4e, h), both within experimental scatter. Impact energy dissipation is similarly concentrated (~4.65–4.71 J, Fig. 4f, i), with the impact-first design exhibiting a small but reproducible increase (~0.06 J, ~1%) relative to the balanced reference. Expressed as relative differences with respect to the balanced formulation, the shifts remain modest yet directionally consistent with the prescribed preferences, supporting that preference guidance primarily refines trade-off ordering within an already optimized performance envelope.

To contextualize these validated formulations within a broader materials landscape, we benchmark their mass-normalized performance against representative natural and engineered materials using Ashby-type maps (Fig. 4j–l; Tables S3 and S4). All four designs occupy a narrow yet high-performance region across the examined property spaces, consistent with convergence toward a Pareto-optimal plateau governed by a shared deformation and energy-dissipation mechanism. In the specific strength–specific fracture toughness map (Fig. 4j), the optimized TICs fall within the advanced composite and bone-like regime, positioned well above conventional polymers and non-technical ceramics. When density is explicitly considered (Fig. 4k), all formulations remain in a lightweight regime (<1000 kg/m3) while maintaining fracture toughness on the order of 10–20 MPa m1/2, highlighting a rare combination of low mass and damage tolerance. A similar clustering is observed in the specific impact force–specific energy dissipation map (Fig. 4l), where the TICs exhibit concurrent enhancement of impact resistance and energy absorption, comparable to or exceeding reported Bouligand and gradient composite architectures.

The close grouping of the four designs across these Ashby maps therefore reflects convergence toward a Pareto-optimal performance plateau governed by a shared deformation and energy-dissipation mechanism, rather than experimental scatter. This convergence underscores the robustness of the preference-guided optimization and provides a reliable basis for Ashby-map–guided materials comparison and selection.

Mechanism analysis
Dynamic stress-adaptive interfacial response
The tight clustering of validated properties and the plateau-like distribution on Ashby maps suggest that the optimized designs share a common, stable dissipation mechanism. We therefore interrogate the dynamic interfacial response of a representative TIC formulation (balanced preference) using in situ tensile SEM and strain-resolved Raman spectroscopy, complemented by impact/ballistic observations. This formulation sustains tensile strains up to ~600% without catastrophic failure (Fig. S13). In drop-weight impact tests, high-speed imaging captures rapid deformation followed by rebound (Fig. S14a and Movie S1), and post-impact inspection shows only shallow indentation without macroscopic cracking. By contrast, control specimens lacking the trabecular-interlocked architecture and adaptive interface design (POK–BF NAFCs) exhibit immediate penetration under identical conditions (Fig. S14b). Under ballistic loading (80 m/s), TICs decelerate the projectile and rebound without perforation (Figure S15a), whereas NAFCs show complete penetration (Fig. S15b and Movie S2). These macroscopic observations motivate the following in situ analyses aimed at resolving how stress is redistributed and dissipated at the crack and interfacial scales.

In situ tensile SEM directly visualizes the crack evolution and load-transfer mechanisms operating in TICs. At an applied strain of ~300%, crack opening is effectively shielded by pronounced fiber-mediated bridging, where multiple polymer ligaments anchored on basalt fibers span the crack faces and remain mechanically engaged (Fig. 5a, yellow arrows). Instead of propagating along a straight path, the crack follows a tortuous and deflected trajectory, indicating crack deflection induced by interfacial heterogeneity and mechanical interlocking. Upon further loading to ~500% strain (Fig. 5b and Movie S3), the bridging mechanism remains operative. The load-bearing ligaments undergo progressive stretching and thinning with delayed rupture, rather than abrupt interfacial debonding or clean fiber pull-out. Crack advance remains gradual and continues along a deflected path, accompanied by a pronounced expansion of the surrounding plastic zone (Fig. 5b, blue arrows). These observations demonstrate that the TICs interface does not fail through brittle separation. Instead, it sustains continuous mechanical engagement across extreme strains through physically entangled, load-bearing ligaments that redistribute stress in the crack-tip region, thereby stabilizing crack growth and enabling sustained energy dissipation.

We next quantify molecular-scale stress accommodation using strain-resolved tensile Raman spectroscopy. The stacked spectra (Fig. 5c) reveal a continuous evolution of the full Raman window with increasing strain. Consistently, the strain–wavenumber heatmap (Fig. 5d) highlights two distinct response regimes: the carbonyl window (1700–1785 cm−1) exhibits strain-dependent but spatially diffuse intensity redistribution, whereas the CH-stretch window (2800–3100 cm−1) progressively develops increased spectral dispersion at higher strain, indicative of enhanced conformational heterogeneity. To quantify the response of polar groups, we track the centroid position of the carbonyl band (Fig. 5e). In TICs, the centroid shift remains relatively small and non-monotonic across the applied strain range, indicating that tensile stress does not continuously accumulate on a narrow subset of carbonyl bonds. Instead, stress is dynamically repartitioned as the physically entangled interface undergoes progressive reconfiguration under load. In contrast, the POK–TPU control (Fig. S17) exhibits a more systematic and monotonic centroid evolution, consistent with a more homogeneous polymer response in which bond-level stress localization can be directly tracked spectroscopically. Chain-level dissipation is further captured by the full width at half maximum (FWHM) of the CH-stretch envelope (Fig. 5f). TICs display a pronounced broadening of the CH envelope at high strain, consistent with the activation of segmental mobility, chain rearrangement, and sliding. By comparison, the control sample shows a substantially weaker broadening response, indicating more limited conformational freedom under tensile loading.

These molecular-scale signatures are fully consistent with the in situ tensile SEM observations of sustained fiber-bridging ligaments and an expanded plastic zone at the crack. Together, the spectroscopic and microstructural evidence demonstrates that the trabecular interlocking architecture enables a dynamically activated, stress-adaptive interfacial response, which progressively converts tensile loading into frictional and entanglement-mediated dissipation pathways. This multiscale stress redistribution mechanism underpins the enhanced crack resistance and impact tolerance observed in TICs. At the molecular level, the stress-adaptive interface originates from thermodynamically driven interdiffusion of flexible elastomer chains into the amorphous regions of the rigid polymer matrix and onto the polar fiber surface during thermal processing. This process generates a physically entangled interphase rather than permanent chemical bonding. Under mechanical loading, the entangled interphase undergoes strain-dependent reconfiguration through chain sliding and rearrangement, enabling continuous stress redistribution and adaptive energy dissipation.

Strain-induced structural evolution across length scales
To contextualize the dynamic interfacial responses revealed above within a broader structural framework, we further examine the strain-induced structural evolution of the balanced POK–TPU–BF TICs across multiple length scales. Cross-sectional SEM images (Fig. 6a) reveal a hierarchical interfacial architecture underlying the exceptional damage tolerance of the balanced POK-TPU-BF TICs. From low to high magnification (100 μm to 5 μm), multiscale features, including fiber bridging, crack deflection, interfacial tearing, and polymer fibrillation, are observed. At higher resolution, elongated polymer ligaments around the BF and deformed matrix zones indicate the formation of a physically entangled interface. Unlike classical fiber pull-out, the interface maintains continuous contact between BF and the surrounding TPU-POK matrix, without evidence of brittle decohesion. The aligned orientation of BF (yellow arrows) further supports efficient stress transfer before chain disentanglement and relaxation. These features confirm that, beyond traditional toughening mechanisms such as frictional pull-out or shear yielding32–34, the stress-adaptive reconfiguration of the entangled interface promotes energy dissipation and enables simultaneous enhancement of strength, ductility, and impact resistance.

At the initial stage (0% strain), 2D SAXS and WAXS patterns display nearly circular, isotropic scattering rings (Fig. 6b, c), reflecting a uniform distribution of oriented nanodomains and chain segments without macroscopic alignment. XRD analysis reveals sharp peaks at 2θ≈20.6° and 29.8°, corresponding to crystalline POK reflections, indicating intact semicrystalline lamellae (Fig. S18). Molecular dynamics (MD) simulations confirm a dense, interpenetrated interfacial network, where TPU chains form physically entangled bridges across the fiber-matrix boundary (Fig. 6f, left). Thermal analysis shows a sharp melting peak at 229.1 °C and high enthalpy, consistent with well-preserved crystalline domains (Fig. 6g).

Upon moderate deformation (50% strain), TICs exhibit a marked stiffening response, as indicated by the steep increase in the tensile stress-strain curve up to approximately 100% strain (Fig. 4a). SAXS patterns transition from isotropic rings to pronounced anisotropic streaks (Fig. 6b), while WAXS profiles develop well-defined equatorial arcs (Fig. 6c), evidencing progressive alignment of amorphous and crystalline phases along the loading axis. Concurrently, XRD peaks broaden and intensify (Fig. S18), reflecting enhanced short-range order due to chain alignment. MD simulations (Fig. 6f, middle) further support the formation of this partially ordered state, showing coordinated chain disentanglement without interfacial failure. The decrease in melting enthalpy observed in thermal analysis (Fig. 6g) indicates partial disruption of crystalline lamellae, likely resulting from interfacial stress and deformation-induced reorganization.

At higher deformation levels (300% strain), the stress-strain curve enters a plateau with a slight upward trend (Fig. 6a), marking a transition from structural reinforcement to energy dissipation. The reduction in SAXS anisotropy (Fig. 6b) and the emergence of diffuse WAXS rings (Fig. 6c) indicate the breakdown of pre-oriented domains and a loss of microstructural coherence. Increased SAXS intensity at low q-values (Fig. 6d) reflects nanoscale reorganization and interfacial expansion, consistent with the development of a more compliant, dissipative interface. WAXS peaks broaden and diminish in intensity (Fig. 6e), while XRD profiles display flattened peaks and elevated baseline scattering (Fig. S18), evidencing the collapse of long-range crystalline order. DSC analysis reveals a partial increase in melting enthalpy, suggesting local chain rearrangement and recrystallization upon stress release. MD simulations reveal strain-dependent interfacial transitions (Fig. 6f). Initially, TPU chains interdiffuse and form stable entanglements with POK at the BF interface. At 100% strain, chain disentanglement initiates, weakening interfacial coupling. By 300%, chains become disordered yet mobile, allowing dynamic rearrangement and stress dissipation. This strain-dependent evolution suggests an adaptive interfacial response in which the dominant role of the interface shifts from efficient load transfer toward enhanced energy dissipation through chain rearrangement and sliding. In contrast to chemically bonded interfaces that often fail through abrupt decohesion, the physically entangled interface accommodates large deformation through sustained chain mobility and interfacial continuity, thereby delaying catastrophic failure and contributing to the observed impact tolerance.

Hierarchical energy dissipation under impact loading
The balanced POK-TPU-BF TICs feature a hierarchical architecture that imparts outstanding resistance to both low- and high-velocity impacts. Three-dimensional reconstruction (Fig. 6h) reveals a continuous, interconnected porous matrix interwoven with aligned BF, forming a trabecular framework reminiscent of biological load-bearing structures. Cross-sectional SEM images show that BF traverse the porous matrix and are embedded within soft-phase ligaments, promoting both mechanical interlocking and localized deformation accommodation.

Under low-velocity impact (drop-weight testing), reversible pore collapse and matrix shear, together with energy dissipation at entangled interfaces via chain sliding and segmental rearrangement, enable distributed damping. This combination minimizes indentation and provides rebound resilience, consistent with viscoelastic recovery.

During high-speed impact (projectile penetration), alternative energy dissipation mechanisms prevail. Finite element (FE) simulations (Fig. 6i) show that strain localizes beneath the impact site, confined within the trabecular interface region. Local densification of the trabecular skeleton, fiber pull-out and bridging, and interfacial chain disentanglement work in concert to dissipate kinetic energy. These mechanisms effectively suppress crack propagation and prevent catastrophic failure. In contrast, non-architectured fiber composites (NAFCs, Fig. S19) develop extended strain bands and experience complete perforation under identical conditions, highlighting the improved damage localization and structural integrity of TICs.

These findings demonstrate that the TICs design synergistically integrates porous load-diffusing scaffolds, crack-arresting fiber networks, and stress-responsive entangled interfaces to achieve dual-mode impact tolerance. This multiscale energy dissipation strategy preserves mechanical integrity even at low densities, offering a generalizable approach for high-performance protective materials.

Versatility of the design strategy across diverse polymer composite systems
Building on the thermodynamic and structural foundation of the TICs concept, we demonstrate its universality across chemically and morphologically diverse polymer composites. The strategy extends to industrially relevant thermoplastics such as POK and polyamide 6 (PA6), reinforced with one-dimensional fibers (carbon fiber, CF; basalt fiber, BF) and two-dimensional graphene nanosheets (G), reflecting applications in aerospace, automotive, and protective materials. All combinations are compatible with standard processing methods.

For comparability, all TICs are fabricated using standardized formulations and processing conditions based on the balanced-performance baseline. Despite variations in matrix chemistry and reinforcement geometry, each system achieves concurrent improvements in strength, fracture toughness, and ductility. The composites use industrially compatible thermoplastic elastomers (TPU, SEBS) without chemical coupling agents or compatibilizers. As shown in Figs. 7a–c and S21, both POK-TPU-CF TICs and POK-SEBS-BF TICs display ductile deformation and enhanced mechanical properties, with strengths up to 250 MPa and fracture toughness approaching 15 MPa m1/2. The approach also applies to PA6-based systems, where introducing SEBS as a soft phase yields PA-SEBS-BF TICs with markedly improved toughness and ductility compared to conventional PA6 composites lacking interfacial design (Supplementary Note 6). The design strategy is further validated in 2D nanosheet-reinforced systems, the results are shown in Fig. 7d–f. Incorporating graphene into POK-TPU-G TICs yields over a threefold increase in fracture toughness compared to control composites lacking interfacial and architectural design.

The broad applicability of the TICs approach derives from its reliance on stress-adaptive interfaces and trabecular interlock architectures, governed by thermodynamic chain dynamics and topological constraints rather than specific chemical interactions. Overall, these findings establish TICs as a scalable, material-agnostic platform for developing lightweight composites with high mechanical resilience, enabling broad implementation across structural applications in aerospace, transportation, and defense.

Methods
Materials
Poly(phenylene ether ketone) (POK, Mw≈3.2×105) was used as the rigid thermoplastic matrix. Thermoplastic polyurethane (TPU, Shore A≈85) and styrene–ethylene–butylene–styrene (SEBS) were employed as soft phases to promote interfacial chain mobility. Basalt fibers (BF, average diameter≈13 μm, length 1–2 cm) and carbon fibers (CF, diameter≈7 μm) were used as one-dimensional reinforcements, while graphene nanosheets (average thickness≈12 nm, lateral size≈5 μm) served as two-dimensional fillers. Polyamide 6 (PA6) was used as an alternative matrix for universality validation. All materials were used as received. Detailed material specifications are provided in Supplementary Note 1.

Composites fabrication
Trabecular interlock composites (TICs) were fabricated via a two-step melt-processing strategy designed to activate thermodynamically driven interfacial entanglement and induce a trabecular architecture. Prior to processing, all polymer components were dried at 80 °C for 12 h under vacuum to remove residual moisture. In the first step, the soft-phase polymer (TPU or SEBS) was melt-blended with the reinforcing phase (BF, CF, or graphene) in a twin-screw extruder at 230 °C, enabling thermal activation of flexible chains and their diffusion onto reinforcement surfaces. In the second step, the thermoplastic matrix (POK or PA6) was introduced and further melt-blended under controlled shear to induce fiber alignment and porous skeleton formation. The extrudates were drawn at a constant speed to promote shear-induced orientation and subsequently cooled to room temperature.

For panel fabrication, extruded strips were stacked in a steel mold and hot-pressed at 230 °C under pressure to obtain composite sheets suitable for mechanical testing and structural characterization. Control composites without interfacial and architectural design (denoted as NAFCs) were prepared under identical thermal conditions but without the sequential soft-phase-assisted interfacial activation. Detailed processing parameters are provided in Supplementary Note 1.

Structural and thermal characterization
Differential scanning calorimetry (DSC) was performed using a TA Instruments Q2000 under nitrogen atmosphere. Samples were heated from 50 °C to 300 °C at a rate of 10 °C min−1, and the second heating cycle was used for analysis. X-ray diffraction (XRD) measurements were conducted on a Rigaku SmartLab diffractometer with Cu Kα radiation (λ=1.5406 Å) over a 2θ range of 5°–40°. Small- and wide-angle X-ray scattering (SAXS/WAXS) experiments were carried out using a XEUSS SAXS/WAXS system to probe nanoscale structural evolution. Fourier-transform infrared (FTIR) spectroscopy was performed in the range of 4000–400 cm−1. Fracture morphologies and microstructural features were examined using scanning electron microscopy (SEM), and three-dimensional porous architectures were visualized by X-ray computed tomography. Detailed characterization procedures are described in Supplementary Note 1.

Mechanical testing
Uniaxial tensile tests were performed on an Instron universal testing machine at a constant crosshead speed of 3 mm min−1. All specimens were tested under identical gauge length and ambient conditions, with at least three specimens per formulation to ensure reproducibility. Fracture toughness was evaluated using single-edge-notched tension tests in accordance with ASTM D5045. Critical stress intensity factors were calculated following standard procedures. Low-velocity impact resistance was assessed using a drop-weight impact tester equipped with a hemispherical striker, while high-speed impact performance was evaluated using a gas-driven projectile system. Detailed test geometries and data reduction methods are provided in Supplementary Note 1.

In situ characterization
In situ tensile SEM experiments were conducted to visualize crack evolution and interfacial deformation during loading. Strain-resolved tensile Raman spectroscopy was employed to quantify molecular-scale stress redistribution, with spectra acquired continuously during deformation. Raman peak positions, centroid shifts, and FWHM values were extracted to analyze strain-dependent interfacial dynamics. Experimental details are provided in Supplementary Notes 4 and 5.

Machine learning workflow
A machine-learning-guided optimization framework combining Pareto set learning (PSL) and active learning (AL) was used to explore the multi-objective formulation space. The initial dataset consisted of experimentally measured formulations spanning controlled ranges of composition and processing parameters. Gaussian process regression models were trained to predict mechanical properties and associated uncertainties, which were subsequently used to guide iterative candidate selection. Equal experimental budgets were enforced for benchmarking against response surface methods. Full algorithmic details and code implementation are available in Supplementary Information and the associated Code Ocean capsule.

Simulation methods
All-atom molecular dynamics simulations were performed to investigate interfacial chain entanglement and deformation mechanisms, while finite element simulations were used to analyze stress localization and impact-induced damage evolution in trabecular architectures. Simulation methodologies and parameters are detailed in Supplementary Notes 2 and 3.

Data availability
The data supporting the findings of this study, including mechanical performance data and datasets used for machine-learning model training and evaluation, are provided in the Supplementary Information and Source Data files associated with this article. All data are available without restriction. The machine-learning codes and workflows used in this study are publicly available via Code Ocean at https://doi.org/10.24433/CO.5966617.v1. All data are available from the corresponding author upon request. Source data are provided with this paper.

Code availability
The complete source code and scripts used for Gaussian Process (GP) modeling, Pareto Set Learning (PSL), and Active Learning (AL) are deposited on Code Ocean under DOI: 10.24433/CO.5966617.v1.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The core value of this paper lies in proposing a composite strategy combining “bone-inspired trabecular interlocking structures + stress-adaptive interfaces,” and integrating it with a machine learning framework to achieve multi-objective optimization. This approach demonstrates innovation and application potential in synergistically enhancing the strength, toughness, and impact resistance of polymer composite materials. Experimental results have preliminarily validated the superior performance of the optimized composite materials, and they exhibit a certain degree of universality across various polymer systems and reinforcing agents, thereby possessing foundational scientific research value.\nCurrently, the manuscript suffers from insufficient evidence (e.g., inadequate explanation of the dynamic mechanisms of the stress-adaptive interface, need to strengthen the sample size and reliability of the machine learning model, and need to supplement experimental replication verification), and it is recommended for major revision before acceptance.\n1. “Cross-validation demonstrates high predictive accuracy (R2 > 0.82)” The evaluation model should not be limited to R2. Please provide other data such as MSE and MAE.\n2. Please provide the relevant code in the supporting documents.\n3. In this section, the author prepared and tested three representative balanced POK-TPU-BF TIC samples, each optimized for different design objectives within the Pareto set: maximum strength, maximum fracture toughness, and maximum impact energy dissipation. Did the author use ML optimization to synthesize materials with the best overall strength, fracture toughness, and impact energy dissipation, since some applications require all three properties to a certain extent?\n4. Where did the initial data used for ML training come from? Was it the initial 30 experimental data sets? It seems that the sample size is small, which may lead to model overfitting. Please provide the specific parameters for ML model training.\n5. Active Learning was only iterated four times (the number of new samples added was not specified), and a “supervolume improvement rate of 0.1%” is insufficient to prove convergence. Additional learning curve analysis (such as changes in model error as the number of samples increases) is required. It is recommended to add independent test set validation (e.g., randomly selecting 10% of formulas not used in training for experimental validation) and compare the optimization efficiency with traditional experimental design methods (e.g., response surface method).\n6. The optimized high-performance region (strength 205-225 MPa, etc.) deviates from the experimentally verified 259 MPa, and the source of the error (such as model prediction error or experimental repeatability issues) has not been explained.\n7. Why are there no significant differences in the strain values of the stress-strain curves in Figure 4a? Most of this degree of strain difference is due to error.\n8. Figure 4d shows the preference vectors λ=[0.032:0.073:0.896], λ=[0.028:0.081:0.891], and λ=[0.034:0.089:0.887], which correspond to the optimal impact energy dissipation, fracture toughness, and strength of TIC, respectively. As is well known, TPU has excellent fracture toughness, BF has excellent strength, and POK has excellent impact energy dissipation. The mechanical test results for λ=[0.034:0.089:0.887] and λ=[0.032:0.073:0.896] are normal and understandable. However, the mechanical results for λ=[0.028:0.081:0.891] appear to be anomalous. Please explain why, among these three materials, the fracture toughness is optimal when the TPU content is reduced.\n9. The images in Figures 4g, h, and i are unclear.\n10. The SAXS and WAXS tests in Figures 5b, c, d, and e can only demonstrate that the crystallinity and crystalline morphology of the polymer change under tensile stress. In fact, most polymer materials undergo changes in SAXS and WAXS even with minor tensile stress. These tests are unrelated to the innovative aspect of the manuscript—biomimetic beams. Additionally, phenomena such as “fiber bridging” and “crack deflection” were only observed via SEM, without quantitative characterization of the dynamic process of interfacial stress response (e.g., interfacial shear strength and changes in chain entanglement density under different strains). Molecular dynamics simulations only depicted the static state of chain entanglement and did not simulate the dynamic evolution under stress loading (e.g., the rate of entanglement breakage, or quantitative calculations of energy dissipation). It is recommended to supplement in situ mechanical testing with spectroscopic characterization (e.g., Raman spectroscopy to track changes in the interfacial chemical environment) or to quantify the strain dependence of interfacial adhesion forces using atomic force microscopy.", "reviewer_2": "In this study, a sequential multi-objective experimental design algorithm was applied to explore the optimal composition of a three-component polymer composite material that forms a bone trabecular interconnection structure, aiming to improve and balance three mechanical properties. As a result, a polymer composite with properties equivalent to or exceeding those of conventional engineering composites and biomaterials was successfully produced. Through comprehensive structural analysis, mechanical testing, and MD simulations, the mechanisms that enable the fabricated material to exhibit high strength, toughness, and impact resistance were elucidated.\nIn recent years, there has been increasing attention to research on the design and development of polymer materials utilizing machine learning, but studies focusing on the mechanical properties of polymer composites relatively few. In this context, the proof of concept demonstrated in this study is significant.\nHowever, the material system under investigation (bone trabecular interconnection structure) is somewhat limited, and the mechanical performance achieved is not yet a drastic improvement over conventional materials. Therefore, it would likely attract more attention if published in a specialized journal in polymer chemistry or materials science, rather than a general journal like Nature Communications, which targets a broader audience.\nThe methodology and the path to the conclusions are deemed rational. I would like to offer the following minor comments mainly regarding the machine learning aspects of the study. I hope these can contribute to improving the paper.\n(1) The input variables for the machine learning models are described as consisting of the composition formula and process parameters. However, any explanation of the process variables is not provided.\n(2) It is stated that the model input includes cubic terms of the composition formula. Is this description correct? If so, please explain why the cubic terms are included without the quadratic terms. Additionally, does the model include cubic interaction terms? A more precise description is required.\n(3) It is mentioned that feature standardization and outlier removal were performed, but the specific procedures are not explained.\n(4) In Fig. S7, the predictive accuracy of the machine learning models is shown. However, instead of using the hold-out method, the author should report the mean and variance of the prediction performance after performing a different split of the training and testing sets.\n(5) In multi-objective optimization, how are the three material properties with different scales adjusted? Without scale transformation, the optimization is generally biased towards the material properties with larger scales.\n(6) Eq S12: Please change the first letter of \"Min\" to lowercase.\n(7) In Table S3, a performance comparison with existing materials, such as engineering composites, is shown, but it is necessary to describe specifically what materials are being referred to.", "reviewer_3": "They propose a universal toughening strategy that integrates a bone-inspired trabecular interlock architecture with a thermodynamically driven, stress-adaptive interface to enable efficient energy dissipation under mechanical loading. However, this article fails to demonstrate the necessity of employing machine learning in its approach. Moreover, it does not adequately clarify the chemical or physical mechanisms underlying its core material design methodology, nor does it provide a clear explanation of the principles governing the dynamic behavior of the proposed system. Therefore, I recommend rejection.\nSpecific Comments:\n1. The chemical or physical mechanism behind the proposed trabecular interlock architecture for interface toughening is not clearly explained.\n2. It is unclear in what way the “adaptive interfaces” are self-adaptive. How are they fundamentally different from traditional chemically constrained interfaces? Please elaborate on the underlying adaptive mechanism.\n3. The correlation between the chemical composition/molecular structure of the adaptive interfaces and their interfacial interactions/dynamics needs to be explicitly established.\n4. The necessity of using machine learning (ML) is not sufficiently justified. Typically, theoretical calculations provide target values for strength and toughness, which are often not met in practice. Please provide direct evidence that ML-based screening of molecular structures, material systems, and processing techniques can effectively bridge this gap and achieve properties closer to the theoretical predictions.\n5. Please provide a detailed description of how the molecular structure design and surface micro/nano-structures of the adaptive interfaces are translated into machine-readable formats (e.g., descriptors, feature vectors) for the ML model.\n6. The sample size used for the machine learning training must be sufficient. Please ensure the dataset adequately considers the variations from different molecular structure designs and micro/nano-structural parameters."}, "author_response": "Point-to-point responses to comments from Reviewers\nReviewer #1: The core value of this paper lies in proposing a composite strategy\ncombining “bone-inspired trabecular interlocking structures + stress-adaptive\ninterfaces,” and integrating it with a machine learning framework to achieve multi-\nobjective optimization. This approach demonstrates innovation and application\npotential in synergistically enhancing the strength, toughness, and impact resistance of\npolymer composite materials. Experimental results have preliminarily validated the\nsuperior performance of the optimized composite materials, and they exhibit a certain\ndegree of universality across various polymer systems and reinforcing agents, thereby\npossessing foundational scientific research value.\nCurrently, the manuscript suffers from insufficient evidence (e.g., inadequate\nexplanation of the dynamic mechanisms of the stress-adaptive interface, need to\nstrengthen the sample size and reliability of the machine learning model, and need to\nsupplement experimental replication verification), and it is recommended for major\nrevision before acceptance.\nComment 1: “Cross-validation demonstrates high predictive accuracy (R2 > 0.82)”\nThe evaluation model should not be limited to R2. Please provide other data such as\nMSE and MAE.\nResponse: We thank the reviewer for this valuable suggestion. In the revised\nmanuscript, we have expanded the cross-validation analysis to include additional error-\nbased metrics. Specifically, alongside the coefficient of determination (R2), we now\n2\nexplicitly report the mean absolute error (MAE) and the mean squared error (MSE) for\nstrength, fracture toughness, and impact energy dissipation. As shown in Fig. S7,\nlearning curves are presented for all three metrics, demonstrating monotonic\nimprovement in predictive accuracy with increasing training set size. The MAE and\nMSE values decrease systematically, while the associated uncertainty bands across\nrepeated cross-validation folds shrink, indicating enhanced generalization and reduced\nsensitivity to data partitioning.\nImportantly, the magnitudes of MAE and MSE remain small relative to the\nexperimental ranges of each property, confirming the absence of outlier-dominated\nprediction errors and supporting the physical reliability of the Gaussian process\nregression models. These additions provide a more complete and quantitative\nassessment of model performance, as requested by the reviewer.\nAccordingly, the corresponding text in main text has been revised as follows:\n“An initial baseline dataset is established using a design-of-experiments (DoE)\nstrategy by selecting 50 representative composite formulations spanning the feasible\ncomposition space. Each formulation is experimentally fabricated and mechanically\ncharacterized to measure tensile strength, fracture toughness, and impact energy\ndissipation. These data are used to train three independent Gaussian process regression\n(GPR) surrogate models.\nAs the training set size increases from 30 to 50, cross-validated predictive\nperformance improves systematically for all three targets, as evidenced by increasing\ncoefficients of determination (R2), decreasing mean absolute errors (MAE), and\n3\ndecreasing mean squared errors (MSE), together with a concurrent reduction in\npredictive uncertainty quantified by the standard deviation across validation folds\n(Supplementary Note 3.1 and Fig. S7). Although toughness-related properties exhibit\ncomparatively higher uncertainty at smaller sample sizes, all surrogate models attain\nstatistically reliable accuracy at 50 samples, thereby establishing a robust foundation\nfor subsequent multi-objective optimization.\nBuilding on the validated GPR surrogates, multi-objective optimization is\nperformed using Pareto set learning (PSL). PSL learns a mapping from preference\nvectors λ, which encode user-defined trade-offs among competing objectives, to\noptimized formulations x by scalarizing the multi-objective problem into a tractable\nsingle-objective form [29]. By sampling 1000 preference vectors, the trained PSL\nmodel generates a diverse and structured set of Pareto-consistent candidate solutions\nspanning a broad range of engineering priorities (Supplementary Note 3.2 and Fig.\nS8).”\nIn addition, the corresponding text in the Supplementary Information has been\nrevised as follows:\n“Note 3. Machine learning (ML) driven multi-objective optimization\nTo efficiently explore the high-dimensional formulation space while minimizing\nexperimental burden, we integrate machine learning (ML) with a structured\nexperimental design strategy. An initial design-of-experiments (DoE) protocol is used\nto select 50 representative composite formulations that uniformly span the admissible\ncomposition and processing space. Each formulation corresponds to a unique\n4\ncombination of polymer constituents and processing parameters.\nAll selected samples are fabricated and experimentally characterized to obtain\nthree target mechanical properties—tensile strength, fracture toughness, and impact\nenergy dissipation—which together define the multi-objective performance landscape\nof the composite system. This experimentally grounded dataset serves as the foundation\nfor surrogate modeling, preference-guided optimization, and active learning–assisted\nexploration described below.\n3.1 Gaussian Process Regression (GPR) modeling\nGaussian process regression (GPR) [4,5] is employed to construct probabilistic\nsurrogate models that map the composite formulation space to the corresponding\nmechanical responses, including strength, fracture toughness, and impact energy\ndissipation. GPR is particularly well suited for this task because it provides flexible,\nnon-parametric function approximation together with calibrated predictive uncertainty.\nThe ability of GPR to quantify both the predictive mean and uncertainty is\nessential in the present context, where the available experimental data are necessarily\nlimited by fabrication and testing costs. Beyond accurate interpolation, the uncertainty\nestimates enable risk-aware decision making and form the basis for subsequent active\nlearning and multi-objective optimization, allowing the model to guide experimental\neffort toward regions of high performance potential and high epistemic uncertainty.\nInput features and target variables. The input features are defined by the\nexperimentally controlled formulation parameters, specifically the weight fractions of\nthermoplastic polyurethane (TPU) and basalt fiber (BF). The content of\n5\npoly(octenylene ketone) (POK) is implicitly determined by mass conservation and is\ntherefore not treated as an independent variable. This formulation ensures that the input\nspace reflects physically admissible composite compositions without introducing\nredundant degrees of freedom. Three mechanical properties are modeled as\nindependent targets: tensile strength, fracture toughness, and impact energy dissipation.\nPrior to model training, all input features are standardized using z-score normalization\nto improve numerical conditioning and to prevent scale imbalance during kernel\noptimization. Each target property is modeled using an independent Gaussian process\nregression (GPR) surrogate. This decoupled modeling strategy avoids imposing\nartificial correlations between objectives and preserves the interpretability of individual\nstructure–property relationships, which is essential for subsequent preference-guided\noptimization and physical analysis of trade-offs among competing mechanical\nresponses.\nKernel design and hyperparameter optimization. Each GPR model employs a\ncomposite kernel of the form:\n𝑘(𝑥,𝑥′) = 𝜎2𝑘 (𝑥,𝑥′)+𝑘 (𝑥,𝑥′) (S9)\n𝑀𝑎𝑡𝑒𝑟𝑛 𝑛𝑜𝑖𝑠𝑒\nwhere 𝑘 (𝑥,𝑥′) is a Matérn kernel with smoothness parameter 𝜈 ∈\n𝑀𝑎𝑡𝑒𝑟𝑛\n{0.5,1.5,2.5} , and 𝑘 (𝑥,𝑥′) is a white-noise kernel accounting for experimental\n𝑛𝑜𝑖𝑠𝑒\nuncertainty.\nThe Matérn kernel is chosen because it allows controlled deviation from overly\nsmooth response surfaces and is therefore well suited to capturing the non-smooth,\nheterogeneous structure–property relationships characteristic of polymer-based\n6\ncomposite systems. The smoothness parameter 𝜈 is selected from the set {0.5,1.5,2.5},\nenabling the model to flexibly span regimes ranging from rough, locally varying\nresponses to moderately smooth trends without enforcing unrealistic differentiability\nassumptions.\nInitial values of the characteristic length scales and noise variance are drawn from\na predefined grid to ensure robust coverage of plausible hyperparameter regimes. The\nfinal kernel hyperparameters are obtained by maximizing the log marginal likelihood\nusing gradient-based optimization with multiple random restarts, mitigating sensitivity\nto local optima and ensuring stable surrogate construction across different target\nproperties.\nOutlier handling and prevention of information leakage. To enhance surrogate\nrobustness while rigorously preventing information leakage, a two-stage outlier\nhandling strategy is adopted. During cross-validation and learning-curve analysis, a\njoint three-sigma criterion is applied exclusively within each training fold, using only\nstatistics computed from the corresponding training data. Validation samples are not\naltered or screened at any stage of this process, ensuring strictly leakage-free\nperformance evaluation. For final model training and deployment, a global joint three-\nsigma filtering across all target properties is applied to the full experimental dataset.\nThis step is performed only after model assessment is complete and serves to construct\na clean, physically consistent training set for the final surrogate models used in\noptimization and active learning.\nModel evaluation and learning-curve analysis. Model performance is evaluated\n7\nusing repeated five-fold cross-validation to ensure robust estimation under limited data\nconditions. Predictive accuracy is quantified using three complementary metrics: the\ncoefficient of determination (R2), mean absolute error (MAE), and mean squared error\n(MSE). While R2 captures the strength of correlation between predicted and\nexperimental values, MAE provides a physically interpretable measure of average\nabsolute deviation, and MSE offers increased sensitivity to occasional large prediction\nerrors. Learning-curve analysis is performed by training the GPR models on\nprogressively larger subsets of the experimental dataset (N = 30–50). For each subset\nsize, cross-validation is repeated multiple times with fixed kernel hyperparameters, and\nthe mean and standard deviation of R2, MAE, and MSE across validation folds are\nreported. When MSE is not directly available from cross-validation outputs, it is\ncomputed as the square of the corresponding root-mean-squared error (RMSE). The\njoint use of these metrics enables a comprehensive assessment of model accuracy,\nrobustness, and susceptibility to outlier-driven errors, providing a reliable basis for\nsubsequent uncertainty-aware optimization.\nFinal model training and deployment. After hyperparameter selection, final\nGPR models are trained on the cleaned dataset using the optimized kernel parameters.\nThe trained models are stored together with the optimized kernel expressions and log\nmarginal likelihood values, enabling full reproducibility and transparent inspection of\nmodel behavior. These surrogate models are subsequently integrated into the multi-\nobjective optimization framework described below.\nResults. Using the initial set of 50 experimental samples, the Gaussian process\n8\nregression (GPR) surrogate models exhibit systematic improvement in predictive\nperformance as the effective training size increases from 30 to 50 (Fig. S7). For all three\ntarget properties, namely strength, fracture toughness, and impact energy dissipation,\nthe cross-validated coefficient of determination (R2) increases monotonically,\naccompanied by concurrent decreases in both the mean absolute error (MAE) and mean\nsquared error (MSE). These consistent trends across complementary performance\nmetrics indicate progressively enhanced generalization accuracy rather than metric-\nspecific artifacts.\nIn parallel, the uncertainty bands (±1 standard deviation across repeated cross-\nvalidation splits) contract steadily with increasing training size, reflecting reduced\nsensitivity to data partitioning and stabilization of the posterior predictions. Toughness-\nrelated properties, particularly impact energy dissipation, exhibit comparatively larger\nuncertainty at smaller sample sizes, consistent with their higher sensitivity to\nmicrostructural heterogeneity and experimental variability. Nevertheless, by 45–50\nsamples, all three models enter a statistically reliable accuracy regime characterized by\nlow absolute errors, reduced variance, and diminishing marginal performance gains.\nTogether, these results demonstrate that the initial DoE dataset provides a robust\nand physically meaningful surrogate foundation for subsequent uncertainty-aware\nmulti-objective optimization, while also quantitatively motivating the use of active\nlearning to further refine model fidelity in performance-critical regions of the design\nspace.\nc. Learning curves for predicting strength, fracture toughness, and impact energy dissipation,\nrespectively, as functions of the training set size (N = 30, 35, 40, 45, and 50). For each target,\nthree complementary cross-validation metrics are reported: the coefficient of determination (R2, top row),\nmean absolute error (MAE, middle row), and mean squared error (MSE, bottom row). Solid lines\nrepresent the mean performance over repeated five-fold cross-validation, and shaded regions\nindicate ±1 standard deviation across validation splits. The monotonic increase in R2 together with\nthe concurrent decrease in MAE and MSE, as well as the progressive narrowing of the uncertainty\nbands, demonstrate systematic improvement in generalization accuracy and stabilization of the\nsurrogate models with increasing training data.”\nComment 2: Please provide the relevant code in the supporting documents.\nResponse: All codes relevant to the machine learning and optimization framework\nhave been made publicly available via a dedicated Code Ocean computational capsule.\nThis approach ensures full computational reproducibility in a controlled execution\nenvironment.\nComment 3: In this section, the author prepared and tested three representative\n10\nbalanced POK-TPU-BF TIC samples, each optimized for different design objectives\nwithin the Pareto set: maximum strength, maximum fracture toughness, and maximum\nimpact energy dissipation. Did the author use ML optimization to synthesize materials\nwith the best overall strength, fracture toughness, and impact energy dissipation, since\nsome applications require all three properties to a certain extent?\nResponse: We thank the Reviewer for raising this practical and important question.\nYes. While this section initially presents three preference-resolved representative\nformulations (strength-first, fracture-toughness–first, and impact-energy–first) to\ndemonstrate that the learned Pareto set can be queried according to different\nengineering priorities, we also used the same AL–PSL framework to extract and\nexperimentally validate an all-round (“best overall”) balanced formulation for\napplications that require all three properties simultaneously.\nImportantly, for a multi-objective system governed by Pareto optimality, there is\ngenerally no single composition that globally maximizes all three objectives at once.\nInstead, the appropriate “best overall” solution corresponds to a Pareto-optimal\ncompromise (equivalently, a formulation that satisfies engineering constraints across\nall three metrics while remaining non-dominated). Therefore, our goal was not to claim\na global triple-maximum, but to deliver a balanced Pareto-optimal design that jointly\nmaintains high values in all three properties.\nAlgorithmic extraction (no subjective choice/no cherry-picking). After completing\nfive active-learning iterations (25 additional experiments, total n = 75), we retrained the\nPSL model on the expanded dataset and generated the updated Pareto-consistent\n11\nsolution set (Fig. 3d–f). We then performed preference-guided extraction by prescribing\nan equal-weight balanced preference vector (λ = [1/3, 1/3, 1/3]) and selecting the\nformulation whose realized preference vector is closest to the target, as defined in\nSupplementary Note 3.5 and summarized in Table S2. This procedure is fully\nalgorithmic (distance minimization in preference space) and avoids any post hoc\nmanual selection.\nExperimentally validated balanced (“best overall”) formulation. The extracted\nbalanced formulation in Table S2 (TPU 3.82 wt%, BF 10.33 wt%, POK 85.85 wt%)\nachieves simultaneously high performance: strength 249.71 MPa, fracture toughness\n14.83 MPa·m1/2, and impact energy dissipation 4.68 J. In addition, the experimentally\nvalidated preference-guided set is tightly clustered in absolute performance (Fig. 4d–\nf), indicating that the AL–PSL search converges to a high-performance basin/plateau\nwhere all three metrics are jointly high, and preference guidance provides fine-grained\nprioritization within this optimized regime rather than trading off one property severely\nto boost another.\nAccordingly, the revised manuscript now explicitly includes and experimentally\nvalidates a balanced, best-overall Pareto-optimal formulation, demonstrating that the\nproposed AL–PSL framework supports both preference-resolved designs and all-round\nmaterial solutions for practical applications requiring concurrent mechanical\nperformance.\nPreference-guided solution extraction\nAfter completion of five active learning iterations (25 additional experiments),\nthe\nPSL model is retrained using the expanded dataset, yielding the optimized Pareto set\nshown in Fig. 3d–f. The resulting Pareto-consistent solutions exhibit a pronounced and\norganized structure across both composition and performance spaces. In the TPU–BF–\nPOK composition space (Fig. 3d), the solutions are confined to a relatively narrow\ncompositional window rather than being dispersed across the full admissible domain,\nrevealing a well-defined high-performance basin.\nA consistent structural organization is observed in the three-dimensional\nperformance space (Fig. 3e), where the Pareto-consistent solutions form a smooth,\nmanifold-like trajectory instead of a diffuse cloud. This behavior indicates that the\nachievable multi-objective trade-offs are constrained by a common structural regime,\nwithin which improvements in tensile strength are accompanied by predictable and\ncoupled variations in fracture toughness and impact energy dissipation, rather than\narbitrary performance fluctuations.\nThis structured Pareto manifold enables preference-guided engineering decision\nmaking (Fig. 3f and Table S2). By specifying target preference vectors, representative\nsolutions corresponding to balanced, strength-prioritized, fracture-toughness–\nprioritized, and impact-energy–prioritized design modes can be directly extracted along\nthe learned manifold. The close agreement between the prescribed preference vectors\nand their matched counterparts confirms that PSL preserves preference information\n13\nduring solution generation, allowing explicit design priorities to be translated into\nconcrete compositions and predicted performance outcomes.\nThese results demonstrate that the combined AL–PSL framework not only\naccelerates the discovery of high-performance Pareto fronts but also constructs a\nmechanistically admissible and queryable design manifold, enabling practical and\ntransparent composite formulation design under competing performance objectives.\nFigure 3 d–f. Pareto set obtained by retraining the PSL model after completion of five AL iterations\n(25 additional experiments, total n=75). d. Approximate Pareto set in the TPU–BF–POK\ncomposition space. Green points denote PSL-generated Pareto-consistent solutions; semi-\ntransparent connections indicate preference-guided associations between neighboring solutions\n(shown as a guide to the eye). The solutions are confined to a relatively narrow compositional basin.\ne. Three-dimensional representation of the Pareto-consistent solutions in the strength–fracture\ntoughness–impact energy dissipation space, revealing a smooth, manifold-like trade-off structure. f.\nProjection of the Pareto set onto the strength–fracture toughness plane with impact energy\ndissipation encoded by color. Highlighted markers indicate representative solutions extracted for\n14\nbalanced, strength-prioritized, fracture-toughness–prioritized, and impact-energy–prioritized design\npreferences. The inset ternary diagram compares prescribed preference vectors (stars) with the\ncorresponding matched preference vectors (circles), demonstrating close alignment between\nintended design priorities and realized PSL solutions.\nExperimental validation of preference-guided Pareto solutions\nTo validate the machine-learning–selected Pareto solutions experimentally, we\nfabricated four representative formulations corresponding to distinct design preferences\n(balanced, strength-first, fracture-first, and impact-first). Across tensile, notched tensile,\nand impact loading, the macroscopic response curves remain highly consistent (Fig.\n4a–c), with comparable yielding, strain hardening, and terminal failure behavior. This\nsimilarity indicates that preference-guided selection operates within a narrowly defined\nhigh-performance basin, rather than shifting the material into qualitatively different\ndeformation regimes.\nQuantitatively, the absolute properties cluster tightly. Tensile strength spans\napproximately 248–250 MPa (maximum deviation ~1.2 MPa relative to the balanced\nformulation, <0.5%, Fig. 4d and 4g), and fracture toughness remains near 14.8–14.9\nMPa·m1/2 (differences ~0.02 MPa·m1/2, Fig. 4e and 4h), both within experimental\nscatter. Impact energy dissipation is similarly concentrated (~4.65–4.71 J, Fig. 4f and\n4i), with the impact-first design exhibiting a small but reproducible increase (~0.06 J,\n~1%) relative to the balanced reference. Expressed as relative differences with respect\nto the balanced formulation, the shifts remain modest yet directionally consistent with\nthe prescribed preferences, supporting that preference guidance primarily refines trade-\noff ordering within an already optimized performance envelope.\nPareto set. a. Tensile stress–strain curves of four representative designs (balanced, strength-first,\nfracture-first, and impact-first; compositions indicated in the legend). b. Force–displacement\nresponses from notched tensile (fracture) tests, with the inset schematically illustrating the loading\nconfiguration. c. Impact force–displacement curves showing the energy-dissipation process. d–f.\nSummary of experimentally measured tensile strength, fracture toughness, and impact energy\ndissipation (mean ± s.d., n = 5 independent specimens per design). The shaded region denotes the\n±1 s.d. envelope of the balanced formulation. g–i. Relative differences (Δ) of the three mechanical\nmetrics with respect to the balanced design. Dashed lines indicate zero difference. The small Δ\nvalues compared with experimental scatter highlight that preference selection reorders trade-offs\nwithin a tightly clustered high-performance region rather than inducing large absolute shifts in\nmechanical performance..\nSupplementary Notes 3.5 Preference-guided solution extraction\nAfter completion of the AL process, the trained PSL model enables direct\nextraction of representative material formulations corresponding to prescribed\npreference vectors λ. These preference vectors encode distinct engineering design\npriorities, including balanced performance, strength-first, fracture-first, and impact-\nfirst designs objectives.\nminimizing the distance between the prescribed preference vector and the realized\npreference vector associated with the predicted performance of candidate solutions.\nThis selection criterion ensures consistency between the intended design priorities and\nthe predicted multi-objective performance, while avoiding manual, subjective, or post\nhoc solution selection.\nThe resulting preference-guided formulations provide a compact and interpretable\nsubset of the learned Pareto set, enabling direct translation of the multi-objective\noptimization results into engineering-relevant material choices. These representative\nsolutions and their corresponding predicted and experimentally validated properties are\nsummarized in Table S2, illustrating how the proposed framework supports\nengineering-driven decision making across different performance trade-offs.”\nTable S2 Engineering-driven material selection via preference-guided Pareto set learning.\ntarget_ target_ target λ_ λ_ TP PO Stren Fract Imp\nDesign Target λ_I BF\nλ_S λ_F _λ_I S F U K gth ure act\nBalanced (equal 0.3 0.3 0.3 3. 10. 85. 249. 14.8\n0.333 0.333 0.333 4.68\npreference) 31 39 31 82 33 85 71 3\n0.8 0.0 0.0 3. 10. 85. 250. 14.8\nStrength-first 0.85 0.075 0.075 4.65\n57 86 56 78 39 83 0 2\n0.0 0.8 0.0 3. 10. 85. 249. 14.8\nFracture-first 0.075 0.85 0.075 4.68\n71 51 78 82 36 82 86 3\n0.0 0.0 0.8 3. 10. 85. 248. 14.8\nImpact-first 0.075 0.075 0.85 4.71\n64 71 66 93 22 85 52 1\nComment 4: Where did the initial data used for ML training come from? Was it the\ninitial 30 experimental data sets? It seems that the sample size is small, which may lead\n17\nto model overfitting. Please provide the specific parameters for ML model training.\nResponse: We thank the Reviewer for this important question regarding data\nprovenance, sample size, and overfitting control. In the revised manuscript, the ML\ntraining dataset and training parameters have been clarified and updated to reflect the\nfinal workflow. Please also refer to our response to Comment 1.\nSource of the initial ML training data and sample size. The initial dataset used to\ntrain the surrogate models is not limited to 30 samples in the revised version. We\nestablish the baseline dataset using a design-of-experiments strategy with n = 50\nexperimentally fabricated and mechanically tested formulations spanning the\nadmissible TPU–BF–POK composition space. These 50 measurements constitute the\ninitial training data for the three Gaussian process regression (GPR) surrogates\n(strength, fracture toughness, and impact energy dissipation). The previously\nmentioned 30-sample setting corresponds to an earlier stage/version and is now used\nonly for learning-curve diagnostics, where we explicitly evaluate how model\nperformance evolves as the training size increases from 30 to 50.\nTo reduce overfitting risk under limited data and to quantify robustness, we use\nmultiple complementary safeguards and evaluations:\n1. Repeated five-fold cross-validation with mean and variance reporting.\nFigure S7 reports the mean and standard deviation of predictive performance across\nrepeated five-fold CV splits, rather than a single hold-out result, addressing split\nsensitivity directly.\n2. Learning-curve analysis. We evaluate performance systematically from N =\n18\n30 to 50 (Fig. S7) and also report the learning behavior across the full AL-updated\ndataset (Fig. S10). The monotonic improvement with increasing N supports improving\ngeneralization rather than overfitting.\n3. Independent hold-out validation. In addition to CV, we include an\nindependent test evaluation (Fig. S11) to verify out-of-sample behavior and check for\nsystematic bias.\n4. Model regularization through probabilistic training. GPR hyperparameters\nare optimized by maximizing the log marginal likelihood with multiple random restarts,\nwhich penalizes overly complex fits in a Bayesian manner.\nML model training parameters. Each target property (strength, fracture toughness,\nand impact energy dissipation) is modeled independently using Gaussian process\nregression. Input features consist of the experimentally controlled composition\nvariables (TPU and BF weight fractions), with POK implicitly constrained by mass\nconservation. All inputs are standardized using z-score normalization. The GPR models\nemploy a Matérn kernel with smoothness parameter ν ∈ {0.5, 1.5, 2.5}, combined\nwith a white-noise kernel to account for experimental uncertainty. Kernel\nhyperparameters, including characteristic length scales and noise levels, are optimized\nby maximizing the log marginal likelihood using gradient-based optimization with\nmultiple random restarts. Leakage-free three-sigma filtering is applied only to training\nfolds during cross-validation, and a global three-sigma criterion is used for final model\ndeployment. Full details are provided in Supplementary Note 3.1.\nThe revisions made in the manuscript are as follows: “\nResults. Please refer to our response to Comment 1.\nSupplementary Note3.3 Active Learning (AL)\nThe evolution of surrogate model performance over five sequential AL iterations\nis summarized by the learning-curve analysis shown in Fig. S10. Within each iteration,\npredictive accuracy improves systematically as the training set size increases, as\nevidenced by monotonic increases in cross-validated R2 and corresponding decreases\nin mean absolute error (MAE) for strength, fracture toughness, and impact energy\ndissipation. This consistent intra-iteration behavior indicates that the surrogate models\nbenefit directly from additional training data and do not exhibit signs of overfitting. In\nparallel, the uncertainty associated with cross-validation progressively decreases with\nincreasing training set size, reflecting reduced sensitivity to data partitioning and\nenhanced model stability. These trends demonstrate well-behaved learning\ncharacteristics of the Gaussian process surrogates within each AL round.\nBeyond these intra-iteration effects, a clear inter-iteration improvement is\nobserved across successive AL rounds. For a given training set size, the baseline\npredictive performance shifts upward from iteration 1 to iteration 5, accompanied by\nan overall reduction in MAE and tighter cross-validation dispersion. This cumulative\nimprovement is particularly pronounced for fracture toughness and impact energy\ndissipation, which display larger variability and broader uncertainty at early stages but\nbecome substantially more stable and accurate in later iterations. In contrast, strength\nprediction remains comparatively robust even in early iterations, consistent with a\n20\nsmoother and less complex response surface with respect to the available compositional\ndescriptors.\nAt larger training set sizes, the learning curves for all three objectives exhibit\nsaturation behavior, characterized by diminishing marginal gains in both R2 and MAE.\nThis plateau indicates that, within the explored design space, the surrogate models\napproach a near-converged performance regime. These results show that the AL\nstrategy not only improves average predictive accuracy but also accelerates the\nstabilization of surrogate models, particularly for objectives with higher intrinsic\nvariability, thereby enabling more efficient allocation of experimental effort toward\nregions of greatest uncertainty.\nFigure S10 Learning curves over five sequential active learning iterations. a–e.\nLearning curves corresponding to active learning iterations 1–5, respectively. Within each iteration,\npredictive performance is evaluated as a function of training set size using cross-validation. Upper\npanels report the coefficient of determination (R2), and lower panels report the mean absolute error\n21\n(MAE), for strength, fracture toughness, and impact energy dissipation. Solid lines denote the mean\ncross-validation performance, while shaded regions represent variability across repeated cross-\nvalidation splits. Increasing training set size consistently improves predictive accuracy and reduces\nuncertainty. Comparison across panels (a–e) shows cumulative improvements in baseline\nperformance and progressive stabilization of learning behavior, with particularly pronounced gains\nfor fracture toughness and impact energy dissipation. At larger training sizes, all objectives exhibit\nsaturation behavior, indicating diminishing marginal gains within the explored design space.\nComplementing the learning-curve analysis, we further assess the final surrogate\nmodels using an independent hold-out test set (Fig. S11). This test set is strictly\nexcluded from both model training and final refitting and is used solely for out-of-\nsample performance evaluation. Across all three target mechanical properties, the\npredicted values show close agreement with experimental measurements, as indicated\nby high coefficients of determination and low prediction errors.\nAmong the three targets, tensile strength is predicted with the highest fidelity (R2\n= 0.959, MAE = 1.364 MPa), reflecting a comparatively smooth and well-structured\nresponse surface with respect to the selected compositional descriptors within the\nexplored domain. Fracture toughness also exhibits excellent predictive accuracy (R2 =\n0.944, MAE = 0.153 MPa·m1/2), demonstrating that the surrogate model reliably\ncaptures variations in toughness despite the increased complexity of this property.\nImpact energy dissipation remains the most challenging target; nevertheless, it is\npredicted with reasonable accuracy (R2 = 0.758, MAE = 0.184 J), consistent with its\nhigher intrinsic variability under the present experimental conditions.\nNo systematic deviation from the parity line is observed across the independent\ntest set, indicating the absence of pronounced prediction bias and supporting robust out-\nof-sample generalization within the explored composition space. Taken together with\nthe progressive improvements observed in the learning-curve analysis, these results\n22\ndemonstrate that the surrogate models achieve reliable and transferable predictive\nperformance, providing a sound foundation for subsequent multi-objective\noptimization and data-driven materials design.\nFigure S11 Independent test validation of surrogate model predictions. Parity plots comparing\npredicted and experimentally measured values for a. strength, b. fracture toughness, and c. impact\nenergy dissipation on an independent test set excluded from model training. The solid lines indicate\nperfect agreement (y = x). The coefficient of determination (R2), mean absolute error (MAE), and\nroot mean squared error (RMSE) are reported for each target. The close alignment of predictions\nwith the parity line demonstrates strong generalization performance of the surrogate models across\nall three mechanical properties.”\nComment 5: Active Learning was only iterated four times (the number of new\nsamples added was not specified), and a “supervolume improvement rate of 0.1%” is\ninsufficient to prove convergence. Additional learning curve analysis (such as changes\nin model error as the number of samples increases) is required. It is recommended to\nadd independent test set validation (e.g., randomly selecting 10% of formulas not used\nin training for experimental validation) and compare the optimization efficiency with\ntraditional experimental design methods (e.g., response surface method).\nResponse: We thank the Reviewer for this constructive comment regarding AL\niteration details, convergence evidence, learning-curve analysis, independent validation,\nand benchmarking against traditional DoE/RSM. We have substantially revised the\n23\nmanuscript and Supplementary Information to address each point with explicit\nexperimental counts and quantitative diagnostics.\nAL iterations and number of new samples. We now clearly state that AL was\nperformed for five iterations, selecting five candidates per iteration using an\nuncertainty-aware acquisition strategy (LCB + HVI), yielding 25 additional\nexperimentally validated formulations. Combined with the baseline dataset (n = 50),\nthe final dataset contains n = 75 experimentally measured formulations used for the\nfinal surrogate training and analysis (main text “Active learning (AL)-guided\noptimization”, Fig. 3 and Fig. S9).\nIn the revision, convergence is supported by multiple complementary lines of\nevidence:\nFrontier expansion behavior in property space: AL-selected formulations\nprogressively populate the upper-right high-performance boundary (strength > ~245\nMPa, fracture toughness > ~14.5 MPa·m1/2) while maintaining high impact energy\ndissipation (>4.6 J), consistent with Pareto-front expansion rather than uniform\nsampling (Fig. 3b and Fig. S9).\nSampling redistribution and refinement: the density evolution in Fig. 3b and\nFig. S9 shows an early boundary-seeking phase followed by later\ndensification/refinement of the Pareto-relevant region, indicating transition toward\npractical convergence rather than continued unconstrained exploration.\nHypervolume saturation: HV analysis shows that >90% of the total normalized\nHV gain is achieved within the first ten additional experiments, driven by a dominant\n24\nearly increase (ΔHV ≈ 0.012) and followed by rapid saturation with marginal gains\nbelow ~0.001 per iteration, demonstrating diminishing returns at later iterations (Fig.\n3c).\nLearning-curve analysis (model error vs. sample size). Following the\nReviewer’s recommendation, we added a learning-curve analysis across five AL\niterations (Fig. S10). Within each iteration, cross-validated predictive accuracy\nimproves monotonically with training-set size (increasing R2 and decreasing MAE) and\nthe cross-validation dispersion decreases, indicating progressive stabilization of the\nsurrogate models. Inter-iteration improvements are particularly evident for fracture\ntoughness and impact energy dissipation, which show larger early uncertainty but\nbecome substantially more stable in later iterations. Please refer to our response to\nComment 4.\nIndependent test-set validation. We further evaluated the final surrogate models\non an independent hold-out test set strictly excluded from training and refitting. The\nparity plots show close agreement between predicted and experimentally measured\nvalues with no systematic bias, confirming robust out-of-sample generalization for all\nthree targets (Fig. S11). Please also refer to the answer in reply to Comment 4.\nBenchmarking against RSM under identical experimental budget. We\nimplemented a controlled RSM baseline using the identical starting point (n = 50) and\nidentical additional budget (n = 25, final n = 75). Under this fair protocol, RSM-selected\nformulations remain largely in intermediate-performance regions and exhibit limited\nPareto-front expansion, whereas the proposed PSL + AL framework consistently drives\n25\nthe Pareto front toward higher-performance regimes and reaches >90% of its terminal\nnormalized HV improvement within the first ten additional experiments. In contrast,\nRSM shows only marginal improvement even after exhausting all 25 additional trials\n(Fig. S12).\nIn summary, we have clarified the AL iteration/sample counts and replaced the\ninsufficient single-threshold argument with a convergent body of evidence—property-\nspace evolution (Fig. 3 and S9), HV saturation (Fig. 3c), learning-curve stabilization\n(Fig. S10), independent test validation (Fig. S11), and a budget-matched RSM\nbenchmark (Fig. S12)—demonstrating both practical convergence and superior sample\nefficiency of the uncertainty-aware AL–PSL strategy. These additions fully address the\nconcerns raised in Comment 5.\nThe relevant text in the revised manuscript is as follows: “\nActive learning (AL)-guided optimization\nTo further improve optimization efficiency and solution quality, we integrate\nactive learning (AL) with multi-objective Bayesian optimization to iteratively refine\nthe learned Pareto set and guide exploration of the trabecular interlocked composites\n(TICs) design space [30, 31]. Specifically, we adopt an acquisition strategy that\ncombines a lower confidence bound (LCB) criterion with hypervolume improvement\n(HVI) to select small batches of informative and non-dominated candidates from the\nPSL-generated solution pool for experimental validation (Supplementary Note 3.3).\nThis uncertainty-aware strategy balances exploitation of high-performing regions with\nexploration of poorly sampled areas, enabling targeted expansion of the experimentally\n26\nvalidated performance envelope under a fixed experimental budget.\nStarting from the baseline dataset (Fig. 3a), the experimentally sampled\nformulations in the strength–fracture toughness space follow an approximately\nmonotonic trend, with impact energy dissipation co-varying with the high-performance\nregion. However, the initial data are concentrated around intermediate strength (≈220–\nperformance boundary sparsely sampled. With successive AL iterations (Fig. 3b and\nFig. S9), newly selected formulations progressively populate the upper-right region of\nthe design space, corresponding to tensile strengths exceeding ~245 MPa and fracture\ntoughness values above ~14.5 MPa·m1/2, while maintaining high impact energy\ndissipation (>4.6 J). Furthermore, the selected candidates are not clustered in regions\nof high sample density but are preferentially drawn toward areas associated with\nelevated predictive uncertainty and potential Pareto-front expansion, consistent with\nthe intended uncertainty-aware acquisition strategy.\nThe efficiency of this frontier expansion is quantified by hypervolume (HV)\nanalysis. More than 90% of the total normalized HV gain is achieved within the first\nten additional experiments, driven by a dominant early HV increase (ΔHV ≈ 0.012),\nfollowed by rapid saturation with marginal gains below 0.001 per iteration. Across five\nAL iterations, learning-curve analysis further reveals systematic improvements in\ncross-validated predictive accuracy and progressive stabilization of the surrogate\nmodels, particularly for fracture toughness and impact energy dissipation (Fig. S10).\nIndependent hold-out testing confirms robust out-of-sample generalization for all three\n27\nproperties, with no evidence of systematic prediction bias (Fig. S11).\nFinally, benchmarking against a classical response surface method (RSM) under\nidentical experimental budgets (50 initial formulations plus 25 additional trials)\ndemonstrates that the combined PSL + AL framework consistently drives the Pareto\nfront toward higher-performance regimes and achieves over 90% of its terminal\nnormalized hypervolume improvement within the first ten additional experiments. In\ncontrast, RSM yields only marginal performance gains even after exhausting the full\nexperimental budget (Fig. S12). These results establish that uncertainty-aware, iterative\nlearning enables substantially higher sample efficiency and reliability than conventional\nquadratic design-of-experiments baselines.\nFigure 3 Active learning–guided Pareto-front expansion. a. Initial experimental dataset (n=50)\n28\nprojected in the strength–fracture toughness space, with impact energy dissipation encoded by color.\nMarginal density profiles summarize the univariate distributions of strength and fracture toughness.\nb. The same property space after successive active learning (AL) iterations, highlighting previously\nselected AL formulations and the current AL batch. Newly selected points preferentially populate\nthe high-strength/high-toughness boundary while maintaining high impact energy dissipation. c.\nNormalized hypervolume (HV) of the non-dominated set as a function of the number of experiments\n(left axis), together with the incremental hypervolume gain (ΔHV) per AL update (right axis). The\nHV exhibits a pronounced early increase followed by saturation, indicating diminishing marginal\ngains after rapid Pareto-front expansion.\nSupplementary Notes 3.3 Active Learning (AL)\nThe evolution of the experimental dataset under the AL strategy in the strength–\nfracture toughness plane, with impact energy dissipation encoded by color are\nillustrated in Figs. 3a, 3b and Fig. S9. The initial experiments (Fig. 3a) define a clear\nperformance manifold in which fracture toughness and strength increase concomitantly\nover a broad range. The highest impact energy dissipation values (≈4.5–5.0 J) are\npredominantly associated with the upper-right corner of the property space. Consistent\nwith this observation, the marginal density distributions shown along the top and right\naxes indicate that the initial sampling is biased toward intermediate fracture toughness\nand strength values, whereas the high-performance extremes remain sparsely populated.\nFollowing the initiation of AL, the newly selected formulations consistently\nconcentrate on this high-performance region (Fig. 3b and S9), indicating that the\nacquisition strategy prioritizes Pareto-front expansion over uniform sampling of\nintermediate-performance regimes. In the first two iterations (Fig. S9b and S9c), the\nselected points advance toward the empirical performance envelope while\nsimultaneously extending coverage along the boundary, which is reflected not only by\nthe scatter distribution but also by a progressive shift of the marginal density profiles\n29\ntoward higher fracture toughness and strength. This redistribution of sampling density\nis essential for constructing a stable and continuous approximation of the Pareto front.\nIn subsequent iterations (Fig. S9d and S9e), sampling is primarily focused on\ndensifying the same region, consistent with a refinement phase in which additional\nexperiments enhance local resolution and confidence near the most decision-relevant\ntrade-offs. At later iterations (Fig. 3b), the finite width of the marginal density curves\nwithin the high-performance domain indicates sustained exploration of multiple non-\ndominated solutions rather than collapse to a single extremum along the Pareto-relevant\nboundary.\nFigure S9 Active learning–driven evolution of the experimental design space in multi-\nobjective composite optimization. a-d. Progressive selection of new formulations over four AL\niterations. In each panel, formulations selected in the current iteration are highlighted with distinct\nmarkers, while those selected in previous iterations are retained for reference. Across successive\niterations, the active learning strategy systematically samples the high-performance region defined\n30\nby concurrently high strength, high fracture toughness, and enhanced impact energy dissipation,\nwhile progressively refining the local resolution of the Pareto-relevant boundary. The initial\nexperimental dataset (n=50) is projected onto the strength–fracture toughness plane, with marker\ncolor encoding the experimentally measured impact energy dissipation. Shaded curves along the top\nand right margins represent kernel density estimates of the marginal distributions of fracture\ntoughness and strength, respectively. These marginal projections provide a statistical view of how\nexperimental samples populate the one-dimensional subspaces of the design space and how their\ndistributions evolve as active learning proceeds.\n3.4 Benchmarking against Response Surface Method (RSM)\nInitial dataset and benchmarking protocol. To rigorously benchmark the\noptimization efficiency of different formulation design strategies, all comparative\nanalyses are conducted under a strictly controlled and identical experimental starting\npoint. An initial dataset comprising 50 experimentally measured composite\nformulations is used as the common baseline for both the classical response surface\nmethod (RSM) and the proposed Pareto set learning combined with active learning\n(PSL + AL). Each formulation is characterized by three target mechanical properties.\nDuring the candidate selection stage, neither approach accesses any information\nbeyond this initial dataset. Both RSM and PSL + AL are required to propose exactly 25\nadditional formulations, resulting in final datasets of identical size (n=75) for direct\ncomparison. By enforcing equal initial data, identical sampling budgets, and a shared\ndesign space, this benchmarking protocol ensures that any observed differences in\noptimization outcomes arise solely from the intrinsic search and decision strategies of\nthe respective methods, rather than from disparities in data availability or experimental\nresources.\nResponse Surface Method (RSM) baseline construction. The RSM is\nimplemented as a conventional design-of-experiments (DoE) baseline using\n31\nexclusively the initial dataset (n = 50). For each target property 𝑓 ∈\n{𝑆𝑡𝑟𝑒𝑛𝑔𝑡ℎ,𝐹𝑟𝑎𝑐𝑡𝑢𝑟𝑒 𝑡𝑜𝑢𝑔ℎ𝑛𝑒𝑠𝑠,𝐼𝑚𝑝𝑎𝑐𝑡 𝑒𝑛𝑒𝑟𝑔𝑦 𝑑𝑖𝑠𝑠𝑖𝑝𝑎𝑡𝑖𝑜𝑛} , an independent\nsecond-order polynomial response surface model is fitted according to:\n𝑓(𝑥) = 𝛽 +∑ 𝛽 𝑥 +∑ 𝛽 𝑥2 +𝑖 < ∑ 𝛽 𝑥 𝑥 (S23)\n0 𝑖 𝑖 𝑖 𝑖 𝑖𝑖 𝑖 𝑖<𝑗 𝑖𝑗 𝑖 𝑗\nwhere 𝑥 denotes the composition variables and 𝛽 are the regression coefficients\nestimated by least squares.\nThe fitted response surfaces are subsequently used to screen the feasible design\nspace and identify candidate formulations based on the predicted multi-objective trade-\nof fs implied by the quadratic models. A total of 25 new formulations are selected in a\nsingle batch from this global screening step. No iterative model updating or explicit\nuncertainty quantification is incorporated during this process, reflecting the standard\nand widely adopted usage of RSM in formulation optimization. The resulting 25\nformulations, together with the initial 50 experimental samples, constitute the RSM-\naugmented dataset (n=75), which is subsequently used for Pareto-front construction and\noptimization-efficiency analysis.\nBenchmarking metrics. To quantitatively compare the two optimization\nstrategies, a dominance gap is defined as the area enclosed between the Pareto-front\nenvelopes obtained by the proposed machine-learning framework (PSL + AL) and the\nRSM baseline, evaluated over the overlapping strength range with respect to the\nsecondary objective under consideration.\nOptimization efficiency is further assessed using the normalized hypervolume\nimprovement as a function of the number of additional experiments 𝑘:\n32\n𝐻𝑉𝑘−𝐻𝑉(0)\n𝛥𝐻𝑉𝑛𝑜𝑟𝑚(𝑘) = (S24)\n𝐻𝑉(𝑚𝑎𝑥)−𝐻𝑉(0)\nwhere 𝐻𝑉(0) denotes the hypervolume of the initial dataset and 𝐻𝑉𝑚𝑎𝑥 is the\nhypervolume obtained after all 25 new experiments are included. This metric provides\na model-independent measure of how efficiently each strategy expands the Pareto front\nas experimental resources increase.\nResults. Under the controlled benchmarking protocol described above, the\noptimization performance of the proposed Pareto set learning combined with active\nlearning (PSL + AL) framework is assessed against a classical response surface method\n(RSM). Both strategies start from the same experimental dataset (n=50) and are\npermitted to propose an identical number of additional formulations (n=25), ensuring a\nfair and unbiased comparison. As shown in Fig. S12a and S12b, formulations suggested\nby RSM remain largely confined to intermediate-performance regions in both the\nstrength–fracture toughness and strength–impact energy dissipation spaces.\nConsequently, the corresponding Pareto fronts exhibit only limited outward expansion\nrelative to the initial dataset. This behavior is consistent with the limited\nrepresentational capacity of quadratic response surface models when approximating\nstrongly nonlinear trade-offs and navigating toward the boundaries of the attainable\nperformance envelope.\nIn contrast, PSL + AL yields Pareto fronts extending into regions that jointly\nexhibit high tensile strength and enhanced fracture toughness or impact energy\ndissipation. Across a broad overlapping strength range, the Pareto fronts obtained by\nPSL + AL systematically dominate those derived from RSM, as quantified by the\n33\ndominance gap highlighted in Fig. S12a and S12b. These results indicate that the\nlearning-guided strategy more effectively identifies candidate formulations that\nadvance strength–toughness and strength–impact trade-offs under an identical\nexperimental budget.\nBeyond Pareto-front geometry, optimization efficiency is further quantified using\nhypervolume-based metrics (Fig. S12c). PSL + AL achieves more than 90% of its\nterminal normalized hypervolume improvement within the first ten additional\nexperiments (k=10) and approaches saturation thereafter, whereas RSM exhibits only\nmarginal hypervolume gains even after completion of all 25 additional trials.\nThis pronounced difference highlights the importance of uncertainty-aware, iterative\nselection in efficiently navigating multi-objective formulation spaces.\nCollectively, these results demonstrate that PSL + AL not only expands the Pareto\nfront toward higher-performance regimes but does so with substantially improved\nsample efficiency, underscoring the practical relevance of learning-guided optimization\nfor composite material design under limited experimental budgets.\nbaseline. a. Pareto fronts in the strength–fracture toughness space obtained using the response\nsurface method (RSM, blue) and the proposed framework combining Pareto set learning with active\nlearning (PSL + AL, orange). Gray circles denote the initial experimental dataset (n=50), and open\nsymbols indicate newly proposed formulations (n=25) from each strategy. The shaded region marks\nthe dominance gap, defined as the area enclosed between the Pareto-front envelopes of PSL + AL\nand RSM within the overlapping strength range. b. Corresponding Pareto fronts in the strength–\nimpact energy dissipation space, showing a similar systematic dominance of PSL + AL over RSM\nunder an identical experimental budget. c. Normalized hypervolume improvement as a function of\nthe number of additional experiments for both strategies. PSL + AL approaches its terminal\nhypervolume value within the first ten additional experiments, whereas RSM exhibits only marginal\nimprovement over the full set of 25 trials, indicating superior sample efficiency for the uncertainty-\naware, iterative learning strategy.”\ndeviates from the experimentally verified 259 MPa, and the source of the error (such as\nmodel prediction error or experimental repeatability issues) has not been explained.\nResponse: We thank the Reviewer for pointing out this important inconsistency.\nWe agree that, as written in the previous version, the stated “optimized high-\nperformance region (strength 205–225 MPa)” could be misinterpreted as the final ML-\noptimal outcome and would appear inconsistent with the experimentally verified high\nstrength (~259 MPa). In the revised manuscript, we have corrected and clarified this\npoint in two ways: (i) we corrected the numerical range to reflect the final AL-updated\nPareto-optimal regime, and (ii) we explain why a single extreme value (e.g., 259 MPa)\nshould not be interpreted as the multi-objective optimum.\n(1) The “205–225 MPa” range is not the final optimized regime and has been\ncorrected. The 205–225 MPa interval corresponded to an early-stage/intermediate-\nperformance portion of the dataset (or a preliminary visualization range) and does not\nrepresent the final AL-updated Pareto-optimal high-performance basin. In the revised\ntext and figures, the optimized region is consistently reported based on the final AL-\nexpanded dataset and the corresponding preference-guided Pareto solutions.\n(2) Our AL–PSL framework performs multi-objective optimization under\nuncertainty-aware criteria (LCB + HVI). Accordingly, the framework does not target\nstrength maximization alone but identifies a Pareto-consistent high-performance basin\nin which strength, fracture toughness, and impact energy dissipation are simultaneously\nhigh and non-dominated. Therefore, a single point exhibiting the highest strength (e.g.,\n36\n~259 MPa) is not necessarily the ML-selected optimum if it is accompanied by weaker\ntoughness/impact performance or is less robust under uncertainty criteria. In this\ncontext, the appropriate “optimized region” is a distribution of Pareto-optimal\ncompromises, not a single-point upper bound.\n(3) Alignment with experimental validation and explanation of error sources. After\nfive AL iterations, the experimentally validated preference-guided formulations (Fig. 4\nand Table S2) cluster tightly in the final optimized basin: strength ≈ 248–250 MPa,\nfracture toughness ≈ 14.8–14.9 MPa·m1/2, and impact energy dissipation ≈ 4.65–4.71 J.\nThe small spread across these independently fabricated specimens (mean ± s.d., n = 5)\ndemonstrates good experimental repeatability and indicates that the identified region is\nphysically realizable and reproducible, rather than an artifact of model error.\nIn addition, we now explicitly support model reliability using (i) learning-curve\nstabilization across AL iterations (Fig. S10) and (ii) independent hold-out validation\nshowing strong out-of-sample agreement without systematic bias (Fig. S11). Together,\nthese results indicate that the final optimized region reflects the multi-objective and\nrobustness-aware definition of optimality adopted in this work, rather than limitations\narising from prediction accuracy.\nIn summary, we have corrected the previously confusing “205–225 MPa”\nstatement and clarified that the ML framework identifies a Pareto-optimal high-\nperformance basin under multi-objective and uncertainty-aware criteria. The\nexperimentally validated solutions fall squarely within this basin with high repeatability,\nand the apparent discrepancy with a single high-strength extreme value is explained by\n37\nthe fact that multi-objective optima are not defined by strength alone.\nComment 7: Why are there no significant differences in the strain values of the\nstress-strain curves in Figure 4a? Most of this degree of strain difference is due to error.\nResponse: We thank the Reviewer for this comment. The absence of pronounced\ndifferences in the strain values among the stress–strain curves in Fig. 4a is intentional\nand physically meaningful, rather than a consequence of experimental error.\nAll four formulations shown in Fig. 4a are selected from the same optimized high-\nperformance Pareto regime identified by the ML-guided framework and share an\nidentical trabecular-interlocked and physically entangled interfacial architecture. This\ncommon structural motif governs the dominant deformation mechanisms under tensile\nloading, including large-scale polymer chain stretching, interfacial sliding, and\nprogressive ligament engagement. As a result, the global deformation capacity (strain-\nto-failure) is primarily controlled by the interface architecture and matrix mobility,\nrather than by the fine compositional variations introduced by preference-guided\nselection.\nImportantly, the ML optimization is designed to reorder trade-offs within an\nalready optimized performance envelope, rather than to alter the fundamental\ndeformation mode. Consequently, preference-guided formulations exhibit tightly\nclustered strain responses, while subtle but reproducible differences are manifested\ninstead in stress levels, energy dissipation, and fracture-related metrics (Fig. 4d–i),\nwhich are more sensitive to the imposed preference vectors.\nobserved strain-to-failure values fall within the expected scatter for highly ductile\npolymer-based composites tested under identical conditions. No systematic shift in\nstrain behavior is observed across formulations beyond this intrinsic variability,\nsupporting that the overlap in Fig. 4a reflects mechanistic consistency and experimental\nreproducibility, rather than measurement error. Details of the revised Fig. 4 are\nprovided in our response to Comment 3.\nComment 8: Figure 4d shows the preference vectors λ=[0.032:0.073:0.896],\nλ=[0.028:0.081:0.891], and λ=[0.034:0.089:0.887], which correspond to the optimal\nimpact energy dissipation, fracture toughness, and strength of TIC, respectively. As is\nwell known, TPU has excellent fracture toughness, BF has excellent strength, and POK\nhas excellent impact energy dissipation. The mechanical test results for\nλ=[0.034:0.089:0.887] and λ=[0.032:0.073:0.896] are normal and understandable.\nHowever, the mechanical results for λ=[0.028:0.081:0.891] appear to be anomalous.\nPlease explain why, among these three materials, the fracture toughness is optimal when\nthe TPU content is reduced.\nResponse: We thank the Reviewer for carefully checking the preference vectors\nand for raising the apparent inconsistency. This concern arose from an earlier\nplotting/labeling version of Fig. 4, in which the displayed λ triplets were not the final,\npreference-resolved solutions used for experimental validation. In the revised\nmanuscript, we have replaced those λ values and updated the entire preference-guided\n39\nselection consistently using the retrained AL–PSL model and the experimentally\nvalidated solutions summarized in Fig. 3f and Table S2.\nIn the revised Fig. 3d-f, Fig. 4 and Table S2, the representative designs are\nextracted using clear engineering targets, namely a balanced design and three single-\npriority designs (strength-first, fracture-first, and impact-first). The target preference\nvectors are explicitly set as [0.333, 0.333, 0.333] (balanced), [0.85, 0.075, 0.075]\n(strength-first), [0.075, 0.85, 0.075] (fracture-first), and [0.075, 0.075, 0.85] (impact-\nfirst). The matched realized λ values closely align with these targets (Fig.3f),\nconfirming that the extracted solutions are preference-consistent after retraining with\nthe expanded dataset.\nImportantly, once the corrected preference-resolved solutions are used, the\n“anomalous” situation described by the Reviewer does not persist. The experimentally\nvalidated fracture-first formulation does not rely on a reduced TPU content. Instead, all\nfour preference-guided formulations fall within a narrow high-performance basin with\nvery similar compositions, with TPU around 3.8–3.9 wt% and BF around 10.2–10.4\nwt%. The measured fracture toughness values also cluster tightly (approximately\n14.81–14.83 MPa·m1/2), indicating that within this already optimized regime, the\npreference operation reorders trade-offs on a Pareto-optimal plateau rather than\ninducing large compositional shifts.\nMore generally, we also clarify in the revised text why “more TPU” does not\nguarantee “higher fracture toughness” once the system is constrained to a Pareto-\noptimal basin. In TICs, fracture toughness is governed by a coupled mechanism,\n40\nincluding (i) trabecular-interlock–stabilized crack deflection and fiber-bridging, and (ii)\na physically entangled interphase that enables crack-tip plastic-zone enlargement and\ndissipative chain rearrangement. Within such a mechanism-dominated basin, small\nvariations in TPU can be offset by compensating changes in BF-mediated bridging,\ninterfacial entanglement topology, and local crystallinity/phase organization. Therefore,\ntoughness optimization is not a monotonic “TPU-only” effect, and the correct\ninterpretation is that toughness emerges from architecture–interface synergy under\ncomposition constraints, not from a single constituent acting independently.\nTo avoid any remaining ambiguity, we have (1) replaced the outdated λ values in\nFig. 3d-f and Fig. 4, (2) provided the exact target λ and matched λ together with the\nextracted compositions and predicted/measured properties in Table S2, and (3) revised\nthe caption/text to make the preference-guided extraction criterion explicit (minimal\ndistance between target and realized preference vectors). This ensures that the\npreference vectors reported in the figure correspond exactly to the formulations that\nwere fabricated and tested. Details of the revised content are provided in our response\nto Comment 3.\nComment 9: The images in Figures 4g, h, and i are unclear.\nResponse: We thank the Reviewer for pointing this out. In the revised manuscript,\nwe have reformatted and enlarged these images and relocated them to the\nSupplementary Information as Figs. S13–15, where each panel is shown at higher\nresolution and with clearer contrast to facilitate visual inspection.\ncamera recordings of the impact experiments. The apparent reduction in sharpness\noriginates from the intrinsic trade-off between temporal resolution and spatial\nresolution in high-speed imaging. To provide full experimental context, the\ncorresponding high-speed videos are now included as Supplementary Videos S1 and\nS2.\nThe revised text reads as follows: “\nFigure S13 Large-strain deformation behavior of TICs under uniaxial tension. Photographs of\na representative TICs specimen captured at 0% strain and 600% engineering strain during\nmonotonic tensile loading. The sample undergoes large, uniform elongation without visible cracking\nor interfacial debonding, indicating excellent stretchability and structural integrity under extreme\ndeformation. The arrows denote the tensile loading direction.\nsequence showing the impact process of a representative TICs sample, including initial contact,\nimpact loading, and impactor rebounding. The sample remains intact after impact, with only a\nshallow surface indentation observed in the post-impact optical image (right). b. Impact sequence\nof a control composite without structural interface design (POK-BF NAFCs) under identical testing\nconditions, where catastrophic perforation and cracking occur upon impact, as evidenced by the\nlarge through-thickness damage region (right). Scale bars: 1 cm. All image frames are representative\nsnapshots extracted from the corresponding high-speed video recordings (see Video S1).\nimages showing the ballistic impact process of a TICs sample, including initial approach, projectile\nsample contact, and clear projectile rebounding after impact. The dashed arrow indicates the\nrebound direction of the projectile, demonstrating efficient momentum reversal and energy\ndissipation without perforation. b. Corresponding impact sequence for a brittle reference sample\nunder identical testing conditions, where the projectile penetrates and passes through the sample,\nindicating insufficient energy dissipation and catastrophic failure. Yellow arrows denote the\ntemporal progression of the impact event. All image frames are representative snapshots extracted\nfrom the corresponding high-speed video recordings (see Video S2).”\nComment 10: The SAXS and WAXS tests in Figures 5b, c, d, and e can only\ndemonstrate that the crystallinity and crystalline morphology of the polymer change\nunder tensile stress. In fact, most polymer materials undergo changes in SAXS and\nWAXS even with minor tensile stress. These tests are unrelated to the innovative aspect\nof the manuscript—biomimetic beams. Additionally, phenomena such as “fiber\nbridging” and “crack deflection” were only observed via SEM, without quantitative\n44\ncharacterization of the dynamic process of interfacial stress response (e.g., interfacial\nshear strength and changes in chain entanglement density under different strains).\nMolecular dynamics simulations only depicted the static state of chain entanglement\nand did not simulate the dynamic evolution under stress loading (e.g., the rate of\nentanglement breakage, or quantitative calculations of energy dissipation). It is\nrecommended to supplement in situ mechanical testing with spectroscopic\ncharacterization (e.g., Raman spectroscopy to track changes in the interfacial chemical\nenvironment) or to quantify the strain dependence of interfacial adhesion forces using\natomic force microscopy.\nResponse: We thank the Reviewer for the rigorous critique. We agree with the\ncentral concern that SAXS/WAXS signatures under tensile loading can occur in many\npolymeric systems and, by themselves, do not prove the novelty of the biomimetic\ntrabecular-interlock concept. We also agree that qualitative SEM snapshots of fiber\nbridging and crack deflection require complementary, strain-resolved quantification of\nthe interfacial stress response. In the revised manuscript, we therefore restructured the\nmechanistic section and added new in situ evidence to directly address these points.\nFirst, we revised how SAXS/WAXS are used and interpreted. In the original\nversion, the SAXS/WAXS discussion could be read as implying that crystallinity\nchanges themselves constitute the key innovation. In the revised version, SAXS/WAXS\nare explicitly repositioned as supporting measurements that track statistically averaged,\nstrain-dependent reorganization of semicrystalline domains and nanophase orientation.\nWe clarified in the text that these scattering changes are not unique to TICs, and we use\n45\nthem only to contextualize the deformation regime transitions and to bridge length\nscales between macroscopic stress–stress–strain behavior and nanoscale structural evolution.\nAccordingly, we revised the wording in the “Strain-induced structural evolution across\nlength scales” subsection so that the innovative claim is not attributed to SAXS/WAXS,\nbut rather to the architecture-enabled interfacial mechanics demonstrated by the new in\nsitu evidence described below.\nSecond, we strengthened the evidence chain for architecture-enabled toughening\nbeyond qualitative SEM. In situ tensile SEM (Fig 5a, 5b, and Video S3) is now\npresented as dynamic crack-evolution evidence that the trabecular-interlocked\narchitecture stabilizes crack growth under extreme strain, rather than as static “fracture\nmorphology.” The revised text emphasizes two strain-resolved observations: persistent\ncrack-wake bridging by load-bearing polymer ligaments anchored on basalt fibers and\nthe concurrent expansion of the plastic zone with tortuous crack advance. These\nfeatures are explicitly linked to the architectural constraint and interlocking-induced\nload transfer, and are discussed as mechanisms that suppress unstable crack opening\nand delay catastrophic interfacial separation.\nThird, to directly address the Reviewer’s request for quantitative characterization\nof the dynamic interfacial stress response, we added strain-resolved in situ Raman\nspectroscopy and integrated it as core mechanistic evidence. This addition directly\nanswers the suggestion to supplement in situ mechanical testing with spectroscopic\ntracking of the interfacial chemical environment. In the revised manuscript (Fig. 5 and\nFigs. S13–S17), we quantify two Raman descriptors as functions of tensile strain. The\n46\ncarbonyl-band centroid in the 1700–1785 cm⁻¹ window is used to track the evolution\nof the local environment and stress repartitioning of polar groups, providing a bond-\nlevel proxy for whether stress accumulates monotonically or is dynamically\nredistributed. The CH-stretch envelope full width at half maximum in the 2800–3100\ncm⁻¹ window is used to quantify conformational dispersion and segmental mobility,\nproviding a chain-level proxy for rearrangement and sliding. Importantly, we\nemphasize in the revised text that these Raman metrics are interpreted together with in\nsitu SEM. The SEM shows sustained mechanical engagement through bridging\nligaments and plastic-zone enlargement, while Raman shows that stress does not\nlocalize monotonically but is dynamically repartitioned and accompanied by enhanced\nchain mobility at high strain. This cross-modal convergence provides direct, strain-\nresolved evidence of the proposed stress-adaptive interfacial response, rather than post\nhoc inference from scattering or fracture surfaces.\nFourth, we added a control comparison to demonstrate that the observed dynamic\nresponse is not a generic polymer effect. The revised manuscript explicitly contrasts\nTICs with control composites lacking the designed architectural and interfacial features.\nIn the control, the Raman descriptors evolve more monotonically and the\nmicrostructural response lacks the same sustained bridging and crack stabilization. This\ncontrol comparison is used to support the causal interpretation that the dynamic stress\naccommodation arises from the physically entangled, architecture-stabilized interface\nrather than from unavoidable crystallinity changes common to polymers.\nFifth, we revised how molecular dynamics simulations are used in the mechanistic\n47\nargument. We agree that static MD snapshots alone cannot quantify the rate of\nentanglement breakage or energy dissipation under load. Therefore, we revised the\nmanuscript to explicitly define the role of MD as providing atomistic insight into the\nformation and topology of the entangled interphase and its strain-dependent\nconfigurational tendency, while avoiding claims that require true non-equilibrium\nmechanical loading simulations. All statements about dynamic stress adaptation and\ndissipation pathways are now anchored primarily in the newly added strain-resolved in\nsitu experiments, with MD used for complementary molecular interpretation.\nFinally, regarding the Reviewer’s suggestion to quantify interfacial adhesion\nforces by AFM, we acknowledge that AFM-based interfacial shear mapping in a porous,\nfiber-reinforced, heterogeneous trabecular architecture is experimentally nontrivial and\nnot directly comparable to the in situ crack-tip loading state. Instead, we addressed the\ncore intent of the suggestion by providing strain-resolved, quantitative proxies for\ninterfacial stress accommodation using Raman spectroscopy and by coupling this with\ndirect crack-evolution visualization by in situ SEM. Together, these additions provide\nthe requested dynamic, mechanistically specific evidence for how the interface\nresponds and adapts under loading.\nIn summary, the revision makes three concrete changes that directly address the\nReviewer’s criticism: the SAXS/WAXS discussion is repositioned as supportive rather\nthan novelty-defining, dynamic crack evolution is strengthened through in situ tensile\nSEM, and quantitative, strain-resolved interfacial stress accommodation is added\nthrough in situ Raman spectroscopy with appropriate control comparisons. The revised\n48\ncontents are incorporated into the “Mechanism analysis” section and Supplementary\nInformation, including the newly added in situ mechanistic evidence (Figs. 5–6, Figs.\nS13–S17, and Video S3).\nThe manuscript has been revised as follows: “\nMechanism analysis\nDynamic stress-adaptive interfacial response\nThe tight clustering of validated properties and the plateau-like distribution on\nAshby maps suggest that the optimized designs share a common, stable dissipation\nmechanism. We therefore interrogate the dynamic interfacial response of a\nrepresentative TIC formulation (balanced preference) using in situ tensile SEM and\nstrain-resolved Raman spectroscopy, complemented by impact/ballistic observations.\nThis formulation sustains tensile strains up to ~600% without catastrophic failure (Fig.\nS13). In drop-weight impact tests, high-speed imaging captures rapid deformation\nfollowed by rebound (Fig. S14a, Video S1), and post-impact inspection shows only\nshallow indentation without macroscopic cracking. By contrast, control specimens\nlacking the trabecular-interlocked architecture and adaptive interface design (POK–BF\nNAFCs) exhibit immediate penetration under identical conditions (Fig. S14b). Under\nballistic loading (80 m s-1), TICs decelerate the projectile and rebound without\nperforation (Fig. S15a), whereas NAFCs show complete penetration (Fig. S15b, Video\nS2). These macroscopic observations motivate the following in situ analyses aimed at\nresolving how stress is redistributed and dissipated at the crack and interfacial scales.\nIn situ tensile SEM directly visualizes the crack evolution and load-transfer\n49\nmechanisms operating in TICs. At an applied strain of ~300%, crack opening is\neffectively shielded by pronounced fiber-mediated bridging, where multiple polymer\nligaments anchored on basalt fibers span the crack faces and remain mechanically\nengaged (Fig. 5a, yellow arrows). Instead of propagating along a straight path, the crack\nfollows a tortuous and deflected trajectory, indicating crack deflection induced by\ninterfacial heterogeneity and mechanical interlocking. Upon further loading to ~500%\nstrain (Fig. 5b and Video S3), the bridging mechanism remains operative. The load-\nbearing ligaments undergo progressive stretching and thinning with delayed rupture,\nrather than abrupt interfacial debonding or clean fiber pull-out. Crack advance remains\ngradual and continues along a deflected path, accompanied by a pronounced expansion\nof the surrounding plastic zone (Fig. 5b, blue arrows). These observations demonstrate\nthat the TICs interface does not fail through brittle separation. Instead, it sustains\ncontinuous mechanical engagement across extreme strains through physically\nentangled, load-bearing ligaments that redistribute stress in the crack-tip region, thereby\nstabilizing crack growth and enabling sustained energy dissipation.\nWe next quantify molecular-scale stress accommodation using strain-resolved\ntensile Raman spectroscopy. The stacked spectra (Fig. 5c) reveal a continuous evolution\nof the full Raman window with increasing strain. Consistently, the strain–wavenumber\nheatmap (Fig. 5d) highlights two distinct response regimes: the carbonyl window\n(1700–1785 cm-1) exhibits strain-dependent but spatially diffuse intensity redistribution,\nwhereas the CH-stretch window (2800–3100 cm-1) progressively develops increased\nspectral dispersion at higher strain, indicative of enhanced conformational\nheterogeneity. To quantify the response of polar groups, we track the centroid position\nof the carbonyl band (Fig. 5e). In TICs, the centroid shift remains relatively small and\nnon-monotonic across the applied strain range, indicating that tensile stress does not\ncontinuously accumulate on a narrow subset of carbonyl bonds. Instead, stress is\ndynamically repartitioned as the physically entangled interface undergoes progressive\nreconfiguration under load. In contrast, the POK–TPU control (Fig. S17) exhibits a\nmore systematic and monotonic centroid evolution, consistent with a more\nhomogeneous polymer response in which bond-level stress localization can be directly\ntracked spectroscopically. Chain-level dissipation is further captured by the full width\nat half maximum (FWHM) of the CH-stretch envelope (Fig. 5f). TICs display a\npronounced broadening of the CH envelope at high strain, consistent with the activation\nof segmental mobility, chain rearrangement, and sliding. By comparison, the control\nsample shows a substantially weaker broadening response, indicating more limited\nconformational freedom under tensile loading.\nThese molecular-scale signatures are fully consistent with the in situ tensile SEM\nobservations of sustained fiber-bridging ligaments and an expanded plastic zone at the\ncrack. Together, the spectroscopic and microstructural evidence demonstrates that the\ntrabecular interlocking architecture enables a dynamically activated, stress-adaptive\ninterfacial response, which progressively converts tensile loading into frictional and\nentanglement-mediated dissipation pathways. This multiscale stress redistribution\nmechanism underpins the enhanced crack resistance and impact tolerance observed in\nTICs.\nRaman spectroscopy. a, b. In situ tensile SEM images of crack evolution in TICs at ~300% and\n~500% tensile strain, respectively. At ~300% strain (a), polymer ligaments anchored on basalt fibers\nbridge the crack faces (yellow arrows), accompanied by crack deflection. At ~500% strain (b), the\nbridging ligaments remain engaged while an expanded plastic zone develops ahead of the crack tip\n(blue arrows), leading to gradual crack advance. c, d. Strain-resolved Raman spectra and\ncorresponding strain–wavenumber heatmap acquired during tensile loading, showing distributed\nspectral changes in the carbonyl region (1700–1785 cm-1) and progressive dispersion in the CH-\nstretch region (2800–3100 cm-1). e, f. Quantitative Raman descriptors showing a small, non-\nmonotonic evolution of the carbonyl-band centroid and a pronounced high-strain increase in the\nCH-stretch envelope full width at half maximum (FWHM) for TICs.\nIn situ tensile deformation is performed inside a scanning electron microscope\nusing a custom micro-mechanical tensile stage equipped with a 4000 N load cell. The\ntensile stage is mounted directly inside the SEM chamber (Fig. S17a), enabling real-\ntime mechanical loading under vacuum conditions. Specimens are clamped between\nthe tensile grips and are loaded under displacement-controlled mode at a constant cross-\nhead speed of 10 µm s-1. SEM imaging is conducted continuously during loading,\nallowing direct observation of crack initiation, crack opening, and subsequent crack\npropagation at the same region of interest (Fig. S17b). All in situ tensile SEM\nexperiments were carried out under identical loading and imaging conditions to ensure\nreproducibility.\nFigure S16 In situ tensile SEM experimental setup. a. Photograph of the custom micro-\nmechanical tensile stage integrated into the SEM chamber, showing the motor-driven loading unit\nand specimen gripping configuration. b. Representative SEM view under vacuum with the specimen\nmounted between the tensile grips, enabling continuous imaging of the same region during tensile\nloading.\nSupplementary Note 5. In situ tensile Raman spectroscopy\nIn situ tensile Raman measurements are performed using a confocal Raman\nmicroscope equipped with a 633 nm excitation laser, a 600 lines mm-1 grating, and a 5×\noperated under displacement-controlled mode at a constant stretching rate of 300 µm s-\n1. Specimens are stretched to predefined strain levels and held during spectral\nacquisition following a step-and-hold protocol. Raman spectra are collected at a fixed\nlocation on the specimen surface with an integration time of 10 s and 5 accumulations\nat each strain level. Autofocus is disabled throughout the measurements to avoid focus-\ninduced spectral artefacts, and the focal position is optimized prior to testing. Spectral\nacquisition covers the range of 800–3200 cm-1, including both the fingerprint region\nand the CH-stretch region. All measurements are conducted under identical optical and\nmechanical conditions to ensure consistency and comparability between samples.\nFor data analysis, all spectra are first corrected for background using polynomial\nbaseline subtraction and normalized to the integrated intensity over the full spectral\nwindow to minimize variations arising from laser power fluctuation and collection\nefficiency. The processed spectra are assembled into a two-dimensional strain–\nwavenumber intensity map to visualize strain-dependent spectral evolution. To enhance\ncontrast and reveal subtle redistribution of spectral intensity across strain, a uniform\ngamma correction is applied to the intensity scale for all datasets.\nQuantitative Raman descriptors are extracted from two predefined spectral\nwindows. The carbonyl region is defined between 1700 and 1785 cm-1, within which\nthe spectral centroid is calculated at each strain level as the intensity-weighted average\nRaman shift. This centroid reflects changes in the local vibrational environment and\nconstraint of carbonyl groups during deformation. The CH-stretch region is defined\nbetween 2800 and 3100 cm-1, and the full width at half maximum (FWHM) of the\ncomposite CH envelope is determined by fitting the intensity profile to a single effective\nband after smoothing. The CH-envelope FWHM serves as a descriptor of\nconformational dispersion and segmental mobility. The strain-dependent carbonyl\ncentroid shifts and CH-envelope FWHM values are plotted as a function of tensile strain\nto generate the quantitative curves shown in Fig. 5e and Fig. 5f. Identical spectral\nwindows, fitting procedures, and processing parameters are applied to both TICs and\nthe POK–TPU control to ensure direct comparability.\nFigure S17 Strain-resolved in situ Raman spectra of POK–TPU sample. a. Offset Raman\nspectra of POK–TPU sample acquired at successive tensile strain levels under step-and-hold loading\nconditions. The color scale indicates the applied tensile strain. b. Corresponding Raman intensity\nheatmap as a function of strain, highlighting the carbonyl (1700–1785 cm-1) and CH-stretch (2800–\n3100 cm-1) regions. Intensities are normalized to facilitate comparison across strain levels.”\nwas applied to explore the optimal composition of a three-component polymer\ncomposite material that forms a bone trabecular interconnection structure, aiming to\nimprove and balance three mechanical properties. As a result, a polymer composite with\nproperties equivalent to or exceeding those of conventional engineering composites and\nbiomaterials was successfully produced. Through comprehensive structural analysis,\nmechanical testing, and MD simulations, the mechanisms that enable the fabricated\nmaterial to exhibit high strength, toughness, and impact resistance were elucidated.\nIn recent years, there has been increasing attention to research on the design and\ndevelopment of polymer materials utilizing machine learning, but studies focusing on\nthe mechanical properties of polymer composites relatively few. In this context, the\nproof of concept demonstrated in this study is significant.\nHowever, the material system under investigation (bone trabecular interconnection\nstructure) is somewhat limited, and the mechanical performance achieved is not yet a\ndrastic improvement over conventional materials. Therefore, it would likely attract\nmore attention if published in a specialized journal in polymer chemistry or materials\nscience, rather than a general journal like Nature Communications, which targets a\nbroader audience.\nThe methodology and the path to the conclusions are deemed rational. I would like to\noffer the following minor comments mainly regarding the machine learning aspects of\nthe study. I hope these can contribute to improving the paper.\nResponse: We thank the Reviewer for the careful evaluation and constructive\n56\ncomments on the scope and significance of this work. While the present study focuses\non a specific polymer composite system with a trabecular interlocked architecture, the\nprimary contribution of this work lies not in a single material formulation, but in the\ngeneralizable integration of multi-objective machine learning, active experimental\ndesign, and mechanism-resolved validation for mechanically demanding polymer\ncomposites.\nImportantly, the proposed ML-guided framework is material-agnostic and can be\nreadily extended to other polymer matrices, reinforcement types, and interfacial\narchitectures, as well as to alternative combinations of mechanical and functional\nobjectives. Moreover, the mechanistic insights revealed here—namely stress-adaptive\ninterfacial reconfiguration and entanglement-mediated energy dissipation—are not\nspecific to the studied system but represent broad physical principles applicable to a\nwide class of composite materials.\nTo further emphasize the general relevance of the approach, we have strengthened\nthe manuscript by clarifying the methodological transferability, expanding\nbenchmarking against representative engineering composites and polymers (Table S3),\nand explicitly discussing the broader implications of the ML-driven design paradigm in\nthe main text and Supplementary Information. We therefore believe that the\ncombination of general methodology, mechanistic insight, and experimentally validated\nmulti-objective optimization makes this work suitable for a broad readership, consistent\nwith the scope of Nature Communications.\nconsisting of the composition formula and process parameters. However, any\nexplanation of the process variables is not provided.\nResponse: We thank the Reviewer for pointing out this ambiguity. We clarify that,\nin the present study, the machine learning models do not include independent\nprocessing parameters as explicit input variables.\nAll machine learning models are trained using composition variables only, namely\nthe weight fractions of TPU and BF, with the POK content implicitly constrained by\nmass conservation. The processing conditions (including mixing procedure, molding\ntemperature, pressure, and curing protocol) are kept strictly identical for all samples\nand therefore do not introduce additional degrees of freedom in the dataset. As a result,\nprocessing parameters are treated as controlled constants rather than model inputs, and\ntheir effects are implicitly embedded in the experimental measurements.\nThe revised text now explicitly states that the input feature space consists solely\nof compositional descriptors, while processing conditions are fixed across all\nexperiments. This clarification ensures consistency between the described model inputs,\nthe experimental design, and the reported results, and avoids introducing ill-defined or\nnon-informative variables into the surrogate models.\nThe corresponding modifications in the Supplementary Information are as follows:\n“ Note 3. Machine learning (ML) driven multi-objective optimization\nTo efficiently explore the high-dimensional formulation space while minimizing\nexperimental burden, we integrate machine learning (ML) with a structured\nexperimental design strategy. An initial design-of-experiments (DoE) protocol is used\nto select 50 representative composite formulations that uniformly span the admissible\ncomposition space. Each formulation corresponds to a unique combination of polymer\nconstituents, while all processing conditions are kept strictly identical across samples\nto eliminate confounding effects from processing variability.”\nComment 2: It is stated that the model input includes cubic terms of the composition\nformula. Is this description correct? If so, please explain why the cubic terms are\nincluded without the quadratic terms. Additionally, does the model include cubic\ninteraction terms? A more precise description is required.\nResponse: We thank the Reviewer for raising this important point and for\nrequesting a more precise description of the input feature construction. We clarify\nbelow and have revised both the main text and Supplementary Information accordingly.\nClarification of polynomial terms. In the revised manuscript, we clarify that the\nsurrogate model does not include only cubic terms. Instead, the input feature space is\nconstructed using a hierarchical polynomial expansion that includes linear, quadratic,\nand cubic terms of the compositional variables. This expansion is applied consistently\nto capture nonlinear composition–property relationships while maintaining model\ninterpretability. Specifically, starting from the independent compositional variables\nTPU and BF (with POK determined by mass conservation), the pure-power terms\ninclude 𝑥 , 𝑥2, and 𝑥3, where 𝑥 denotes the standardized composition variables.\n𝑖 𝑖 𝑖 𝑖\nImportantly, cubic terms are included in addition to, not instead of, quadratic terms. The\n59\nearlier wording in the manuscript may have unintentionally suggested that cubic terms\nwere used alone, which we have now corrected.\nOn cubic interaction terms. In addition to pure-power terms, we include interaction\nfeatures up to total degree three in the two independent composition variables, TPU and\nBF. Specifically, the feature set contains the second-order interaction (𝑥 𝑥 ) and\n𝑇𝑃𝑈 𝐵𝐹\nthe third-order interaction terms (𝑥𝑇𝑃𝑈𝑥 and 𝑥 𝑥2 ). POK is determined by mass\n2 𝐵𝐹 𝑇𝑃𝑈 𝐵𝐹\nconservation ( 𝑥 = 100−𝑥 −𝑥 ) and is therefore not treated as an\n𝑃𝑂𝐾 𝑇𝑃𝑈 𝐵𝐹\nindependent input. As a result, no three-variable interaction term of the form\n𝑥 𝑥 𝑥 is introduced in the model.\n𝑇𝑃𝑈 𝐵𝐹 𝑃𝑂𝐾\nRationale for including cubic univariate terms. Although Gaussian process\nregression (GPR) kernels are inherently capable of representing nonlinear functions,\nthe inclusion of explicit polynomial features provides a more physically interpretable\nand structured representation of curvature in the composition–property relationship. In\nparticular, cubic univariate terms allow the surrogate to capture asymmetric response\nand saturation behavior commonly observed in polymer–fiber composite systems, such\nas non-monotonic effects arising from TPU–BF interfacial interactions. These\npolynomial features are introduced in addition to, rather than in place of, the GPR\nkernel and serve as a structured inductive bias that enhances surrogate expressiveness\nin data-scarce regimes, without introducing high-order cross-interaction terms or\nrelying on excessively flexible kernel constructions. Potential overfitting is mitigated\nthrough z-score standardization of all inputs and kernel hyperparameter optimization\nvia marginal likelihood maximization with multiple restarts. Model generalization is\n60\nfurther assessed using repeated cross-validation and learning-curve analysis, which\nshow consistent improvement in predictive accuracy with increasing training size (Figs.\nS7 and S10), together with independent hold-out validation of the final surrogate\nmodels after active learning (Fig. S11).\nThe corresponding modifications in the Supplementary Information are as follows:\n“Supplementary Note 3.1 Gaussian Process Regression (GPR) modeling\nInput features and target variables. The input features are defined by the\nexperimentally controlled formulation parameters, specifically the weight fractions of\nthermoplastic polyurethane (TPU) and basalt fiber (BF). The content of\npoly(octenylene ketone) (POK) is implicitly determined by mass conservation and is\ntherefore not treated as an independent variable. This formulation ensures that the input\nspace remains physically admissible without introducing redundant degrees of freedom.\nTo enhance surrogate expressiveness under limited data conditions, the standardized\ncompositional variables are augmented using a polynomial feature expansion up to\nthird order, including linear, quadratic, and cubic pure-power terms as well as\ninteraction terms up to third order. This structured feature construction enables explicit\nrepresentation of curvature and compositional coupling effects, while the nonlinear\nmapping between inputs and properties is subsequently learned through the GPR kernel.\nThree mechanical properties are modeled as independent targets, namely tensile\nstrength, fracture toughness, and impact energy dissipation. Each target is described\nusing a separate GPR surrogate model. This decoupled modeling strategy avoids\nimposing artificial correlations between objectives and preserves the interpretability of\n61\nindividual structure–property relationships, which is critical for preference-guided\noptimization and for physically meaningful analysis of trade-offs among competing\nmechanical responses.”\nComment 3: It is mentioned that feature standardization and outlier removal were\nperformed, but the specific procedures are not explained.\nResponse: We thank the Reviewer for requesting clarification. We have now\nprovided a precise and explicit description of both feature standardization and outlier\nhandling procedures in the Methods and Supplementary Information.\nFeature standardization. All input features (TPU and BF contents) are standardized\nusing z-score normalization prior to model training. During cross-validation and\nlearning-curve analysis, standardization is performed strictly within each training fold.\nThe mean and standard deviation are computed exclusively from the training subset of\neach split and then applied to the corresponding validation subset. At no point are\nvalidation or test samples used to estimate standardization parameters, thereby fully\npreventing information leakage and ensuring unbiased evaluation of predictive\nperformance.\nOutlier filtering (per-target ±3σ exclusion). Outlier handling is performed using a\nper-target three-sigma (±3σ) criterion applied independently to each of the three\nexperimentally measured properties (tensile strength, fracture toughness, and impact\nenergy dissipation). For a given training dataset, the mean and standard deviation are\ncomputed separately for each target. A sample is excluded from training if any one of\n62\nits target values lies outside ±3σ of the corresponding training-set distribution. No\nmultivariate or distance-based outlier criterion is used.\nTwo leakage-free stages are employed:\n1. Cross-validation and learning-curve analysis. During cross-validation and\nlearning-curve construction, outlier filtering is applied only to the training subset of\neach split, with all ±3σ thresholds determined solely from the training targets.\nValidation samples are never used to define thresholds and are left completely unaltered\nthroughout model evaluation.\n2. Final model training for optimization. After hyperparameter selection and\nmodel assessment are completed, a global per-target ±3σ filtering is applied to the final\ntraining dataset used to fit the surrogate models for prediction and multi-objective\noptimization. For independent evaluation, the hold-out test set is defined in advance\nand remains strictly excluded from standardization parameter estimation, outlier-\nthreshold determination, and model selection.\nThis two-stage, leakage-free preprocessing strategy improves robustness against\noccasional experimental anomalies while preserving a strict separation between\ntraining, validation, and test data. The procedures are now explicitly documented in\nSupplementary Note 3.1, and their effectiveness is reflected in the stabilized learning\nbehavior observed in the learning curves (Figs. S7 and S10) and in the strong\ngeneralization performance on the independent hold-out test set (Fig. S11).\nComment 4: In Fig. S7, the predictive accuracy of the machine learning models is\n63\nshown. However, instead of using the hold-out method, the author should report the\nmean and variance of the prediction performance after performing a different split of\nthe training and testing sets.\nResponse: We thank the Reviewer for this constructive suggestion. Following this\ncomment, the evaluation protocol and the corresponding figures have been revised and\nexpanded.\nFirst, Fig. S7 has been updated to report predictive performance obtained from\nrepeated five-fold cross-validation with multiple randomized data splits, rather than\nfrom a single hold-out partition. For each training set size, the coefficient of\ndetermination (R2), mean absolute error (MAE), and mean squared error (MSE) are\ncomputed on the validation folds and averaged across different splits. The associated\nvariability is explicitly quantified by the standard deviation across cross-validation runs,\nwhich is shown as shaded uncertainty bands in Fig. S7. The revised figure caption and\nSupplementary Note 3.1 now clearly state this procedure.\nIn addition, to further assess model stability and learning behavior beyond a single\ntraining stage, we have added Fig. S10, which presents learning curves over five\nsequential active learning iterations. In Fig. S10, predictive performance is evaluated\nas a function of training set size within each iteration using cross-validation, and both\nthe mean performance and its variance across different data splits are reported. This\nanalysis demonstrates not only monotonic improvement with increasing data but also\nprogressive stabilization of model performance as active learning proceeds.\nTogether, the revised Fig. S7 and the newly added Fig. S10 provide a statistically\n64\nrobust evaluation of predictive accuracy by reporting the mean and variance of model\nperformance under different training–testing splits, thereby fully addressing the\nReviewer’s concern. For completeness, an independent hold-out test set is still used\nseparately to assess final out-of-sample generalization of the surrogate models (Fig.\nS11), complementing the cross-validation–based analyses.\nThe corresponding modifications are as follows: “\nSupplementary Note 3.1 Gaussian Process Regression (GPR) modeling\nResults. Using the initial set of 50 experimental samples, the Gaussian process\nregression (GPR) surrogate models exhibit systematic improvement in predictive\nperformance as the effective training size increases from 30 to 50 (Fig. S7). For all three\ntarget properties, namely strength, fracture toughness, and impact energy dissipation,\nthe cross-validated coefficient of determination (R2) increases monotonically,\naccompanied by concurrent decreases in both the mean absolute error (MAE) and mean\nsquared error (MSE). These consistent trends across complementary performance\nmetrics indicate progressively enhanced generalization accuracy rather than metric-\nspecific artifacts.\nIn parallel, the uncertainty bands (±1 standard deviation across repeated cross-\nvalidation splits) contract steadily with increasing training size, reflecting reduced\nsensitivity to data partitioning and stabilization of the posterior predictions. Toughness-\nrelated properties, particularly impact energy dissipation, exhibit comparatively larger\nuncertainty at smaller sample sizes, consistent with their higher sensitivity to\nmicrostructural heterogeneity and experimental variability. Nevertheless, by 45–50\nsamples, all three models enter a statistically reliable accuracy regime characterized by\nlow absolute errors, reduced variance, and diminishing marginal performance gains.\nTogether, these results demonstrate that the initial DoE dataset provides a robust\nand physically meaningful surrogate foundation for subsequent uncertainty-aware\nmulti-objective optimization, while also quantitatively motivating the use of active\nlearning to further refine model fidelity in performance-critical regions of the design\nspace.\nFigure S7 Learning curve analysis of Gaussian process regression (GPR)\nsurrogate models. a–c. Learning curves for predicting strength, fracture toughness,\nand impact energy dissipation, respectively, as functions of the training set size (N = 30,\n35, 40, 45, and 50). For each target, three complementary cross-validation metrics are\nreported: the coefficient of determination (R2, top row), mean absolute error (MAE,\nmiddle row), and mean squared error (MSE, bottom row). Solid lines represent the\nmean performance over repeated five-fold cross-validation, and shaded regions indicate\n±1 standard deviation across validation splits. The monotonic increase in R2 together\nwith the concurrent decrease in MAE and MSE, as well as the progressive narrowing\nof the uncertainty bands, demonstrate systematic improvement in generalization\naccuracy and stabilization of the surrogate models with increasing training data.\nSupplementary Note3.3 Active Learning (AL)\nis summarized by the learning-curve analysis shown in Fig. S10. Within each iteration,\npredictive accuracy improves systematically as the training set size increases, as\nevidenced by monotonic increases in cross-validated R2 and corresponding decreases\nin mean absolute error (MAE) for strength, fracture toughness, and impact energy\ndissipation. This consistent intra-iteration behavior indicates that the surrogate models\nbenefit directly from additional training data and do not exhibit signs of overfitting. In\nparallel, the uncertainty associated with cross-validation progressively decreases with\nincreasing training set size, reflecting reduced sensitivity to data partitioning and\nenhanced model stability. These trends demonstrate well-behaved learning\ncharacteristics of the Gaussian process surrogates within each AL round.\nBeyond these intra-iteration effects, a clear inter-iteration improvement is\nobserved across successive AL rounds. For a given training set size, the baseline\npredictive performance shifts upward from iteration 1 to iteration 5, accompanied by\nan overall reduction in MAE and tighter cross-validation dispersion. This cumulative\nimprovement is particularly pronounced for fracture toughness and impact energy\ndissipation, which display larger variability and broader uncertainty at early stages but\nbecome substantially more stable and accurate in later iterations. In contrast, strength\nprediction remains comparatively robust even in early iterations, consistent with a\n67\nsmoother and less complex response surface with respect to the available compositional\ndescriptors.\nAt larger training set sizes, the learning curves for all three objectives exhibit\nsaturation behavior, characterized by diminishing marginal gains in both R2 and MAE.\nThis plateau indicates that, within the explored design space, the surrogate models\napproach a near-converged performance regime. These results show that the AL\nstrategy not only improves average predictive accuracy but also accelerates the\nstabilization of surrogate models, particularly for objectives with higher intrinsic\nvariability, thereby enabling more efficient allocation of experimental effort toward\nregions of greatest uncertainty.\nFigure S10 Learning curves over five sequential active learning iterations. a–e. Learning curves\ncorresponding to active learning iterations 1–5, respectively. Within each iteration, predictive\nperformance is evaluated as a function of training set size using cross-validation. Upper panels\nreport the coefficient of determination (R2), and lower panels report the mean absolute error (MAE),\nfor strength, fracture toughness, and impact energy dissipation. Solid lines denote the mean cross-\nvalidation performance, while shaded regions represent variability across repeated cross-validation\nsplits. Increasing training set size consistently improves predictive accuracy and reduces uncertainty.\nComparison across panels (a–e) shows cumulative improvements in baseline performance and\n68\nprogressive stabilization of learning behavior, with particularly pronounced gains for fracture\ntoughness and impact energy dissipation. At larger training sizes, all objectives exhibit saturation\nbehavior, indicating diminishing marginal gains within the explored design space.\nComplementing the learning-curve analysis, we further assess the final surrogate\nmodels using an independent hold-out test set (Fig. S11). This test set is strictly\nexcluded from both model training and final refitting and is used solely for out-of-\nsample performance evaluation. Across all three target mechanical properties, the\npredicted values show close agreement with experimental measurements, as indicated\nby high coefficients of determination and low prediction errors.\nAmong the three targets, tensile strength is predicted with the highest fidelity (R2\n= 0.959, MAE = 1.364 MPa), reflecting a comparatively smooth and well-structured\nresponse surface with respect to the selected compositional descriptors within the\nexplored domain. Fracture toughness also exhibits excellent predictive accuracy (R2 =\n0.944, MAE = 0.153 MPa·m1/2), demonstrating that the surrogate model reliably\ncaptures variations in toughness despite the increased complexity of this property.\nImpact energy dissipation remains the most challenging target; nevertheless, it is\npredicted with reasonable accuracy (R2 = 0.758, MAE = 0.184 J), consistent with its\nhigher intrinsic variability under the present experimental conditions.\nNo systematic deviation from the parity line is observed across the independent\ntest set, indicating the absence of pronounced prediction bias and supporting robust out-\nof-sample generalization within the explored composition space. Taken together with\nthe progressive improvements observed in the learning-curve analysis, these results\n69\ndemonstrate that the surrogate models achieve reliable and transferable predictive\nperformance, providing a sound foundation for subsequent multi-objective\noptimization and data-driven materials design.\nFigure S11 Independent test validation of surrogate model predictions. Parity plots comparing\npredicted and experimentally measured values for a. strength, b. fracture toughness, and c. impact\nenergy dissipation on an independent test set excluded from model training. The solid lines indicate\nperfect agreement (y = x). The coefficient of determination (R2), mean absolute error (MAE), and\nroot mean squared error (RMSE) are reported for each target. The close alignment of predictions\nwith the parity line demonstrates strong generalization performance of the surrogate models across\nall three mechanical properties.”\nComment 5: In multi-objective optimization, how are the three material properties\nwith different scales adjusted? Without scale transformation, the optimization is\ngenerally biased towards the material properties with larger scales.\nResponse: We thank the Reviewer for raising this important point regarding scale\nbias in multi-objective optimization. We confirm that explicit scale normalization is\napplied to ensure that no objective dominates the optimization due to differences in\nnumerical magnitude.\nIn the present work, all target properties are normalized prior to multi-objective\noptimization. Specifically, for each mechanical property (strength, fracture toughness,\nand impact energy dissipation), the Gaussian process surrogate predictions are\ntransformed to a dimensionless normalized scale using min–max normalization based\non the predicted property ranges within the feasible design space. This normalization\n70\nmaps each objective to a comparable range and ensures balanced contribution during\nscalarization and Pareto-front construction.\nThis normalized representation is consistently used throughout all subsequent\noptimization steps, including the construction of scalarized objectives in Pareto set\nlearning (PSL), the evaluation of hypervolume and hypervolume improvement (HVI),\nand active learning–based candidate selection. As a result, the optimization process is\ninvariant to the absolute numerical scales of individual properties and is guided solely\nby the relative trade-offs encoded by the preference vectors. This normalization\nprevents bias toward objectives with larger raw magnitudes, such as tensile strength,\nand enables fair exploration of the multi-objective performance landscape. To avoid\nambiguity, we clarify in the Methods and Supplementary Information that all objectives\nare transformed to a common normalized scale prior to scalarization and Pareto analysis.\nComment 6: Eq S12: Please change the first letter of \"Min\" to lowercase.\nResponse: We thank the Reviewer for pointing this out. The first letter of “Min”\nin Eq. S12 has been corrected to lowercase (“min”) in the revised Supplementary\nInformation.\nComment 7: In Table S3, a performance comparison with existing materials, such as\nengineering composites, is shown, but it is necessary to describe specifically what\nmaterials are being referred to.\nResponse: We thank the Reviewer for this helpful suggestion. In response, we\n71\nhave revised Table S3 and its caption to explicitly specify the material categories and\nrepresentative systems included in the comparison.\nSpecifically, the category of engineering composites now explicitly refers to\ncontinuous fiber–reinforced and advanced composite systems, including continuous\ncarbon fiber–reinforced polymers (CFRP), continuous glass fiber–reinforced polymers\n(GFRP), Kevlar-reinforced composites, natural fiber composites (e.g., hemp, flax, and\njute–reinforced polymers), and representative nanocomposites incorporating nanoscale\nreinforcements such as carbon nanotubes (CNTs) and graphene.\nIn addition, the category of engineering polymers is now clearly defined to include\ncommonly used structural and engineering plastics, such as acrylonitrile butadiene\nstyrene (ABS), polycarbonate (PC), polyamide (PA), polyether ether ketone (PEEK),\npolyketone (POK), polyethylene (PE), poly(methyl methacrylate) (PMMA),\npolystyrene (PS), polytetrafluoroethylene (PTFE), polyvinyl chloride (PVC),\npolyoxymethylene (POM), polyester, epoxy resins, and polylactic acid (PLA).\nThese clarifications have been incorporated directly into Table S3 and the\ncorresponding Supplementary Information to ensure that the basis of comparison is\ntransparent and unambiguous.\nReviewer #3: They propose a universal toughening strategy that integrates a bone-\ninspired trabecular interlock architecture with a thermodynamically driven, stress-\nadaptive interface to enable efficient energy dissipation under mechanical loading.\nin its approach. Moreover, it does not adequately clarify the chemical or physical\nmechanisms underlying its core material design methodology, nor does it provide a\nclear explanation of the principles governing the dynamic behavior of the proposed\nsystem. Therefore, I recommend rejection.\nResponse: We thank the Reviewer for the careful reading of our manuscript and\nfor raising important concerns regarding the mechanistic clarity and the role of machine\nlearning. In response, we have substantially revised the manuscript and performed a\nseries of new experiments that directly address these points. In particular, we have\nadded in situ tensile scanning electron microscopy (SEM) and strain-resolved in situ\nRaman spectroscopy, which provide direct, real-time, and multiscale experimental\nevidence for the proposed stress-adaptive interfacial mechanism. These additions\nfundamentally strengthen the mechanistic foundation of the work and clarify the\nnecessity of the ML framework.\nSpecific Comments:\nComment 1: The chemical or physical mechanism behind the proposed trabecular\ninterlock architecture for interface toughening is not clearly explained.\nResponse: We thank the Reviewer for pointing out the need to clarify the chemical\nand physical mechanisms underlying the trabecular interlock architecture. In the\nrevised manuscript, we have substantially strengthened the mechanistic explanation by\nexplicitly linking architecture, interfacial physics, and deformation dynamics,\n73\nsupported by newly added in situ experiments and multiscale characterization.\nFirst, the trabecular interlock architecture provides a geometric and mechanical\nframework for stress redistribution. The TICs architecture consists of a continuous,\nthree-dimensionally interconnected porous skeleton reinforced by aligned fibers\nembedded within a soft–hard polymer matrix. This trabecular geometry is not\nintroduced merely as a structural motif, but functions to (i) diffuse applied stress over\na large interfacial area, (ii) promote crack deflection and tortuosity, and (iii) stabilize\nfiber-mediated bridging during deformation. Cross-sectional SEM and 3D\nreconstruction (Fig. 6a and 6h) directly show that fibers traverse the porous network\nand remain mechanically engaged across pores, forming load-bearing paths analogous\nto trabecular struts in bone.\nSecond, the architecture enables the formation of a physically entangled, non-\ncovalent interfacial interphase. Unlike conventional fiber composites that rely on\nchemical coupling or rigid interphases, the trabecular architecture is deliberately\ncombined with a thermodynamically driven interfacial design. During thermal\nprocessing, flexible TPU chains interdiffuse into the semicrystalline POK matrix and\nalong the fiber surface, lowering interfacial enthalpy and promoting physical\nentanglement rather than covalent bonding. This mechanism is supported by DSC,\nXRD, and FTIR analyses (Fig. S5), which show reduced crystallinity and\nconformational rearrangement without the formation of new chemical bonds, and by\nmolecular dynamics simulations (Fig. S6), which visualize the emergence of a dense,\ninterpenetrated interfacial network.\ninterface coupling under load. Newly added in situ tensile SEM and strain-resolved\nRaman spectroscopy (Fig. 5, Figs. S13–S17) directly demonstrate how the trabecular\ninterlock architecture and entangled interface operate together during deformation. At\nmoderate strain, crack growth is shielded by fiber-anchored polymer ligaments and\ncrack deflection within the trabecular network. At higher strain, these ligaments\nundergo progressive stretching and delayed rupture, while the surrounding interface\nreorganizes through chain sliding and rearrangement, enabling sustained energy\ndissipation rather than brittle failure. This behavior is fundamentally different from\nclassical pull-out or interfacial debonding mechanisms.\nFinally, the physical mechanism explains the observed macroscopic performance\nplateau. Because stress is redistributed across the trabecular network and dissipated\nthrough adaptive interfacial rearrangement, TICs avoid the typical strength–toughness\ntrade-off. This leads to convergence toward a Pareto-optimal performance plateau,\nexperimentally validated through preference-guided testing (Fig. 4) and Ashby-type\nbenchmarking. The clustered high-performance region is therefore a direct\nmanifestation of the underlying architecture-enabled interfacial toughening mechanism,\nrather than an empirical coincidence.\nIn summary, the trabecular interlock architecture toughens the interface by\ngeometrically stabilizing crack growth and mechanically activating a\nthermodynamically entangled interphase, which together enable multiscale stress\nredistribution and energy dissipation. To address this comment explicitly, we have\n75\nsubstantially strengthened the mechanistic analysis by adding a dedicated subsection\nentitled “Dynamic stress-adaptive interfacial response”. We integrated in situ tensile\nscanning electron microscopy (SEM) and strain-resolved in situ Raman spectroscopy\nas central mechanistic evidence (Fig. 5, Figs. S13–S17 and Video S3), providing direct,\nreal-time insight into crack evolution, interfacial engagement, and molecular-scale\nstress redistribution.\nThe manuscript has been revised as follows: “\nMechanism analysis\nDynamic stress-adaptive interfacial response\nThe tight clustering of validated properties and the plateau-like distribution on\nAshby maps suggest that the optimized designs share a common, stable dissipation\nmechanism. We therefore interrogate the dynamic interfacial response of a\nrepresentative TIC formulation (balanced preference) using in situ tensile SEM and\nstrain-resolved Raman spectroscopy, complemented by impact/ballistic observations.\nThis formulation sustains tensile strains up to ~600% without catastrophic failure (Fig.\nS13). In drop-weight impact tests, high-speed imaging captures rapid deformation\nfollowed by rebound (Fig. S14a, Video S1), and post-impact inspection shows only\nshallow indentation without macroscopic cracking. By contrast, control specimens\nlacking the trabecular-interlocked architecture and adaptive interface design (POK–BF\nNAFCs) exhibit immediate penetration under identical conditions (Fig. S14b). Under\nballistic loading (80 m s-1), TICs decelerate the projectile and rebound without\nperforation (Fig. S15a), whereas NAFCs show complete penetration (Fig. S15b, Video\nresolving how stress is redistributed and dissipated at the crack and interfacial scales.\nIn situ tensile SEM directly visualizes the crack evolution and load-transfer\nmechanisms operating in TICs. At an applied strain of ~300%, crack opening is\neffectively shielded by pronounced fiber-mediated bridging, where multiple polymer\nligaments anchored on basalt fibers span the crack faces and remain mechanically\nengaged (Fig. 5a, yellow arrows). Instead of propagating along a straight path, the crack\nfollows a tortuous and deflected trajectory, indicating crack deflection induced by\ninterfacial heterogeneity and mechanical interlocking. Upon further loading to ~500%\nstrain (Fig. 5b and Video S3), the bridging mechanism remains operative. The load-\nbearing ligaments undergo progressive stretching and thinning with delayed rupture,\nrather than abrupt interfacial debonding or clean fiber pull-out. Crack advance remains\ngradual and continues along a deflected path, accompanied by a pronounced expansion\nof the surrounding plastic zone (Fig. 5b, blue arrows). These observations demonstrate\nthat the TICs interface does not fail through brittle separation. Instead, it sustains\ncontinuous mechanical engagement across extreme strains through physically\nentangled, load-bearing ligaments that redistribute stress in the crack-tip region, thereby\nstabilizing crack growth and enabling sustained energy dissipation.\nWe next quantify molecular-scale stress accommodation using strain-resolved\ntensile Raman spectroscopy. The stacked spectra (Fig. 5c) reveal a continuous evolution\nof the full Raman window with increasing strain. Consistently, the strain–wavenumber\nheatmap (Fig. 5d) highlights two distinct response regimes: the carbonyl window\n(1700–1785 cm-1) exhibits strain-dependent but spatially diffuse intensity redistribution,\nwhereas the CH-stretch window (2800–3100 cm-1) progressively develops increased\nspectral dispersion at higher strain, indicative of enhanced conformational\nheterogeneity. To quantify the response of polar groups, we track the centroid position\nof the carbonyl band (Fig. 5e). In TICs, the centroid shift remains relatively small and\nnon-monotonic across the applied strain range, indicating that tensile stress does not\ncontinuously accumulate on a narrow subset of carbonyl bonds. Instead, stress is\ndynamically repartitioned as the physically entangled interface undergoes progressive\nreconfiguration under load. In contrast, the POK–TPU control (Fig. S17) exhibits a\nmore systematic and monotonic centroid evolution, consistent with a more\nhomogeneous polymer response in which bond-level stress localization can be directly\ntracked spectroscopically. Chain-level dissipation is further captured by the full width\nat half maximum (FWHM) of the CH-stretch envelope (Fig. 5f). TICs display a\npronounced broadening of the CH envelope at high strain, consistent with the activation\nof segmental mobility, chain rearrangement, and sliding. By comparison, the control\nsample shows a substantially weaker broadening response, indicating more limited\nconformational freedom under tensile loading.\nThese molecular-scale signatures are fully consistent with the in situ tensile SEM\nobservations of sustained fiber-bridging ligaments and an expanded plastic zone at the\ncrack. Together, the spectroscopic and microstructural evidence demonstrates that the\ntrabecular interlocking architecture enables a dynamically activated, stress-adaptive\ninterfacial response, which progressively converts tensile loading into frictional and\nentanglement-mediated dissipation pathways. This multiscale stress redistribution\nmechanism underpins the enhanced crack resistance and impact tolerance observed in\nTICs.\nFigure 5 Dynamic stress-adaptive interfacial response revealed by in situ tensile SEM and\nRaman spectroscopy. a, b. In situ tensile SEM images of crack evolution in TICs at ~300% and\n~500% tensile strain, respectively. At ~300% strain (a), polymer ligaments anchored on basalt fibers\nbridge the crack faces (yellow arrows), accompanied by crack deflection. At ~500% strain (b), the\nbridging ligaments remain engaged while an expanded plastic zone develops ahead of the crack tip\n(blue arrows), leading to gradual crack advance. c, d. Strain-resolved Raman spectra and\ncorresponding strain–wavenumber heatmap acquired during tensile loading, showing distributed\nspectral changes in the carbonyl region (1700–1785 cm-1) and progressive dispersion in the CH-\nstretch region (2800–3100 cm-1). e, f. Quantitative Raman descriptors showing a small, non-\nmonotonic evolution of the carbonyl-band centroid and a pronounced high-strain increase in the\nCH-stretch envelope full width at half maximum (FWHM) for TICs.\nIn situ tensile deformation is performed inside a scanning electron microscope\nusing a custom micro-mechanical tensile stage equipped with a 4000 N load cell. The\ntensile stage is mounted directly inside the SEM chamber (Fig. S17a), enabling real-\ntime mechanical loading under vacuum conditions. Specimens are clamped between\nthe tensile grips and are loaded under displacement-controlled mode at a constant cross-\nhead speed of 10 µm s-1. SEM imaging is conducted continuously during loading,\nallowing direct observation of crack initiation, crack opening, and subsequent crack\npropagation at the same region of interest (Fig. S17b). All in situ tensile SEM\nexperiments were carried out under identical loading and imaging conditions to ensure\nreproducibility.\nFigure S16 In situ tensile SEM experimental setup. a. Photograph of the custom micro-\nmechanical tensile stage integrated into the SEM chamber, showing the motor-driven loading unit\nand specimen gripping configuration. b. Representative SEM view under vacuum with the specimen\nmounted between the tensile grips, enabling continuous imaging of the same region during tensile\nloading.\nSupplementary Note 5. In situ tensile Raman spectroscopy\nIn situ tensile Raman measurements are performed using a confocal Raman\nmicroscope equipped with a 633 nm excitation laser, a 600 lines mm-1 grating, and a 5×\noperated under displacement-controlled mode at a constant stretching rate of 300 µm s-\n1. Specimens are stretched to predefined strain levels and held during spectral\nacquisition following a step-and-hold protocol. Raman spectra are collected at a fixed\nlocation on the specimen surface with an integration time of 10 s and 5 accumulations\nat each strain level. Autofocus is disabled throughout the measurements to avoid focus-\ninduced spectral artefacts, and the focal position is optimized prior to testing. Spectral\nacquisition covers the range of 800–3200 cm-1, including both the fingerprint region\nand the CH-stretch region. All measurements are conducted under identical optical and\nmechanical conditions to ensure consistency and comparability between samples.\nFor data analysis, all spectra are first corrected for background using polynomial\nbaseline subtraction and normalized to the integrated intensity over the full spectral\nwindow to minimize variations arising from laser power fluctuation and collection\nefficiency. The processed spectra are assembled into a two-dimensional strain–\nwavenumber intensity map to visualize strain-dependent spectral evolution. To enhance\ncontrast and reveal subtle redistribution of spectral intensity across strain, a uniform\ngamma correction is applied to the intensity scale for all datasets.\nQuantitative Raman descriptors are extracted from two predefined spectral\nwindows. The carbonyl region is defined between 1700 and 1785 cm-1, within which\nthe spectral centroid is calculated at each strain level as the intensity-weighted average\nRaman shift. This centroid reflects changes in the local vibrational environment and\nconstraint of carbonyl groups during deformation. The CH-stretch region is defined\nbetween 2800 and 3100 cm-1, and the full width at half maximum (FWHM) of the\ncomposite CH envelope is determined by fitting the intensity profile to a single effective\nband after smoothing. The CH-envelope FWHM serves as a descriptor of\nconformational dispersion and segmental mobility. The strain-dependent carbonyl\ncentroid shifts and CH-envelope FWHM values are plotted as a function of tensile strain\nto generate the quantitative curves shown in Fig. 5e and Fig. 5f. Identical spectral\nwindows, fitting procedures, and processing parameters are applied to both TICs and\nthe POK–TPU control to ensure direct comparability.\nFigure S17 Strain-resolved in situ Raman spectra of POK–TPU sample. a. Offset Raman\nspectra of POK–TPU sample acquired at successive tensile strain levels under step-and-hold loading\nconditions. The color scale indicates the applied tensile strain. b. Corresponding Raman intensity\nheatmap as a function of strain, highlighting the carbonyl (1700–1785 cm-1) and CH-stretch (2800–\n3100 cm-1) regions. Intensities are normalized to facilitate comparison across strain levels.”\nare they fundamentally different from traditional chemically constrained interfaces?\nPlease elaborate on the underlying adaptive mechanism.\nResponse: We thank the Reviewer for this important question regarding the nature\nof the “self-adaptive” interfaces. In the revised manuscript, we clarify that the\nadaptivity of the interfaces in trabecular interlocked composites (TICs) is mechanical\nand thermodynamic in origin, and is fundamentally distinct from conventional\nchemically constrained (covalently bonded or compatibilized) interfaces.\nFirst, the defining distinction lies in the interfacial constraint mechanism.\nTraditional chemically constrained interfaces rely on permanent covalent bonds, strong\nspecific interactions, or compatibilizers to enforce adhesion. While such interfaces can\nprovide high initial stiffness, they impose fixed topological constraints and typically\nfail through abrupt debonding or brittle fracture once the local stress exceeds the bond\nstrength. In contrast, the interfaces in TICs are physically entangled rather than\nchemically fixed, formed through thermally activated chain interdiffusion and\ntopological interlocking between TPU, POK, and the fiber surface. No new covalent\nbonds are introduced, as confirmed by FTIR and XRD analyses (Fig. S5), indicating\nthat the interface remains non-covalent and dynamically reconfigurable.\nSecond, the “self-adaptive” behavior arises from a strain-dependent interfacial\ntransition. As demonstrated by newly added in situ tensile SEM, strain-resolved Raman\nspectroscopy, and SAXS/WAXS measurements (Figs. 5–6 and Video S3), the\n83\nentangled interface exhibits a continuous functional shift under loading. At low strain,\nthe interfacial network remains highly connected, enabling efficient load transfer and\ncrack-tip shielding through fiber-mediated bridging. With increasing strain, partial\nchain disentanglement and rearrangement occur, allowing stress redistribution and\nenergy dissipation via chain sliding, segmental mobility, and interfacial expansion. This\ntransition is gradual and reversible, rather than catastrophic, which is the hallmark of\nthe adaptive response. Please also refer to the answer in reply to Comment 1.\nThird, the adaptivity is thermodynamically driven rather than kinetically locked.\nThe interface is designed to reduce interfacial enthalpy (ΔH) during processing,\nlowering the energetic barrier for chain mobility and entanglement formation under the\nGibbs free energy criterion (ΔG = ΔH − TΔS). As a result, the interface remains capable\nof reorganization when subjected to mechanical work. Molecular dynamics simulations\n(Fig. S6) directly visualize this behavior, showing progressive interdiffusion and strain-\ndependent evolution of entanglement density rather than bond rupture.\nFinally, the macroscopic consequence of this self-adaptive interface is sustained\nmechanical integrity across extreme deformation. Unlike chemically constrained\ninterfaces that typically exhibit strength–toughness trade-offs, TICs maintain\nsimultaneous high strength, fracture toughness, and impact resistance. This is\nexperimentally validated by the convergence of preference-guided designs into a\nnarrow Pareto-optimal performance plateau (Fig. 4) and by the absence of abrupt failure\nin both quasi-static and high-rate impact tests (Figs. S14–S15).\nIn summary, the interfaces are termed self-adaptive because their load-bearing role\n84\nevolves continuously with deformation, transitioning from stress transfer to energy\ndissipation through thermodynamically enabled, physically entangled chain dynamics.\nThis adaptive mechanism is fundamentally different from traditional chemically\nconstrained interfaces and is now explicitly supported by multiscale experimental and\nsimulation evidence added in the revised manuscript.\nComment 3: The correlation between the chemical composition/molecular structure of\nthe adaptive interfaces and their interfacial interactions/dynamics needs to be explicitly\nestablished.\nResponse: We thank the Reviewer for requesting a clearer and explicit link\nbetween chemical composition, molecular structure, and interfacial dynamics. In the\nrevised manuscript, we now establish this correlation as a closed, multiscale causal\nchain supported by thermodynamic rationale, direct characterization, molecular\nsimulations, and newly added in situ evidence (Results and Discussion,\nSupplementary Note 2, Figs. 5–6, Video S3, Figs. S5–S6 and Figs. S13–S17).\nFirst, chemical composition defines the interfacial molecular structure. The\nadaptive interface is intentionally engineered in the ternary POK–TPU–BF system,\nwhere semicrystalline POK provides a rigid load-bearing framework, TPU contributes\nhighly mobile soft segments, and basalt fibers create a mechanically interlocking\nscaffold. This compositional combination is designed to promote interfacial chain\ninterdiffusion during thermal processing by reducing the interfacial enthalpy barrier,\nconsistent with the Gibbs free energy criterion (ΔG = ΔH − TΔS) discussed in the\n85\nrevised “Biomimetic trabecular architecture with thermodynamically adaptive\ninterfaces” section and Supplementary Note 2.\nSecond, molecular structure and interfacial interactions are experimentally\nverified to be entanglement-mediated rather than chemically bonded. Thermal and\nstructural characterizations (Fig. S5) show a reduced melting enthalpy and disrupted\ncrystalline order in the ternary composites compared with binary controls, indicating\nenhanced chain intermixing and interphase formation. FTIR reveals no signatures of\nnew covalent bonds, while peak broadening/shifts support non-covalent interfacial\ninteractions governed by chain rearrangement. Molecular dynamics simulations (Figs.\nS6 and 6f) visualize TPU chain interpenetration into POK and along the BF surface,\naccompanied by an increase in interfacial entanglement density, confirming the\nformation of a physically interpenetrated and mechanically engaged interphase.\nThird, interfacial dynamics under load are directly captured by newly added in situ\nexperiments. We added in situ tensile SEM and strain-resolved Raman spectroscopy\n(Fig. 5 and Figs. S13–S17), together with strain-dependent SAXS/WAXS signatures\n(Fig. 6), to explicitly demonstrate that the entangled interface undergoes strain-\ndependent reconfiguration. In situ SEM shows persistent fiber-bridging ligaments and\ncrack deflection without brittle decohesion, indicating sustained mechanical\nengagement. Raman metrics capture the transition from load transfer to dissipation: the\ncarbonyl response remains small and non-monotonic (distributed bond-level stress\nwithout progressive localization), while CH-stretch broadening at high strain indicates\nactivated segmental mobility and chain rearrangement. This dynamic response is absent\n86\nor substantially weaker in control composites lacking the designed interphase,\nconfirming its compositional and molecular origin.\nFinally, macroscopic consequences follow from this composition-driven\ninterfacial dynamics. The shared adaptive interfacial mechanism explains why the\nexperimentally validated preference-guided formulations converge into a high-\nperformance Pareto plateau (Fig. 4) and cluster tightly in Ashby-type benchmarking\n(Fig. 4j–l), reflecting a common deformation and energy-dissipation regime rather than\nexperimental scatter.\nTaken together, the revised manuscript now provides an explicit, experimentally\ngrounded linkage from composition → molecular interphase structure →\nentanglement-mediated interfacial interactions → strain-adaptive dynamics →\nmacroscopic strength/toughness/impact performance, which directly addresses the\nReviewer’s concern.\nThe revised contents are provided in the Results and Discussion (“Mechanism\nanalysis”), Supplementary Note 4 and 5, and the newly added in situ mechanistic\nevidence (Figs. 5–6, Figs. S5–S6, S13–S17, and Video S3).\nComment 4: The necessity of using machine learning (ML) is not sufficiently justified.\nTypically, theoretical calculations provide target values for strength and toughness,\nwhich are often not met in practice. Please provide direct evidence that ML-based\nscreening of molecular structures, material systems, and processing techniques can\neffectively bridge this gap and achieve properties closer to the theoretical predictions.\nthe necessity of machine learning in the present study. We would like to clarify that ML\nis not introduced as a replacement for physical theory, nor as a direct predictor of\ntheoretical limits, but rather as a practical tool for navigating a strongly coupled,\nexperimentally constrained, multi-objective design space that cannot be efficiently\nexplored by theory or intuition alone.\nFirst, although theoretical models can estimate upper bounds for individual\nproperties such as strength or toughness, they do not prescribe realizable compositions\nand processing windows that simultaneously satisfy multiple, often competing,\nobjectives under real manufacturing constraints. In the present system, strength,\nfracture toughness, and impact energy dissipation are strongly coupled through\ncomposition, interfacial structure, and processing-induced morphology. The feasible\ndesign space is therefore highly nonconvex and sparsely sampled by experiments,\nmaking exhaustive experimental exploration impractical.\nSecond, the role of ML in this work is to learn an experimentally grounded\nsurrogate of the structure–property landscape from real data and to explicitly resolve\ntrade-offs among competing objectives. By combining Gaussian process regression,\nPareto set learning, and active learning, the ML framework identifies Pareto-consistent\nregions that are both high-performing and experimentally accessible, rather than merely\nextrapolating toward theoretical optima. This is directly evidenced by the rapid\nhypervolume expansion, accelerated convergence, and superior optimization efficiency\ncompared with a classical response surface method under identical experimental\n88\nbudgets (Fig. S12).\nThird, the ML-guided optimization does not merely select extreme points for\nsingle properties. Instead, it reveals a structured Pareto manifold, within which multiple\npreference-guided solutions cluster into a high-performance plateau. Experimental\nvalidation (Fig. 4) confirms that these solutions achieve simultaneously high strength\n(~250 MPa), high fracture toughness (~15 MPa·m1/2), and high impact energy\ndissipation (~4.7 J), which a combination that is rarely realized in polymer composites\nusing conventional trial-and-error or theory-driven approaches. Importantly, these\nresults are obtained using industrially compatible materials and processing routes,\nunderscoring their practical relevance.\nFinally, while the universality across different polymer matrices and\nreinforcements (Fig. 7 and Fig. S21) arises from the interfacial and architectural design\nprinciples, ML plays an indispensable role in discovering, quantifying, and validating\nthese principles in a representative system. Establishing such a mechanistically\nadmissible and experimentally verifiable design pathway would be infeasible through\ntheoretical calculations alone, which lack direct access to processing-induced disorder,\ninterfacial entanglement, and kinetic constraints, or through brute-force\nexperimentation, which is prohibitively costly in a multi-dimensional design space.\nIn summary, ML is essential in this work not because it predicts theoretical limits,\nbut because it bridges the gap between physical mechanisms and experimentally\nrealizable, multi-objective optimal materials, enabling efficient discovery of design\nrules that can subsequently be generalized and validated across material systems.\ndesign and surface micro/nano-structures of the adaptive interfaces are translated into\nmachine-readable formats (e.g., descriptors, feature vectors) for the ML model.\nResponse: We thank the Reviewer for raising this important point. We clarify that,\nin the present work, machine learning does not directly encode atomistic structures or\nmicro/nano-morphologies as explicit descriptors. Instead, the molecular and interfacial\ndesign principles are implicitly embedded through experimentally controllable\nformulation variables, which act as physically grounded proxies for interfacial structure\nand dynamics.\nSpecifically, the ML model operates at the formulation level, where the input\nfeature vector is constructed from the experimentally defined composition variables,\nnamely the weight fractions of thermoplastic polyurethane (TPU) and basalt fiber (BF),\nwith the poly(octenylene ketone) (POK) content determined implicitly by mass\nconservation. These formulation variables uniquely determine the molecular mobility\ncontrast, entanglement density, and interfacial topology that emerge during thermal\nprocessing.\nThe physical rationale for this representation is as follows. The molecular structure\ndesign of the adaptive interface is governed by the relative fraction of soft, highly\nmobile TPU chains and rigid, semicrystalline POK chains, together with the presence\nand spatial distribution of BF. Increasing TPU content systematically enhances chain\nmobility, interdiffusion length scales, and entanglement probability at the interface,\n90\nwhile BF content controls the density of mechanically interlocking sites and the\ntrabecular load-transfer framework. These effects are experimentally verified by DSC,\nXRD, FTIR, SAXS/WAXS, in situ SEM, in situ Raman spectroscopy, and molecular\ndynamics simulations (Figs. 5–6, Figs. S5–S6 and Figs. S13–S18).\nTo enable the ML model to capture nonlinear coupling between these composition-\ncontrolled interfacial mechanisms and macroscopic mechanical performance, the\nformulation variables are expanded using a polynomial feature construction up to third\norder, including linear, quadratic, cubic, and interaction terms (Supplementary Note\n3.1). This feature expansion provides a compact, interpretable representation of\ncurvature and coupling effects (e.g., TPU–BF synergy) without introducing ad hoc\nstructural descriptors.\nImportantly, micro- and nano-structural features such as interfacial entanglement\ndensity, ligament formation, crack-bridging morphology, and stress-adaptive\nrearrangement are not assumed a priori. Instead, they emerge naturally from the\nformulation variables through processing, and their influence on performance is learned\nindirectly by the ML model via experimentally measured mechanical responses. In this\nsense, the ML framework learns a coarse-grained surrogate of the structure–property\nrelationship, where formulation parameters serve as machine-readable encodings of\ncomplex interfacial physics.\nWe emphasize that this design choice is intentional. Direct encoding of high-\ndimensional microstructural descriptors would require extensive in situ or tomographic\ndata for every sample, which is experimentally impractical and would severely limit\n91\ndataset size. By contrast, the formulation-based descriptor strategy enables statistically\nrobust learning under data-scarce conditions, while retaining a clear physical\ninterpretation validated by independent mechanistic experiments.\nTo avoid ambiguity, we have revised the manuscript and Supplementary\nInformation to explicitly state that:\n1. ML inputs correspond to formulation-level descriptors rather than explicit\nmicrostructural images or atomistic features;\n2. Molecular and interfacial structures are encoded implicitly through\ncomposition-controlled thermodynamic and kinetic effects;\n3. The physical meaning of these descriptors is validated independently through\nin situ multiscale experiments and simulations.\nThe relevant clarifications are provided in Methods (GPR modeling section),\nSupplementary Note 3.1, and are mechanistically supported by Figs. 5–6, Figs. S5–\nS6, and Figs. S13–S18.\nComment 6: The sample size used for the machine learning training must be sufficient.\nPlease ensure the dataset adequately considers the variations from different molecular\nstructure designs and micro/nano-structural parameters.\nResponse: We thank the Reviewer for raising the important issue of dataset\nsufficiency. We clarify that the adequacy of the dataset in this work is ensured not by\nsheer sample count, but by a deliberate combination of controlled formulation space\ncoverage, uncertainty-aware learning, and extensive experimental validation of\n92\nemergent micro/nano-structural variations.\nFirst, the machine learning models are trained on experimentally measured\nformulations, starting from an initial design-of-experiments (DoE) dataset of 50\ncompositions that uniformly span the admissible formulation space defined by TPU and\nBF contents, with POK constrained by mass conservation. This initial dataset already\ncaptures a broad range of molecular mobility contrasts, interfacial entanglement\npropensities, and reinforcement densities, which are the dominant drivers of interfacial\nstructure and dynamics in this system.\nSecond, rather than relying on a single static dataset, we employ a combination of\nPareto set learning (PSL) and active learning (AL) to adaptively expand the dataset in\ndirections that are maximally informative for the multi-objective performance\nlandscape. Over five AL iterations, an additional 25 formulations are selected based on\npredictive uncertainty and Pareto relevance, yielding a final dataset of 75\nexperimentally validated samples. This sequential design ensures that sampling is\npreferentially directed toward regions associated with strong nonlinearity and high\nperformance, as well as regions indicative of emerging structural changes, rather than\nredundantly populating low-information areas of the design space.\nThird, dataset sufficiency is quantitatively assessed through repeated cross-\nvalidation, learning-curve analysis, and independent hold-out testing. As shown in Fig.\nS7 and Fig. S10, predictive accuracy improves monotonically with increasing training\nsize, accompanied by systematic reductions in MAE, MSE, and variance across\nvalidation splits. Independent test validation (Fig. S11) further confirms the absence of\n93\nsystematic bias and demonstrates robust out-of-sample generalization. These results\ncollectively indicate that the models approach a data-saturated regime within the\nexplored formulation space.\nFourth, while micro- and nano-structural parameters are not used as explicit ML\ninputs, their variation is inherently encoded through formulation-controlled\nthermodynamic and kinetic effects. This is experimentally confirmed by extensive\ncharacterization, including in situ tensile SEM, strain-resolved Raman spectroscopy,\nSAXS/WAXS, XRD, DSC, and molecular dynamics simulations (Figs. 5–6, Figs. S5–\nS6, Figs. S13–S18). These results demonstrate that variations in composition\nsystematically translate into distinct interfacial morphologies, entanglement states, and\ndynamic responses under load, which are then reflected in the measured mechanical\nproperties learned by the ML model.\nFinally, we emphasize that the scope of the ML model is intentionally restricted to\na well-defined material system and processing window, within which the dataset is\ndemonstrably sufficient. The purpose of ML here is not to enumerate all possible\nmicrostructures, but to efficiently resolve the coupled composition–interface–property\nlandscape and identify Pareto-optimal regimes under realistic experimental constraints.\nTo clarify these points, we have revised the manuscript to explicitly describe the\ndataset construction, active learning protocol, and validation procedures in the Methods\nand Supplementary Information, with supporting evidence provided in Figs. S7, Figs.\nS10–S12, and Supplementary Notes 3–5.\nReviewer #1 (Remarks to the Author):\nThe manuscript has been revised well. The current draft could be accepted.\nResponse: We thank the Reviewer for the careful evaluation of the revised\nmanuscript and for confirming that all comments and requested revisions have been\naddressed appropriately and sufficiently. We greatly appreciate the Reviewer’s\nconstructive feedback throughout the review process.\nReviewer #2 (Remarks to the Author):\nAll comments and requested revisions have been addressed appropriately and\nsufficiently.\nResponse: We sincerely thank the Reviewer for the positive assessment of our\nrevised manuscript. We appreciate the recognition that the revisions have adequately\naddressed the previous concerns, and we are grateful for the Reviewer’s support for\nacceptance of the current version.\nReviewer #3 (Remarks to the Author):\nThe authors have carefully revised the manuscript in response to the reviewers'\ncomments and have also added more data samples for machine learning. The\nmanuscript is now acceptable, with only one minor suggestion remaining: the molecular\nstructure of the stress-adaptive interface and its formation mechanism need to be\nclarified.\nReviewer #3 (Remarks on code availability):\nThe authors have carefully revised the manuscript in response to the reviewers'\ncomments and have also added more data samples for machine learning. The\nmanuscript is now acceptable, with only one minor suggestion remaining: the molecular\nstructure of the stress-adaptive interface and its formation mechanism need to be\nclarified.\nResponse: We thank the Reviewer for the positive assessment and for the\nremaining minor suggestion. In response to the Reviewer’s suggestion, we have\nclarified the molecular structure and formation mechanism of the stress-adaptive\ninterface at the end of the “Dynamic stress-adaptive interfacial response” subsection,\nexplicitly defining it as a thermodynamically driven, physically entangled interphase\nwith strain-dependent reconfiguration under load.\nThe revised text as follows:\n“At the molecular level, the stress-adaptive interface originates from\nthermodynamically driven interdiffusion of flexible elastomer chains into the\namorphous regions of the rigid polymer matrix and onto the polar fiber surface during\nthermal processing. This process generates a physically entangled interphase rather than\npermanent chemical bonding. Under mechanical loading, the entangled interphase\nundergoes strain-dependent reconfiguration through chain sliding and rearrangement,\nenabling continuous stress redistribution and adaptive energy dissipation.”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript has been revised well. The current draft could be accepted.", "reviewer_2": "All comments and requested revisions have been addressed appropriately and sufficiently.", "reviewer_3": "The authors have carefully revised the manuscript in response to the reviewers' comments and have also added more data samples for machine learning. The manuscript is now acceptable, with only one minor suggestion remaining: the molecular structure of the stress-adaptive interface and its formation mechanism need to be clarified."}, "author_response": "The authors have carefully revised the manuscript in response to the reviewers' comments and have also added more data samples for machine learning. The manuscript is now acceptable, with only one minor suggestion remaining: the molecular structure of the stress-adaptive interface and its formation mechanism need to be clarified.\n\nResponse: We thank the Reviewer for the positive assessment and for the\nremaining minor suggestion. In response to the Reviewer’s suggestion, we have\nclarified the molecular structure and formation mechanism of the stress-adaptive\ninterface at the end of the “Dynamic stress-adaptive interfacial response” subsection,\nexplicitly defining it as a thermodynamically driven, physically entangled interphase\nwith strain-dependent reconfiguration under load.\nThe revised text as follows:\n“At the molecular level, the stress-adaptive interface originates from\nthermodynamically driven interdiffusion of flexible elastomer chains into the\namorphous regions of the rigid polymer matrix and onto the polar fiber surface during\nthermal processing. This process generates a physically entangled interphase rather than\npermanent chemical bonding. Under mechanical loading, the entangled interphase\nundergoes strain-dependent reconfiguration through chain sliding and rearrangement,\nenabling continuous stress redistribution and adaptive energy dissipation.”"}]

Decision letter