The Monte Carlo approach to model and predict the melting point of imidazolium ionic liquids using hybrid optimal descriptors.

Ionic liquids (ILs) have captured intensive attention owing to their unique properties such as high thermal stability, negligible vapour pressure, high dissolution capacity and high ionic conductivity as well as their wide applications in various scientific fields including organic synthesis, catalysis, and industrial extraction processes. Many applications of ionic liquids (ILs) rely on the melting point (T m). Therefore, in the present manuscript, the melting points of imidazolium ILs are studied employing a quantitative structure-property relationship (QSPR) approach to develop a model for predicting the melting points of a data set of imidazolium ILs. The Monte Carlo algorithm of CORAL software is applied to build up a robust QSPR model to calculate the values T m of 353 imidazolium ILs. Using a combination of SMILES and hydrogen-suppressed molecular graphs (HSGs), the hybrid optimal descriptor is computed and used to generate the QSPR models. Internal and external validation parameters are also employed to evaluate the predictability and reliability of the QSPR model. Four splits are prepared from the dataset and each split is randomly distributed into four sets i.e. training set (≈33%), invisible training set (≈31%), calibration set (≈16%) and validation set (≈20%). In QSPR modelling, the numerical values of various statistical features of the validation sets such as R Validation 2, Q Validation 2, and IICValidation are found to be in the range of 0.7846-0.8535, 0.7687-0.8423 and 0.7424-0.8982, respectively. For mechanistic interpretation, the structural attributes which are responsible for the increase/decrease of T m are also extracted. This journal is © The Royal Society of Chemistry.

Introduction

Generally, ionic liquids (ILs) are a class of room temperature molten salts with melting points near or below 100 °C and these are made up of organic cations and organic/inorganic anions. ILs have been widely studied because of their potential applications in academia and industry, such as organic synthesis, electrochemistry, nanotechnology, catalysis and industrial extraction processes.[1-6] Melting point (Tm) is an important physical property of ILs as it determines their physical state and is also used in the calculation of physicochemical properties such as aqueous solubility, vapour pressure and phase equilibrium properties.[7,8] Melting points (Tm) are especially important for ILs since they have a wide liquidus range due to their low melting points and high thermal decomposition temperature (Td) and their solubility in water or organic solvents is highly associated with their Tm.[9,10]

The melting point of organic molten salts based on the imidazolium cation is low and many of them are liquid at ambient temperature. Imidazolium ILs have attracted a lot of attention because of their solvating characteristics, high conductivity, and a broad range of electrochemical stability as well as their wide applications in different areas including catalysis, supercapacitors and photovoltaics.[11-13] The development of task-specific ILs is feasible by in-depth standardization of the anions and cations. On the other hand, methods related to the design and development of ionic liquids is laborious and expensive. So, using computer-assisted approaches to develop ILs may save a lot of time and money. Hence, the development of the quantitative structure–property relationship (QSPR) models for the target endpoint i.e. melting point (Tm) of the ionic liquids has met with due success.

