Development of Prediction Models on Base-Catalyzed Hydrolysis

Subscriber access provided by - Access paid by the | UCSB Libraries

Environmental Modeling

Development of Prediction Models on Base-Catalyzed Hydrolysis Kinetics of Phthalate Esters with Density Functional Theory Calculation Tong Xu, Jingwen Chen, Zhongyu Wang, Weihao Tang, Deming Xia, Zhiqiang Fu, and Hong-bin Xie Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.9b00574 • Publication Date (Web): 08 Apr 2019 Downloaded from http://pubs.acs.org on April 8, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

Environmental Science & Technology

1

Development of Prediction Models on Base-Catalyzed Hydrolysis Kinetics of

2

Phthalate Esters with Density Functional Theory Calculation

3

Tong Xu, Jingwen Chen,* Zhongyu Wang, Weihao Tang, Deming Xia, Zhiqiang Fu,

4

Hongbin Xie

5

Key Laboratory of Industrial Ecology and Environmental Engineering (Ministry of

6

Education), School of Environmental Science and Technology, Dalian University of

7

Technology, Dalian 116024, China Table of Contents Graphic

8

9 10

ABSTRACT

11

Many phthalate esters (PAEs) are chemicals of high production volume and of

12

toxicological concern. Second-order rate constant for base-catalyzed hydrolysis (kB) is a key

13

parameter for assessing environmental persistence of PAEs. However, kB for most PAEs are

14

lacking, and experimental determination of kB encounters various difficulties. Herein, density

15

functional theory (DFT) methods were selected by comparing empirical kB values of 5 PAEs

16

and 5 carboxylic acid esters with the DFT calculated ones. Results indicate that PAEs with

17

cyclic side chains are more vulnerable to base-catalyzed hydrolysis than PAEs with linear

18

alkyl side chains, followed by PAEs with branched alkyl side chains. By combining

19

experimental and DFT calculated second-order rate constants for base-catalyzed hydrolysis

ACS Paragon Plus Environment

1


Page 2 of 32

20

of one side chain in PAEs (kB_side chain), quantitative structure-activity relationship (QSAR)

21

models were developed. The models can differentiate PAEs with departure of leaving-group

22

(or nucleophilic attack of OH‾) as the rate-determining step in the hydrolysis, and estimate

23

kB values, which provides a promising way to predict hydrolysis kinetics of PAEs. Half-lives

24

of the investigated PAEs were calculated and vary from 0.001 hours to 558 years (pH = 7 ~

25

9), further illustrating the necessity of prediction models for hydrolysis kinetics in assessing

26

environmental persistence of chemicals.

27 28

INTRODUCTION

29

Phthalate esters (PAEs) are extensively used as additives in plasticizers, pesticides,

30

cosmetics and personal care products.1-8 PAEs are not chemically bonded to products, so they

31

can be easily released into the environment during manufacturing, application, and disposal

32

of relevant products.9,10 They have been detected in various aquatic environments, including

33

groundwater, surface water, wastewater and seawater, with concentrations up to hundreds

34

μg·L-1.11-15 Several PAEs are suspected endocrine disruptors,4,5,16 and have drawn regulatory

35

concerns due to their pervasive pollution and human exposure.5,17-21 The United States

36

Environmental Protection Agency (US EPA) has listed 6 PAEs [di(2-ethylhexhyl) phthalate

37

(DEHP), dimethyl phthalate (DMP), diethyl phthalate (DEP), butyl-benzyl phthalate (BBP),

38

di-n-butyl phthalate (DBP), di-n-octyl phthalate (DNOP)] as priority pollutants.22,23

39

Specifically, DEHP is categorized as substances that have known or presumed carcinogenic

40

potential within the globally harmonized system of classification and labeling of

41

chemicals.8,24 Hence, it is of importance to understand the environmental behavior of PAEs

42

for assessment of their persistence and ecological risk.


2

Page 3 of 32


