Thermal Conductivity Estimation of Diverse Liquid Aliphatic Oxygen-Containing Organic Compounds Using the Quantitative Structure-Property Relationship Method.

Thermal conductivity is an essential thermodynamic data in chemical engineering applications. Liquid aliphatic oxygen-containing organic compounds are important organic intermediates and raw materials. As a result, estimating thermal conductivity of liquid aliphatic oxygen-containing organic compounds is of significance in industry production. In this study, the genetic function approximation method was applied to screen descriptors and develop a 6-descriptor linear quantitative structure-property relationship model. The entire data set of these compounds covering 1064 thermal conductivity values was divided into 694-member training set, 298-member test set, and 72-member prediction set. The average absolute relative deviation of the training set, test set, and prediction set were 4.14, 4.41, and 4.16%, respectively. Model validation and Y-randomization test proved that the developed model has goodness-of-fit, predictive power, and robustness. In addition, the applicability domain of the developed model was visualized by the Williams plot. This study can provide a convenient method to estimate the thermal conductivity for researchers in chemical engineering production.

Introduction

Thermal conductivity acts as an important physical property that can reflect the ability of heat transfer of a substance. In chemical industry, petroleum industry, and energy engineering, thermal conductivity is one of the essential basic thermodynamic data for heat transfer design.

From the middle of the 18th century, many researchers have made a lot of exploration on thermal conductivity measuring methods.[1−6] However, these methods are time-consuming, costly, and technique limited, and sometimes, there is a certain error in the measured thermal conductivity value. The main reason is that the heat loss caused by convection and radiation cannot be controlled during the measurement. As a result, developing a theoretical method for estimating thermal conductivity is of theoretical significance.

Some estimation methods have been proposed to predict thermal conductivity of aliphatic oxygen-containing organic compounds, which were based on the liquid molecular motion theory model and knowledge of the liquid heat conduction mechanism. However, it should be aware that these methods need be supported by empirical theory and experimental data. As a result, these methods are empirical with big errors, generally between 2 and 10%, as shown in Table .

Table 1

Results of Aliphatic Oxygen-Containing Organic Compound Calculations

	Rodenbush et al.		Nagvekar and Daubert		Baroncini et al.
classes	no. of data	MRDa/%	no. of data	MRD/%	no. of data	MRD/%
alcohols	267	2.7	634	6.3	592	7.7
acids	68	3.2	787	6.5	236	4.1
esters	92	2.8	243	9.7	197	5.7
ketones	48	3.2	68	6.3	72	8.3
ethers	18	3.1	75	7.3	60	3.2
aldehydes	28	2.1	43	8.3	44	6.4
refs	(7)		(8)		(9)

MRD stands for mean relative deviation.

The quantitative structure–property/activity relationship (QSPR/QSAR) was defined as a mathematical relationship connecting the chemical structure with compound properties in a quantitative manner.[10] It can reliably predict the physicochemical, biological, and pharmacological properties of compounds from the molecular structures of compounds. This method became prevalent when it is not possible to obtain accurate values in the experiment limited to economy, time, or technology.[11−13] As for thermal conductivity, few works are reported with QSPR. Results of some works are listed in Table .

Table 2

Retrospect of Works for Thermal Conductivity Prediction

methods	authors	compound class	parameters	N	R2	results	refs
MLR	Gao and Cao	alkanes	4	155	0.9510	sa= 0.0033	(14)
MLR	Kauffman and Jurs	organic solvents	9	213	0.953	RMSEb= 0.0136	(15)
GFAc	Khajeh and Modarress	alcohols	5	116	0.9438	RMSE = 0.0474	(16)
MLR	Liu et al.	alcohols	4	139	0.9738	RMSE = 0.0029	(17)
GFA	Liu et al.	alkyl halides	6	410	0.9745	RMSE = 0.0035	(18)

s stands for the standard deviation.

RMSE stands for root-mean-square error.

GFA stands for genetic function approximation.

