Prediction of Second-Order Rate Constants of Sulfate Radical with Aromatic Contaminants Using Quantitative Structure-Activity Relationship Model

Abstract

Predicting the second-order rate constants between aromatic contaminants and a sulfate radical ($k_{{\mathrm{SO}}_4^{•-}}$) is vital for the screening of pollutants resistant to sulfate radical-based advanced oxidation processes. In this study, a quantitative structure-activity relationship (QSAR) model was developed to predict the values for aromatic contaminants. The relationship between $\log k_{{\mathrm{SO}}_4^{•-}}$ and three molecular descriptors (electron density, steric energy, and ratio between oxygen atoms and carbon atoms) was built through multiple linear regression. The goodness-of-fit, robustness, and predictive ability of the model were characterized statistically with indicators showing that the model was reliable and applicable. Electron density was found to be the most influential descriptor that contributed the most to $\log k_{{\mathrm{SO}}_4^{•-}}$. All data points fell within the applicability domain, and no outliers existed in the training set. The comparison with other models indicates that the QSAR model performs well in elucidating the mechanism of the reaction between aromatic compounds and sulfate radicals.

1. Introduction

Sulfate radical-based advanced oxidation processes (SR-AOPs) are considered as a promising technology for the treatment of wastewater with recalcitrant organic contaminants [1]. With a high standard reduction potential comparable to that of HO^• (E°(${\mathrm{SO}}_4^{•-}$/ ${\mathrm{SO}}_4^{2-}$) = 2.60 V vs. E°(HO^•/OH⁻) = 2.80 V), the sulfate radical (${\mathrm{SO}}_4^{•-}$) is highly reactive and capable of mineralizing recalcitrant organic contaminants [2]. Unlike HO^•, which reacts unselectively with organic compounds via hydrogen abstraction or addition, ${\mathrm{SO}}_4^{•-}$ prefers to attack organics rich in electron moieties through the single-electron transfer reaction [3]. Therefore, ${\mathrm{SO}}_4^{•-}$ is less affected by the complex water matrix in real wastewaters than HO^•. Moreover, ${\mathrm{SO}}_4^{•-}$ has a much longer lifetime than HO^•, making its reaction with pollutants more efficient [4]. Peroxydisulfate and peroxymonosulfate have been frequently used as the precursors of ${\mathrm{SO}}_4^{•-}$, which are more stable and safer for transportation than the commonly used HO^• precursor H₂O₂ [5]. The advantages of ${\mathrm{SO}}_4^{•-}$ make SR-AOPs a potential substitute to HO^•-based AOPs.

Aromatic compounds are widely produced and used. Some synthetic aromatic compounds such as pesticides and pharmaceuticals are resistant to conventional wastewater treatment processes and discharged directly into the aquatic environment [6]. They constitute a major group of organic pollutants in the aquatic ecosystem and pose threats to human health due to their potential carcinogenicity [7,8]. SR-AOPs have been extensively studied for the removal of recalcitrant aromatic contaminants [9]. The second-order rate constants of ${\mathrm{SO}}_4^{•-}$ with organic compounds ($k_{{\mathrm{SO}}_4^{•-}}$) in aqueous solution range from 10⁵ to 10⁹ M⁻¹ s⁻¹ [10]. Considering that the scavenging of ${\mathrm{SO}}_4^{•-}$ by inorganic anions and natural organic matters in the water matrix will inhibit the degradation of target contaminants [11], only contaminants with higher $k_{{\mathrm{SO}}_4^{•-}}$ than those compounds in the water matrix could be removed efficiently in real waters. Currently, the available $k_{{\mathrm{SO}}_4^{•-}}$ values are quite limited and obtained mainly through experiments [12]. As there are thousands of aromatic contaminants, experimental measurement of all the $k_{{\mathrm{SO}}_4^{•-}}$ values would be laborious, costly, and time consuming. Therefore, it is necessary to develop an alternative method to estimate the $k_{{\mathrm{SO}}_4^{•-}}$ values as they are imperative to the assessment of whether a target compound is suitable for SR-AOPs treatment.