43

Due to the hydrolytic functional groups of carboxylic ester bonds, hydrolysis is assumed

44

to be a transformation pathway for some PAEs in the aquatic environment. PAE hydrolysis

45

products can have varied toxicological effects and higher environmental persistence than the

46

parent compounds.25 For example, mono (2-ethylhexhyl) phthalate (MEHP), a hydrolysis

47

product of DEHP, exhibits stronger peroxisome proliferator-activated receptor γ (PPARγ)

48

binding potential than DEHP.25 Therefore, it is of importance to investigate hydrolysis of

49

PAEs.

50

Previous studies indicate that for all the hydrolyzable compounds, reactions with OH‾

51

(base catalysis) is important even at conditions of pH < 7, and acid catalysis is relevant only

52

for compounds showing rather slow hydrolysis kinetics.26 PAEs are also typically stable

53

under acid and neutral conditions and tend to be hydrolyzed into phthalate mono-esters under

54

alkaline conditions.26-31 The hydrolysis half-lives (t1/2) of PAEs can be calculated from

55

second-order rate constant for base-catalyzed hydrolysis (kB), t1/2 = ln2/(kB·[OH¯]), vary

56

from seconds to years.28,29 Therefore, kB is a key parameter for assessing environmental

57

persistence of PAEs.32-35 Based on a survey of published literature and public available

58

databases, it was found there were only five PAEs that have experimental kB values, which

59

is quite insufficient relative to hundreds of commercially produced PAEs.28,36 Therefore, it

60

is necessary to determine kB for PAEs.

61

Experimental determination of kB values for PAEs is generally laborious and time-

62

consuming, due to possible long t1/2 for some PAEs, unavailability of authentic chemical

63

standards for some PAEs and analytical difficulty of PAEs caused by their structural

64

similarity and widespread pollution. Density functional theory (DFT) calculation, may

65

overcome the limitations encountered in the experimental determination.33,37-40 Zhang et al.


3


Page 4 of 32

66

employed DFT to investigate hydrolysis pathways and kinetics for cephradine, and the

67

hydrolysis products and rates predicted by DFT were consistent with experimental

68

observations.37 Sviatenko et al. employed DFT to unveil the mechanism of alkaline

69

hydrolysis and degradation rates of an energetic material octahydro-1,3,5,7-tetranitro-

70

1,3,5,7-tetrazocine, and found DFT predicted products corresponded to experimental

71

species.33 Therefore, it is feasible to investigate PAE hydrolysis pathways and kinetics by

72

DFT calculation.

73

Nevertheless, it is not pragmatic to calculate the kB values for all PAE molecules either,

74

due to the high cost of DFT computation. Hence, it is necessary to develop high-throughput

75

methods, such as quantitative structure-activity relationship (QSAR) models. Based on

76

multiple linear regression (MLR) and artificial neural network analysis, Wang et al.

77

developed QSAR models for estimating hydrolysis kinetics of disinfection byproducts.41

78

Actually, the Hammett equation or linear free energy relationships (LFERs) that lay a

79

theoretical foundation for the QSAR methodology, have already illustrated the influence of

80

ring substituents on the base-catalyzed hydrolysis of substituted benzoic acid esters in 1937.26

81

It deserves mentioning that the EPI SuiteTM which has been developed by the US EPA since

82

1980s for estimating chemical properties for risk assessment of chemicals, predicts kB values

83

of esters based on the 2D key structure R1-C(=O)-OR2, where R represents fragment

84

substituents of compounds.42-44 However, many fragment substituents in PAE molecules are

85

not covered in the EPI SuiteTM. Therefore, there is no available QSAR model that can be

86

employed for predicting kB values of PAEs.45

87

In this study, 5 PAEs with experimental kB values were adopted as reference compounds.

88

Based on the kB values of the references, appropriate DFT method was evaluated and selected


4

Page 5 of 32


89