Aliphatic alcohols, ethers, aldehydes, ketones, acids, and esters are diverse important chemical products and intermediates, which are widely used as industrial raw materials, lubricants, solvents, medicines, daily necessities, food additives, and so forth. Thermal conductivity is an important thermodynamic property for these organic compounds in industrial application. Therefore, how to accurately predict the thermal conductivity of these diverse aliphatic oxygen-containing organic compounds has important significance for engineering application.

Structures of these compounds have some certain similarities, which satisfies the basic conditions for the QSPR study.[19−21] This work intends to (1) collect thermal conductivity data of aliphatic alcohols, ethers, aldehydes, ketones, acids, and esters; (2) extract molecular descriptors from molecules; (3) develop a QSPR model and validate the model; and (4) predict the thermal conductivity of diverse aliphatic oxygen-containing organic compounds using the model.

Results and Discussion

Using the GFA (100 initial population size, 500 the maximum generations, 10% the mutation probability, and the smooth parameter of LOF α = 0.5), a 6-descriptor linear QSPR model was developed:where SM2_B(s) is the spectral moment of order 2 from the Burden matrix weighted by I-state, SIC0 is the structural information content index (neighborhood symmetry of 0-order), IC1 is information content index (neighborhood symmetry of 1-order), T is the temperature, Eta_F is the eta functionality index, and MATS2m is Moran autocorrelation of lag 2 weighted by mass. The regression statistics of these involved parameters are given in Table .

Table 3

Regression Statistics of Parameters Involved in the QSPR Model

parameters	type	coefficients	standardized coefficients	t	p	VIF
SM2_B(s)	2D matrix-based descriptors	0.0711	1.0443	47.33	0.000	3.88
SIC0	Information indices	0.2960	0.9626	58.79	0.000	2.14
IC1	Information indices	–0.0511	–0.5876	–36.67	0.000	2.04
T	Temperature	–0.0002	–0.5853	–49.16	0.000	1.16
Eta_F	ETA indices	–0.0127	–0.5579	–27.77	0.000	3.22

Model Validation

As can be seen from the absolute values of standardized coefficients in Table , the most significant influence of descriptors on thermal conductivity in turn is SM2_B(s), SIC0, IC1, T, Eta_F, and MATS2m. The variance inflation factor (VIF) values of these descriptors lower than 10 together with the intercorrelation matrix (Table ) manifest the absence of multicollinearities among the descriptors.[22] In addition, the confidence level p values of the descriptors all far less than 0.0001 and F = 1214.35 prove the statistical significances of the descriptors and the developed model, respectively.

Table 4

Correlation Matrix of the Involved Descriptors

	SM2_B(s)	SIC0	IC1	T	Eta_F	MATS2m
SM2_B(s)	1.000
SIC0	–0.366	1.000
IC1	0.114	0.575	1.000
T	0.206	–0.277	–0.250	1.000
Eta_F	0.810	–0.231	0.093	0.106	1.000
MATS2m	0.119	0.193	0.349	–0.168	0.239	1.000

The squared correlation coefficient of training set Rtr2 = 0.9138 and root-mean-square error (RMSE) = 0.0067 indicated that the developed model was acceptable with goodness-of-fit. Using this model to predict the test set and the prediction set, the squared correlation coefficients were 0.8922 and 0.8816, which suggested that the developed model had predictive ability. The correlation between experimental thermal conductivity and calculated values was plotted in Figure (detailed information are given in Supporting Information, Table S2) which shows that the experiment thermal conductivity were highly in agreement with the calculated values.

Figure 1

Experimental thermal conductivity vs calculated values.

In addition, predicted thermal conductivity values of different chemical classes versus temperature are shown in Figure . It can be found that the tendency of predicted thermal conductivity variation with temperature is basically the same with the experimental thermal conductivity variation with temperature. As for the leave-one-out cross validation, QCV2 = 0.9118 and RMSECV = 0.0072 manifested the stability of the developed model.

Figure 2