Quantitative structure-activity relationship (QSAR) analysis focuses on correlating the molecular descriptors calculated from molecular structures with the activities of the corresponding chemicals (e.g., toxicity and biodegradability) [13,14,15]. QSAR models are mainly used to predict a particular physical or chemical property of a chemical compound and interpret the mechanism behind the prediction [16]. Numerous QSAR models have been developed to predict the rate constants of reactive species such as HO^•, O₃, singlet oxygen, and hydrated electron with organic contaminants [17,18,19,20]. However, QSAR models for ${\mathrm{SO}}_4^{•-}$ are quite limited. Xiao et al. [21] found that $k_{{\mathrm{SO}}_4^{•-}}$ was negatively related to the ratio of oxygen atoms to carbon atoms and the energy gap between the lowest occupied molecule orbital and highest occupied molecule orbital. Unfortunately, their model did not contain a variable positively correlated with the $k_{{\mathrm{SO}}_4^{•-}}$ values; thus, it is unable to tell what property of the chemical compound would facilitate its reaction with ${\mathrm{SO}}_4^{•-}$. Ye et al. [12] developed a model that linked ln$k_{{\mathrm{SO}}_4^{•-}}$ with 32 frequencies of structural fragment. However, too many independent variables in the model may cause overfitting problems and reduce the predictive ability of the model [22]. Much more effort could be made to improve the QSAR models for $k_{{\mathrm{SO}}_4^{•-}}$ prediction due to a lack of relevant studies currently.

In this study, a QSAR model for predicting the $k_{{\mathrm{SO}}_4^{•-}}$ of aromatic compounds was developed and validated by using the multiple linear regression (MLR) method. Two new descriptors (electron density and steric energy) were incorporated into the model to bring new insight into the ${\mathrm{SO}}_4^{•-}$ reaction with aromatic compounds. Following the guidance of The Organization for Economic Cooperation and Development (OECD) for QSAR development [23], statistical characteristics of the developed model were analyzed, the mechanism behind the model was interpreted, and the applicability domain of the model was assessed. Finally, a comparison with previous models was made. The results of this work would help to judge whether an aromatic compound is suitable for SR-AOPs treatment.

2. Materials and Methods

2.1. Dataset