Quantitative structure–property relationships (QSPR) is a statistical approach to link the physical characteristics of a chemical substance to its chemical structure and this approach is also used to predict the property of an unknown compound. The CORAL (CORrelation And Logic) programme (available at http://www.insilico.eu/coral) has been recommended as a tool for doing QSPR analysis on a variety of endpoints.[14-19] The simplified molecular input line-entry system (SMILES) notations of the chemical structures are used to compute the descriptor correlation weight (DCW) in the CORAL software using Monte Carlo optimization.[20-23] In recent times, many publications utilized the ‘index of ideality of correlation (IIC)’ as a unique criterion to construct the best predictive QSPR models.[24-27]

The goal of the present study is to construct the QSPR models based on a hybrid optimal descriptor obtained from SMILES and a hydrogen-suppressed graph (HSG) to predict the Tm for a dataset of 353 imidazolium ILs. The index of ideality of correlation (IIC) is assessed as a criterion of predictive potential of the QSPR models of Tm. The accuracy and reliability of the constructed QSPR model are rigorously examined utilizing four random splits.

Method

Data

The experimental data on the Tm of 353 imidazolium ILs were obtained from Zhang et al.[28] The range of melting point temperature (Tm) values was from 180.65 to 541.15 K. The molecular structure of all imidazolium ILs was sketched using the free software BIOVIA draw 2019 and displayed by the SMILES notation. The compound IDs, SMILES codes, and corresponding experimental and predicted Tm are given in the ESI in Table S1.† Four splits were prepared from the 353Tm data and each split was randomly distributed into four sets i.e. training set (≈33%), invisible training set (≈31%), calibration set (≈16%) and validation set (≈20%). Each set was assigned a specific task in the development of the QSPR models and was well cited in the literature.[29-31]

Hybrid optimal descriptor

The CORAL platform can provide three optimal descriptors: graph-based, SMILES-based, and hybrid, computed with both the graph and the SMILES. Further, three molecular graph-based descriptors can also be created using the CORAL platform: HSG (hydrogen-suppressed graph), HFG (hydrogen-filled graph) and GAO (graph of atomic orbitals). Earlier studies indicate that QSPR models constructed with hybrid optimal descriptors are more reliable than models constructed by only using SMILES or molecular graphs.[32,33] In the present QSPR modelling, a hybrid descriptor, a combination of SMILES and HSG, is applied.

The Tm of imidazolium ILs, computed by the Monte Carlo optimization of CORAL software, is represented in the following equation:here, C0, C1 and DCW (descriptor of correlation weights) are the regression coefficients, the slope assessed by the least-squares method and the optimal descriptor calculated by molecular features extracted from HSG and SMILES notations, respectively. The T* denotes the threshold value and N* denotes the number of epochs for the Monte Carlo optimization.

So, the hybrid optimal descriptor employed here is represented by the following equation

The SMILES and graph-based optimal descriptors are computed using the following equations:

The notation S, SS, and SSS are local SMILES features interpreted as one-, two-, and three-symbols of SMILES notations; NOSP and HALO are global SMILES attributes which are computed according to the absence or presence N, O, S, P and halogens; the BOND is a mathematical function that describes the existence or absence of double (), triple (#), or stereochemical bonds (@ or @@); PAIR is the pairing of NOSP, BOND, and HALO; HARD represents the existence or absence of all of the above-mentioned characteristics (BOND, NOSP, and HALO) of the molecular structure; Cmax, Nmax and Omax and Smax are the maximum numbers of rings, nitrogen, oxygen, and atoms present in a molecule structure.

In eqn (4), pt2 and pt3 are the numbers of paths length 2 and 3 starting from a given vertex in a graph, respectively; S2 is the valence shell of the second orders; nn is the nearest neighbours; C5 and C6 are codes of rings (five-member and six-member rings, with the data on presence or absence of heteroatoms, aromaticity, and the total number of given rings in the molecule).[6,15,34]

The Monte Carlo optimization procedure can be performed based on two types of target functions: the balance of the correlation without IIC (TF1) and the balance of correlation with IIC (TF2). The mathematical equations of the two target functions are as the following:RTRN and RiTRN are the correlation coefficients between optimal descriptor and observed melting point for the training and invisible-training sets, respectively.

Here, Monte Carlo optimization based on IIC (TF2) is used to generate QSPR models. However, in eqn (6), the weight of IIC (WIIC) is an empirical coefficient (here WIIC = 0.2) and IICCAL is the index of ideality of correlation for calibration-set.

The following mathematical equation is applied to compute the IIC:[35-38]RCAL is the correlation coefficient between experimental values and predicted values of endpoint for the calibration set. The negative and positive mean absolute errors are shown with −MAE and +MAE, which are calculated using the following equations:

The ‘k’ is the index (1, 2, …. N) and the observed and calculated are related to values of the endpoint.

The numerical data on the above equations can be achieved with the CORAL software.

Domain of applicability

The applicability domain (AD) is a hypothetical chemical space area that encompasses both the model descriptors and the predicted response. The AD of substances is employed to estimate the ambiguity in the prediction of a given chemical based on how close it is to the substances employed to create the model in the building of a QSPR model. Because it is impossible to determine the whole world of compounds employing a particular QSPR model, the prediction of a modelled response utilizing QSPR is only meaningful if the molecule being predicted falls inside the AD of the model.[39]

In the QSPR models developed by CORAL software, the applicability domain is computed by the allocation of SMILES attributes in the training, invisible training, and calibration sets and is defined as ‘DefectA’.[40,41]DefectPTRN(A) and PCAL(A) are the probability of an attribute ‘A’ in the training and the calibration sets; NTRN(A) and NCAL(A) are frequencies of A in the training and calibration sets, respectively.

The statistical defect (D) can be defined as the sum of statistical defects of all attributes present in the SMILES notation.

NA is the number of active SMILES attributes for the given compounds.

In CORAL, a substance is an outlier if inequality 13 is fulfilled:D is an average of statistical defect for the dataset of the training set.

Validation of the model

It is important to validate the predictive potential of a created QSPR model. As per OECD guidelines, validation is defined as “the process by which the reliability and relevance of a particular approach, method, process, or assessment is established for a defined purpose”.[30,42]

Here, to determine the robustness, reliability and the predictive capability of the QSPR models for Tm of ILs three strategies were used: (i) internal validation or cross-validation; (ii) external validation and (iii) Y-scrambling or data randomization.[43,44] The mathematical relationship of different validation parameters employed herein is given in Table 1. Finally, IIC is employed to judge better models.[45,46]

The mathematical relationship of validation parameters used for the predictive potential of QSPR models

The criterion of the predictive potential	References
	22
	47
	48
	48
	48
	22
	49
	22
	22

Results and discussion

QSPR modelling for Tm

The QSPR models for the prediction of the Tm of different imidazolium ILs were constructed based on the hybrid optimal descriptor and Monte Carlo optimization using target function TF2 with the IIC. These QSPR models are demonstrated by the following equations:

Split 1

Split 2

Split 3

Split 4

The statistical parameters of constructed QSPR models based on eqn (14)–(17) for four random splits are summarized in Table 2.

The summary of statistical quality and criteria of predictability of the QSPR models

Split	Set	n	R 2	CCC	IIC	Q 2	Q F 1 2	Q F 2 2	Q F 3 2	R m 2		r̄ m 2	Δrm2	Y-r	S	MAE	F
1	Training	113	0.7862	0.8803	0.6548	0.7802					0.7848				31.3	24.6	408
	Invisible training	112	0.7864	0.8865	0.7868	0.7780					0.7830				23.5	18.8	405
	Calibration	55	0.8196	0.9029	0.9053	0.8079	0.8112	0.8103	0.3341		0.8070				24.4	18.8	241
	Validation	54	0.8204	0.8954	0.8972	0.8060				0.8086		0.7419	0.1335	0.0279	28.0	22.8	238
2	Training	116	0.8023	0.8903	0.7278	0.7958					0.7998				26.4	20.6	463
	Invisible training	111	0.8334	0.8859	0.6057	0.8277					0.8299				26.4	20.4	545
	Calibration	57	0.8256	0.9071	0.9005	0.8137	0.8163	0.8144	0.8301		0.8136				24.5	20.2	260
	Validation	50	0.8535	0.9133	0.8982	0.8423				0.8271		0.7889	0.0764	0.0215	24.7	20.7	280
3	Training	109	0.8116	0.8960	0.7922	0.8052					0.8064				24.9	19.2	461
	Invisible training	107	0.8226	0.8900	0.8274	0.8155					0.8195				27.9	22.5	487
	Calibration	62	0.7809	0.8687	0.8837	0.7665	0.7287	0.7267	0.6810		0.7747				33.3	26.1	214
	Validation	56	0.7846	0.8818	0.7784	0.7687				0.6838		0.6965	0.0256	0.0218	27.5	22.3	197
4	Training	118	0.8232	0.9031	0.8195	0.8183					0.8188				25.2	18.8	540
	Invisible training	107	0.8551	0.9038	0.7369	0.8503					0.8471				23.1	17.6	620
	Calibration	62	0.8177	0.8952	0.9042	0.8035	0.8224	0.8154	0.7975		0.8093				26.8	22.0	269
	Validation	47	0.8323	0.8986	0.7424	0.8163				0.7888		0.7077	0.1623	0.0182	23.9	17.3	223

Taking into account the statistical criteria of the various validation parameters, all the constructed QSPR models were statistically good and had robust predictions. Furthermore, all of the generated QSPR models matched the MAE requirements of the error-based measure.[50] Also, the Y-randomization test was carried out to evaluate the robustness of constructed models and the numerical value of was found more than 0.5 for all models. The numerical value of the determination coefficient (RValidation2 = 0.8535) of the validation set for split 2 was found highest, so this model was nominated as the dominant model (Table 2). Plots of predicted Tmversus experimental Tm of four models developed by TF2 are shown in Fig. 1A. The plots of residual Tmversus predicted Tm for all subsets of all splits are represented in Fig. 1B. The dispersion of residual Tm was found to be near the horizontal line centred around zero and confirmed that all developed QSPR models were well fitted.

Fig. 1

Experimental Tmversus predicted Tm values (A) and residual of Tmversus predicted Tm (B) for four QSPR models constructed by TF2.

In the CORAL-QSPR model constructed utilizing the Monte Carlo method, the applicability domain was employed to identify outliers. The average of defect SMILES was 1.39684 for split 2 and imidazolium ILs appeared into the domain of applicability if defect SMILES < 2.79367. The number of outliers present in the QSPR models constructed by TF2 was 9, 6, 8, and 8 for the splits 1, 2, 3 and 4 respectively.

Mechanistic interpretation

One of the OECD principles is the mechanistic interpretation of the generated QSPR model which means molecular structure information may be retrieved from such models. Because established optimal chemical descriptors may be connected with suitable chemical pieces, all models created using the Monte-Carlo technique follow the concept of mechanistic interpretation. In CORAL QSPR, the structural attributes (SA) are used for mechanistic interpretation. The comparable statistical reliability of the model should be theoretically achieved in multiple rounds of Monte Carlo optimization. So, four classes of the attributes can be classified depending upon the numerical data for correlation weights (CWs) in three or more independent optimization runs: (i) if CW(SA) is positive in all runs then these attributes are classified as a promoter of Tm increase; (ii) if CW(SA) is negative in all runs then these attributes are classified as a promoter of Tm decrease; (iii) if CW(SA) is both positive and negative in all runs then these attributes are classified as undefined; and (iv) if CW(SA) = 0 these attributes are kept in rare or blocked category. The structural attributes as a promoter of increase and decrease extracted from split 2 (the best QSPR model) are listed in Table 3.

The list of the promoter of increase/decrease of Tm extracted from split 2 using TF2

Type of descriptors	No.	SAk	CWs			NT, a	NiT, b	NC	Defect [SAk]	Comments
			Probe 1	Probe 2	Probe 3
Promoter of increase
Graph-based descriptors	1	VS2-C…5…	1.34960	2.43334	1.31153	108	99	50	0.0003	Valence shell of the second order for aliphatic carbon atom equal to 5
	2	PT3-C…5…	0.10558	1.51087	0.92938	99	87	41	0.0006	The presence of the path of length 3 equal to 5 for a carbon atom
	3	PT2-C…4…	0.24752	2.00169	0.35048	72	72	37	0.0003	The presence of the path of length 2 equal to 4 for a carbon atom
SMILES based descriptors	1	1………..	7.29691	5.95405	6.34387	113	112	55	0.0000	Presence of a cyclic ring
	2	c………..	1.71525	2.42427	0.64969	109	102	55	0.0002	Presence of aromatic carbon
	3	n………..	2.38850	1.33337	0.67039	108	102	55	0.0002	Presence of aromatic nitrogen
	4	c…(…….	2.49408	0.42892	2.49631	101	91	48	0.0002	Branching at an aromatic carbon
	5	(…C…(…	1.27092	2.30761	2.86614	88	89	39	0.0006	Combination of aliphatic carbon with two branching
	6	c…(…C…	1.33043	0.11257	0.06315	83	80	39	0.0003	Aromatic carbon joined by branching with the aliphatic carbon atom
	7	[…1…….	0.89010	1.76115	0.62621	73	74	32	0.0006	Presence of branching connected to the ring
	8	n…(…….	0.94747	1.97741	2.03406	71	60	32	0.0005	Presence of aromatic nitrogen and branching
Promoter of decrease
Graph-based descriptors	1	VS2-F…6…	−0.31025	−0.58746	−0.85323	67	70	30	0.0005	Valence shell of second-order equal to 6 for a fluorine atom
	2	PT2-C…1…	−0.76524	−0.99261	−0.74397	85	81	45	0.0004	The presence of the path of length 2 equal to 1 for a carbon atom
SMILES based descriptors	1	BOND10000000	−2.42991	−1.67430	−2.03277	60	58	30	0.0001	Presence of double bonds and absence of triple and stereochemical bonds
	2	[…S…….	−0.94276	−0.55305	−0.90267	35	37	13	0.0015	Combination of branching and aliphatic sulphur
	3	[…C…….	−1.90870	−0.15225	−0.77063	76	67	32	0.0009	Presence of branching connected to aliphatic carbon

According to the outcomes mentioned in Table 3, the graph-based descriptors as promoters of Tm increase for ILs are: (i) VS2-C…5…, PT3-C…5…, PT2-C…4… and the graph-based descriptors as promoters of decrease are VS2-F…6…, PT2-C…1….

In the same way, some SMILES based descriptors as Tm enhancer are: 1……….., c……….., n……….., c…(……., (…C…(…, c…(…C…, […1……., n…(……., and some promoter of Tm decrease SMILES based attributes are: BOND10000000, […C……., […S……..

Comparison with the previous report

Data of QSPR models reported in the literature is shown in Table 4. These results are commonly based on the use of quantitative QSPR, group contribution methods (GCM), artificial neural network (ANN) model clustering techniques and CODESSA program.[9,51-55] Huo et al.[52] developed a QSPR model (R2 = 0.8984) based on the group contribution method to predict the melting points of 190 imidazolium and benzimidazolium ILs. J. A. Cerecedo-Cordoba et al.[51] used QSPR techniques based on clustering methods for the prediction of melting points using two data sets of 281 (R2 = 0.78)and 134 (R2 = 0.93) imidazole ionic liquids. A. R. Katritzky et al.[54] reported QSPR models based on CODESSA program to predict the melting points of 104 imidazolium bromides and 45 benzimidazolium bromides. They could not develop a unified correlation for the investigated ILs. Instead, they suggested several QSPRs models. They divided 149 substituted ILs on the basis of the N-substituents into four subsets: subset A having 57 compounds, subset B having 29 compounds, subset C having 18 compounds and subset D having 45 benzimidazolium bromides. The numerical value of determination coefficient were (i) set A, R2 = 0.7442, (ii) set B, R2 = 0.7517, (iii) set C, R2 = 0.943, and (iv) set D, R2 = 0.68992. López-Martin et al.[55] developed a QSPR model (R2 = 0.869) for prediction of the Tm of 84 imidazolium-based ILs. They used the molecular descriptors generated by CODESSA and DRAGON software.

The comparison between some of the previous models and the present study for the prediction of Tm of imidazolium ILsa

Descriptor type	Feature selection method	Machin learning method	Data set size		R 2		RMSD		Ref.
			Training	Test	Training	Test	Training	Test
CODESSA	BMLR	MLR	16	3	0.90	0.9815	19.2	13.2	9
			25	4	0.92	0.8622	15.2	29.1
PaDEL-descriptor	Tree feature selection	MLR	291	—	0.78	—	18.2	—	51
Group contribution descriptors	—	Group contribution method	190	—	0.90	—	28.2	—	52
Artificial neural networks	Multilayer perceptron network (MLP)	ANN	97	—	0.99	—	—	—	53
CODESSA	BMLR	MLR	57	—	0.74	—	29.2	—	54
			25	—	0.75	—	14.5	—
			18	—	0.94	—	17.7	—
			45	—	0.69	—	20.0	—
Dragon and CODESSA	—	PLS	22	—	0.95	—	—	—	55
	—		62	—	0.87	—	—	—
Materials Studio	Genetic algorithm	MLR	50	10	0.88	0.74		29.9	56
		BA-ANN	50	10	0.91	0.95		12.2
CORAL	Monte-Carlo	LR	226	109	0.83	0.85	26.0	24.7	This work

BMLR: best multilinear regression method, PLS: partial least squares, MLR: multiple linear regression, ANN: artificial neural network.

The comparison depicts that the R2 of the calculated CORAL models are qualitatively comparable or even better than most of the other models suggested in the previous report for Tm of imidazolium ILs. However, previously reported models were developed utilizing a very small dataset, whereas the present dataset contains 353 imidazolium ILs, which is more than any reported work. In the present study. The numerical value of MAE and RMSE of the constructed QSPR models is smaller than all the reported QSPR models. Another superiority of the present method is that CORAL software does not require any physicochemical parameters or 3D optimization structure, as well as does not require the previous calculations of the chemicals descriptors for the creation of the QSPR models. Generally, the comparative data demonstrated the supremacy of the built models over the earlier published models.

Conclusion

For the prediction of the Tm of 353 imidazolium ILs, QSPR models were developed using a hybrid descriptor based on SMILES and HSG attributes. Using the balance of correlation method with the index of ideality correlation (IIC), reliable QSPR models for the Tm of 353 imidazolium ILs were successfully constructed. The IIC improves the predictive potential of described models. Four models were developed from four random splits and all the models were found statistically satisfactory and robust. The numerical value of the correlation coefficient (R2 = 0.8535) of the validation set of split 2 was highest than the other models and therefore, it was designated as the key model. The authenticity and robustness of the developed models were predicted by various statistical parameters such as R2, CCC, IIC, Q2, Q2, Q2, Rm2, etc. The present hybrid QSPR models were more robust and predictive than models reported in the literature. The structural attributes acting as promoters of Tm increase or decrease were also recognized.

Conflicts of interest

The authors declare no conflicts of interest.