Development of Prediction Models for the Self-Accelerating Decomposition Temperature of Organic Peroxides

Thermal risk assessment is very important in the primary stages of chemical compound development. In this study, a model to estimate the self-accelerated decomposition temperature of organic peroxides was developed. The structural information of compounds was used to calculate descriptors, on which partial least-squares (PLS) regression and support vector regression were applied for temperature prediction. Molecular mechanics and density functional theory calculations were performed before descriptor calculations, for structure optimization, using a genetic algorithm for variable selection. Structure optimization and variable selection immensely improved the prediction accuracy. Thus, a PLS model, with R2 = 0.95, root mean square error = 5.1 °C, and mean absolute error = 4.0 °C, exhibiting higher accuracy than existing self-accelerating decomposition temperature prediction models, was constructed.


■ INTRODUCTION
Thermal risk assessment is extremely crucial in the development of chemical compounds. Self-accelerating decomposition temperature (SADT) is a key parameter characterizing the thermal risk of organic peroxides. It is the lowest temperature for self-accelerating decomposition in organic peroxides and selfreactive substances (used in transportation packaging). Thus, it determines the optimum temperature-control to avoid thermal hazards during material storage and transport. 1 Several experimental methods measure SADT for thermal risk assessment; however, the associated cost, risk, and chemicals make early-stage evaluation very difficult. Therefore, it is beneficial to develop a simple and high-accuracy SADT prediction method.
Previous studies have proposed quantitative structure− property relationship models to predict SADT. Wang et al. 2 performed density functional theory (DFT) calculations [6-31G(d)/B3LYP] in Gaussian 09 to obtain descriptors. Geometrical descriptors (bond length, bond angle, and dipole moment) and quantum chemical descriptors [highest occupied molecular orbital (HOMO)/ lowest unoccupied MO (LUMO), bond dissociation energy] were used to construct prediction models with multiple linear regression (MLR) and support vector regression (SVR) to estimate SADT. He et al. 3 used the semiempirical molecular orbital technique (AM1) for preprocessing before geometry optimization and frequency calculations. Descriptors (excluding quantum chemical descriptors) were calculated in DRAGON 6.0, and the genetic algorithm (GA) was applied for variable selection, followed by MLR and SVR construction, to estimate SADT. The first and second methods suffered from the limitations of high computational load and low accuracy, respectively.
In this study, a model with high accuracy and low calculation load was developed to estimate the SADT of organic peroxides. Descriptors were calculated using the optimal molecular conformation, determined by molecular mechanics (MM) and DFT calculations. GA was used for variable selection, followed by the application of partial least-squares (PLS) regression, as a linear regression method, and SVR, as a nonlinear regression method, to predict SADT. Prediction models, including and excluding structural optimization and variable selection, were developed and analyzed.
■ METHODS Data Set. The data set included 65 organic compounds, with 90.14−571.00 molecular weights, and −5.0 to 196.5°C SADTs, obtained from the literature, 2−4 determined using different calorimetric methods (TG-DSC and C80). However, as previously reported, SADT is independent of the determination method used. 4 Compounds included commonly used organic peroxides, such as dialkyl peroxide, diacyl peroxide, hydroperoxide, peroxyester, ketone peroxide, peroxy carbonate, and  Table 1. The data set was divided into two subsets: a training set (52 samples) and a test set (13 samples), following the method reported in a previous publication. 3 Geometry Optimization. Molecular structures of the 65 compounds were prepared in the molfile format. Some descriptors were molecular-structure-dependent, making geometry optimization very important. Auto geometry optimization on Avogadro, 5 a molecular modeling software for quantum chemical calculations, was performed for each mol file. The basic MM potential energy function includes bonded terms (for covalently bonded atomic interactions) and nonbonded terms (for long-range electrostatic and van der Waals forces). Here, UFF (universal force field) was used to improve bond lengths and bond angles to obtain a minimum-energy conformation. Prior MM calculations improved DFT convergence, with less computational load. DFT calculations using Firefly (PC GAMESS), 6 via MoCalc2012, 7 optimized the molecular structures. "Geometry optimization" job type, using hybrid density functional B3LYP with 6-31G basis set, was used. B3LYP hybrid functional incorporates approximations to the exchange−correlation energy functional in DFT (combination of exact exchange from Hartree−Fock theory, and from other sources). Polarization-incorporated 6-31G, 6-31G(d), is commonly used for organic compounds. Here, 6-31G was chosen to reduce the computational load and verify the effect of d-orbitals.
Descriptor Calculation and Selection. After structural optimization, 5666 and 552 molecular descriptors of organic peroxides were calculated by alvaDesc 2.0.8 8 and CODESSA 3, 9 respectively, most of which were irrelevant to this study. These methods calculate molecular descriptors and fingerprints from structural information. Descriptors with small standard deviations and strong multicollinearity were eliminated. CODESSA was used to calculate the quantum chemical descriptors, such as HOMO/LUMO, and GA was used to find a descriptor set for model construction.
GA, applied as a variable selection method, is an iterative procedure that continuously improves the fitness function, from which the fitness score (indicating probability of descriptor set selection) is calculated. An initial population of descriptors (a few hundred sets) was selected at random or heuristically. Each iteration step calculated and assigned a fitness value to the descriptor sets, and proportional probabilities were used to select a new descriptor population. This selection procedure cannot independently generate a new point in the search space; thus, crossover and mutation were additionally used by GA to generate new descriptor sets. Here, the coefficient of determination (R 2 ), after fivefold cross validation in PLS and SVR modeling, labeled GA-PLS and GA-SVR, respectively, was used as the fitness function. Number of components in the population was 100, crossover probability was 0.5, mutation probability was 0.2, and maximum number of generations was 200.
Regression and Validation. Here, PLS and SVR were used to develop the prediction models. PLS is a statistical linear regression method used to find fundamental relations between explanatory variables (X) and response variables (Y). It is widely used when the number of explanatory variables is significantly larger than the number of samples. Decompositions of X and Y, as shown below, were constructed to maximize the covariance between the latent variables (T) and response variables (Y).
where A is the number of latent variables, t a is the ath latent variable, p a is the ath loading, q a is the weight on ath latent variable, and E and f are the error terms that cannot be explained by X and Y. p a and q a were calculated to minimize the sum of squares of errors. The number of latent variables with the highest R 2 , obtained via fivefold cross validation, was used in the prediction model. Variable importance in projection (VIP) 10 scores of each descriptor were calculated to identify variables that contributed significantly to the prediction. VIP scores, as shown below, are defined for each X variable and j, as the sum of latent variables of its PLS-weight value (w ij ), weighted by the percentage of explained Y variance.
where h is the number of latent variables, w i is the weight vector, and R 2 (y, t a ) is the percentage of explained Y variance. SVR is also a regression method, performing linear and nonlinear regression using the kernel trick, implicitly mapping inputs into high-dimensional feature spaces. If x (i) is the explanatory variable for the ith sample, the response variable, f(x (i) ), is expressed as follows where b is a constant term, w is the weight vector, and k is the number of dimensions. The error function is expressed as follows where C is the regularization parameter, n is the number of training samples, and ε specifies the epsilon tube within which no penalty is associated with the loss function. The RBF kernel, a kernel function used in machine learning, shown below, was used. x x where γ is the RBF kernel parameter. SVR includes three hyperparameters (C, ε, and γ) to be provided before model construction.
Model development and descriptor selection were performed using Python 3.7. Scikit-learn, 11 a machine learning library in Python, was used for the PLS and SVR calculations. Grid-SearchCV, 12 a parameter estimator in Python, was used for hyperparameter optimization. Fast optimization of hyperparameters was implemented following the procedure adopted by a previous publication. 13 Hyperparameter optimization and descriptor selection were performed only on the training data set and then evaluated the performance of the model for the test data set. Model accuracy was evaluated on the common statistical parameters: root mean square error (RMSE), mean absolute error (MAE), and R 2 . These parameters were calculated as follows Additionally, prediction accuracies of existing models were compared ( Table 2).  Figure 2. The statistical parameters of the model for the training set were as follows: R 2 = 0.99, RMSE = 0.98, and MAE = 0.76 for PLS, and R 2 = 0.99, RMSE = 2.15, and MAE = 0.34 for SVR, and those for the test set were as follows: R 2 = 0.82, RMSE = 9.80, and MAE = 7.70 for PLS, and R 2 = 0.77, RMSE = 11.07, and MAE = 8.68 for SVR. Prediction performance was significantly higher than case 1 and comparable to existing models. Prediction accuracy of both the models could be improved by DFT calculations, and the addition of quantum chemical descriptors, before model building. Changing the number of latent variables in PLS, and hyperparameter values in SVR, did not improve the prediction accuracy, similar to case 1.
Case 3 (Geometry Optimization: MM/DFT, Variable Selection: GA). In this case, in addition to preprocessing (similar to case 2), the descriptors were selected using GA-PLS and GA-SVR, before model building. After variable selection, the fitness function of GA was improved from R 2 = 0.69 to R 2 = 0.91 for PLS and from R 2 = 0.38, to R 2 = 0.90 for SVR. The model was developed by PLS with 11 latent variables and SVR (C = 466, ε = 0.02343, and γ = 0.00000714). A parity plot of actual values versus calculated values is shown in Figure 3. The statistical parameters of the model for the training set were as follows: R 2 = 0.99, RMSE = 1.57, and MAE = 1.14 for PLS, and R 2 = 0.99, RMSE = 3.69, and MAE = 1.78 for SVR, and those for the test set were as follows: R 2 = 0.95, RMSE = 5.11, and MAE = 4.03 for PLS, and R 2 = 0.91, RMSE = 6.87, and MAE = 5.15 for SVR. Prediction performance dramatically improved compared to   case 2, and the prediction accuracy was also better than the existing models. Thus, it can be said that the prediction accuracy could be improved by the appropriate selection of variables before model building.
Comparison of Prediction Accuracy between PLS and SVR. Comparison between PLS and SVR prediction accuracies are listed in Table 3. "Descriptors" line shows the change in the number of descriptors due to preprocessing, and "variable selection" line shows the number of descriptors after applying GA. The prediction accuracy was good overall, with a higher value for PLS than SVR, for the data set used in this study. Geometry optimization using MM/DFT calculations, addition of quantum chemical descriptors, and variable selection using GA influenced prediction accuracy.
Comparison of Prediction Accuracy with Existing Models. The model in the literature 2 had small RMSE and high prediction accuracy; however, the computational load was very high due to DFT calculation. On the other hand, the model in the literature, 3 which used the semiempirical molecular orbital method (AM1), had relatively larger RMSE, although the computational load was low.
The prediction accuracy of the proposed model was significantly improved by the addition of quantum chemical descriptors and variable selection using GA. Changing the basis set for DFT calculation to 6-31G could reduce the computational load while maintaining the prediction accuracy.
The addition of quantum chemical descriptors and appropriate optimization of molecular conformation before calculating the descriptors were important for a model with high accuracy. However, improvement in prediction accuracy reached a ceiling at some point, so improvement in prediction accuracy should be balanced with a computational load to create an effective model (Table 4).
Descriptors with High Impact on Prediction Accuracy. Top 15 descriptors with highest VIP scores are shown in Table  5. Various descriptors related to oxygen bonding (bond order, valence, and charge), and quantum chemical descriptors (LUMO and repulsion/attraction energy), are included in the table. The result shows that these descriptors obviously influenced the prediction accuracy of SADT.
Analyzing Effects of Preprocessing Including and Excluding DFT/GA. In no. 53, 54, 55, 59, 64, and 65, the prediction accuracy was improved by performing the DFT calculations before descriptor calculations, but in no. 57, 58, and 60, preprocessing by DFT calculations did not significantly improve the prediction accuracy. Additionally, no. 56, 62, and 63 showed high prediction accuracy from the beginning. In all the cases, except no. 62, variable selection via GA improved the prediction accuracy, but its overall influence on predictionaccuracy improvement was smaller than the effect of DFT optimization (Figures 4).
No. 53, 54, 56, and 58 were compared as representatives ( Figures 5−8). Conformation (bond length and angle) near the O−O bond of each molecule significantly changed after DFT calculation, as presented in Table 6. No. 53 and 54 exhibited relatively larger bond length changes than no. 56 and 58 due to the greater difference between the initial and optimal states for no. 53 and 54 than between no. 56 and 58. No. 56 exhibited a low prediction error for case 1, despite only 18 iterations, which could be due to the optimal initial conformation, where only MM calculations were performed. A large number of iterations yielded a high prediction accuracy, with some exceptions. Prediction error for no. 58 could not be reduced significantly.  This could be due to the DFT calculation errors because 6-31G, instead of 6-31G(d), was used as the basis function. Thus, appropriate optimization of molecular conformation and addition of quantum chemical descriptors immensely influence SADT prediction (Tables 7−10).
Double Cross-Validation. Here, following a previous study, the data set was divided into training and test data to conduct holdout validation of the prediction accuracy. The temperature range of the test data was 15−70°C, and it was unclear whether the model was applicable to a wider temperature range. To verify the extrapolation (generalization performance) for a high temperature range, using a small data set, double cross-validation was conducted using case 3 data, variables, and model building conditions. In inner crossvalidation, hyperparameters were determined via fivefold cross-validation, whereas in outer cross-validation, leave-one-    (Figure 9).

■ CONCLUSIONS
In this study, we constructed a model to estimate the SADT of organic peroxides. PLS regression and SVR were applied on the descriptors calculated using the structural information of the        compounds, to predict the SADTs. MM and DFT calculations were performed before calculating the descriptors, and GA was used for variable selection. In DFT calculation, B3LYP with the 6-31G basis set was used instead of 6-31G(d), significantly improving prediction accuracy and reducing computational load. Thus, a model with higher accuracy than the existing SADT prediction models was developed. Appropriate preprocessing and variable selection were important for a model with high accuracy, and optimizing compound conformation before descriptor calculations improved SADT prediction accuracy. However, the improvement in prediction accuracy should be balanced with a computational load to create an effective model.
In the future, application of machine learning models other than SVR, preprocessing methods for descriptor calculations, and descriptor selection with high contribution to SADT prediction, to further improve prediction accuracy, will be investigated.