to calculate kB values of other PAEs. By combination of the experimental kB values and those

90

calculated by DFT, QSAR models were developed for predicting kB of PAEs, which promises

91

to fill the data gap for kB of PAEs efficiently.

92 93

COMPUTATIONAL DETAILS

94

Model Chemicals. By searching literature and database such as SciFinder46 and PubChem,47

95

25 common PAEs with definite molecular structures and different side chains were chosen

96

as model compounds (Figure 1), excluding those with unclear molecular structures (i.e.,

97

mixtures, mixture of isomers). The side chains are linear or branched in most cases, and their

98

structural diversity can facilitate QSAR modeling and understanding the effects of side

99

chains on PAE hydrolysis kinetics. As shown in Table S1 of the supporting information (SI),

100

five PAEs (DMP, DEP, DBP, DIBP and DEHP) have experimental kB values,28,36 and five

101

carboxylic acid esters (formic acid methyl ester, acetic acid methyl ester, acetic acid ethyl

102

ester, acetic acid butyl ester and acetic acid propyl ester) have experimentally derived

103

activation energy (Ea) values.48,49 These ten compounds were selected as reference

104

compounds. Ea (kJ·mol-1) relates with experimentally observed kB [kB_obs, (s-1·mol-1·L)]

105

through the following Arrhenius equation:

106

kB_obs  A  exp(

 Ea ) R T

(1)

107

where A (s-1·mol-1·L) is the pre-exponential factor or frequency factor, R (8.314 J·mol-1·K-1)

108

is molar gas constant, T (K) is temperature. It deserves mentioning that one kB value of DEHP

109

reported in 1980 by a EPA laboratory at Athens, Georgia, was (1.1 ± 0.1) × 10-4 s-1·mol-1·L

110

at 30 C.28 In 2017, the laboratory reported a renewed value for kB value of DEHP at 25 C,

111

5.09 × 101 s-1·mol-1·L.36 Therefore, the updated value was adopted in the current study.


5


O

O O O

O Dimethyl phthalate DMP

O O O

O Diethyl phthalate DEP

O

O O O

O

O O O O

O Diisodecyl phthalate DIDP

O

O

O O O O

O

O O O

Dicyclohexyl phthalate DCP O

O O

O Octyl n-decyl phthalate ODP

O O O

Ditridecyl phthalate DTDP

O

Butyl decyl phthalate BDP

Diisononyl phthalate DINP

O O O

Diisotridecyl phthalate DIUP

O O

O 2-octanol butyl and octyl phthalate, OBOP

O O O

O

Diundecyl phthalate DUP

O

Di-n-octyl phthalate DNOP

O O O

O O

O O O

Di(2-propylheptyl) phthalate DPHP

Diisoheptyl phthalate DIHP O

O O O

O O O

Di-n-hexyl phthalate DNHP

O

O O O

O

O

Bis(4-methylpentyl) phthalate, BMPP

O O O Di-n-butyl phthalate DBP

O O

O

Diisooctyl phthalate DIOP

O

O O O

O O

O

O Diallyl phthalate DAP

O

Di-n-pentyl phthalate DNPP

O O O

O

O

Di(2-ethylhexyl) phthalate DEHP

O Di-n-propyl phthalate DPP

O O

O Diisobutyl phthalate DIBP

O O O

Page 6 of 32

O Butyl cyclohexyl phthalate BCP

O O O Butyl benzyl phthalate BBP

112 113

Figure 1. Molecular structures of PAEs under study (The PAEs marked in green are

114

reference compounds.)

115

DFT calculation. The flexible side chains can lead to PAEs with various conformations.

116

Therefore, it is necessary to identify the optimal conformation for DFT calculations. The

117

search for global minimum of various conformations was carried out by TURBOMOLE

118

program,50 as detailed in the SI. The conformation with the lowest energy was chosen for

119

further calculations.

120

Gaussian 09 program51 was adopted for the DFT calculation. Since the B3LYP52 and M06-

