Predicting Thermal Decomposition Temperature of Binary Imidazolium Ionic Liquid Mixtures from Molecular Structures.

Ionic liquids (ILs) have been regarded as "designer solvents" because of their satisfactory physicochemical properties. The 5% onset decomposition temperature (T d,5%onset) is one of the most conservative but reliable indicators for characterizing the possible fire hazard of engineered ILs. This study is devoted to develop a quantitative structure-property relationship model for predicting the T d,5%onset of binary imidazolium IL mixtures. Both in silico design and data analysis descriptors and norm index were employed to encode the structural characteristics of binary IL mixtures. The subset of optimal descriptors was screened by combining the genetic algorithm with the multiple linear regression method. The resulting optimal prediction model was a four-variable multiple linear equation, with the average absolute error (AAE) for the external test set being 12.673 K. The results of rigorous model validations also demonstrated satisfactory model robustness and predictivity. The present study would provide a new reliable approach for predicting the thermal stability of binary IL mixtures.

Introduction

By definition, ionic liquids (ILs) are novel green solvents composed of organic cations and organic or inorganic anions, which are liquid at room temperature. They are labeled as “designer solvents” because of their excellent physical and chemical properties such as high conductivity, low melting point, remarkable thermal stability, and good solubility.[1−3] ILs have achieved a rapid development and are widely applied in various fields such as catalysis,[4] carbon capture,[5] electrolytes in batteries,[6] and pharmaceutical processing.[7]

However, it is still difficult to find the ideal IL that satisfies the deeper demand of desired properties. Since each new IL requires a complete property analysis and registration, it can be costly and time-consuming to study and develop a new IL. A feasible and alternative solution is the IL mixtures. IL mixtures, which are also known as double salt ILs (DSILs),[8] make it possible to further fine-tune their properties. DSILs often exhibit different physicochemical properties from pure ILs since new ionic associations generate. By changing the combinations of ILs, the desired IL mixtures can be designed and synthesized for certain properties. A great deal of research has been published on the applications of IL mixtures in multiple fields, e.g., electrochemistry,[9,10] carbon dioxide capture,[11,12] catalysis,[13] and extraction.[14]

Since many industrial processes generate substantial heat or need to be carried out under high-temperature conditions, high requirements are put forward for the thermal stability of ILs. Although the ILs are almost nonflammable, there still exists thermal decomposition. During high-temperature industrial processes, the ILs may change their structures and thermally decompose, which leads to unwanted byproducts and unexpected accidents.[15] Therefore, the thermal decomposition temperature (Td) is commonly referred to characterize the thermal hazards of ILs.[16−19] Furthermore, compared to the traditional onset decomposition temperature (Td,onset), the temperature at 5% decomposition (Td,5%onset) is considered as a more realistic and conservative indicator to determine the maximum permissible working temperature during the operational process.[20,21]

Extensive research on Td,5%onset of IL mixtures has been performed.[12,22−24] All of them tested the thermal stability by experimental studies, which requires a lot of manpower and material resources. Furthermore, for those samples with poisonous and harmful substances, there remains difficulties and hazards in measurement. Currently, there is a trend that the quantitative structure–property relationship (QSPR) method is widely employed to predict the physicochemical properties of chemicals. Compared with the experimental method, the QSPR is apparently a more convenient and efficient approach.

In the past few years, considerable efforts have been made to predict the Td of pure ILs and IL composites. Lazzús established a theoretical prediction model for the Td of 198 ILs based on a group contribution method.[25] The data set contains 27 cationic groups and 31 anionic groups, and the resulting model provided an alternative approach to predict the Td of ILs. Venkatraman and Alsberg employed partial least squares and random forests (RF) to develop models for predicting Td of 995 diverse ILs that consisted of 461 cations and 119 anions.[26] Electronic, thermodynamic, and geometrical descriptors calculated from semiempirical PM6 were selected, and the RF model showed better predictive performances. Zhao et al. established a QSPR model to describe Td values for a total of 168 ILs.[27] They used the DRAGON program (version 6) to calculate descriptors and employed genetic algorithm with the multiple linear regression (GA-MLR) to filter the most relevant variables. Both the structure for single ions and the cation–anion interactions were considered, and the results showed satisfactory predictability. Zeeshan et al. combined 29 diverse imidazolium ILs with two diverse metal–organic frameworks (MOFs) and proposed QSPR models based on the MLR method to predict Td values.[28] The results demonstrated satisfactory prediction of Td values of IL/MOF composites, and they found that ILs with the [NTf2]− anion possess higher thermal stability in the composites. Recently, Duan et al. developed the QSPR model for 158 different imidazolium ILs via norm index descriptors.[29] The descriptors they employed not only characterize the structures of anions and cations but also describe the interaction between anions and cations.