Predicted thermal conductivity values of different chemical classes vs temperature.

The Y-randomization test was carried out to testify the absence of chance correlation in this study. Random shuffles of the thermal conductivity were conducted ten times, and the results are shown in Table . The low Rrand2 and Qrand2 values lower than the original model suggests that good stability and predictive performance of the model are not based on the chance correlation.[23] In addition, CRp2 equals to 0.9094, which is close to the value of Rtr2, manifesting the robustness of the developed model.

Table 5

R2 and Q2 Values after Several Y-Randomization Tests

iteration	R2	Q2
1	0.01	0.01
2	0.01	0.01
3	0.01	0.01
4	0.01	0.02
5	0.01	0.02
6	0.01	0.01
7	0.01	0.01
8	0.01	0.01
9	0.02	0.00
10	0.00	0.02

For the training set, test set, and prediction set, the average absolute relative deviation (AARD) are 4.14, 4.41, and 4.16%, respectively, resulting in an AARD of 4.23% for the whole dataset, and for the entire data set, the squared correlation coefficient and the RMSE are 0.9038 and 0.0067, respectively. The calculated statistical metrics recommended by Tropsha et al.[24,25] for the test set are as follows, which testify the predictive ability of the developed modelor

According to the ref (26), both experimental and calculated endpoint values were taken in the log scale, and MAE values were calculated from long transformed values. In this study, the MAE, training set range, and σ of the multiple linear regression (MLR) model are 0.0174, 7.820 and 0.0113, respectively. As a result, it can be found that two conditions are satisfying the “GOOD model” criteria as the following

By using the “Prediction Reliability Indicator” tool, it can be found that none of the external compounds is “bad”. Detailed information can be seen in Supporting Information (Table S2).

Based on these results, it can be concluded that the GFA–MLR model has high accuracy, robustness, and good predictive ability.

Applicability Domain

As shown in Figure , the blue vertical line stands for the threshold value h* and the two horizontal red lines stand for ±3 standard deviation units. The AD is located in the region of 0 ≤ h ≤ 0.0303 and −3 ≤ R ≤ 3. It can be found that a majority of data points fell within the AD, which further indicates the robustness and validity of the developed model.

Figure 3

AD of the developed model.

Meanwhile, data points (thermal conductivity of formic acid at 300, 310, 320, 330, 340, 350, 360, and 370 K, and thermal conductivity of isoheptanol at 230, 240, 250, and 260 K) are located in the region of h > 0.0303 and −3 ≤ R ≤ 3, which are so-called as X outliers (good influence points). These compounds were well predicted by the developed model, which can stabilize it and make it more precise.[27] There is no Y outlier, which proves the reliability of the developed model.

Descriptor Interpretation

As discussed previously, the most significant influence of descriptors on thermal conductivity is SM2_B(s), SIC0, IC1, T, Eta_F, and MATS2m, in turn.

SM2_B(s)[28] is one of the 2D matrix-based descriptors. It can be interpreted as a spectral moment of order 2 from the Burden matrix weighted by I-state, which was proposed to evaluate the spectral moment from the Burden matrix and is expressed aswhere nSK represents the number of graph vertices and λ represents the eigenvalues of the Burden matrix. The Burden matrix is the augmented adjacency matrix derived from the H-depleted molecular graph. The adjacency matrix can describe the relationship of each vertex. Thus, SM2_B(s) can be used to investigate the effect of different molecule group interaction on thermal conductivity. The positive correlation coefficient for SM2_B(s) indicates that the stronger the molecule group interaction is, the larger the thermal conductivity is.

IC1(28,29) is the information content index (neighborhood symmetry of 1-order), one of the information indices. The general form of the neighborhood information content (ICk) is calculated as followswhere g runs over the equivalence classes, A represents the cardinality of the gth equivalence class, and nAT represents the number of molecule atoms. IC1 can be used to evaluate and describe the effect of molecules size, shape, and branching information on thermal conductivity. From the negative correlation coefficient, it can be found that the more the molecular branches are, the shorter the main chain is, and the lower the thermal conductivity is.