The $k_{{\mathrm{SO}}_4^{•-}}$ values (M⁻¹s⁻¹) of 88 aromatic compounds were collected from published literature, and the corresponding references are cited in Table S1. Five descriptors, i.e., electron density, steric energy (kcal/mol), the ratio between the number of oxygen atoms and carbon atoms, volume of the molecule (ų), and octanol–water partitioning coefficient were calculated to represent the physical properties of selected chemicals. The electron density (E) represents the probability of an electron appearing in a specific space around an atom or molecule. The electron density on each atom of an aromatic compound was calculated by Gaussian 16 Rev. A.03 with HF method at 6-31G level [24]. The highest electron density on the benzene ring was selected as the descriptor E. Steric energy (S) is the sum of energies that resulted from bonded and non-bonded energies within a molecule, which reflects the energy due to the geometry of a molecule [25]. The calculation of S was run on the Chem3D 20.1.1 with the MM2 Dynamics method. The ratio between the number of oxygen atoms and carbon atoms (O/C) of a molecule was found to be negatively correlated with $k_{{\mathrm{SO}}_4^{•-}}$ [21]. The molecular volume (V) and octanol–water partition coefficient (logP) of each compound were obtained from Molinspiration Cheminformatics free web services (https://www.molinspiration.com, accessed date: 25 June 2021). The $\log k_{{\mathrm{SO}}_4^{•-}}$ and descriptor values of selected aromatic compounds are listed in Table S1. The 88 aromatic compounds were randomly divided into a training set and validation set with a ratio of 3:1. The data in the training set were used to develop the QSAR model, while the data in the validation set were used to verify the predictive ability of the developed model. MATLAB R2021b was used for the generation of the QSAR model and its validation.

2.2. QSAR Model Development and Characterization

The QSAR model in this study depicts the linear relationship between $\log k_{S{O}_4^{•-}}$ and the molecular descriptors, as shown in Equation (1). A stepwise multiple linear regression method was used to determine the significant descriptors in the model. A model with all five descriptors was built first. Then, the p-value of each descriptor was calculated. The descriptors with a p-value lower than the significance level (0.05) were considered as insignificant and then excluded from the model one by one until all the descriptors left were significant.

$$\log k_{{\mathrm{SO}}_4^{•-}}={\beta }_0+{\beta }_{1}E+{\beta }_{2}S+{\beta }_3\left(\mathrm{O}/\mathrm{C}\right)+{\beta }_4\log V+{\beta }_5\mathrm{logP}$$

(1)

The degree of multicollinearity among the descriptors was reflected by the variance inflation factor (VIF), which was calculated by Equation (2) [26]. The $r_i^{2 }$ is the coefficient of determination when conducting multiple linear regression between the ith descriptor and all the other descriptors. A value of 10 was suggested as the threshold of VIF, above which the multicollinearity is considered as severe [27]. The goodness-of-fit of the model was assessed by $R^2$, adjusted $R^2$($R_{\mathrm{adj}}^2$) and the root mean square error (RMSE). $R_{\mathrm{adj}}^2$ and RMSE were calculated by Equations (3) and (4), respectively, where n is the number of compounds in the training set, k is the number of descriptors, y_i is the experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ of the ith compound, and $\widehat{y_{\mathrm{I}}}$ is the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$. Leave-one-out cross-validation was conducted to estimate the robustness of the developed model. For all compounds in the training set, one compound was removed from the training set each time to test the model, and the rest of the compounds were used to train the model [28]. The indicator of leave-one-out cross-validation $Q_{\mathrm{LOO}}^2$ was calculated by Equation (5), where $\widehat{y_{\mathrm{I}}}$ is the log$k_{{\mathrm{SO}}_4^{•-}}$ predicted from model with the ith compound removed from the training set, and $\overline{y_t}$ is the average experimental log$k_{{\mathrm{SO}}_4^{•-}}$ of the training set [28]. $Q_{\mathrm{LOO}}^2$ measures the robustness of the developed model [29].

$$\mathrm{VIF}=\frac{1}{1-r_i^2}$$

(2)

$$R_{\mathrm{adj}}^2=1-\left(\frac{n-1}{n-k-1}\right)\left(1-R^2\right)$$

(3)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{n}}$$

(4)

$$Q_{\mathrm{LOO}}^2=1-\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y_t})}^2}$$

(5)

2.3. Validation of the Model

The validation set was used to test the predictive ability of the model developed from the training set. The indicator $Q_{\mathrm{ext}}^2$ was defined as Equation (6), where $y_j$ is the experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ of the jth compound in the validation set, $\widehat{y_{\mathrm{j}}}$ is the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ of the jth compound in the validation set, and $\overline{y_t}$ is the average $\log k_{{\mathrm{SO}}_4^{•-}}$ of all compounds in the training set. A higher $Q_{\mathrm{ext}}^2$ value means a better prediction by the model, and 0.5 was suggested as the threshold [30,31]. The external RMSE was calculated by Equation (7), which gauges the deviation of predicted values from the experimental ones.

$$Q_{\mathrm{ext}}^2=1-\frac{{\sum }_{j=1}^{test}{(y_j-\widehat{y_{\mathrm{j}}})}^2}{{\sum }_{j=1}^{test}{(y_j-\overline{y_t})}^2}$$

(6)

$${\mathrm{RMSE}}_{\mathrm{ext}}=\sqrt{\frac{{\sum }_{j=1}^{test}{(y_j-\widehat{y_{\mathrm{j}}})}^2}{n}}$$

(7)

Moreover, the conditions below should be met so that the predictive ability of the model could be considered satisfactory [30]:

$$R_{\mathrm{ext}}^2>0.6$$

(8)

$$\frac{R_{\mathrm{ext}}^2-R_0^2}{R_{\mathrm{ext}}^2}<0.1 \mathrm{or} \frac{R_{\mathrm{ext}}^2-R_0^{\prime 2}}{R_{\mathrm{ext}}^2}<0.1$$

(9)

$$0.85\le k \mathrm{and}{k}^{\prime }\le 1.15$$

(10)

where $R_{\mathrm{ext}}$ is the correlation coefficient between the predicted and experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ values; $R_0^2$ is the coefficient of determination when the experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ is regressed against the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ with the fitting curve through the origin and k is the corresponding slope; $R_0^{\prime 2}$ is the coefficient of determination when the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ is regressed against the experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ with the fitting curve through the origin, and $k^{\prime }$ is the corresponding slope.

A y-randomization test was used to verify the robustness of the developed model by measuring the degree of chance correlation between $\log k_{{\mathrm{SO}}_4^{•-}}$ and the descriptors [32]: the randomly shuffled $\log k_{{\mathrm{SO}}_4^{•-}}$ values in the training dataset were regressed against the fixed descriptors for several times to generate new MLR models, of which the R² shall be poor. The y-randomization test was conducted with Scikit-learn 1.0.2 [33] and the code was used to calculate the R² for each shuffle was listed in Text S1.

2.4. Relative Contribution of Each Descriptor

The relative contribution of each descriptor to the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ was estimated with Equation (11), where E_i, S_i, and (O/C)_i are the descriptors of the ith compound, and D_i is one of the three descriptors in the model. The analysis of relative contribution would help to elaborate on the role of each descriptor in the model.

$$\mathrm{Relative}\mathrm{contribution}(\%)=\frac{\left|{\beta }_i{D}_i\right|}{{\beta }_0+|{\beta }_1{E}_i\left|+\left|{\beta }_2{S}_i\right|+|{\beta }_3{\left(\mathrm{O}/\mathrm{C}\right)}_i\right|}\times 100$$

(11)

2.5. Applicability Domain

The applicability domain (AD) is used to define a region of chemicals with specific structures where the model could make an accurate prediction [29]. The visualization of AD was fulfilled by the Williams plot, of which the X-axis refers to hat values (leverages), and the Y-axis refers to standardized residuals (δ) [34]. The hat value was calculated with Equation (12), where h_i is the hat value of the ith compound, X is the n × k descriptor matrix containing all compounds in the training set, and x_i is the 1 × k descriptor vector of the ith compound (n is the number of compounds in the training set and k is the number of descriptors) [35]. The critical hat value (h*) was calculated with Equation (13) [36]. When h_i < h* and δ < 3, the predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ of the ith compound is considered reliable [21]. The standardized residual of the ith compound was calculated with Equation (14) [37].

$$h_i=x_i{(X^{T}X)}^{-1}{x}_i^T$$

(12)

$$h^{*}=\frac{3\left(k+1\right)}{n}$$

(13)

$${\delta }_i=\frac{y_i-\widehat{y_i}}{\sqrt{\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{n-k-1}}}$$

(14)

3. Results

3.1. Selection of Significant Descriptors

Backward stepwise regression was performed on the training set (n = 66) to select significant descriptors for the model. As shown in

Table 1. Selection of significant descriptors for the QSAR model by backward stepwise regression (n = 66).

Steps	Descriptors	Coefficient	t-Statistic	p-Value	${\mathit{R}}^2$	${\mathit{R}}_{\mathbf{adj}}^2$	RMSE	Decision
1	E	1.5755	5.7113	0.0000	0.685	0.659	0.219	Exclude logP
	S	0.0042	6.5814	0.0000
	O/C	−0.1924	−0.5424	0.5895
	LogV	−0.2572	−1.2209	0.2269
	LogP	0.0374	1.4067	0.1647
	Constant	1.8972
