Efficient Chemical Space Exploration Using Active Learning Based on Marginalized Graph Kernel: an Application for Predicting the Thermodynamic Properties of Alkanes with Molecular Simulation

The AL algorithm used in this work is shown in Algorithm 1. At the beginning, the selected set $S$ (two molecules) are randomly selected from the target CCS $M$, where the remaining molecules constitute the pool set $P$ to be selected during the AL process. In each iteration of AL, a Gaussian process regression (GPR) [rasmussen_gaussian_2006] model is trained using the selected set $S$ and used to predict the uncertainties of $P$

where $K_{PP}$ is the kernel matrix of $P$, $K_{SS}$ is the kernel matrix of $S$, $K_{SP}$ is the kernel matrix between $S$ and $P$. The molecule $\alpha$ with the highest predicted uncertainty is then added to the selected set $S$. In practice, the pool set $P$ may contain a large number of molecules, so it is too time-consuming to calculate the predicted uncertainty of all molecules in the pool set. We accelerated the active learning process by randomly selecting $N_t$ molecules from $P$ for prediction, discarding all data points with prediction uncertainty less than the threshold $U_t$, and then selecting the one with the highest prediction uncertainty to add to $S$. AL iterations are repeated until the pool set $P$ is empty (all predicted uncertainties of P are less than the threshold value $U_t$). The AL algorithm described above is called explorative AL because the criterion for selecting the most valuable molecules is the predicted uncertainty, i.e., the AL algorithm keeps exploring the CCS by finding the molecule in $P$ that is least similar to $S$. Therefore, $S$ is a representative subset of the CCS $M$ at the end of AL. Every molecule in $S$ always has similar molecules in S, which guarantees that an ML model trained using $S$ as a training set can make accurate predictions for all molecules in $P$. Furthermore, the molecules in $S$ are not similar to each other, which indicates the number of molecules in $S$, as well as the cost of obtaining the target properties of these molecules, are minimized.

We chose GPR-MGK as the surrogate ML model to compute the predicted uncertainty as shown in equation 1. The MGK was used to compute the kernel matrix describing the similarity between molecules by comparing the similarity between simultaneous random walk paths on a pair of graphs. We used the MGK proposed in our previous work [xiang_predicting_2021], which is designed to predict the thermodynamic properties of molecular liquids. The formula of MGK is expressed as:

The bond kernel $K_e(\cdot ,\cdot )$ is defined in the same way as atom kernel. In addition, the MGK is normalized to improve prediction performance

The atom features, bond features, and hyperparameters were tuned by minimizing the cross-validation error in the critical temperature data set. For the details about MGK, we refer the reader to ref [xiang_predicting_2021].

Here, we highlight the advantages of GPR-MGK in AL. The first advantage of GPR-MGK is its reliable uncertainty quantification, as demonstrated in our previous studies [xiang_predicting_2021, xiang_comparative_2021]. The quality of the surrogate model’s predicted uncertainty is critical to the AL algorithm because it guarantees that the selected sample is truly different from the molecules in $S$, resulting in maximum information gain. The second advantage is the compatibility with high-throughput simulations or experiments, i.e., data selection can be performed prior to data generation. According to our previous studies, the hyperparameters optimized using the critical temperature data set work well in various data sets such as density, heat capacity, viscosity, etc[xiang_predicting_2021]. The transferability of hyperparameters allows us to use the same hyperparameters during the AL process. According to equation 1, the predicted uncertainties of GPR only depend on the kernel matrices, which are computed from the molecular graphs, and their target properties of the molecules in S are not required. Therefore, the selected subset $S$ can be obtained prior to data generation using MS. The logic behind this is transfer learning [cai_transfer_2020], where the hyperparameters are the “knowledge” learned from critical temperature data sets that can be transferred to other data sets. In contrast, with graph neural networks or random forests as surrogate models, each AL iteration must wait for the MS of the newly added molecules in the previous iteration to complete before proceeding. Batched active learning that add multiple samples in each AL iteration is a common solution for this compatibility [zhang_bayesian_2019, graff_accelerating_2021], but the performance suffers severely due to the high similarity between the samples added in each batch.

We used the TEAM-AMBER force field [gong_temperature_2018] and HT-MS developed by Gong et al. [gong_predicting_2018, cao_high-throughput_2017] to generate the liquid density, heat capacity, and heat of vaporization (HOV) of alkane liquids. Liquid density was computed by averaging the density of the trajectories. Heat capacity was computed by combining the intramolecular and intermolecular contributions