121

2X53 functional have been previously employed to investigate hydrolysis mechanisms and

122

kinetics with satisfactory performance,33,37 they were initially tested by predicting the kB

123

values for the reference compounds. The integral equation formalism of the polarized

124

continuum model (IEFPCM) was selected to simulate the water solvent effect. As the atomic

125

radii defining the solute cavity can influence calculation results,54 DEP was selected as a

126

reference to evaluate the reliability of the IEFPCM with three different solute cavities [the


6

Page 7 of 32


127

universal force field (UFF),55 Pauling56 and Bondi57]. The effect of explicit water molecules

128

on hydrolysis of PAEs was also evaluated with DEP as a case.

129

Geometries of reactants (R), transition states (TS), intermediates (IM) and products (P)

130

were optimized at the IEFPCM/B3LYP/6-311++G(2d,2p) level. Frequency analysis was

131

performed at the same level to ensure that TSs have only one imaginary frequency, and other

132

geometries have real frequencies. Intrinsic reaction coordinate (IRC) calculations were

133

performed to confirm that TSs connect the reactants and expected products. Single point

134

energy was calculated at the IEFPCM/B3LYP/6-311++G(3df,2pd) level.

135

The Gibbs free energy of reaction (∆G) for all hydrolysis pathways was calculated by the

136

energy difference between R and P, which can give insights into the spontaneity of assumed

137

hydrolysis pathways. For thermodynamic favorable hydrolysis pathways, difference of Gibbs

138

free energy between R and TS of rate-determining steps was calculated as the Gibbs free

139

energy of activation (∆G‡). Details for ∆G‡ and Ea calculations are provided in the SI. The

140

aqueous molar standard state was considered for all the reactions,33,37 as detailed in the SI.

141

Transition state theory (TST) was employed to calculate the second-order rate constant for

142

the base-catalyzed hydrolysis of one side chain in PAEs (kB_side chain):

143

kB_side chain   rd  (

kb  T G ‡ )  exp( ) h R T

(2)

144

where κrd is the transmission coefficient, kb (J·K-1) is the Boltzmann’s constant, h (J·s) is the

145

Planck constant, T was set as 298.15K. For PAEs with two same side chains, the DFT

146

calculated kB-DFT was calculated by kB-DFT = 2kB_side chain, with the assumption that these two

147

side chains have same probabilities to be hydrolyzed and have same kB_side chain values. For

148

PAEs with two different side chains, kB-DFT was obtained by summing up kB_side chain of the

149

two different side chains (kB-DFT = kB_side chain 1 + kB_side chain 2), as these two side chains can be


7


Page 8 of 32

150

hydrolyzed with different probabilities and kB_side chain values.

151

QSAR Modeling. Both the experimentally determined and DFT-calculated second-order

152

base-catalyzed hydrolysis rate constants were combined in the QSAR modeling. For the

153

PAEs with two same side chains and with experimental kB-exp values, the endpoint for QSAR

154

modeling was log(kB-exp/2). For the PAEs with DFT-derived kB-DFT values, the endpoint was

155

logkB_side chain. The endpoint values were randomly split into training and validation sets with

156

a ratio of 5:1.

157

As quantum chemical descriptors have clear physicochemical definitions58,59 and Dragon

158

molecular descriptors60 can describe molecular structural fragments, these two categories of

159

molecular structural descriptors were considered for the modeling. The quantum chemical

160

descriptors were calculated at the IEFPCM/B3LYP/6-311++G(2d,2p) level52,54 with the

161

Gaussian 09 software package.51 Based on the molecular structures optimized by Gaussian

162

09, Dragon molecular descriptors were calculated with Dragon (Ver. 6.0).60

163

Specially, quantum chemical descriptors and Dragon descriptors for sub-structures of PAE

164

molecules were calculated. When considering a specific ester bond in one side chain of PAE

165