2	E	1.6799	6.2729	0.0000	0.675	0.653	0.221	Exclude logV
	S	0.0039	6.4515	0.0000
	O/C	−0.5938	−2.7933	0.0070
	LogV	−0.1412	−0.7225	0.4727
	Constant	1.2381
3	E	1.6882	6.3344	0.0000	0.672	0.656	0.220	Accept E, S and O/C as the significant descriptors
	S	0.0037	6.7496	0.0000
	O/C	−0.6035	−2.8559	0.0058
	Constant	0.8868

, when all five descriptors were included in the model, three descriptors (O/C, logV, and logP) were found to be insignificant, as their p-values were above the significance level of 0.05.

After excluding logP from the model, O/C became significant while logV was still insignificant. The decrease in $R^2$ and $R_{\mathrm{adj}}^2$ and the increase in RMSE were results from reducing the overfitting by taking out logP. Finally, the logV was removed from the model, and the descriptors left were all significant. Meanwhile, $R_{\mathrm{adj}}^2$ was increased and RMSE was decreased, indicating that the overfitting was reduced with an unnecessary descriptor removed. The formula of the model is shown below:

$$\begin{array}{c}\log k_{{\mathrm{SO}}_4^{•-}}=0.8868+1.6882E+0.0037S-0.6035\left(\mathrm{O}/\mathrm{C}\right)\\ n=66, R^2=0.672, R_{\mathrm{adj}}^2=0.656,\mathrm{RMSE}=0.220,F=42.30, p<0.0001\end{array}$$

(15)

3.2. Exclusion of Outliers

The analysis of residuals of the training set could help to identify outliers that affected the goodness-of-fit of the model. As shown in Figure 1, the compounds whose residual confidence intervals did not include zero point were classified as outliers and colored in red. After removing the identified outliers, the MLR would be repeatedly conducted to find new outliers until all the residual confidence intervals included the zero point.

After four rounds of MLR, five compounds in the training set were removed to improve the goodness-of-fit of the model. The final formula of the developed model is shown in Equation (16). By removing five outliers, $R^2$ was increased from 0.672 to 0.748, and RMSE was decreased from 0.220 to 0.193, indicating that the goodness-of-fit was significantly improved.

$$\begin{array}{c}\log k_{{\mathrm{SO}}_4^{•-}}=0.4641+1.7751E+0.0036S-0.5836\left(\mathrm{O}/\mathrm{C}\right)\\ R^2=0.748, R_{\mathrm{adj}}^2=0.735,\mathrm{RMSE}=0.193, F=56.43, p<0.0001\end{array}$$

(16)

3.3. Statistical Characteristics of the Developed Model

The statistical indicators of the model before and after excluding the five outliers are listed in

Table 2. Statistical characterization of the model before and after excluding the outliers.

Model	Goodness-of-Fit			Robustness	Predictive Ability			${\mathit{R}}_{\mathbf{0}}^{\mathbf{2}}$	${\mathit{R}}_{\mathbf{0}}^{\mathbf{\prime }\mathbf{2}}$	k	k′	$\frac{{\mathit{R}}_{\mathbf{ext}}^{\mathbf{2}}\mathbf{-}{\mathit{R}}_{\mathbf{0}}^{\mathbf{2}}}{{\mathit{R}}_{\mathbf{ext}}^{\mathbf{2}}}$	$\frac{{\mathit{R}}_{\mathbf{ext}}^{\mathbf{2}}\mathbf{-}{\mathit{R}}_{\mathbf{0}}^{\mathbf{\prime }\mathbf{2}}}{{\mathit{R}}_{\mathbf{ext}}^{\mathbf{2}}}$
	${\mathit{R}}^2$	${\mathit{R}}_{\mathbf{adj}}^2$	RMSE	${\mathit{Q}}_{\mathbf{LOO}}^2$	${\mathit{Q}}_{\mathbf{ext}}^2$	RMSEext	${\mathit{R}}_{\mathbf{ext}}^2$