However, with more research on IL mixtures, there is still no prediction model for the Td,5%onset of IL mixtures because of the lack of experimental data. In this study, the data set was selected from previous publications, aiming to develop a simple and reliable model to predict the decomposition temperatures of imidazolium IL mixtures via the QSPR method. In this study, descriptors for both structural characteristics of ions and cation–anion interactions were calculated and combined together to represent the characteristics of the selected imidazolium IL mixtures. In particular, in silico design and data analysis (ISIDA), which represents the molecular structure based on a finite number of topological fragments,[30,31] was employed to explore the relationship between target properties and ionic structures. Besides, the norm index descriptor was utilized to describe the interaction between cations and anions.[29,32,33] The QSPR modeling was carried out in this study for a total of 31 different imidazolium IL mixtures, and various rigorous model validation strategies were then performed to ensure the model performance. Through our work, it was expected to provide an efficient approach to evaluate the thermal stability of IL mixtures and provide some guidance for experimental design and mixing of ILs.

Results and Discussion

Results of Prediction

The most optimal subset was selected after the process of GA-MLR. In this study, three structural descriptors with “augmented atoms” type “IIAB” (including a cation descriptor and two anion descriptors) and an interaction descriptor were obtained to encode the characteristics of imidazolium IL mixtures. To analyze the contribution of substructures, the representations of these selected descriptors are listed in Table . The optimal MLR model is shown as follows:where R2 represents the squared correlation coefficient, Q2LOO represents the value of leave-one-out (LOO) cross-validation (CV), SE represents the standard error of the model, n represents the number of mixtures in the training set, and F represents the Fischer F-ratio.

Table 1

Descriptors Selected for the Present Model for the Prediction of IL Mixturesa

structure descriptor	representation	ionic type	mixing rule	interaction descriptor	representation
	C(−N′)	cation			normin1(CM2,1,2)
	F(−C′)	anion		/	/
	F(−C′)	anion		/	/

Here, N, C, and F represent the kinds of atoms and “–” represent the single bond, with those outside of brackets denoting augmented atoms and those in brackets denoting their neighbor atoms and bonds

Then, the predictive ability of the model was tested by comparing the predicted and observed values, which is shown in Figure . The predicted values are presented in the Supporting Information (Table S1). The plot clearly depicts that the majority of data points are located in the diagonal line, which indicates a good correlation between the predicted and observed values of IL mixtures. Table summarizes the main statistical parameters of the MLR model. The result shows that the Td,5%onset values of IL mixtures can be effectively predicted by the model equation with an AAE of 7.294 K. In addition, the predicted percentage error of 31 samples is shown in Figure . The distribution showed the specific number of samples in each error interval. The obtained minimum relative error is 0.14%, the maximum relative error is 5.91%, and the average relative error is 1.93%.

Figure 1

Plot of the predicted vs observed values of the Td,5%onset of binary IL mixtures.

Table 2

Main Statistical Parameters of the Obtained MLR Model

statistical parameters	training set	test set
R2	0.969	0.923
Q2LOO	0.949
Q2EXT		0.924
RMSE	9.463	14.152
AAE	7.294 K	12.673 K
n	24	7

Figure 2

Percent errors obtained by the model and the number of mixtures in each range.

Model Stability Validation and Result Analysis

For the present model, the Y-randomization test was employed as a validation approach to determine model stability. It was carried out 100 times for the training set samples in this paper. The visual summary of the Y-randomization test is performed by demonstrating the relationship between the R2 values of randomized models and the model frequency (Figure ). The obtained maximum, minimum, and average R2 values of the randomized model are 0.465, 0.0153, and 0.168, respectively, while the value of standard deviation (SD) is 0.110. Since the difference of R2 between the original model and the mean highest random is much greater than 3SD, the result indicates a negligible chance correlation between the randomly shuffled models and the original model, which emphasizes the unique connection between the developed model and the selected variable set.

Figure 3

Histogram of R2 of randomization vs frequency of occurrence of the randomized models.

Furthermore, the predicted residuals and observed values are analyzed in Figure . Evidently from the diagram, all predicted residuals are distributed on both sides of the zero baseline uniformly and randomly, which indicates that no systematic errors occurred during the establishment of the MLR model.

Figure 4

Plots of the residual vs the observed Td,5%onset values of binary IL mixtures.

With the rigorous validation of multiple strategies, it can be concluded that the obtained model shows excellent robustness and predictability. The QSPR model derived from molecular structural descriptors can efficiently and conveniently predict the Td,5%onset values of IL mixtures.

Applicability Domain of the Present Model

The domain application of the QSPR model is defined as a chemical structure space, which is characterized by the properties of the compounds in the training set. The Williams plot is employed to analyze the domain application and is shown in Figure . The applicability domain (AD) is depicted as a squared area within ±3 standard deviations and a leverage threshold h* of 0.625.

