Open Access Article
Seong-Hoon Jang
*ab,
Di Zhang
ac,
Xue Jia
a,
Hung Ba Tran
a,
Linda Zhang
ac,
Ryuhei Sato
d,
Yusuke Hashimoto
c,
Yusuke Ohashi
e,
Toyoto Sato
e,
Kiyoe Konno
f,
Shin-ichi Orimo
*ae and
Hao Li
*a
aAdvanced Institute for Materials Research (WPI-AIMR), Tohoku University, Sendai 980-8577, Japan. E-mail: jang.seonghoon.b4@tohoku.ac.jp; shin-ichi.orimo.a6@tohoku.ac.jp; li.hao.b8@tohoku.ac.jp
bUnprecedented-scale Data Analytics Center, Tohoku University, Sendai 980-8578, Japan
cFrontier Research Institute for Interdisciplinary Sciences (FRIS), Tohoku University, Sendai 980-8577, Japan
dDepartment of Materials Engineering, The University of Tokyo, Tokyo 113-8656, Japan
eInstitute for Materials Research, Tohoku University, Sendai, 980-8577, Japan
fInstitute of Fluid Science, Tohoku University, Sendai, 980-8577, Japan
First published on 25th May 2026
Hydrogen is a promising energy carrier, yet its practical deployment is limited by the lack of storage materials that simultaneously achieve high storage capacity (w) and practical equilibrium pressure at room temperature (Peq,RT). Interstitial metal hydrides offer fast kinetics and favorable thermodynamics (high Peq,RT) but suffer from intrinsically low w. Here, we establish a physically interpretable, data-driven framework to uncover descriptor–property relationships in interstitial hydrides using a curated database of pressure-composition-temperature measurements (Digital Hydrogen Platform, DigHyd) and white-box symbolic regression. Strikingly, the analysis reveals a clear separation of governing mechanisms, in which w is governed by geometric and lattice conditions, captured by the average atomic radius (〈rM〉) and average thermal conductivity (〈κ〉), with an optimal regime of 〈rM〉 ∼ 1.47
Å and relatively low 〈κ〉. In contrast, Peq,RT is governed by elastic properties, captured by the average shear modulus (〈G〉) and average Poisson's ratio (〈ν〉), reflecting the role of lattice rigidity and mechanical compliance. These relationships are translated into compositional optimization pathways that follow the descriptor trends above, enabling the design of candidate materials with enhanced w under practical equilibrium conditions (Peq,RT ∼ 0.1 MPa). This work establishes a general, interpretable strategy for physics-informed design of energy materials systems.
Despite decades of extensive investigation, the compositional landscape of hydride-forming alloys remains far from fully explored. While a vast number of binary and multicomponent systems are theoretically accessible, only a limited fraction has been experimentally synthesized and evaluated. This challenge is further exacerbated by the absence of predictive frameworks that are both quantitatively accurate and physically interpretable, which hinders rational materials design. Although recent machine learning approaches have shown promise in accelerating property prediction,17–19 they often rely on relatively small or inconsistently curated datasets and employ black-box models that provide limited insight into the underlying physicochemical mechanisms. To overcome these limitations, we previously developed the Digital Hydrogen Platform (DigHyd: https://www.dighyd.org), a curated database constructed through large-scale extraction of experimental pressure-composition-temperature (PCT) data from the literature.20,21 Building on this dataset, symbolic regression was performed using a white-box modeling approach, enabling the construction of explicit relationships between materials descriptors and key hydrogen storage metrics, namely gravimetric capacity (w) and equilibrium pressure at room temperature (Peq,RT).22,23 This approach identified a compact set of physically meaningful descriptors, including contributions from atomic mass, electronic structure, and packing characteristics, which govern hydrogen storage behavior. In particular, systems containing light elements were found to favor higher hydrogen capacity, whereas electronic and structural descriptors play a dominant role in determining equilibrium pressure. Within this descriptor space, beryllium (Be)-containing alloys emerged as promising candidates; however, their practical applicability is severely limited by toxicity and associated handling constraints.24
These considerations motivate a shift toward more practical materials systems, as schematically illustrated in Fig. 1. Fig. 1a shows the construction of a curated database of PCT measurements, DigHyd, which serves as the foundation for this study. Building on this dataset, Fig. 1b presents the development of interpretable, white-box symbolic regression models that link materials descriptors to key hydrogen storage properties, w and Peq,RT. Focusing on interstitial hydrides, we aim to identify the key descriptors governing both w and Peq,RT by extracting physically meaningful descriptors from the symbolic models (see Fig. 1c). Finally, as shown in Fig. 1d, these descriptor–property relationships are translated into materials design guidelines that enable enhanced w under practical operating conditions. In particular, we target equilibrium pressures near ambient conditions, Peq,RT ∼ 0.1 MPa, corresponding in this study to the window of −1.5 log10[MPa] < log10Peq,RT < −0.5 log10[MPa].
Fig. 2a shows the distribution of data points in the two-dimensional w − Peq,RT materials map for cases where multi-temperature PCT data are available. Most structural classes with different parent structures, such as Laves (C14), LaNi5, LaMgNi4, and TiFe, are located in the low-w region (w < 2.5%), whereas the BCC class extends into the higher-capacity region (2 < w < 5%). However, BCC metals and alloys generally exhibit two distinct plateaus in their PCT curves, with the lower plateau typically located far below atmospheric pressure at room temperature, thereby limiting the practical accessibility of their full w.25 Fig. 2b presents ndata for each reported class. Notably, three classes, Laves (C14), LaNi5, and BCC, out of 14 classes account for ndata = 198, corresponding to 66.2% of the dataset with multi-temperature PCT measurements (ndata = 299), which is close to the Pareto principle.
For symbolic regression modeling, features were extracted as candidate descriptors for w and Peq,RT. The full list of features and their denotations is provided in Table 1; hereafter, the denotations are omitted for simplicity. A total of 57 features were constructed, comprising 18 chemical, physical, and structural features applied to their averaged values ⋯, standard deviations σ(⋯), and skewness r(⋯), along with three additional compositional features (RTr, RTr(IV), and RTr+RE). For example, for a compound AaBbCc, the averaged atomic mass is defined as M = (aMa + bMb + cMc)/(a + b + c), where Mx denotes the atomic mass of element x. Structural descriptors such as Ω, Ωσ, and Vp require the definition of coordination polyhedra XYn (X, Y = A, B, and C crystallographic sites; X-centered). The maximum cutoff distance for metal–metal pairs was set to 3.5 Å, corresponding to the shoulder of the first peak in the radial distribution function averaged across crystal structures (see the section “Averaged Radial Distribution Function of Metal-alloy Parent Structures and the Construction of XYn Polyhedra” in the SI). In addition, the heatmap of Pearson correlation coefficients (rcol) among all feature pairs is provided in the section “Pearson Correlation Heatmap” in the SI.
| Properties | Description | Unit |
|---|---|---|
| χ–χH | Difference in electronegativity between metal atom and hydrogen | — |
| M | Atomic mass | g mol−1 |
| nve | The number of valence electrons. For s/p-block and d/f-block, electrons in the shells of nps and npp, and in those of nps and (np − 1)d are counted, respectively; np is the highest principal quantum number | — |
| rM | Metallic radius | Å |
| ρ | Density | g cm−3 |
| ρmol | Molar density, defined as ρ/M | mol cm−3 |
| ηfM | Metallic filling rate (per unit volume), defined as where Navo is the Avogadro's constant |
— |
| B | Bulk modulus | GPa |
| G | Shear modulus | GPa |
| ν | Poisson's ratio | — |
| κ | Thermal conductivity | W m−1 K−1 |
| α | Thermal expansion coefficient (linear not volumetric) | K−1 |
| θD | Debye temperature | K |
| χm | Molar susceptibility | m3 mol−1 |
| nc | Coordination number to neighbor metal atoms within 3.5 Å for each crystallographic site (A, B, and C). For example, in the compound LaMgNi4, A, B, and C correspond to La, Mg, and Ni, respectively | — |
| Ω | Solid angle per face constructed in polyhedra of XYn (X, Y = A, B, and C crystallographic sites; X-centered; and X–Y distance shorter than 3.5 Å). Estimated from the parent structures listed in Fig. 2b | Radian |
| Ωσ | Standard deviation across Ω for each crystallographic site (A, B, and C). Estimated from the parent structures listed in Fig. 2b | Radian |
| Vp | Volume of polyhedra of XYn for each crystallographic site (A, B, and C). Estimated from the parent structures listed in Fig. 2b | Å3 |
| ⋯ | Average over the constituent metal ions | Same with the unit of ⋯ |
| σ(⋯) | Standard deviation over the constituent metal ions | Same with the unit of ⋯ |
| r(⋯) | Skewness over the constituent metal ions | — |
| RTr | Compositional fraction of transition metal elements | — |
| RTr(IV) | Compositional fraction of the fourth-row transition metal elements (Sc, …, Zn) | — |
| RTr+RE | Compositional fraction of transition and rare-earth metal elements (La, …, Lu) | — |
In Fig. 3b, the PDP alone suggests that smaller 〈rM〉, which is often associated with lower atomic mass, favors higher w. However, in real materials, excessively small 〈rM〉 can instead induce steric constraints, requiring lattice expansion to accommodate interstitial hydrogens. As a result, an optimal value emerges around 〈rM〉 ∼ 1.47 Å, whereas larger values lead to expanded interstitial sites that are unable to effectively stabilize hydrogen, likely due to geometric mismatch between the interstitial cage size and the effective size required for hydrogen accommodation. Notably, in the BCC structure, 〈rM〉 = 1.47 Å yields a geometric maximum hard-sphere radius of 0.43 Å for interstitial hydrogen in tetrahedral sites. In Fig. 3c, the PDP suggests that lower thermal conductivity 〈κ〉, which is linked to electronic structure near the Fermi level and the strength of metallic bonding, and typically associated with softer lattice structures, correlates with higher w. Here, 〈κ〉 is interpreted not as an isolated causal factor, but as an effective descriptor reflecting coupled electronic-structure and lattice-related characteristics of the host metal framework. In particular, softer lattice environments can more readily accommodate lattice expansion upon hydrogen insertion, thereby facilitating higher hydrogen uptake. In practice, however, κ exhibits a more nuanced influence. In contrast, for 〈κ〉 > 110 W m−1 K−1, w shows little variation, indicating that the effect of κ becomes saturated in this regime. Taken together, these results indicate that w is maximized when 〈rM〉 is tuned to approximately 1.47 Å in conjunction with relatively low 〈κ〉.
Fig. 3d demonstrates the strong predictive performance of the model for log10Peq,RT, achieving R2 = 0.889, RMSE = 0.527 log10[MPa], and MAE = 0.400 log10[MPa] over the entire dataset (ndata = 299). Fig. 3e and f present the PDPs for two most important descriptors for log10Peq,RT. In Fig. 3e, both the PDP and experimental and regressed data points show good agreement. A higher 〈G〉 indicates a more rigid lattice, which increases the elastic energy penalty associated with interstitial hydrogen insertion and thereby destabilizes hydride formation, leading to higher Peq,RT. In Fig. 3f, the PDP suggests that the influence of 〈ν〉 on Peq,RT is relatively weak. In real materials, however, Peq,RT decreases rapidly when ν is below approximately 0.3 and becomes nearly constant above this threshold. This behavior suggests that materials with higher 〈ν〉 are more mechanically compliant and can more readily accommodate lattice deformation during hydride formation, consistent with the negative correlation between 〈G〉 and 〈v〉 shown in the section “Pearson Correlation Heatmap” in the SI. Taken together, these results indicate that lattice rigidity plays a central role in governing equilibrium pressure, with stiffer crystal structures leading to higher Peq,RT. To further assess the robustness of the model across different structural classes, a class-resolved sensitivity analysis was performed (see the section “Class-resolved Model Performance and Sensitivity Analysis” in the SI), confirming that the predictive performance is not dominated by any single class.
![]() | ||
| Fig. 4 Materials design pathways toward high hydrogen capacity under practical equilibrium pressure, mapped onto the 〈κ〉-〈rM〉 plot. Optimization trajectories for representative parent structure classes: (a) BCC,26 (b) Laves (C14),27 (c) LaMgNi4,28 (d) La2MgNi9 (PuNi3),29 (e) TiFe,30 (f) LaNi5,31 (g) Laves (C15),32 and (h) Ce2Ni7.33 Each panel shows the compositional optimization pathway guided by the symbolic models within the target window −1.5 log10[MPa] < log10Peq,RT < −0.5 log10[MPa]. The marker size and color, as well as the line color at each step, are scaled by w. Initial and optimized compositions are represented by descriptor-performance vectors p = (〈rM〉/Å, 〈κ〉/(W m−1 K−1), w/%), together with intermediate substitution steps. The red rectangule in each panel indicates a shared design target identified by symbolic models, where 〈rM〉 converges to approximately 1.47 Å with reduced 〈κ〉. In most cases, 〈κ〉 either decreases or remains nearly constant, consistent with adherence to the design rule of 〈rM〉 → ∼1.47 Å. The parent crystal structures are also visually represented. The color gradient is used as a visual guide to indicate the overall direction of materials design trends, rather than representing a quantitative colormap. The reliability of the optimized candidate compositions was further evaluated using ensemble uncertainty metrics (see the section “Uncertainty Assessment of Optimized Candidate Materials” in the SI), indicating that the proposed candidates lie within a low-uncertainty prediction regime. These compositions should be interpreted as physically guided design candidates rather than directly synthesizable materials, and further experimental validation will be required. | ||
As an illustrative example, Fig. 4a shows the optimization of the BCC-type alloy TiVNbCrNi2, initially characterized by a descriptor-performance vector p = (〈rM〉/Å, 〈κ〉/(W m−1 K−1), w/%) = (1.42, 63.8, 3.5).26 The optimization proceeds through a multi-step compositional pathway: first, 90% of Ni is replaced with V to yield TiVNbCrV1.8Ni0.2; next, 10% of Cr is substituted with Tm, resulting in TiVNbCr0.9Tm0.1V1.8Ni0.2; finally, 10% of Ti is replaced with Nb, leading to the composition Ti0.9Nb0.1VNbCr0.9Tm0.1V1.8Ni0.2. This sequence yields the optimized composition Ti0.9Nb1.1V2.8Cr0.9Tm0.1Ni0.2, with final p = (1.49, 45.1, 7.57). Similar symbolic-model-guided optimization systematically shifts compositions toward enhanced w across other structure types27–33 by modulating 〈rM〉 and 〈κ〉 (Fig. 4b–h). Here, all descriptors are updated self-consistently with composition and are not fixed to their average values, in contrast to the PDP analysis in Fig. 3.
Across all optimized cases, two consistent design principles emerge in most cases, represented by red rectangules in Fig. 4a–h. First, 〈rM〉 converges toward an optimal value of approximately 1.47 Å, represented by dashed magenta lines. Second, 〈κ〉 either decreases or remains nearly unchanged during optimization. These trends directly reflect the descriptor–property relationships identified in Fig. 3b and c, demonstrating that the symbolic models capture transferable design rules for achieving high w under practical Peq,RT conditions. Although the predicted values of w themselves may be subject to uncertainty due to their extrapolative nature, the proposed optimization pathways provide physically grounded and clear design directions, based on the white-box symbolic models. Also, this design pathway can be further enhanced through integration with a closed-loop discovery framework,34,35 in which symbolic models are iteratively refined using synthesis and measurement data guided by the optimization trajectories, thereby enabling more accurate, real-world-relevant predictions.
In contrast, Peq,RT is primarily governed by elastic properties, described by 〈G〉 and 〈ν〉. A higher 〈G〉 increases the elastic penalty for hydrogen insertion, leading to higher equilibrium pressure, whereas larger ν reflects greater mechanical compliance and lowers Peq,RT. Thus, Peq,RT is fundamentally controlled by the elastic response of the host lattice. These descriptors provide complementary representations of lattice rigidity and deformability rather than independent contributions. This also can be understood in the context of the van't Hoff equation,
, where Peq,T is the equilibrium pressure at temperature T, and ΔH and ΔS are the enthalpy and entropy changes associated with hydride formation. For many metal-hydrogen systems, ΔS is largely dominated by the loss of gas-phase entropy when molecular hydrogen is incorporated as atomic hydrogen into the lattice. Consequently, variations in Peq,RT are primarily governed by changes in ΔH. In this context, 〈G〉 and 〈ν〉 reflecting the elastic energy penalty is consistent with a prior study demonstrating that lattice strain has been shown to modify the effective ΔH and shift Peq,T over several orders of magnitude by altering the elastic energy associated with hydrogen insertion.37 Furthermore, 〈G〉 has recently been reported as a key descriptor governing ΔH in the formation of BCC high-entropy alloy monohydrides, based on a machine learning study.38
Thus, we emphasize that these relationships are consistent mainly with established physical understanding, and the present analysis provides a quantitative framework to support and organize these insights. Here, it should be noted that the present descriptor framework is most applicable to systems undergoing isostructural hydrogenation. In materials with non-isostructural phase transitions, additional thermodynamic and structural factors may become dominant.
Taken together, these results establish a clear separation of roles between descriptors governing capacity and thermodynamics. While w is controlled by geometric-electronic factors, Peq,RT is dictated by elastic constraints. This separation provides a physically transparent and actionable framework for materials design under practical operating conditions. Importantly, these descriptor relationships are consistently observed across different structure types, indicating that they define a transferable design space. Compositions across multiple classes converge toward 〈rM〉 ∼ 1.47 Å with reduced 〈κ〉, demonstrating the robustness of the identified design principles. The present framework is expected to be most reliable within the probed compositional and structural domains of interstitial hydrides. Application beyond this space may require incorporation of additional descriptors to capture different bonding regimes, phase transformations, or thermodynamic contributions not present in interstitial hydrides.
These findings should also be understood in the broader context of descriptor discovery across hydrogen storage materials. In our previous study, we showed that across a global materials space, including both interstitial and saline hydrides, w and Peq,RT are strongly governed by elemental electronegativity, and can therefore impose an intrinsic trade-off between capacity and equilibrium pressure.19 This indicates that, at the global scale, electronegativity acts as a unifying descriptor controlling hydrogen storage behavior across fundamentally different bonding types. In contrast, the present work focuses on the local materials space of interstitial hydrides, where more specific descriptors emerge, reflecting finer structural and mechanical effects. Thus, while global electronegativity-driven trends may introduce trade-offs, the local descriptor framework identified here provides a route to navigate and potentially mitigate these trade-offs for simultaneous optimization of high w and practical Peq,RT. Together, these results highlight a “hierarchical” descriptor framework, in which global trends are governed by simple chemical parameters, while local optimization requires more detailed, physically specific descriptors.
More broadly, the present results demonstrate that symbolic regression can resolve coupled physical mechanisms into interpretable descriptor relationships, enabling a transition from empirical exploration toward mechanism-informed materials design. This approach provides a generalizable strategy for descriptor-based materials design in hydrogen-related systems. Furthermore, such a framework is expected to be directly applicable to saline hydrides and, more broadly, to systems beyond hydrogen storage, including hydride-based solid-state electrolytes for batteries.39
000 curated data entries on hydrogen storage properties and thermodynamic parameters. Importantly, the database spans a wide range of material classes beyond conventional metal hydrides, including interstitial hydrides, complex hydrides, ionic (saline) hydrides, multi-component and destabilized systems, as well as porous materials such as metal–organic frameworks. This breadth enables systematic analysis across fundamentally different hydrogen storage mechanisms. The curated data include gravimetric hydrogen capacity (w), as well as the enthalpy (ΔH) and entropy (ΔS) changes associated with hydrogenation reactions, primarily defined as
.
:
2. Given the set of candidate descriptors listed in Table 1, the number of possible model combinations including interaction terms is given as 1.95 × 10463 for the 109 default scalar transforms provided in the Regressor module.
Model construction was therefore performed in a staged manner. Starting from models with 20 active independent variables, the number of active variables was iteratively reduced one by one (20 → 19 → 18 → ⋯ → 1), while progressively allowing more complex interaction terms (i.e., increasing interaction depth). At each stage, candidate models were evaluated, and the model with the highest R2 over the entire dataset was selected.
All models were generated under the so-called full Fisher condition, whereby the p-values of both the overall model (F-test) and all individual coefficients (t-tests) were required to be less than 0.05. For each target metric (log10(w/M) and log10Peq,RT), 10 symbolic models were constructed and subsequently combined into a single “stacking-ensembled” model. Key descriptors were then identified by examining features that consistently appeared across the 10 symbolic models, as well as by evaluating the magnitude of their contributions through z-scored coefficients. In particular, descriptors associated with terms exhibiting large average z-scored coefficient magnitudes across the ensemble were considered to play dominant roles in determining the target properties.
To ensure model reliability, the standard deviation across the regressed values in the ensemble was constrained to be below 0.5 log10[% mol g−1] and 0.7 log10[MPa] for log10(w/M) and log10Peq,RT, respectively. Predictions exceeding these thresholds were discarded. In the 5-fold benchmark tests, the discard ratios were approximately 1% and 10% for log10(w/M) and log10Peq,RT, respectively. Notably, when models were constructed using the entire dataset, no data points were discarded, as all predictions satisfied the ensemble consistency criteria. Representative examples of symbolic model ensembles, their error distributions, and comparisons across different target metrics (log10w and log10(wM)) are provided in the section “Symbolic Models: Formulation and Error Analysis” in the SI.
Compositional modifications were performed through controlled substitution operations while preserving the original stoichiometry. Specifically, full substitution, 90% substitution, and 10% substitution of constituent metal elements were allowed. Candidate substitution elements were required to satisfy two criteria relative to the original element: the differences in electronegativity and metallic radius must not exceed 0.5 and 0.5 Å, respectively.
Supplementary information (SI): Radial distribution function averaged over metal-alloy parent structures for hydrogen storage and the construction of XYn polyhedra, Pearson correlation heatmap, benchmarking of regression models, class-resolved model performance and sensitivity analysis, uncertainty assessment of optimized candidate materials, and symbolic models: formulation and error analysis. See DOI: https://doi.org/10.1039/d6sc03089k.
| This journal is © The Royal Society of Chemistry 2026 |