Including outliers	0.672	0.656	0.220	0.617	0.605	0.282	0.632	0.066	0.610	1.010	0.990	0.896	0.035
Excluding outliers	0.748	0.735	0.193	0.694	0.603	0.289	0.648	0.110	0.624	1.013	0.987	0.830	0.037

. After excluding the outliers from the training set, the goodness-of-fit and robustness of the model were both improved. Even though $Q_{\mathrm{ext}}^2$ and RMSE_ext were slightly decreased, $R_{\mathrm{ext}}^2$ was increased. Therefore, the predictive ability of the model remained stable after the exclusion of outliers.

The conditions listed in inequalities 8–10 were all met, showing a satisfactory predictive ability of the model. The VIF of E, S, and O/C were 1.48, 1.11, and 1.32, respectively. Therefore, the multicollinearity among the variables was negligible, as the VIF values were well below the threshold of 10 and close to 1. As shown in Figure 2, the values of descriptors were far away from the flat surfaces, which represent the fitting results of E, S, and O/C by MLR when the other two descriptors acted as independent variables. It confirms that there is little multicollinearity among E, S, and O/C, suggesting the model has high stability [38].

The predicted $\log k_{{\mathrm{SO}}_4^{•-}}$ versus the experimental $\log k_{{\mathrm{SO}}_4^{•-}}$ is shown in Figure 3. The prediction made by the QSAR model agreed quite well with the experimental results, showing a high predictive ability of the developed model. The result of the y-randomization test is shown in Figure S1. The poor R² values (<0.16) for the MLR models developed from shuffled $\log k_{{\mathrm{SO}}_4^{•-}}$ against E, S, and O/C indicate that there is no chance correlation. Therefore, the MLR model developed (Equation (16)) can be trusted to predict $\log k_{{\mathrm{SO}}_4^{•-}}$ values from new descriptors.

4. Discussion

4.1. Interpretation of the Model

The coefficient of E is 1.7751, as shown in Equation (16), which is much larger than the other two coefficients, indicating that the change of $\log k_{{\mathrm{SO}}_4^{•-}}$ per unit E is the highest. Therefore, E is the most influential factor in the model. As shown in Figure 4, E also contributed the most to $\log k_{{\mathrm{SO}}_4^{•-}}$. ${\mathrm{SO}}_4^{•-}$ reacts with aromatic compounds mainly through single electron transfer (SET) from the benzene ring to the radical. As the occurrence of the SET process requires the electron extraction from the nucleophile [39], the electron density of the benzene ring may play an essential role in the electron transfer rate between ${\mathrm{SO}}_4^{•-}$ and aromatic compounds. It was assumed in this study that the highest electron density on the carbon atom of the benzene ring was positively correlated with $\log k_{{\mathrm{SO}}_4^{•-}}$; i.e., a higher electron density would lead to a faster reaction rate, which is consistent with the result of the model, as the coefficient of E is positive.

As shown in Figure 5a, R² of $\log k_{{\mathrm{SO}}_4^{•-}}$ vs. E is 0.455, suggesting a moderate correlation between the rate constant and electron density. However, Luo et al. [40] observed no clear relationship between the Gibbs free energy of SET ($\Delta G_{SET}^{°}$, which represents the reactivity of ${\mathrm{SO}}_4^{•-}$ toward aromatic contaminants) and the total electrostatic potential charge on the benzene ring (ESP, which represents the electron density on the benzene ring). The disparity might be due to the different ways in calculating electron density. In this study, the highest rather than the total electron density on the carbon atom in the benzene ring was used as it was assumed that ${\mathrm{SO}}_4^{•-}$ might preferentially attack the carbon atom in the benzene ring with the highest electron density.