SIC0 (structural information content index (neighborhood symmetry of 0-order)) is a type of information indices. It is calculated as the information content IC0 normalized form to delete the graph size influence

It is a topological descriptor that can be used to measure the degree of atom diversity in a molecule and describe its shape. As for aliphatic oxygen-containing organic compounds, deleting the influence of graph size, it can be regarded as a measurement of number of oxygen atoms. The greater the number of oxygen atoms is, the more the charge transfer is, and the larger the thermal conductivity is. Thus, the SIC0 is positively related with thermal conductivity.

In addition, T (temperature) also has a remarkable effect on thermal conductivity. The positive value of coefficient means that T is negatively related to thermal conductivity. As the temperature rises, the liquid molecules move more chaotic, and the molecular directional movement from the high-energy region to the low-energy region weakens, which leads to the decrease of heat conduction.[30]

Eta_F is the abbreviation of eta functionality index, a type of ETA indices, which was proposed to evaluate the molecule functionality, here quantifying the presence of heteroatoms and multiple bonds.[28,31] For aliphatic oxygen-containing organic compounds, the larger the Eta_F value is, the more molecular branches are, and the shorter the main chain is. Accordingly, heat transfer along the main chain is reduced, and so is thermal conductivity. For example, the Eta_F values of n-butanol, isobutyl alcohol, and tertbutyl alcohol are 0.669, 0.722, and 0.772, respectively, and at 340 K, the thermal conductivity values of them are 0.143, 0.128, and 0.11.

The Moran autocorrelation of lag 2 weighted by mass MATS2m is one of the 2D autocorrelations. The general form of Moran autocorrelation is represented by using Moran coefficient to the molecular graphwhere w represents any atomic property, w̅ represents the w average value on the molecule, δ represents the Kronecker delta, and Δ represents the sum of the Kronecker deltas. The Moran coefficient usually takes value in the interval [−1, +1]. Positive spatial autocorrelation corresponds to positive coefficient values, whereas negative spatial autocorrelation corresponds to negative coefficient values. The molecules with positive autocorrelation have more branch chains than molecules with linear chain, of which the main chain is shorter and thermal conductivity value reduces. As a result, the MATS2m is negatively related with thermal conductivity.

All these six parameters have significant contribution to the thermal conductivity of oxygen-containing organic compounds.

Comparison between the Previous Studies and the Present Study

Comparisons with previous studies are shown in Table . Compared with the previous studies, the present study has been further improved. Under the premise of ensuring high fitting ability and predictive power, the present study has fully validated the developed models and analyzed the applicability domains (ADs).

Table 6

Comparisons with Previous Studies

author	methods	no	compounds class	N	R2	results
Gao and Cao[14]	MLR	155	alkanes	4	0.9510	s = 0.0033
Kauffman and Jurs[15]	MLR	213	organic solvents	9	0.953	RMSEP = 0.0136
Khajeh and Modarress[16]	GFA–MLR	116	alcohols	5	0.9521	RMSEP = 0.0474
Liu et al.	MLR	139	alcohols	4	0.9738	RMSEP = 0.0029
this study	GFA–MLR	1064	aliphatic oxygen-containing organic compounds	6	0.8922	AARD % = 4.41%
						RMSEtest = 0.0071

Conclusions

In the present study, a 6-descriptor linear QSPR model was developed for predicting the thermal conductivity of diverse aliphatic oxygen-containing organic compounds (aliphatic alcohols, ethers, aldehydes, ketones, acids, and esters). The molecular descriptors were calculated based on the optimized structures of aliphatic oxygen-containing organic compounds and were screened using the genetic function approximation (GFA) method. The predictive ability, model validation, and the AD of the developed model indicate that the GFA–MLR model has good predictive reliability and robustness. The GFA–MLR model provides a transparent output, with an AARD of 4.23% for the whole data set. Meanwhile, it can be concluded that molecule group interaction, molecular branches, number of oxygen atoms, temperature, and main chain of molecules have great effects on the thermal conductivity of aliphatic oxygen-containing organic compounds. The present study is of great significance not only by providing a robust model for predicting thermal conductivity of diverse aliphatic oxygen-containing organic compounds but also by shedding light on other thermodynamic data estimation.