Figure 5

Williams plot describing the AD of the QSPR model for the Td,5%onset of the binary IL mixtures (h* = 0.625).

As can be easily observed in Figure , all predicted values of IL mixtures are inside the area. Since all IL mixtures are within the residual range, the developed model is considered reliable and can satisfactorily predict Td,5%onset of imidazolium IL mixtures within the application ranges.

Conclusions

In this work, the QSPR approach was employed to determine the relationship between the thermal stability and molecular structures of binary imidazolium IL mixtures. After carrying out internal and external validations, the results showed that the predictive model has been successfully developed with satisfactory model performance. To sum up, this work has the following main findings:

The relationship between the thermal stability of binary IL imidazolium mixtures and their molecular structures was investigated in this study, for the first time, and the corresponding prediction model was proposed to predict the Td,5%onset of the IL mixtures.

The norm index descriptors combined with ISIDA descriptors were employed to extensively and effectively characterize both molecular structure features and interactions of binary IL mixtures, which lead to the resulting QSPR model with satisfactory prediction performance.

The developed model was conceptually simple, convenient for application, and would also be reasonably expected to reliably predict the decomposition temperature of binary imidazolium IL mixtures with satisfactory robustness and predictivity, from only their molecular structures. Additionally, this study could also provide some guidance for synthesizing or designing IL mixtures for safety purposes.

Materials and Methods

Data Set

The observed 5% onset decomposition temperatures of binary imidazolium IL mixtures were selected from previous publications.[12,23] The data set consists of 31 data points containing 4 pure imidazolium ILs and 27 mixtures of binary imidazolium ILs. The Td,5%onset values of these mixtures ranged from 402 to 597 K. The details of 31 binary IL mixtures are listed in Table S2 of the Supporting Information.

Calculation of Molecular Descriptors

Calculation of Descriptors for Ions

The ISIDA fragment descriptor, which is also realized as the substructure molecular fragment (SMF) descriptor, is defined as the number of topological fragments of the molecular structure. There are two main classes of molecular subgraphs: (I) “sequences” and (II) “augmented atoms.” The “sequences” means the shortest path of two atoms connected by chemical bonds. There are three types of “sequences,” which are performed as atom types (A), bond types (B), and atom and bond types (AB). For each type of sequence, the minimal and maximal lengths of the path are defined as nmin and nmax, respectively. In this study, the length of the sequence corresponds to nmin = 2 and nmax = 8. The selected atoms with their environment are defined as “augmented atoms,” which also includes neighbor atom and bond types (AB), atom types (A), or bond types (B).

Since ILs are composed of anions and cations, it is necessary to calculate descriptors for cations and anions separately. The detailed descriptor calculation process is performed as follows: first, the 2D chemical structures of each cation and anion were drawn using MarvinSketch (version 15.6.29.0) and optimized based on the clean in 2D method. Structure data files (SDF) including two cations and three anions were then, respectively, imported into EdiSDF, a program of the ISIDA system, to generate SDF records of cations and anions. The described SMF descriptors were generated by the SMF module of the ISIDA system.

The mixture descriptor for the IL mixture was developed by combining the descriptors of diverse cations and anions separately. The 12 proposed formulas that take into account their respective mole fractions were used to calculate mixed-cation descriptors and mixed-anion descriptors[34] (Table ). Then, descriptors for constituted species were concatenated to describe the structure of the multiple-ion system.[35]

Table 3

List of Mixture Formulasa

no.	formula
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)

D is the descriptor of mixture; x1 and x2 are the molar fractions of components 1 and 2; respectively; d1 and d2 are ISIDA descriptor values for components 1 and 2, respectively; is the absolute value of the difference between x1 and x2; and is the absolute value of the difference between d1 and d2.

Calculation of Interaction Descriptors

Since ILs consist of cations and anions, their thermal stability and the character of the cations and anions are intimately connected. Thus, not only the structural characteristics of ions should be taken into account, but also the contribution of interaction is indispensable. The molecular structure can be intuitively represented by the SMF module, while the interactions between the cations and anions cannot be represented because of the limitation of the software. Considering the diverse associations that possibly exist in a multiple-ion system of IL mixtures, the descriptors based on the norm index were proposed to describe each cation–anion interaction of imidazolium IL mixtures in this work.

First, the 3D structure of each cation and anion was achieved by Chemdraw (version 14), with the optimization based on the MM2 module (the class 1 Allinger molecular mechanics program). For further optimization, the Gaussian (version GaussView 6.0.16) was employed to carry out density functional theory (DFT) M06-2X functional calculation on the basis of 6-311+G (d,p). Then, the optimized molecular structures were characterized by the matrix norm index.[33]