The coefficient of S is 0.0036, of which the absolute value is the least among the three coefficients. Therefore, S is the least influential factor in the model. As shown in Figure 4, the relative contribution of S was comparable to that of O/C for most of the aromatic compounds selected, except for those with the benzenesulfonamide moiety, which shows a much higher contribution of S than that of O/C. The steric energy is a sum of energies from bond stretching, bending, torsion, Van der Waals, and electrostatic interactions within a molecule [41]. The lowest energy conformation of a molecule is most favored, and it is achieved when the steric energy is minimized [25]. As the reactivity of a molecule is affected by its geometry, there should be a certain kind of relationship between the steric energy and the reactivity of a molecule. The positive correlation between S and $\log k_{{\mathrm{SO}}_4^{•-}}$ indicates that higher steric energy tends to make the aromatic molecule more reactive to ${\mathrm{SO}}_4^{•-}$. As shown in Figure 5b, R² of the linear regression between S and $\log k_{{\mathrm{SO}}_4^{•-}}$ is 0.503, which is even higher than that between E and $\log k_{{\mathrm{SO}}_4^{•-}}$. However, the data points in Figure 5b were closely gathered, while in Figure 5a, the data points were scattered around the regression line. Therefore, a higher R² of $\log k_{{\mathrm{SO}}_4^{•-}}$ vs. S does not mean a higher correlation between $\log k_{{\mathrm{SO}}_4^{•-}}$ and S.

The negative coefficient of O/C (–0.5836) indicates that the increase in oxygen content would inhibit the reactivity of aromatic compounds toward ${\mathrm{SO}}_4^{•-}$. Oxygen atoms, especially those attached directly to the benzene ring, have a strong ability to withdraw electrons from the benzene ring. As a result, the electron density on the benzene ring is reduced, and the electron transfer process is hindered. Xiao et al. [21] also reported that $\log k_{{\mathrm{SO}}_4^{•-}}$ was negatively correlated with O/C, and the correlation was quite strong as R² of $\log k_{{\mathrm{SO}}_4^{•-}}\mathrm{vs}.$ O/C was 0.621. They suggested that an increase in O/C would decrease the number of H atoms, thus weakening the H abstraction by ${\mathrm{SO}}_4^{•-}$, which would slow down the reaction. However, no clear linear relationship was observed between $\log k_{{\mathrm{SO}}_4^{•-}}$ and O/C in this study (Figure 5c). There were only six compounds containing no oxygen atom in the training set (65 compounds in total) of Xiao et al. [21], while there were 11 in the training set of this study (61 compounds in total). The difference in training set composition might lead to the disparity in the importance of O/C to the models.

4.2. Applicability Domain

The standardized residual depicts the difference between the experimental and predicted results. As shown in Figure 6, all the δ values were within ±3, indicating no outlier existed. The h_i value represents how far the x_i value of the ith compound deviates from the average X value of all the compounds [42]. There was only one compound (gallate ion) in the training set with an h_i larger than h*, which could be called a “good high leverage point”, as it made the model more stable and accurate [42]. It was suggested that such a compound would have excessive influence during the model development process [34]. The developed model exhibits good extrapolating ability, as all data points from the validation set were within the AD, suggesting that the $\log k_{{\mathrm{SO}}_4^{•-}}$ of chemicals with structures similar to those in the training set may be reliably predicted.

4.3. Comparison with Other Models

Studies regarding the QSAR model for $k_{{\mathrm{SO}}_4^{•-}}$ prediction are quite limited, and only two relevant literatures could be found. The comparison between previous models and the one developed in this study is listed in

Table 3. A comparison of QSAR models for $k_{{\mathrm{SO}}_4^{•-}}$ prediction.

Reference	Model Type	n*	Molecular Descriptors	Goodness-of-Fit			Robustness	Predictive Ability			AD
				${\mathit{R}}^2$	${\mathit{R}}_{\mathbf{adj}}^2$	RMSE	${\mathit{Q}}_{\mathbf{LOO}}^2$	${\mathit{Q}}_{\mathbf{ext}}^2$	RMSEext	${\mathit{R}}_{\mathbf{ext}}^2$