Materials and Methods

Data set Collection

Thermal conductivity of the substance is mainly affected by the molecular structure and the temperature.[32] Therefore, in this study, at first, 992 thermal conductivity data of liquid aliphatic alcohols, ethers, aldehydes, ketones, acids, and esters at different temperatures were collected from the Handbook of Thermal Conductivity of Liquids and Gases.[33] These compounds covered 112 aliphatic oxygen-containing compounds. There are 290 thermal conductivity values of alcohols (29%), 78 thermal conductivity values of ethers (8%), 100 thermal conductivity values of aldehydes (10%), 112 thermal conductivity values of ketones (11%), 140 thermal conductivity values of acids (14%), and 272 thermal conductivity values of esters (28%). Thermal conductivity values of different chemical classes versus temperature are shown in Figure .

Figure 4

Thermal conductivity values of different chemical classes vs temperature.

In this study, Kennard Stone algorithm was performed to split the data set, which was developed based on the different Euclidean distance among all data. At first, the samples with larger difference were selected into the training set. Then, the remaining samples which are closer to the training data enter the test set. Under the circumstance, all representative samples can be divided into the training set. To some extent, the unevenness of the training set can be avoided.[34] In the end, the entire data set was divided into 694-member training set and 298-member test set with Kennard–Stone algorithm. The QSPR model would be developed based on the training set, and the test set was used to evaluate the performance and robustness of the developed model.[35]

In addition, to further testify the predictive power of the developed model, in our laboratory, 72 experimental thermal conductivity data of 15 compounds were measured as the prediction set. These data were measured using the DRE-2A thermal conductivity instrument with nonequilibrium state transient hot-wire method.[17] The transient hot wire method is a method for measuring the thermal conductivity of liquids and has the characteristics of high speed and high precision. According to the relationship between temperature and time, thermal conductivity was calculated automatically by the instrument. Each sample was measured three times, and the final data is an averaged value. The error of measurement is <3%.

The detailed information of the involved compounds and the whole data set which contains 1064 thermal conductivity data are shown in Supporting Information (Tables S1 and S2).

Molecular Descriptors