The calculation process of distance matrices and property matrices is shown as follows:

n is the path between atoms i and j.

is the Euclidean spatial distance between atoms i and j.where AW, Ra, EN, MW, AQ, and BD are atom weight, atom radius, electronegativity, molecular weight, atom charge, and branching degree, respectively.

Then, the distance matrices and the property matrices are combined as CM matrices, which are calculated as follows:where m = 1, 2, 3, 4, and 5, which correspond to matrices MP1 to MP5. Furthermore, the three norm indexes are defined in eqs –14:where p and q represent the number of rows and columns of matrix CM, respectively. is the eigenvalue of the matrix, and the is defined as the Hermite matrix of the matrix.

Considering that the thermal stability is related to the various interactions of the binary IL system, the descriptors for multiple cation–anion interactions are further proposed based on the original calculation for single ILs,[29,33] which are defined as follows:

In this study, all the four interaction descriptors , , , and were proposed to represent the cation–anion interactions individually, and then all those descriptors were concatenated to describe the characteristics of multiple-interaction system, where w = 1, 2, and 3, corresponding to norm(CM,1), norm(CM,2), and norm(CM,3), and CF and AF refer to the mass fractions of the cations and anions, respectively. The , and , refer to the norm indexes of matrices calculated from the cations and anions, individually.

Descriptor Selection and Model Development

Based on the above programs, descriptors for IL mixtures were obtained. Determining the optimal descriptors that are most relevant to the target attribute is a critical step in QSPR modeling. The GA has been proved to be an efficient method applied for variable selection. In this paper, the optimal subset of descriptors that accurately characterize structural features corresponding to their target properties was selected by combining the GA-MLR method. The selection process of GA-MLR was implemented by Materials Studio (version 8.0).

Model Validation

To ensure the model reliability of resulted QSPR models, model validation can be absolutely necessary. Multiple verification strategies were employed for both internal and external validation in the present work.

The squared correlation coefficient (R2) is commonly used to measure the fit between the observed and predicted values. The model has good fit if the R2 value is greater than 0.6. The root-mean-square error (RMSE), the SE, and the AAE were also presented to evaluate the predictive capability.[36]where y is defined as the observed Td,5%onset value, y0 represents the predicted Td,5%onset value, and n is the number of the IL mixtures.

The CV is considered as an approach of great use to verify the robustness of the resulting model. The internal predictive ability was tested by the CV result (Q2) in internal validation. The most economical LOO CV (Q2LOO) was utilized in this work,[37] which is calculated as follows:where y, y0, and ¯Y represent, respectively, the observed, predicted, and mean observed Td,5%onset values of the IL mixtures in the training set.

To determine the generalizability and the prediction ability of the resulting QSPR model, external validation is indispensable. The data set was divided into a training set and an external test set. The training set is used for descriptor selection and model development, while the test set is for model external validation. There are three main strategies commonly applied for data set partition of mixtures, including “Points out,’, ‘‘Compounds out,’’ and ‘‘Mixtures out.’’ The “Points out” strategy was employed in this work.

In this study, the data set was randomly divided into a training set with 24 samples (77.4% of the data set) and a test set (22.6% of the data set) with seven samples based on the “Points out” strategy. Then, the squared correlation coefficient for external validation (Q2ext) was calculated to indicate the predictive capability of the established QSPR model. The calculation process is as follows:where y and y0 represent the observed and predicted Td,5%onset values of the IL mixtures in the test set, and represents the mean observed Td,5%onset values of the IL mixtures in the training set.

Furthermore, the Y-randomization test was applied to verify the stability of the prediction model. It randomly scrambles the dependent variables and keeps independent variables unchanged. Then, the performance of the presented model is compared with the original model. The operation is repeated 50–100 times through the program. If the randomly shuffled models present much lower R2 values than the original model, it can be considered that no chance correlation exists in the modeling process.

AD

An explicitly defined AD is required for any resulting QSPR model according to OECD principle 3.[38] The AD is a theoretical area defined by the molecular similarity of the training set and determines if the developed model can reliably predict the valid range of new chemicals. The popular Williams plot is applied as a practical tool for the definition of the AD. It is presented as a two-dimensional scatter plot determining whether the chemicals of the developed model are located in or out of the AD.

The structural similarity between the sample and the training set is represented by the leverage value h, which is defined aswhere represents the descriptor row-vector of descriptors for the ith sample, and X represents the matrix of descriptors for all samples in the training set. The warning leverage value h* is defined as follows:where k represents the number of selected descriptors in the developed model, and m represents the number of samples in the training set. If the leverage value (h) of the sample is higher than the warning leverage value (h*), the predicted result of the compound is considered out of the valid range of application and the prediction is considered unreliable.