molecules, the whole molecular structure of PAEs can be expressed as R1-(C=O)-OR2, where

166

OR2 stands for the leaving-group. Quantum chemical descriptors for R1-H, HO-R2, and whole

167

PAE molecules were calculated. Details on the calculation of the molecular descriptors, a

168

full list of the molecular descriptors, and values of some descriptors are provided in the SI.

169

All the descriptor values were standardized to eliminate the influence of different data range

170

for different descriptors. MLR was used for establishing QSAR models.

171 172

RESULTS AND DISCUSSION


8

Page 9 of 32


173

Selection of DFT Methods. As shown in Figure 2, hydrolysis of the PAEs proceeds

174

through a pre-reactive complex (IM1), a TS1 formed by nucleophilic attack of OH‾ on the C

175

atom of the ester bond, an intermediate (IM2) and TS2 formed by departure of the leaving-

176

group from the ester bond, and finally leads to a phthalate mono-ester.26 The free energy

177

surfaces for DMP and DEP are shown in Figure 3. For DBP, DIBP and DEHP, the free energy

178

surfaces are shown in Figure S6. It can be observed there are two types of rate-determining

179

steps. For DMP and DEHP, the rate-determining step is TS1. For DEP, DBP and DIBP, the

180

rate-determining step is TS2.

181 182

Figure 2. Schematic mechanism of base-catalyzed hydrolysis of PAEs

183 184

Figure 3. Schematic free energy surfaces for base-catalyzed hydrolysis of DMP and DEP

185

calculated at the IEFPCM/B3LYP/6-311++G(2d,2p) level [The total energy of the reactants

186

and OH‾ is set to zero (reference state). The symbols “R, IM1, TS, IM2, P” refer to reactants,

187

pre-reactive complexes, transition states, intermediates and products involved in the reaction,


9


188

Page 10 of 32

respectively]

189 190

Tables S2 and S3 list the calculated kB-DFT values of the reference PAEs. It can be seen that

191

the logkB-DFT values calculated at the B3LYP/6-311++G(2d,2p) level are closer to the

192

experimental (logkB-exp) values than the M062X/6-311++G(2d,2p) method. The logkB-DFT

193

values calculated with IEFPCM(UFF) as solvation models are closer to the logkB-exp values

194

than those with IEFPCM(Pauling) and IEFPCM(Bondi).

195

As can be seen from Figure S1(A), the logkB-DFT values calculated with the

196

IEFPCM(UFF)/B3LYP/6-311++G(2d,2p) method linearly correlate with the logkB-exp values,

197

and the slope is close to 1. Significantly, for DEHP, the DFT calculated logkB-DFT value is

198

1.77 (mol-1·L·s-1), which is close to the corresponding experimental value [1.71 (mol-1·L·s-

199

1)36].

200

the selected method linearly correlate with the experimentally derived values.48,49 Thus, it

201

can be concluded that the IEFPCM(UFF)/B3LYP/6-311++G(2d,2p) method is suitable for

202

estimating base-catalyzed hydrolysis kinetics of PAEs.

Figure S1(B) also shows the Ea values of the five carboxylic acid esters calculated by

203

It deserves mentioning that for the five carboxylic acid esters, the slope of the regression

204

line (2.38) in Figure S1(B) deviates from 1, implying Ea-DFT overestimated Ea-exp. The

205

deviation of the slope can be caused by fortuitous errors in the Ea-exp values. The Ea-exp values

206

for methyl, ethyl, propyl and butyl acetates were determined in aqueous acetone solutions

207

(62% acetone in water, 25 C), while the Ea-exp value for methyl formate was obtained in pure

208

water (25 C).49 It was concluded that Ea decreases with increasing the concentration of

209

acetone in water,49 which may explain that the Ea-DFT values are higher than the Ea-exp values

210

for methyl, ethyl, propyl and butyl acetates.


10

Page 11 of 32


211