Xiao et al. (2015)	MLR	65	Ratio of oxygen atoms to carbon; ELUMO and EHOMO energy gap	0.866	-	-	0.86	0.89	-	0.89	All but one compound from the validation set was outside the AD
Ye et al. (2017)	MLR and ANN	75	32 molecular fragment descriptors	0.88 (MLR); 0.99 (ANN)	-	-	-	-	-	0.62 (MLR); 0.42 (ANN)	-
This study	MLR	61	Electron density, steric energy, and ratio of oxygen atoms to carbon	0.748	0.735	0.193	0.694	0.603	0.289	0.648	All data points of the validation set fell within the AD

* The number of compounds in the training set.

.

Xiao et al. [21] developed the first QSAR model for $k_{{\mathrm{SO}}_4^{•-}}$ prediction using the MLR method. Their model has a better goodness-of-fit, robustness, and predictive ability than the one in this study. However, the ratio of oxygen atoms to carbon atoms was the dominant descriptor in their model, which would limit the prediction of $k_{{\mathrm{SO}}_4^{•-}}$ for compounds without oxygen atoms. Moreover, both descriptors in the model have negative coefficients, which means that the model could only reveal the factors that reduce the reactivity of compounds toward ${\mathrm{SO}}_4^{•-}$. In this study, the model consists of two descriptors with positive coefficients and one descriptor with a negative coefficient. Therefore, both enhancing and inhibitory factors to the reaction between aromatic compounds and ${\mathrm{SO}}_4^{•-}$ are explained by the model. Ye et al. [12] developed two models with 32 molecular fragment descriptors by using MLR and artificial neural network (ANN), respectively. They found that the ANN model showed much better goodness-of-fit but much lower predictive ability than MLR. Considering that the number of descriptors was close to the number of compounds in the training set, overfitting would weaken the predictive ability of the model developed by Ye et al. [12], which was also reflected by the low $R_{\mathrm{ext}}^2$ (0.62 for MLR and 0.42 for ANN). In addition, the descriptors in their model stand for the frequencies of structural fragments, which possess no physical meaning. Therefore, the mechanism behind the model was difficult to explain. In contrast, our model has successfully reduced the overfitting problem by excluding insignificant descriptors, and the physical properties of each descriptor help to elucidate the mechanism behind the reaction between aromatic compounds and ${\mathrm{SO}}_4^{•-}$.

5. Conclusions

In this study, a QSAR model was developed and validated to predict the second-order rate constants between aromatic compounds and ${\mathrm{SO}}_4^{•-}$. The stepwise MLR was used to exclude insignificant descriptors, and the final model was composed of E (electron density), S (steric energy), and O/C (number of oxygen atoms vs. carbon atoms). Residual analysis was applied to remove outliers from the training set, and the goodness-of-fit of the model was improved. The statistical indicators for goodness-of-fit ($R_{\mathrm{adj}}^2$= 0.735), robustness ($Q_{\mathrm{LOO}}^2$ = 0.694), and predictive ability ($Q_{\mathrm{ext}}^2$ = 0.603 and $R_{\mathrm{ext}}^2$ = 0.648) suggest that the model is satisfactory and applicable. E was the most influential descriptor and contributed most to the $\log k_{{\mathrm{SO}}_4^{•-}}$. The positive coefficient of E suggests that higher electron density on the benzene ring could enhance its reactivity to ${\mathrm{SO}}_4^{•-}$, which is consistent with the assumption that the reaction between ${\mathrm{SO}}_4^{•-}$ and aromatic compounds was achieved mainly through single electron transfer from the benzene ring to the radical. S is positively correlated with $\log k_{{\mathrm{SO}}_4^{•-}}$, suggesting that higher steric energy could improve the reactivity of aromatic compounds via affecting its geometry. An increase in O/C would decrease $\log k_{{\mathrm{SO}}_4^{•-}}$ as oxygen is electron withdrawing, so that the electron density on the benzene ring could be depleted. All data points in the validation set fell into the applicability domain, suggesting that the model is suitable for aromatic compounds with various functional groups. The comparison with other models shows that the QSAR model developed performed better in mechanism elucidation and overfitting reduction.