Veerupaksh
Singla
,
Qiyuan
Zhao
and
Brett M.
Savoie
*
Davidson School of Chemical Engineering, Purdue University, West Lafayette, IN 47906, USA. E-mail: bsavoie@purdue.edu
First published on 1st July 2024
The absence of computational methods to predict stressor-specific degradation susceptibilities represents a significant and costly challenge to the introduction of new materials into applications. Here, a machine-learning framework is developed that predicts stressor-specific stability scores from computationally generated reaction data. The thermal degradation of alkanes was studied as an exemplary system to demonstrate the approach. The half-lives of ∼32k alkanes were simulated under pyrolysis conditions using 59 model reactions. Using a hinge-loss function, these half-life data were used to train machine learning models to predict a scalar representing the relative stability based only on the molecular graph. These models were successful in transferability case studies using distinct training and testing splits to recapitulate known stability trends with respect to the degree of branching and alkane size. Even the simplest models showed excellent performance in these case studies, demonstrating the relative ease with which thermal stability can be learned. The stability score is also shown to be useful in a design study, where it is used as part of the objective function of a genetic algorithm to guide the search for more stable species. This work provides a framework for converting kinetic reaction data into stability scores that provide actionable design information and opens avenues for exploring more complex chemistries and stressors.
Despite progress in developing methods to simulate or empirically predict other functional properties, contemporary stability characterization is overwhelmingly done experimentally via empirical make-and-break testing. Although this can reliably characterize a particular material, it is costly in terms of time and material and crucially only delivers information at the end of the design process. Even when such testing is done, generating mechanistic information requires additional analytical monitoring such as differential scanning calorimetry, thermogravimetric analysis, spectroscopy (XRD, XPS, IR, UV-Vis, and NMR), mass spectrometry, cyclic voltammetry, or high-pressure liquid chromatography, amongst others.22–27 To the extent that stability properties are considered during design, it is typically done only indirectly by limiting the scope of materials to established chemistries with known stability or by accounting for a small number of acute susceptibilities that are domain specific.16,21,27–29 Computational contributions at the design stage are primarily through thermodynamic energy estimates. Although this is a reasonable starting point, it provides no information about the kinetics and susceptibility of a material to degradation via specific stressors.30–34 Famously, diamond is thermodynamically unstable relative to graphite, but its thermal decomposition is not kinetically favorable except at temperatures higher than 2000 K.35
Recent advances in reaction prediction create the possibility that first-principles reaction data might provide a basis for estimating the degradation kinetics of prospective chemistries.36,37 Automatic reaction network exploration and transition state characterization algorithms have become more efficient and accurate.38–41 Although such characterization techniques still remain too expensive to be used in high-throughput applications for virtual screening and comparative down selection of prospective materials, they can potentially serve as the basis for generating kinetic data relevant to broad classes of degradation reactions. Critically, material stability is mainly determined by the kinetics of the first irreversible degradation reaction and thus does not necessarily require a full elaboration of a many-step degradation network. Thus, although data scarcity is no longer a fundamental obstacle, the challenge remains in converting this information into a form that generalizes to similar chemistries to avoid the direct simulation of every new material.
Here, we address the absence of inexpensive kinetic measurements of chemical stability by investigating whether ML can be used to turn simulation based half-life data into a scalar stressor-specific stability score (Fig. 1). The thermal stability of alkanes is taken as an exemplary problem and several model architectures are tested to learn a scalar “stability score” that allows for the comparison of the relative thermal stability of alkanes. Although one of the simplest organic classes, the resulting scores compress information from theoretical reaction networks spanning billions of reactions. Remarkably, the relative stability is shown to be easily learned by all ML architectures, suggesting the feasibility of greatly generalizing this approach to other stressors and broader chemical classes. The utility of the derived stability score for materials design is finally demonstrated by performing a case study using a genetic algorithm to discover thermally stable isomers.
The reaction networks for low temperature (i.e., ≤700 K) alkane pyrolysis consist of five elementary reaction steps (Fig. 1a(i)).42 The initiation reaction involves homolytic C–C bond cleavage, which produces two alkyl radicals (C–H bond cleavage is not preferred at such temperatures). Once an alkyl radical has been generated, it can undergo three types of reactions: one bimolecular reaction, hydrogen abstraction (used interchangeably with H-abstraction in the text), and two unimolecular reactions, isomerization and radical decomposition. The hydrogen abstraction reaction occurs between an alkane and an alkyl radical, where an intermolecular hydrogen radical transfer shifts the alkane and alkyl. In contrast, in the isomerization reaction, an intramolecular hydrogen radical transfer within the alkyl radical generates isomeric alkyl radicals. The transition states of isomerization reactions have ring-like structures. Based on whether the hydrogen is abstracted from the 4th, 5th, 6th, or the 7th nearest neighboring carbon to the radical carbon, isomerization reactions have been named 1–4, 1–5, 1–6, or 1–7 isomerization (or iso, used interchangeably through the text). Alkyl radicals can also decompose into alkenes and smaller alkyl radicals through β scission, which is referred to as the radical decomposition reaction. As a termination step, two alkyl radicals can recombine to form stable alkanes via radical recombination. To prevent unlimited network expansion, only radical recombination reactions that produce alkanes with an equal or fewer number of carbon atoms than the starting alkane are considered in the pyrolysis network generation process. Species-specific activation energies for these reactions were calculated from a tractable number of model reactions (MRs) generated from a graph-truncation approach. The reaction graph truncation consists of conserving the reacting atoms and their first bonded neighboring atoms to define the reaction neighborhood. Applying this truncation to all acyclic alkane pyrolysis reactions results in 59 unique MRs. The activation barriers of these MRs were then used in place of larger reactions of the same class whenever they occur in the network (Fig. 1(ii) and ESI Section 1.2 for further details†). The median free energy errors associated with the model reaction approximation were found to be within 3 kcal mol−1 for the 55 model reactions that were benchmarked.
A pruning procedure was used to assemble effective reaction networks and determine accurate half-lives for each alkane (Fig. 1a). The size of the untruncated pyrolysis network grows exponentially with respect to size, making the direct simulation out to terminal products impractical. For example, the reaction network of n-C15H32 is composed of over a billion reactions. However, the downstream reactions from the parent alkane have little impact on the half-life, which allows these networks to be pruned depth-wise based on their simulated flux, while still accurately calculating the half-life (see ESI Section 1.3†). The network pruning approach was benchmarked against complete reaction networks for all structural isomers of alkanes from ethane to decane (i.e., the 149 systems for which complete simulations out to terminal products were practical). Using a minimum relative flux threshold of 10−9, a robust fit was observed (Fig. 2a) with a reduction of up to 108 times in terms of the number of unique reactions required to include in the reaction network, while still obtaining an accurate half-life (Fig. 2b). The growth of unique reactions vs. alkane length also drops to a linear trend compared to the exponential trend for complete networks. The pruned reaction networks exhibit the largest reduction in H-abstraction reactions among the different mechanisms in the full pyrolysis networks (Fig. S2b†).
The half-lives of 32421 acyclic alkanes, i.e., all structural isomers up to C17H36, were simulated with Cantera,43 based on their depth-wise pruned networks (see ESI Section 1.3†). Although the absolute accuracy of these half-lives is not at issue in the current study, it is still necessary to establish that qualitatively correct stability trends are reflected in the data for later testing of whether various ML architectures can generalize the implicit stability trends. Experimental kinetic data compiled by Sundaram and Froment for the pyrolysis of ethane, propane, n-butane, and isobutane provide several points of comparison (Fig. 3).44 Additionally, these experimental data and pyrolysis mechanisms were consistent across the literature.45,46 These data were used to replace the quantum chemistry based reaction data in the pyrolysis networks of ethane, propane, n-butane, and isobutane (ESI Section 3†) to resimulate the half-life for comparison with the purely computational results (Fig. 3a). For simple alkanes, both the experimental and simulated half-lives recapitulate established trends that increasing the alkane length (comparing ethane, propane, and n-butane) and branching (comparing n-butane and isobutane) reduce half-lives. Comparisons to experimental kinetic data for more complex alkanes were, however, not possible in the absence of systematic experimental literature. Additionally, sensitivity analysis was performed on both the experimental and computational kinetic data with 1000 simulations done for each reaction network, where random uncertainties from a normal distribution centered at zero and a standard deviation of 3 kcal mol−1 were introduced in all activation barriers. The overlapping shaded region in Fig. 3a implies that the experimental and computational half-lives are within this uncertainty range.
Comparing the half-life data with the experimental heats of formation at 298 K (ΔHf), a frequently used stability surrogate, illustrates the qualitative error of neglecting kinetics in assessing stability. A larger exothermic ΔHf corresponds to higher thermodynamic stability, and a lower half-life corresponds to lower kinetic stability. The thermodynamic data compiled by Pedley, Naylor, and Klein (PNK data)47 show that thermodynamic stability increases with the alkane length, while kinetic stability decreases (Fig. 3b). The thermodynamic trend also misses the half-life discontinuity between n-butane and n-pentane. Comparing ΔHf and the half-life of linear (n-alkanes) and branched alkanes reveals that the thermodynamic trend also qualitatively misrepresents the effect of branching (Fig. 2c). While branching marginally increases thermodynamic stability, it decreases kinetic stability and induces a broad range of half-lives depending on the position and degree of branching. Hence, ΔHf lacks the information required to predict degradation kinetics and gets structure-stability trends qualitatively wrong.
Pairwise accuracy was used as a metric to quantify and compare the model predictions. The training and the testing datasets were converted to unique molecule pairs and the pairwise stability trends of the ground truth kinetic data were compared with the pairwise stability trends predicted by models. The accuracy was then measured as the percent of pairs where the model predicts the same stability trend as the ground truth. Several training and testing splits were generated to test the transferability of the stability score under different scenarios (Fig. 4). The consistent accuracy of >90%, even for the relatively simple MLP architecture, illustrates the relative ease of the learning task.
Using a random training and testing split, the elementary MLP and message-passing Chemprop architectures show comparable performance (Fig. 4a(i)). Although random splitting ensures that the testing set is structurally independent of the training data, the diversity of the training data is sufficient to guarantee that the branching and size distributions of the training and testing data obtained from the random splitting are similar. To carry out a more rigorous test of transferability, four other case studies were performed using distinct train and test splits. In the case study definitions, the backbone refers to the longest chain in the alkane structure. Core branches refer to the number of branches originating from the backbone, whereas total branches refer to all branches and sub-branches from the backbone. Total branches ≤6 put alkanes with at most six total branches in the training set and the rest in the testing set, resulting in a 90:10 train:test split (Fig. 4a(ii)). Core branches ≤4 refer to alkanes with at most four core branches being included in the training set, resulting in a 61:39 train:test split (Fig. 4a(iii)). Backbone ≤10 includes alkanes with a backbone at most ten carbons long as the training set and the rest as the testing set, resulting in a 68:32 train:test split (Fig. 4a(iv)). Length ≤ 10 puts alkanes with at most 16 carbons in the training set and the C17 alkanes in the testing set, resulting in a 47:53 train:test split (Fig. 4a(v)).
The >98% accuracy on splits (iv) and (v) demonstrates the ability of the models to extrapolate to longer alkanes. The >90% accuracy on splits (ii) and (iii) shows that while not as good as extrapolating on length, the models perform reasonably well when extrapolating to highly branched alkanes. Chemprop performs better than the Morgan fingerprint-based MLP model in most of the cases, which is reasonable given that a message-passing architecture is more complex than the MLP fingerprints. However, the fact that the MLP model surpasses 90% accuracy in all the cases, is within 1% of Chemprop accuracy in two instances, and even surpasses Chemprop in one case illustrates that predicting relative stability is a reasonably simple learning problem. The learning curves for all models across all the split types are present in ESI Section 4.†
We also compared the n-alkane stability scores predicted by both the models in Fig. 4b to test their ability to recapitulate the jump in half-life observed between n-butane and n-pentane under pyrolysis conditions. Although the absolute values of the two models cannot be compared, they both predict the stability jump between n-butane and n-pentane (Fig. 3b). This jump, which was unseen during training, demonstrates that the models are learning the relative stability associated with different chemical features rather than just memorizing the training data. These results lend credence to our initial assumption that given consistently generated kinetic data for the stressor-specific degradation reactions, it is possible to learn the inherent relative kinetics, which opens up the avenue to apply this strategy to more complex systems.
Three types of operations were applied as part of the genetic search for stable structures: growth, deletion, and insertion of methyl groups (Fig. 5a). Crossovers were not considered because for alkanes they are essentially the same as a couple of generations of mutations. Growth involves adding a methyl group in place of a hydrogen on the alkane, deletion involves removing a primary or secondary carbon atom from the alkane and then generating a smaller alkane by maintaining hybridization using the removal or addition of hydrogens, and insertion involves adding carbon between a C–C bond and balancing the charge with hydrogens.
The stability scores were used to guide genetic searches to find the most stable C17H36 isomer out of the 24894 possible structures. The genetic searches were initialized with 100 randomly selected C17H36 alkanes selected from the least stable 5000 C17H36 structures, and every generation maintained a population of 100 with the constraint that all molecules in every generation should have 17 C atoms. Elitist selection was used, wherein 20% of the most stable predicted molecules from each generation were propagated to the next generation, and the 50% most stable predicted molecules were used as parents for further mutations to generate children. The remaining 80 molecules to be propagated to the next generation were randomly selected from this group of children using probabilities weighted by the predicted stability scores.
For both the models (MLP and Chemprop), 20 independent genetic runs were conducted and from each run the most stable predicted molecules from every generation were obtained. The median stability of the generated molecules predicted to be the most stable was plotted along with the zeroth and fourth quartile stability scores. The limited variation illustrates that irrespective of starting alkanes, all runs quickly arrive at similar structures (Fig. 5b). Within 30 generations, both models converge with respect to discovering alkanes that maximize their respective scores. The MLP and Chemprop models both show a general trend of predicting less-branched alkanes to be more stable. The searches also consistently converge at late stages on species capped with tert-butyl and isopropyl moieties over species capped with unsubstituted CH3 units. Indeed, a comparison of stability scores for linear alkanes of a given C-value with and without tert-butyl capping shows that those with this termination are slightly preferred as termini, although non-terminal branches in general are disfavored. There are insufficient experimental data from the pyrolysis literature to validate this directly, but it stands as a prediction of the case study.
There are several avenues for further improving this work. The current work focuses on the thermal stability of alkanes as an entry point for a broader range of degradation susceptibilities that are relevant to materials development. We found that both machine learning architectures were easily able to learn the thermal stability structure–function relationships such that the stability scores could be reliably used in design tasks. Nevertheless, the generalization to other chemistries and stressors is both obvious and necessary. The ease with which the present relationships were learned suggests that the broad data generation performed here so that the scaffolded training:testing case studies could be performed will not be necessary when considering other chemistries and stressors. As such, we do not anticipate fundamental data generation obstacles to extending this approach to redox stability (e.g., by considering reactions involving anions and cations of neutral parents) or other chemistries beyond alkanes. We further showed that the kinetic data contain material stability information and that this relative information can be learned by neural networks and even be extrapolated to unseen data for material design. This means that if properly defined, kinetics can be learned in a way useful to computational material design and virtual screening.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4dd00036f |
This journal is © The Royal Society of Chemistry 2024 |