Effects of explicit water molecules on the hydrolysis kinetics were investigated with DEP

212

as a case. Due to the constraint in computational time, 1-3 water molecules were considered.

213

The ΔG‡ values of TS1 and TS2 with consideration of explicit water molecules are listed in

214

Table S4. It can be seen that the ΔG‡ values show a slight decrease with the involvement of

215

0-3 explicit water molecules. Based on the kB-exp values, the experimental ΔG‡ values and

216

errors (differences between the experimental and DFT calculated ΔG‡ values) were

217

calculated (Table S5). It can also be seen that consideration of explicit water molecules did

218

not lead to a trend to decrease the errors. The water molecules were not actively involved in

219

the nucleophilic attack of OH‾ on the ester bonds and the departure of the leaving-groups. It

220

can thus be concluded that explicit water molecules may not have significant influences on

221

the calculated ΔG‡ values, with consideration of the errors. Therefore, only the IEFPCM(UFF)

222

implicit solvation model was employed for calculation of logkB-DFT for the other PAEs so as

223

to save computational time.

224

Hydrolysis pathways of PAEs. The rate-determining steps of the other PAEs were

225

investigated and the ∆G‡ values of TS1 and TS2 are listed in Table S6. For DTDP, BDP_C9

226

side chain, BBP_C7 side chain and DEHP, nucleophilic attack of OH‾ on the ester bonds is

227

the rate-determining step. For the others, departure of leaving-groups is the rate-determining

228

step. It deserves mentioning that for both DMP (the PAE with the shortest side chians) and

229

DTDP (the PAE with the longest side chains), the rate-determining step is nucleophilic attack

230

of OH‾. Therefore, charge distribution in the active sites26 plays a more important role in

231

determining the rate-determining step than steric hindrance for the situation.

232

As shown in Figure 2, the leaving-groups may dissociate from the parents by two ways,

233

protonated (R-OH) or unprotonated (R-O‾). According to the IRC calculation, all the PAEs


11


Page 12 of 32

234

except for OBOP and BMPP follow the pathway of proton transfer and have the protonated

235

(R-OH) leaving-groups. Since hydrogen-bond interactions can influence the process of

236

proton transfer, TS2 of OBOP and BMPP with inclusion of 1 - 2 explicit water molecules

237

were further calculated. When explicit water molecules were added into the system, OBOP

238

and BMPP follow the pathway of proton transfer and have protonated leaving-groups. Duarte

239

et al. investigated mechanistic details of phosphate monoester hydrolysis, also found the

240

hydrolysis follows proton transfer concerted with cleavage of the ester bond, and indicated

241

that the nature of leaving-groups plays an important role in this situation.38,61

242

Kinetics of Base-Catalyzed Hydrolysis. The DFT calculated kB_side chain values for the

243

PAEs under study are listed in Table 1. Note that the kB-DFT values can be calculated as 2kB_side

244

chain

245

different side chains. EPI SuiteTM predicts kB values of PAEs based on group/fragment

246

substituent values generally.43 However, EPI SuiteTM cannot predict kB values of DIOP,

247

OBOP, DINP, DIUP, DIHP and BMPP, due to its limit of the applicability domain. EPI

248

SuiteTM gives a same kB value for some PAEs. For example, the kB value predicted by EPI

249

SuiteTM is 1.427 × 10-2 (mol-1·L·s-1) for DNOP, DTDP, ODP and DUP; 2.058 × 10-2 (mol-

250

1·L·s-1)

251

DPP and DNPP. The DFT and kinetics calculation methods reported in this study can well

252

characterize the variability of the base-catalyzed hydrolysis kinetics of the PAEs.

for PAEs with two same side chains or as (kB_side chain 1 + kB_side chain 2) for those with two

for DIBP, DEHP and DPHP; and 3.204 × 10-2 (mol-1·L·s-1) for DBP, DNHP, DIDP,

253