All molecule structures were constructed in GaussView graphical interface software package.[36] Then, the structures were output as the Gaussian input files. After modifying the input files, using the Keywords “Opt freq B3LYP/6-31G(d)”, the molecule structures were optimized with the Gaussian 09W.[37] Confirming that there is no virtual frequency in the convergence of the optimization result, the optimized structures were saved as .sdf format files. In the end, 4885 molecular descriptors for each optimized molecule were calculated in Dragon 6.0 software[38] (https://chm.kode-solutions.net/products_dragon.php) covering most of various theoretical approaches. The list of descriptors includes the simplest atom types, functional groups, fragment counts, topological and geometrical descriptors, three-dimensional descriptors, and so forth.

The calculated descriptors with constant and near-constant values, descriptors with at least one missing value (some descriptors of some molecules can’t be calculated), and descriptors with pair correlation larger than or equal to 0.90 were excluded in Dragon 6.0. As a result, 363 descriptors were retained.

Descriptor Screening and Model Development

Generally, the multiple linear regression method[39] is applied to describe the linear quantitative relationship between a dependent variable and multiple independent variables. It is used commonly in QSPR research on account of its simplicity, easy interpretability, and transparency.[40,41]

The GFA method can simulate biological evolution to generate statistical models, which was derived from the combination of genetic algorithm and multivariate adaptive regression splines.[42] In the present study, this method was used to screen molecular descriptors and provide several MLR models for selecting the best regression model. A flowchart of the GFA process can be seen in Supporting Information (Figure S1). Model development was carried out in the QSAR module of Material Studio 8.0 software. At first, multiple equations are randomly established on the basis of the population of descriptors. From the developed equations, the “parents” are selected according to the probability ratio of fitting, and the next generation equation is generated. After a crossover of each generation, mutations were made.[43] These operations will be repeated continuously for the specified number of iterations unless the convergence criteria are met. The fitting criteria for a GFA model can be evaluated using different scoring functions during evolution. In the present study, lack-of-fit was used as the scoring function because it can determine the most proper number of variables, avoid overfitting, and make the smoothness of the fitting under control.[44,45]

Model validation was performed to evaluate the goodness-of-fit statistic, robustness, and predictive power of the developed model.

Internal Validation

The role of the internal validation is to monitor the accuracy of the developed model and confirm the presence of overfitting. Generally, the commonly applied parameters to measure the goodness-of-fit statistic is the squared correlation coefficient (Rtr2), the RMSE, and AARD (%) of the training set. Moreover, for MLR model, there must be the test of multicollinearities among the descriptors to reduce the redundant parameters. Thus, the VIF was calculated. If VIF values of each descriptor in the developed model are lower than 10, it can be concluded that no multicollinearity exist among the descriptors.[22] In addition, leave-one-out cross validation (LOO cross validation) was applied to further evaluate the possibility of overfitting with QCV2, RMSECV.

External Validation

External validation was performed to assess the generalization ability of the developed model. Validation metrics such as Rtest2, RMSEtest, AARD %,[46,47]Qext-F12, Qect-F22, Qect-F32, rm2, and CCC[48] were applied to evaluate the model predicting performance. Tropsha et al.[24,25] proposed that whether the developed QSPR model successful depends on the following criteria of the test set

The parameters involved in above conditions are presented clearly in the Reference of Tropsha et al.[49]

In addition, in order to further illustrate the stability and robustness of the developed model, Roy et al.[26] proposed the MAE based criteria. If the MAE of the developed model satisfies the below conditions,[50] it can be concluded that the model is acceptablewhere the σ value represents the standard deviation of the absolute error values for the test set data.

Furthermore, in this study, the quality of predictions for the external compounds was checked with the following tool: “prediction reliability indicator”, which is proved very valid for multiple linear regression models.[50] This tool categorized the quality of predictions for the test set into three groups (good, moderate, and bad) based on absolute prediction errors.

Y-Randomization Test

In addition, Y-Randomization was performed to confirm the presence of chance correlation between the dependent variables and the independent variables.[24,52] For Y-Randomization test, low Rrand2 and Qrand2 compared to the original R2 and Q2 of the resulting model are expected.[53] Meanwhile, the CRp2 was applied to evaluate the chance correlation, which is expected to be close to the value of Rtr2.[54] The Y-randomization test was performed in the program package MLR Y-randomization 1.2 (http://dtclab.webs.com/software-tools).

It is essential to define a domain of applicability for the developed model, namely, the AD. The AD is determined by the response (property) of the compounds and the descriptors in the developed QSPR model.[55] The Williams plot with standardized cross validated residuals (R) versus leverage (Hat diagonal) values (h) is highly recommended, which is widely used for visualizing the AD.[56−58] It should be noted that if the h value of a compound in the training set is greater than the threshold value h* (h* = 3p/n, where p is the number of model variables plus one and n is the number of the training set data), the structure of this compound reinforces the developed model.[10,59] If the most data points are located in the region of 0 ≤ h ≤ h* and −3 ≤ R ≤ 3, the developed model can be considered statistically acceptable and valid.

If cross validated standardized residuals of a compounds is greater than three standard deviation units (R > 3σ) while the leverage value is lower than the threshold value (h < h*), this compound could be judged as the response outlier (Y outlier). If the leverage value of a compound is greater than threshold value (h > h*) while the value of R is lower than the 3 standard deviation units, this compound could be judged as the structurally influential compounds (X outlier).