where $C_{trans}$ and $C_{rot}$ are the ideal translational and rotational contributions computed through the equipartition theorem, $C_{vib}$ is the vibrational contribution computed using the hindered-rotor model [ayala_identification_1998], and $U_{inter}$ is the intermolecular potential. In literature [sun_compass_1998, wang_application_2011], HOV was broadly computed as

Gong et al. found that equation 7 systematically overestimates the HOV of alkanes, thus proposing an empirical correction

where $n_C$ is the number of carbon atoms, and the coefficient $-\frac{1}{15}$ was obtained by fitting the experimental data.

For each molecule, the simulated temperature was set to 16 data points ranging from 0.4 $T_c$ to 0.9 $T_c$, and the simulated pressure was set to 1 bar. The critical temperature $T_c$ was predicted by a D-MPNN [yang_analyzing_2019] trained on experimental data from the NIST database [diky_thermodata_2016, noauthor_knovel_nodate]. With the SMILES string, temperature, and pressure as inputs, the HS-MS module automatically generated all input files, submitted the jobs to the HPC clusters, analyzed the simulation trajectory, and collected the results. OpenBabel [oboyle_open_2011] and Packmol [martinez_packmol_2009] were used to generate input 3D molecular configurations. The quantum chemistry calculations were performed using Gaussian 16 [g16]. The molecular dynamics simulations were performed using GROMACS [hess_gromacs_2008]. We developed a series of quality control tests, and simulations that fail were discarded or extended until they pass the test: (1) Check the distribution of the kinetic energy trajectory via the Kolmogorov-Smirnov test. The null hypothesis is the distribution satisfies the Maxwell-Boltzmann distribution, and the simulation is extended if p-value less than 0.01 [merz_testing_2018]. (2) If the density is less than 50 $kg\cdot m^{-3}$, the simulation is considered as vapor and discarded. (3) If the diffusion coefficient is less than ${10}^{-8}c{m}^2\cdot s^{-1}$, the simulated system is considered solid and discarded. (4) The trajectories of temperature, pressure, density, and potential are detected using the method proposed by Chodera to determine the convergence and production interval [chodera_simple_2016]. (5) If the density, heat capacity, and HOV are not monotonic functions of temperature, the simulations are discarded. (6) The density, heat capacity, and HOV are fitted as quadratic functions of temperature. If the fitting score $R^2$ is less than 0.98, the simulations are discarded. For more details on the HT-MS, we refer the reader to ref [gong_predicting_2018].

After the MS was completed, an MS database $\{(x_i,T_i,{\rho }_i,C_{p,i},H_{vap,i}),i=1,2,\cdots \}$ was generated, where $x_i\in S$ is a molecule selected via AL, $T_i$ is the temperature, $\rho$ is the density, $C_{p,i}$ is heat capacity and $H_{vap,i}$ is HOV. Pressure is not considered because the pressure was set to 1 bar for all simulations. Three D-MPNNs were trained using $\{x,T\}$ as input and $\rho ,C_P,H_{vap}$ as output respectively. The density, heat capacity, and HOV of all molecules in $M$ were predicted using the D-MPNNs.

D-MPNN [yang_analyzing_2019] is a variant of message passing neural network (MPNN) [gilmer_neural_2017] which predicts molecular properties with molecular graphs and RDKit features [noauthor_descriptor_2021] as input. D-MPNN contains two phases: message passing and readout. In the message-passing phase, messages are passed through bonds, and the hidden states of atoms and bonds are updated according to the incoming messages. The hidden states of atoms are aggregated (sum or average) to obtain a feature vector of the molecular graph. The molecular representation is the concatenation between the feature vector learned via message passing and the RDKit feature vector [noauthor_descriptor_2021]. The readout phase is a fully connected neural network that connects the learned molecular representation from message passing and target property. We used Chemprop [noauthor_message_2021] library for training and prediction of D-MPNN with default settings except for followings: 100 epochs for training; aggregation function was “mean” for density and HOV, “sum” for heat capacity. In order to incorporate temperature into the model input, we concatenated molecular representation with temperature between the message-passing and readout phases.

We chose D-MPNN because its computational cost is lower than GPR-MGK. Our previous benchmark study shows that the predictions of D-MPNN are correlated to that of the GPR-MGK at a per-sample level [xiang_comparative_2021]. Therefore, it is reliable to first use GP-MGK AL to select samples to construct the training set, and then use D-MPNN to make predictions.

We used alkane with 4 to 19 heavy atoms as target CCS. The SMILES strings for these molecules were generated by an atom-addition algorithm. Alkanes with N heavy atoms were generated by adding a carbon atom to the alkanes with N-1 heavy atoms. Duplicate SMILES strings were discarded by canonicalization using RDKit[landrum_rdkit_2013]. In total, we collected 251728 distinct molecules.