Table 1. DFT calculated kB_side chain values for PAEs at the IEFPCM/B3LYP/6-

254

311++G(2d,2p) level Name

kB-side chain (mol-1·L·s-1)

Name

kB-side chain (mol-1·L·s-1)

DMP

6.24 × 102

DINP

8.91 × 10-1


12

Page 13 of 32


DEP

7.84 × 10-2

DIDP

5.00 × 10-2

DPP

8.05 × 10-3

DUP

3.44 × 100

DAP

5.27 × 10-2

DIUP

1.66 × 10-1

DBP

1.84 × 10-1

DTDP

7.50 × 103

DIBP

1.97 × 10-4

DCP

1.48 × 10-2

DNPP

3.24 × 10-2

BDP_C7

1.24 × 10-1

BMPP

6.28 × 10-2

BDP_C9

3.11 × 100

DNHP

8.46 × 10-2

ODP_C7

1.97 × 100

DIHP

2.77 × 10-2

ODP_C9

1.16 × 100

DEHP

2.95 × 101

BCP_C7

2.39 × 10-2

DIOP

1.65 × 10-2

BCP_C9

5.32 × 10-2

DPHP

6.84 × 10-4

BBP_C7

3.15 × 100

DNOP

2.05 × 10-1

BBP_C9

5.83 × 100

255 256

In natural aquatic environments, pH is normally between 5 and 9.62 pH of seawater is

257

normally between 8.0 and 8.4, while organic activity or other causes can result in temporary

258

pH values from 6 to 10 locally.62 Based on the kB-DFT values, t1/2 of PAEs at different pH was

259

calculated from t1/2 = ln2/(kB-DFT∙[OH‾]). At pH = 7, the t1/2 values range from 0.13 hours to

260

558 years, and for most PAEs are several years (Table S7). At pH = 9, the t1/2 values range

261

from 0.001 hours to 5.6 years, and for most PAEs are several days. At pH = 9, the t1/2 values

262

of DIBP and DPHP are 5.58 and 1.61 years. Thus, DIBP and DPHP are resistant to base-

263

catalyzed hydrolysis. The side chains for the molecular structures of DIBP and DPHP are

264

branched alkyls. It can be speculated that the base-catalyzed hydrolysis rates of PAEs with


13


265

Page 14 of 32

branched alkyl side chains are slow.

266

Effects of Side Chains on PAE Hydrolysis. The favorable hydrolysis pathways were

267

investigated for the PAEs with two different side chains. The free energy surfaces for BCP,

268

BBP, BDP and ODP are presented in Figure 4. For BCP, ∆G‡ for the C8 side chain (80.34

269

kJ·mol-1) is lower than that for the C7 side chain (82.33 kJ·mol-1). For BBP, ∆G‡ for the C9

270

side chain (68.69 kJ·mol-1) is also lower than that for the C7 side chain (70.22 kJ·mol-1).

271

Thus, it can be concluded that the alkaline hydrolysis is more favorable for the cyclic side

272

chains than for the linear side chains. Su et al. experimentally investigated alkaline hydrolysis

273

of organophosphate (OP) triesters, and also found that the OP triesters with alkyl moieties

274

are more stable than OP triesters with aryl moieties.35

275 276

Figure 4. Schematic free energy surfaces for base-catalyzed hydrolysis of the PAEs

277

calculated at the IEFPCM/B3LYP/6-311++G(2d,2p) level [The total energy of the reactants

278

and OH‾ is set to zero (reference state). The symbols “R, IM1, TS, IM2, P” refer to reactants,

279

pre-reactive complexes, transition states, intermediates and products involved in the reaction,


14

Page 15 of 32

280


respectively]

281

As shown in Figure 4, for BDP, the nucleophilic attack of OH‾ is easier to occur at the C7

282

side chain (∆G‡ = 78.24 kJ·mol-1) with butyl as the leaving-group, but ∆G‡ for the hydrolysis

283

pathway of the C7 side chain is higher than that for the C9 side chain (70.26 kJ·mol-1) with

284

decyl as the leaving-group. For ODP, the nucleophilic attack is easier to occur at the C9 side

285

chain (∆G‡ = 72.70 kJ·mol-1) with decyl as the leaving-group, while ∆G‡ for the hydrolysis

286

pathway of the C9 side chain is slightly higher than that for the C7 side chain (71.37 kJ·mol-

287

1)

288

hydrolysis pathways do not show a generic preference over the side chains. OBOP also has

289

two different side chains. However, the ΔG values for the two hydrolysis pathways are > 0.

290

Thus, base-catalyzed hydrolysis of OBOP is not thermodynamically feasible and was not

291

considered.

with octyl as the leaving-group. Thus, for BDP and ODP with two linear side chains, the

The leaving-groups of both DBP and DIBP have four C atoms. DIBP [kB-DFT = 3.94 × 10-

292 293

4

(mol-1·L·s-1)] with branched side chains is more stable than DBP [kB-DFT = 3.69 × 10-1 (mol-

294

1·L·s-1)]

295

other PAEs, i.e. the stability of the PAEs has the order: BMPP > DNHP (6 C atoms), DIOP >

296

DNOP (8 C atoms), DPHP > ODP_C9 > BDP_C9 (10 C atoms), DIUP > DTDP (13 C atoms).

297

It can be concluded that the PAEs with branched alkyl side chains are more resistant to base-

298

catalyzed hydrolysis than the PAEs with linear alkyl side chains, and the PAEs with cyclic

299

side chains tend to have faster base-catalyzed hydrolysis kinetics than the PAEs with linear

300

alkyl side chains. The conclusion can be explained with goodness or electronegativity of

301

leaving-groups. For example, phenylmethanol has smaller acidity constant (pKa) and larger

302

electronegativity, while alkyl alcohols have larger pKa and smaller electronegativity.26,34 In

with linear side chains (Table S8). The similar trend can also be observed for the


15


Page 16 of 32

303

principle, the smaller is pKa of leaving alcohols (R-OH), the larger is electronegativity of R-

304

OH, and the easier is the dissociation of R-O‾ from PAEs.26,34

305

QSAR Modeling. The DFT results indicate nucleophilic attack of OH‾ on the ester bonds

306

(TS1) is the rate-determining step for 5 cases (DMP, DEHP, DTDP, BDP_C9 side chain and

307

BBP_C7 side chain), while departure of leaving-groups (TS2) is the rate-determining step for

308

the other PAEs. To develop QSAR models with a clear mechanism to predict logkB_side chain,

309

a discriminant model was firstly developed. Since the rate-determining step can be

310

distinguished by the differences between ΔG‡TS1 (energy difference between TS1 and

311

reactants) and ΔG‡TS2 (energy difference between TS2 and reactants), a model for predicting

312

ΔG‡TS2 – ΔG‡TS1 was developed: ΔG‡TS2 – ΔG‡TS1 = 26.4D1 – 16.1D2 + 11.8D3 – 8.43

313

(3)

314

ntra = 23, R2tra = 0.891, Q2LOO = 0.841; next = 5, R2ext = 0.945, Q2ext = 0.848

315

where D1 ~ D3 represent the different molecular structural descriptors listed in Table 2; ntra

316

and next represent the number of the endpoint values in the training and validation datasets,

317

respectively; R2 represents the determination coefficient; Q2LOO represents the leave-one-out

318

cross-validated determination coefficient; and Q2ext represents the external validation

319

coefficient.

320

Table 2. Molecular structural descriptors involved in the models, their variable inflation

321

factor (VIF), t-statistics, and significance level p values. Descriptors D1

Depiction

VIF

Hirshfeld charge on the O atom of leaving-groups. 1.01

t

p

12.3

Development of Prediction Models on Base-Catalyzed Hydrolysis

Recommend Documents