The target CCS contains 251728 molecules. We performed a three-stage AL with uncertainty thresholds $U_t=0.5,0.4,0.3$, and three subsets $S_{U_t=0.5}$ (313 molecules), $S_{U_t=0.4}$ (599 molecules), $S_{U_t=0.3}$ (1679 molecules) were obtained, respectively. The AL set with a high uncertainty threshold is the subset of the AL set with a low uncertainty threshold: $S_{U_t=0.5}\subset S_{U_t=0.4}\subset S_{U_t=0.3}$.

The distribution of the CCS and AL sets along the number of heavy atoms is shown in Figure 2. The number of molecules in the CCS grows exponentially with heavy atoms, while the number of molecules selected by AL grows much more slowly. This suggests that the distribution of molecules in CCS becomes denser as the number of heavy atoms increases, so the proportion of molecules selected by AL becomes less and less. Therefore, it is promising to explore a huge CCS with affordable cost by only focusing on the small subset selected by AL. As elucidated below, the ML model trained with only hundreds of molecules selected by AL predicts the entire CCS with extremely high accuracy.

The densities, heat capacities, and HOV of the AL-selected molecules were obtained through HT-MS, and D-MPNNs were trained for each property and each AL set respectively. In order to evaluate the prediction accuracy of the ALMS framework, a test set $S_{test}$ (3000 molecules) were randomly selected from the CCS excluding the AL sets.

For different AL sets as training sets and different target properties, the predictions of D-MPNNs versus the simulation data of the test set are compared in Figure 3. For $S_{U_t=0.5}$, the root mean square errors (RMSEs) of the predictions are 6.2 $kg\cdot m^{-3}$ for density, 9.7 $J\cdot mo{l}^{-1}\cdot K^{-1}$ for heat capacity and 0.94 $kJ\cdot mo{l}^{-1}$ for HOV. Despite the high accuracy, $R^2>0.99$, continued improvement in prediction accuracy was observed as $U_t$ decreased to 0.4, 0.3, and the training set increased from 313 to 599, 1679 molecules.

The uncertainty threshold $U_t$ in AL is a parameter that balances the computational cost and prediction accuracy. The smaller the $U_t$, the more molecules are simulated, and the larger the computational cost, the more accurate the ALMS prediction. For the CCS of alkane considered in this work, setting $U_t$ to 0.5 is accurate enough for practical consideration. Further reduction of $U_t$ requires several times the computational cost, but the improvement in model prediction accuracy is relatively small.

To evaluate the power of the explorative AL algorithm introduced in this study, we conducted a controlled experiment. 313 molecules (the same as $S_{U_t=0.5}$) are randomly selected from $S_{t}est$ as the training set of D-MPNN, and the union of the remaining molecules and $S_{U_t=0.5}$ is used as the test set. Although this is not a truly random sample, we believe it has no effect on our conclusions, as only $\frac{313}{251728}=0.667\%$ of the molecules is ignored. For all target properties, AL outperformed random sampling significantly as indicated in Figure 4. The drawback of random sampling is that it lacks the sampling of sparse regions in the CCS. In our study case, an aspect easily understood is that alkanes with heavy atomic numbers from four to six were all selected by AL, because these molecules are small and they sparsely distributed on the CCS.

It is essential to evaluate our model by comparing the predictions with experimental data, which examines the combinatory prediction errors due to both machine learning and MS. From the NIST database collected by knovel [diky_thermodata_2016, noauthor_knovel_nodate], experimental data on the density, heat capacity, and HOV of 240, 240, and 235 alkanes were collected, respectively. For different AL sets as training sets and different target properties, the predictions of D-MPNNs versus the experimental data of the test set are compared in Figure 5. For $S_{U_t=0.5}$, The RMSEs and R2 are 19.3 $kg\cdot m^{-3}$ and 0.947 for density, 10.2 $J\cdot mo{l}^{-1}\cdot K^{-1}$ and 0.996 for heat capacity and 2.6 $kJ\cdot mo{l}^{-1}$ and 0.958 for HOV. The RMSEs for density and HOV are about twice as high as when the test set was simulated data (Figure 3), indicating that the density and HOV calculated by MS, although accurate, still did not reach the level of experimental reproducibility. On the contrary, the predictions on heat capacity reach the level of experimental reproducibility. Another important phenomenon is that the prediction accuracy of $S_{U_t=0.5}$, $S_{U_t=0.4}$, $S_{U_t=0.3}$ is almost the same, which indicates that simulating more molecules does not improve the prediction accuracy. This is because, simulating more molecules can only reduce the error of the ML but the error of MS dominates. Therefore, only a very small number of molecules need to be simulated to reach the upper limit of the ALMS predictions.