Prediction of Michael-Type Acceptor Reactivity ... - ACS Publications

Sep 30, 2010 - Mark T. D. Cronin,† and Gerrit Schüürmann*,‡,§. School of ... for toxicity, Michael acceptors, covers a wide range of reactivity...
0 downloads 0 Views 373KB Size
1576

Chem. Res. Toxicol. 2010, 23, 1576–1585

Prediction of Michael-Type Acceptor Reactivity toward Glutathione Johannes A. H. Schwo¨bel,† Dominik Wondrousch,‡,§ Yana K. Koleva,† Judith C. Madden,† Mark T. D. Cronin,† and Gerrit Schu¨u¨rmann*,‡,§ School of Pharmacy and Chemistry, LiVerpool John Moores UniVersity, Byrom Street, LiVerpool, L3 3AF, England, UFZ Department of Ecological Chemistry, Helmholtz Centre for EnVironmental Research, Permoserstr. 15, 04318 Leipzig, Germany, and Institute for Organic Chemistry, Technical UniVersity Bergakademie Freiberg, Leipziger Strasse 29, 09596 Freiberg, Germany ReceiVed May 19, 2010

A model has been developed to predict the kinetic rate constants (kGSH) of R,β-unsaturated Michael acceptor compounds for their reaction with glutathione (GSH). The model uses the local charge-limited electrophilicity index ωq [Wondrousch, D., et al. (2010) J. Phys. Chem. Lett. 1, 1605-1610] at the β-carbon atom as a descriptor of reactivity, a descriptor for resonance stabilization of the transition state, and one for steric hindrance at the reaction sites involved. Overall, the Michael addition model performs well (r2 ) 0.91; rms ) 0.34). It includes various classes of compounds with double and triple bonds, linear and cyclic systems, and compounds with and without substituents in the R-position. Comparison of experimental and predicted rate constants demonstrates even better performance of the model for individual classes of compounds (e.g., for aldehydes, r2 ) 0.97 and rms ) 0.15; for ketones, r2 ) 0.95 and rms ) 0.35). The model also allows for the prediction of the RC50 values from the Schultz chemoassay, the accuracy being close to the interlaboratory experimental error. Furthermore, kGSH and associated RC50 values can be predicted in cases where experimental measurements are not possible or restricted, for example, because of low solubility or high volatility. The model has the potential to provide information to assist in the assessment and categorization of toxicants and in the application of integrated testing strategies. Introduction Electrophilic reactivity is important for understanding and describing covalent interactions between xenobiotic toxicants and biological macromolecules (1). A growing part of developing alternatives for the toxicological assessment of chemicals, including in chemico (experimental measurement of reactivity) and in silico (computational) approaches, involves the understanding of reactive toxicity. It has long been understood that compounds capable of covalent interactions with biological macromolecules are liable to elicit toxic effects (2). For instance, there have been recent improvements in our understanding of the reactions involved between xenobiotics and immunoproteins that result in allergic contact dermatitis (skin sensitization) (3, 4) as well as for the irreversible reaction with fish gill membranes (5) that will promote acute toxicity in excess of narcosis. Many of these mechanisms of action are electrophilic in nature. There are many further examples of electrophilic toxicity, e.g., respiratory sensitization, liver toxicity, skin irritation, etc. (6). There are a number of reactive mechanisms that can be established from organic chemistry. One of the most important for toxicity, Michael acceptors, covers a wide range of reactivity, from being so rapid that they need special stopped-flow equipment so that they may be followed kinetically to requiring elevated temperatures to proceed with any form of reaction. In this context, many R,β-unsaturated compounds acting via Michael-type reactions are considered to be particularly reactive * To whom correspondence should be addressed. E-mail: gerrit. [email protected]. Telephone: +49-341-235-1262. Fax: +49-341-2351785. † Liverpool John Moores University. ‡ UFZ Department of Ecological Chemistry. § Technical University Bergakademie Freiberg.

Figure 1. Michael addition domain. R,β-unsaturated carbonyl or carboxyl groups with double or triple bonds, including cyclic systems (in addition, electron-withdrawing groups for “Michael-type” electrophiles might include pyridino heterocyclic groups and, not considered in this paper, CN, SO2R, and NO2 groups) (4).

(7). They are capable of forming irreversible bonds with biological macromolecules, such as proteins or DNA molecules (8). This covalent reactivity is known to result in a number of toxic effects, including elevated acute (aquatic, terrestrial, and human) toxicity, allergenic reactions, mutagenicity, and carcinogenicity. Compounds with an R,β-unsaturated carbonyl or carboxyl fragment are very important industrial substances, for example, in the manufacture of polymers, textiles, or auxiliary materials in medicine (9, 10). The types of chemical fragments associated with the Michael acceptor mechanism define its domain. Michael acceptors are characterized as having carbon-carbon double (CdC) or triple (CtC) bonds with electron-withdrawing substituents (see Figure 1). This results in a polarizable electron density at the π-bond, where the β-carbon atom is positively polarized and becomes the preferred site of an attack for soft nucleophiles, for example, glutathione (GSH) (7). The tripeptide glutathione (GSH, L-γ-glutamyl-L-cysteinylglycine) is one of the most widely used nucleophilic reference

10.1021/tx100172x  2010 American Chemical Society Published on Web 09/30/2010

Michael Addition ReactiVity

Figure 2. Structures of methyl propiolate (1), methyl acrylate (7), methyl crotonate (50), and methyl methacrylate (52) in order of decreasing reactivity. Compound numbers from Table 1.

molecules in reactivity assays. Direct GSH depletion can be measured in a pure chemical reactivity assay. GSH is reactive toward soft electrophiles, for example, polarized R,β-unsaturated compounds, which act predominantly as Michael-type acceptors. Soft molecules are readily polarizable. In contrast, Schiff base formers react with GSH only in the absence of amino groups (11). The thiol group of GSH is taken as a surrogate for protein thiol groups as a typical site of covalent attack of soft electrophiles. Moreover, GSH itself is the most prevalent cellular thiol and the most abundant low-molecular weight peptide in cells (12). It protects cells by detoxifying electrophilic compounds and acts as an antioxidant. The concentration of GSH is depleted during attack by electrophilic compounds, commonly by alkylation mechanisms. A high GSH depletion rate makes other endogenous thiol groups susceptible to attack, especially soft cysteine SH moieties. Therefore, if the GSH concentration falls below a critical level in a cell, this can cause an accumulation of damage (13). Chemical reactivity assays have the potential to provide information to assist in the assessment of toxicants. Their use is foreseen in integrated testing strategies (ITS) (14), for instance, in complying with the needs of the European Union’s REACH regulation (15). Various studies have shown good qualitative relationships between GSH reactivity of chemicals acting as Michael-type acceptors and biological end points, e.g., skin sensitization (12), acute fish toxicity (13), or toxicity to Tetrahymena pyriformis (16). For example, the differential toxicity of methacrylate, acrylate, and crotonate can be explained by their varying reactivity, although the structures are very similar. Methacrylate is poorly reactive; acrylate is very reactive, and crotonate lies somewhere between them (structures shown in Figure 2) (17). There is current interest in using in silico techniques to predict potential in chemico reactivity (14). This is an attractive approach, especially in cases in which experimental measurements are not possible or restricted, for example, because of low solubility or high volatility. The aim of this investigation was to construct a quantitative in silico model for kinetic rate constants of the Michael addition reaction. As experimental reference values, data from chemical reactivity assays for the reaction of GSH with R,β-unsaturated xenobiotics in an aqueous phosphate buffer solution were taken. Apart from reaction rate constants (16, 18–20), results from the 2 h spectrophotometric assay developed by Schultz and coworkers were used as reference values for parametrization, providing values of the 50% reaction concentration (RC50) as a quantification of the electrophilic reactivity toward GSH (11, 17). Recently, the advantages and disadvantages of the Schultz assay have been discussed in comparison to a kinetic GSH chemoassay that also accounts for oxidative GSSG formation (16). Moreover, a test set of RC50 values could be reproduced with an interlaboratory accuracy of ∼0.3 log unit. However, the static RC50 assay appears to be limited by water solubility and a lack of sensitivity to the most reactive chemicals. Schultz et al. incorporated some compounds with limited water solubility into the assay by applying a water/methanol mixture instead of an aqueous solution (21). It should be noted that water solubility limits the maximum amount available for reaction. Thus,

Chem. Res. Toxicol., Vol. 23, No. 10, 2010 1577

compounds with similar kinetic rate constants, but different water solubilities, would have different RC50 values. Therefore, only fully soluble compounds were considered with regard to the RC50 values used within this publication. The model developed here considers both electronic and steric effects, with the descriptors based on easily computed ground-state properties, performed by quantum chemical calculations.

Materials and Methods Data Set. Experimental kinetic rate constants have been collected for the reaction between R,β-unsaturated compounds and glutathione (GSH) from various literature sources. Data for 66 Michael-type acceptor compounds are included in the data set. Of these, secondorder kinetic rate constants (kGSH) were known for 45 compounds, and for 21 additional compounds RC50 values from the Schultz assay could be converted into kGSH values (see below). The compounds are listed in Table 1. The data set contains R,β-unsaturated compounds with a double bond or a triple bond, cyclic systems, and those with or without substituents at the R-carbon. The compound classes included aldehydes, ketones, amides, and esters. The kinetic GSH assay (16) for measuring rate constants for electrophilic Michael acceptor-type compounds was performed in an aqueous phosphate buffer solution at 25 °C and pH 7.4. At certain time intervals, the concentration of free thiol groups was measured. Free thiol groups were quantified by UV-vis spectroscopy at 412 nm after reaction with the chromophore 5,5′-dithiobis(2nitrobenzoic acid) (DTNB). However, for some compounds in the data set, slightly different experimental conditions were used. Those values were corrected to correspond with the standard assay, using the following equations. For any experimental pHexp value other than pH 7.4, we used the following equation, with the pKa value for GSH being 8.56 (22), to convert the experimental rate constant k into the value at pH 7.4:

(

log k(pH 7.4) ) log k(pHexp) - log

1 + 10pKa-pHexp 1 + 10pKa-7.4

)

(1)

To account for different temperature conditions (Texp between 20 and 37 °C), corrections to 25 °C (298 K) are made by the application of the Arrhenius equation (23):

log k(25 °C) ) log k(Texp) +

(

1 EA 1 1 ln 10 R Texp 298 K

)

(2)

The activation energy EA can be calculated if the variation of experimental rate constant k is known for the same compound at two different temperatures (T1 and T2):

EA ) ln 10 × R × [log k(T1) - log k(T2)] ×

(

1 1 T2 T1

)

-1

(3)

The RC50 value is the millimolar concentration of electrophile that gives a defined half-life (t1/2) for the nucleophile. In other words, the reaction time is fixed and the concentration changes in a way that half of the reaction is completed at the fixed time. If secondorder kinetics are assumed, the following equation can be used to convert RC50 values into rate constants (16):

[

log k ) log

1 0 t1/2(RC50 - cNu )

(

× ln 2 -

0 cNu RC50

)]

(4)

In eq 4, c0Nu is the initial concentration of GSH. In the Schultz assay, 0 cNu is 1.4 mM and the half-life t1/2 is 120 min. Note that we did

1578

Chem. Res. Toxicol., Vol. 23, No. 10, 2010

Schwo¨bel et al.

Table 1. Experimental Rate Constants (kGSH in M-1 min-1 at 25 °C and pH 7.4) and Predicted Rate Constants (eq 8, model B, parameters listed in Table 2) for Adduct Formation with GSH and RC50 Values (in millimolar) of Schultz’s Chemoassay 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

name

CAS

log kGSHa (exptl)

Rem.b

log kGSH (predicted)

log RC50a

log(S/RC50)c

solubled

e

922-67-8 1679-36-3 623-47-2 1629-58-9 625-33-2 1629-60-3 96-33-3 930-30-3 140-88-5 2497-21-4 925-60-0 4312-99-6 1119-44-4 4643-27-0 141-32-2 1120-73-6 1669-44-9 107-02-8 818-61-1 123-73-9 1576-87-0 6728-26-3 2548-87-0 78-94-4 2177-18-6 930-68-7 10477-47-1 764-39-6 100-69-6 505-57-7 5362-56-1 142-83-6 5910-85-0 999-55-3 14861-06-4 104-55-2 868-77-9 1423-60-5 2579-22-8 1817-57-8 23326-27-4 4341-76-8 111-12-6 16205-90-6 141-79-7 565-62-8 14309-57-0 5392-40-5 2758-18-1 623-43-8 923-26-2 80-62-6 623-70-1 97-63-2 6622-76-0 18829-55-5 18829-56-6 5910-87-2 557-48-2 25152-84-5 97-86-9 96-05-9 106-63-8 2499-95-8 25584-83-2 79-06-1

2.07 1.90 2.02 3.10 1.43 3.07 1.06 1.41 1.03 1.38 1.01 3.03 1.10 1.42 0.93 -0.70 1.06 2.12 1.43 1.68 1.43 1.21 1.12 2.14 1.95 1.24 1.54 0.89 0.32 0.90 0.74 0.59 0.72 0.23 0.14 0.78 -0.76 2.66 2.66 2.48 0.68 0.67 1.43 0.62 -0.68 -0.11 1.03 -0.23 -1.13 -0.79 -0.03 -1.14 -0.79 -1.24 -2.15 1.23 1.07 0.71 1.30 0.87 -0.44 -0.04 1.23 0.90 1.43 -0.88

kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kRC50 kexp kexp kexp kexp(pH) kexp kexp kexp(pH) kexp kexp kexp kexp kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(pH) kexp(T) kexp(T)

2.12 1.87 1.93 2.89 1.26 2.91 1.11 1.44 0.92 0.92 0.89 2.93 0.77 1.07 0.86 -0.47 0.47 1.97 1.54 1.39 1.19 1.11 1.10 2.87 1.40 1.19 1.66 1.04 -0.02 1.05 0.93 0.85 0.90 1.09 0.40 1.38 -0.33 2.60 2.60 2.60 1.02 0.87 0.83 0.74 -0.15 -0.29 0.72 -0.17 -0.52 -1.03 0.18 -0.73 -1.17 -0.88 -2.33 1.14 1.15 0.81 1.18 0.87 -0.82 -0.71 0.99 0.83 1.00 0.06

-0.85 -0.92 -0.89 -1.29 -0.68 -1.28 -0.30 -0.24 -0.28 -0.47 -0.31 -0.96 -0.37 -0.66 -0.26 0.92 -0.20 -1.07 -0.57 -0.66 -0.48 -0.40 -0.55 -1.05 -0.96 -0.43 -0.68 -0.11 0.45 -0.12 0.04 0.18 0.06 0.54 0.63 0.00 1.52 -1.24 -1.60 -1.30 0.10 0.11 -0.59 -0.02 1.50 1.00 -0.36 0.31 1.43 1.51 1.85 1.85 1.54 NRf NRf

3.95 3.57 3.47 3.71 3.42 3.17 3.06 2.83 2.46 2.43 1.81 1.81 1.81 1.58 1.45 1.35 1.12 4.65 4.51 4.07 2.73 2.13 1.24 3.99 3.12 3.01 2.87 2.36 1.97 1.85 1.75 1.75 1.35 1.10 1.07 1.03 1.37 4.98 3.32 2.98 2.44 1.90 1.04 1.00 0.97 0.91 0.76 0.63 0.59 0.58 0.58 0.33 0.03 -

yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes yes partly partly partly partly partly partly partly partly partly no no

methyl propiolate 3-hexyn-2-onee ethyl propiolatee 1-penten-3-one 3-penten-2-one 1-hexen-3-one methyl acrylate 2-cyclopenten-1-one ethyl acrylate 4-hexen-3-one n-propyl acrylate 1-octen-3-one 3-hepten-2-one 2-octen-4-one n-butyl acrylate 2-methyl-2-cyclopenten-1-one 3-octen-2-one acrolein 2-hydroxyethyl acrylate (E)-crotonaldehyde trans-2-pentenal trans-2-hexenal trans-2-octenal methyl vinyl ketone vinyl acrylate 2-cyclohexen-1-one propargyl acrylate 2-pentenal 2-vinylpyridine 2-hexenal 4-methyl-2-pentenal 2,4-hexadienal 2,4-heptadienal allyl acrylate vinyl crotonate cinnamic aldehyde 2-hydroxyethyl methacrylate 3-butyn-2-onee phenylpropargyl aldehydee 4-phenyl-3-butyn-2-onee methyl tetrolatee ethyl 2-butynoatee methyl 2-octynoatee ethyl hex-2-ynoatee mesityl oxide 3-methyl-3-penten-2-one 3-nonen-2-one 3,7-dimethyl-2,6-octadienal (Citral) 3-methyl-2-cyclopenten-1-one methyl crotonate 2-hydroxypropyl methacrylate methyl methacrylate ethyl crotonate ethyl methacrylate methyl tiglate trans-2-heptenal trans-2-nonenal trans,trans-2,4-nonadienal trans-2-cis-6-nonadienal trans,trans-2,4-decadienal isobutyl methacrylate allyl methacrylate isobutyl acrylate hexyl acrylate hydroxyl propyl acrylate acrylamide

a Values taken from literature (2, 7, 11–13, 16, 17, 20, 21, 46–49). b kexp is rate constant kGSH for standard conditions (pH 7.4 the kGSH recalculated for pH 7.4 via eq 1. kexp(T) is the kGSH recalculated for 25 °C via eq 2. kRC50 is the kGSH based on the RC50 available, experimental solubility values (S, millimolar units) taken from the ChemIDPlus (50) database; otherwise calculated d Fully soluble in Schultz’s (11) RC50 assay [log(S/RC50) g 1]. e Acetylenic Michael acceptor, calculated with model B2; all other with model B1. f Not reactive at the limit of saturation.

not include RC50 values if the water solubility was lower than, or on the same order of magnitude as, the particular RC50 value.

and 25 °C). kexp(pH) is value (eq 3). c Where with EPISUITE (51). compounds calculated

Computational Details. The geometries of all molecules were optimized at the density functional theory B3LYP/6-31G** level,

Michael Addition ReactiVity

Chem. Res. Toxicol., Vol. 23, No. 10, 2010 1579

Figure 3. Michael-type addition. Reaction mechanism of an R,β-unsaturated compound in the reaction with glutathione (GSH) in an aqueous phosphate buffer solution: (a) GSH with the Michael-type acceptor, (b) carbanionic transition state (see Figure 4), (c) enolic intermediate, and (d) peptide adduct complex.

Figure 4. Resonance stabilization of the Michael-type addition transition state and solvation-mediated participation of the R-carbon atom in the transition-state (TS) structure (actually taking place in an aqueous phosphate buffer solution). Dashed lines indicate developing or vanishing bonds of the resonance structures representing the TS, and dotted lines represent hydrogen bonds between the TS and solvating H2O and H3O+.

employing Gaussian 03 (24). For all compounds, frequency analysis was conducted to confirm that correct ground-state geometries had been obtained. Natural bond orbital (NBO) analysis methods were applied by means of the NBO 3.1 implementation in Gaussian 03 (25, 26). For the calculation of the local electrophilicity and effective electron donor and acceptor energies [originally introduced by Klamt (27), and subsequently reformulated (28–30)] and optimization of the nonlinear parameters, we used an in-house. NET program. Parameters for steric accessibility were calculated with the freely available MSMS software (31), using published Connolly radii (32). Multilinear regression analysis was performed using the Minitab 15.1 statistical software. Michael Addition Reactivity Model. The aim of the method described is to predict rate constants of the Michael addition reaction, based on ground-state properties of a single molecule. To estimate those rate constants in a physically reliable way, we took the following considerations into account. In the chemical reaction between the peptide GSH and an R,βunsaturated electrophilic compound, it is probable that multiple steps are involved, all of which have an influence on the overall rate constant of Michael addition reactions. One important step is the formation of an adduct complex. GSH binds to the β-carbon atom of an R,β-unsaturated compound (see Figure 3a). Therefore, a sound reactivity descriptor is required that models the ability to accept electrons at the β-carbon atom reliably. The reactive species of GSH is mostly probably the negatively charged thiolate ion GS-, keeping in mind that the additional lone pair and excess charge increase the nucleophilic power significantly. Note that under experimental conditions the rate constant increases with higher pH values, reflecting the greater amount of thiolate as a reactive species. A common method applied in predictive toxicology to describe the reactivity of an electrophilic compound is the molecular electrophilicity index ω (33). It can be calculated by the following equation:

ω)

EHOMO2 + 2EHOMOELUMO + ELUMO2 µ2 ) 2η 4(ELUMO - EHOMO)

(5)

where µ is the molecular chemical potential and η the molecular hardness, both of which can be replaced by the energies of the

highest occupied molecular orbital, EHOMO, and the lowest unoccupied molecular orbital, ELUMO, via application of Koopmans’ theorem (34). The descriptor ω has been widely used. Domingo et al. (35, 36) calculated the electrophilicity index for a series of activated ethylenes. Enoch et al. (37, 38) described a method to rank the reactivity of Michael acceptors by the electrophilicity index, using this ranking for a read-across approach to predict skin sensitization and respiratory sensitization capabilities of compounds. In a related study, the electrophilic index was used to model the cytotoxicity of a series of sugars acting via Michael addition (39, 40). (Note that some definitions of hardness lack the factor of 1 /2; therefore, versions of electrophilicity with a prefactor of 8 in the denominator can be found in the literature, instead of 4.) As a site-specific reactivity descriptor is desirable in this study for the prediction of reactivity, a local counterpart of the electrophilicity index was applied. In principle, this can be provided by multiplying the molecular electrophilicity index ω with an atomcentered Fukui index f +, as described by Chattaraj et al. (34). This is a method of obtaining the most reactive site within a molecule or describing trends in a qualitative manner. Because of the dependence on the population analysis scheme used for the calculation of the Fukui index f + (41), we have followed a different route. We have defined the local charge-limited electrophilicity index ωq(r) at an atomic reaction site r in the following way (41):

ωq(r) ) [EQocc(qocc, r)2 + 2EQocc(qocc, r)EQvac(qvac, r) + EQvac(qvac, r)2]/{4[EQvac(qvac, r) - EQocc(qocc, r)]}

(6)

where EQocc(qocc,r) is the effective (charge-limited) electron donor energy and EQvac(qvac,r) the analogous electron acceptor energy. Moreover, qocc and qvac represent the (typically fractional) charges, the release (qocc) or acceptance (qvac) of which is characterized by EQocc and EQvac, respectively. Both functions are site-specific, located at the β-carbon. These, and related, descriptors have been used in earlier studies to predict rate constants of indirect photolysis by •OH radicals (29, 42), and the ability to form hydrogen bonds (28, 30, 43). The exact definitions for the electron donor and acceptor energies are given below. Furthermore, we have defined q the maximum electrophilicity index ωmax that accounts for the empirical observation that beyond a certain (calculated) level of electrophilicity, a further increase does not translate into increased (experimental) rate constants. This finding suggests further exploration of the underlying reason is needed in future studies. In a previous publication, it was shown that the site-specific definition of eq 6 captures intrinsic reactivity better, when compared to the use of global frontier orbitals (41). In particular, local electrophilicity alone could well explain intrinsic reactivity within compound classes. However, an indicator variable was needed to distinguish between ketones and esters. This indicator variable should account for the difference in reactivity caused by the -COR and -COOR activating groups, respectively (i.e., contributions not being reflected by the particular atomic orbitals leading to ωq), and can be explained by resonance stabilization. As shown below, the indicator variable can be traced back to reflecting differences in a distinct (calculated) bond energy that is now used instead as a model parameter. We have assumed that resonance stabilization of the

1580

Chem. Res. Toxicol., Vol. 23, No. 10, 2010

Schwo¨bel et al.

(nominally) carbanionic transition state plays an important role in the reaction pathway of formation of the GSH-electrophile complex (shown in Figure 3, the reaction from a to b should involve the rate-limiting step). The actual reaction takes place in water where the reactive nucleophile is presumed to be the thiolate ion, and the transition state (TS) to be stabilized through solvation in the aqueous buffer solution (see below). Apart from solvation, the TS appears to be stabilized electronically through resonance between a reactantlike and product-like mesomeric form as shown in Figure 4. It follows that the closer in energy both TS resonance forms are, the better the TS resonance stabilization should be, and hence the higher the reaction rate. In the reactant-like TS resonance form (Figure 4a), there is, formally, a single bond between the ipso carbon atom (the carbonyl or carboxyl C atom) and the R-carbon atom. In the product-like (right-hand) resonance form (Figure 4b), this bond has turned into a double bond. Note further that the latter resonance form is likely to have a higher energy than the former, considering the analogy in structural change between the reactant-like and product-like resonance form and keto-enol tautomerism. We now assume specifically that the higher the energy of the original σ-bond, the closer in energy both resonance forms become, and thus the easier it is to convert the single bond to a double bond. Note that in the representation of local orbital interaction, the initial C-C σ-bond becomes stabilized upon formation of the additional π-bond (that is higher in energy than the σ-bond), and this energetic lowering is larger if its original energy level is higher. Moreover, any solutionphase proton with vacant orbitals solvating the reaction site region of CR (see Figure 4) and Cβ might interact better with the higher energy levels of the occupied orbital of the C-C σ-bond. The energy level Eσ(C-C) of this σ-bond can be calculated by Weinhold’s natural bond orbital (NBO) analysis (25, 26). Accordingly, we have quantified the original C-C σ-bond energy through respective NBO calculations in the ground-state geometry of the compounds. For the different activating groups, -COR (ketone), -COOR (ester), -CHO (aldehyde), -CONH2 (acrylate), and -C5H4N (pyridine moiety), the average NBO energies Eσ(CR-Cβ) are -17.83, -18.56, -18.29, -17.69, and -17.60 eV, respectively, showing only small variations typically of Pocc(i, r) + pi(r)

qocc - Pocc(i, r) if Pocc(i, r) e qocc e Pocc(i, r) + pi(r)

0 wk(qvac,r) analogous

if qocc < Pocc(i, r) (10)

The electron population pi of an occupied molecular orbital i or the electron space pk available in an unoccupied molecular orbital k is quantified by taking twice the sum of the square of the coefficients cµi (cσk) of the LCAO-MO approach (linear combination of atomic orbitals to molecular orbitals) at an atomic center r:

Michael Addition ReactiVity

pi(r) ) 2

∑ (cµi)2;

Chem. Res. Toxicol., Vol. 23, No. 10, 2010 1581

pk(r) ) 2

µ(r)

∑ (cFk)2

(11)

F(r)

The partial charge of selected occupied or vacant MOs i, k at an atomic center r is i+1

Pocc(i, r) )



k-1

pj(r);

Pvac(k, r) )

j)HOMO



pl(r)

l)LUMO

(12) Results and Discussion The purpose of this investigation was to develop and validate a model based on local molecular parameters for the reactivity of toxicologically important compounds able to react by the Michael addition mechanism. Model Calibration. The Michael addition reactivity model has been applied and tested for the compounds listed in Table 1. The model was calibrated initially for all 66 compounds, where experimental kinetic rate constants for GSH depletion are known or can be calculated from their RC50 values. 2-Vinylpyridine was modeled in the same way as the carbonyl and carboxyl compounds, with the nitrogen as the electronwithdrawing group and the vinyl substituent representing the R- and β-carbon atoms. Because the descriptor ωq contains three nonlinear parameters (namely, the limiting charge amount qocc of the descriptor EQocc for occupied MOs, the equivalent limiting charge amount qvac q of the descriptor EQvac for virtual MOs, and the ωmax value), those were obtained in an iterative way by searching for a minimum of the root-mean-square errors (rms) between experimental and model values (eq 6). In each iteration step, the q were varied slightly, nonlinear parameters qocc, qvac, and ωmax followed by a subsequent multilinear calibration of the variables a, b, c, d, and e. On the basis of the optimized parameters, a regression analysis was performed using the Minitab statistical software. The parameters for four different models are listed in Table 2. The first model (model A) includes only compounds with double bonds between the R-carbon and the β-carbon and is based only on measured rate constants (kexp). The second model (model B1) for double bonds includes both measured rate constants and calculated rate constants based on measured RC50 values in the calibration (kexp and kRC50). The third model (model B2) is intended for compounds with triple bonds. Finally, the fourth model (model C) contains the whole data set of 66 compounds. The regression coefficients in Table 2 follow the expected trend. As indicated by the positive sign for a, reactivity in terms

of the rate constants increases with increasing electrophilicity ωq. The positive sign for b can be interpreted as follows: The higher the energy for the σ-bond between the R-carbon and the ipso carbon, the more likely it is to convert into the π-bond between both atoms stabilizing the transition state. Moreover, any proton with vacant orbitals, which eventually will be added to the R-carbon, might interact better with the higher energy levels of the occupied orbitals at this site. The positive values for the regression coefficients c and d can also be interpreted. The modeled rate constant increases with greater steric accessibility of the reaction site, defined by the steric effects of both the R-carbon and the β-carbon. Although the principle mechanism of Michael addition remains the same, it seems reasonable to parametrize compounds with triple bonds differently, related to slightly different energy levels of the participating MOs. It should be noted that for compounds with triple bonds only the variation in the local electrophilicity index is statistically significant, again proving that the local electrophilicity index is a parameter capable of predicting intrinsic reactivity. It is certain that the parameters accounting for the resonance stabilization of the transition state or steric hindrance could improve the model further, but because of limited experimental data for triple-bond compounds in the data set, it would not be legitimate. It should be noted that model B2, derived from a calibration with 10 compounds, employs one regression coefficient (a), the regression constant (e), and three nonlinear parameters (qocc, q ). As such, model B2 could be called a tentative qvac, and ωmax model, acknowledging that the corresponding parameter:compound ratio of 1:2.5 is below the generally recommended minimum of 1:5. Accordingly, the combined model B is also tentative with regard to the subset of triple-bond compounds. Model Performance. Applying the parameters listed in Table 2 leads to a good correlation between calculated and experimental values, as shown in Figure 5. Information about the data source is also given (e.g., rate constants kexp from the standard assay, kexp from a different assay with varying pH values or temperatures, and rate constants kRC50 from the RC50 assay). Here, model B1 was used for compounds with double bonds and model B2 for compounds with triple bonds. In the following, the combination of model B1 and model B2 is denoted as model B. As one can see from Figure 5, the performance of the model is surprisingly similar for the different data sources. The statistical parameters for the combined overall model B are as follows: number of compounds (n), 66; squared correlation coefficient (r2), 0.91; root-mean-square error (rms), 0.34. Various types of compounds are included in the calibration; the model contains R,β-unsaturated compounds and cyclic systems with and without substituents at the R-carbon. The

Table 2. Variables for the Michael Acceptor Model (eq 8), Fitted with Experimental Rate Constants (kexp) and Calculated Rate Constants, Based on RC50 Values (kRC50) model A

model B1

data set n CR-Cβ

parameter

kexp 42 double bond

kexp and kRC50 56 double bond

model B2

a (eV-1) b (eV-1) c d e (constant)

5.29 ((0.31) 1.46 ((0.20) 2.83 ((0.59) 3.66 ((0.35) 2.55

Multilinear Parameters 4.90 ((0.30) 1.56 ((0.19) 3.17 ((0.59) 3.67 ((0.34) 5.35

-4.88

q ωmax qocc (e) qvac (e)

3.69 1.0 0.5

Nonlinear Parameters 3.69 1.0 0.5

3.33 0.6 0.3

kexp and kRC50 10 triple bond 2.24 ((0.25)

model C kexp and kRC50 66 all compounds 2.69 ((0.26) 0.41 ((0.18) 5.33 ((0.77) 3.56 ((0.49) -8.62 3.39 0.5 0.3

1582

Chem. Res. Toxicol., Vol. 23, No. 10, 2010

Schwo¨bel et al. Table 4. Statistical Performance of Model B for Individual Compound Classes of r,β-Unsaturated Moleculesa compound class

data set

n

r2

rms

me

bias

aldehydes

kexp and kRC50 kexp kexp and kRC50 kexp kexp and kRC50 kexp

18 11 18 15 28 18

0.90 0.97 0.92 0.95 0.91 0.95

0.21 0.15 0.37 0.35 0.35 0.29

0.16 0.11 0.29 0.29 0.28 0.23

0.04 -0.07 -0.05 -0.10 -0.02 -0.10

ketones esters

a No statistics given for amides and heteroaromatics, because there is only one compound in the data set, particularly.

Figure 5. Plot of predicted and experimental rate constants for GSH adduct formation (using model B1 from Table 2 for R,β-unsaturated compounds with double bonds and model B2 for compounds with triple bonds).

compound classes include aldehydes, ketones, amides, esters, and heteroaromatics. This shows good overall model quality; the experimental values have a range of ∼6 orders of magnitude, while the rms error is approximately 1/3 of an order of magnitude (∼7%), which can be considered to be within the range of the experimental error. More details are listed in Table 3. As expected, model A (only compounds with measured kinetic rate constants kexp) is slightly superior to model B1, and the latter does also include values based on Schultz’s static chemoassay (r2 ) 0.93 vs 0.90; rms ) 0.31 vs 0.35). Nevertheless, improvement of the chemical domain vindicates taking values from the static chemoassay. The performance of model C, with all compounds included in the calibration and thus allocating the same regression coefficients to both double-bond and triple-bond compounds, is reasonable and demonstrates the general validity of the assumptions made (r2 ) 0.79; rms ) 0.52), but it is significantly inferior to that of model B. Therefore, separating compounds with double from those with triple bonds is justified. The statistics of model B for individual compound classes are listed in Table 4. The results for R,β-unsaturated aldehydes are excellent, with an r2 of 0.97 and an rms of 0.15 when the predicted values were compared to the measured rate constants. In addition, the model performs well for both ketones and esters (r2 ) 0.95 and 0.95, and rms ) 0.35 and 0.29, respectively). As one can see in Figure 6, there are no significant outliers to the model, neither an individual compound nor a whole compound class. Table 4 confirms that there is no significant bias for individual compound classes. In general, the model fares better with measured, rather than calculated, rate constants (e.g., for esters, r2 ) 0.95 vs 0.91 and rms ) 0.29 vs 0.35). In conclusion, the good model performance, with respect to both the overall model

Figure 6. Plot of predicted vs experimental rate constants for GSH adduct formation for different compound classes (using model B).

and individual compound classes, warrants the recommendation of combined model B. On the basis of structural similarity, the kinetics of the three compounds, ethyl 2-hexynoate (C3H7CtCCO2Et), ethyl 2-butynoate (MeCtCCO2Et), and methyl 2-octynoate (C5H11Ct CCO2Me), should all be very similar, yet the experimental result for methyl 2-octynoate (log kGSH ) 1.43) is up to 6 times more reactive than the former two compounds (0.62 and 0.67, respectively). The predicted results for those three compounds, however, are very close to each other (0.74, 0.87, and 0.83). This might be an indication that the experimental result for the latter is inaccurate and overestimates the reactivity of this compound. If global frontier orbitals were used for the calculation of electrophilicity index ω (eq 5), instead of site-specific ones, ωq(r) (eq 6), the performance of model B worsened (r2 ) 0.85 vs 0.90; rms ) 0.44 vs 0.35). In accordance with the purpose of this study, the site-specific definition is advantageous when considering molecules with more than one possible reaction site, for example, as in ethylene glycol dimethacrylate or in β-damascone (results not shown). Model Validation. To assess the prediction performance, leave-50%-out cross validation was applied. Each data set (of

Table 3. Statistical Performance of Models A-C from Table 2a model A B1 B2b B (B1 and B2) C

data set kexp kexp kexp kexp kexp

and and and and

kRC50 kRC50 kRC50 kRC50

CR-Cβ bond type

n

r2

rms

qCV2

rmsCV

biasCV

double double triple all all

42 56 10 66 66

0.93 0.90 0.91 0.91 0.79

0.31 0.35 0.25 0.34 0.52

0.91 0.86 0.87 0.87 0.76

0.37 0.41 0.30 0.40 0.55

0.00 -0.01 -0.01 -0.01 0.00

a The statistical parameters are as follows: n, number of compounds; r2, Pearson squared correlation coefficient (calibration); rms, root-mean-square error (calibration); qCV2, predictive squared correlation coefficient (cross validation); rmsCV, root-mean-square error (cross validation); biasCV, bias (cross validation). b The number of compounds (n ) 10) is limited compared to the number of parameters (one linear parameter, the regression constant, and three nonlinear parameters); therefore, nonlinear parameters listed in Table 2 were not recalibrated in cross validation.

Michael Addition ReactiVity

Chem. Res. Toxicol., Vol. 23, No. 10, 2010 1583

models A, B1, B2, and C) was divided into two subsets, both covering (almost) the whole target value range. This was achieved in the following way. The compounds were ordered according to their experimental log kGSH values, followed by their alternate allocation to the first and second subset, respectively. Subsequently, calibration of eq 8 was performed individually for both subsets. Finally, the trained model of one subset was used to predict all values of the second subset and vice versa. The predictive squared correlation coefficients of the cross validation (qCV2 in Table 3) were found to be close to particular r2 values of the calibration, and root-mean-square errors of the cross validation (rmsCV) are close to rms values of the calibration. Moreover, the calibration results of the training subsets are similar to those for the complete data set [for example, r2 ) 0.96 (first subset of model A; n ) 21), r2 ) 0.92 (second subset; n ) 21), and r2 ) 0.93 (data set of model A; n ) 42); rms ) 0.27, 0.30, and 0.31, respectively]; we can therefore conclude that the model appears to be statistically robust, reflecting a sufficiently balanced composition of the data set. Individual predictive squared correlation coefficients (q2, as defined in ref 45) of particular test subsets indicate a good prediction quality [for example, q2 ) 0.93 (first subset of model A, calibrated by the second subset), and q2 ) 0.87 (second subset, calibrated by the first subset)], which are close to r2 values of corresponding training subsets. We also tested the stability of the model when excluding methyl tiglate from the calibration set, because its experimental value (kGSH ) -2.15) is out of the range of the experimental values of the other compounds in the data set (from -1.24 to 3.10). Neither the r2 values [0.93 (n ) 41) vs 0.93 (data set of model A; n ) 42)] nor the rms values (0.31 vs 0.31) changed following the exclusion of this compound, and the regression coefficients remained within the calculated error range {a ) 5.28 (n ) 41) [5.29 (n ) 42)]; b ) 1.48 (1.46); c ) 2.95 (2.83); d ) 3.70 (3.66); e ) 2.73 (2.55)}. Therefore, it is justified to include methyl tiglate in the calibration set. These findings demonstrate that our model appears to be robust with a good predictive power. Prediction of RC50 Values. Considering eq 4 (16), the model should be able to predict the RC50 values reported from the Schultz chemoassay, for compounds in the intermediate reactivity range. If we apply model B to 44 compounds, whose experimental RC50 values are known, the following regression equation is obtained:

log RC50 ) -0.67((0.06) × log kGSH(modeled) + 0.49((0.09)

(13)

For the sake of completeness, if experimental rate constants from Table 1 were used (instead of modeled ones), the regression equation is

log RC50 ) -0.68((0.03) × log kGSH(experimental) + 0.51((0.05) (14) Equation 14 is very similar to the model equation. The rootmean-square error of eq 13 is 0.32, which is close to the interlaboratory experimental error of 0.3 (16). Estimated (eq 13) and experimental RC50 values are shown in Figure 7. Note that the slope and intercept of eq 13 (and of eq 14) deviate from those of a previously published regression equation that had been confined to R,β-unsaturated carbonyls in the intermediate reactivity range (16); here, all fully soluble compounds have been included in the regression.

Figure 7. Predicted RC50 values for GSH adduct formation (eq 13). The 10 partly soluble compounds (gray) were not included in the calibration.

In addition, 10 additional compounds, which are assumed to be only partly soluble, are shown in Figure 7. In more detail, the solubility is below or on the same order of magnitude as the RC50 value, given in Table 1. Therefore, their predicted RC50 values are generally lower than those published. In turn, their intrinsic chemical reactivity will be greater than that indicated by the published RC50 values, because of their limited solubility. This shows one limitation of the static chemoassay, which is obviously highly dependent on the solubility of a compound. Kinetic assays are, by definition, not affected by this limitation. It should be noted that measured kinetic rate constants are given for those compounds, which are only partly soluble in the RC50 assay, and the calculated rate constants are predicted correctly by kinetic model B. This is one further advantage of the application of the in silico model presented. It allows for the estimation of chemical (and biochemical) potency of a compound, even if the experimental conditions are limited, for example, related to low solubility or high volatility.

Conclusions We have introduced a model to predict the reactivity, in terms of GSH adduct formation, of R,β-unsaturated compounds. A quantum chemical model for various compound classes assumed to react by Michael addition has been shown to work in a quantitative manner. The model has good statistical performance and is mechanistically interpretable. Furthermore, the model is based on easily computed ground-state properties, which are less computationally demanding than, for example, quantum chemical transition-state calculations. Starting from the initial model A that utilized measured rate constants, our inclusion of recalculated RC50 values from the Schultz assay confined to compounds in the intermediate reactivity range led simultaneously to good model performance and to an enlargement of the chemical model domain. This should be valid for a wide range of Michael acceptor compounds with a single reaction site. With regard to the best performance and inclusion of the full chemical domain, refined model B is recommended. This is a combination of model B1 for compounds with double bonds and model B2 for those with triple bonds. Despite the poorer performance of model C, in which compounds with double and triple bonds between the R-carbon and the β-carbon are calibrated together, it is still a reasonable predictive approach.

1584

Chem. Res. Toxicol., Vol. 23, No. 10, 2010

All subdomains for compounds capable of acting directly by Michael addition, as proposed by Schultz et al. (21), with one single reaction site are covered by the model. However, the following three subdomains are not yet included: first, quinones (with more than one reaction site), which are known to be very highly reactive; second, pre-Michael acceptors that require initial air oxidization of the compounds; and third, pro-Michael acceptors, for which metabolic activation to the actual Michael acceptor is requried. Because the model is based on local, sitespecific properties, further extensions to molecules with more than one reaction site should be possible. An indication of Michael addition reactivity for the latter two subdomains could, in principle, be obtained by calculations of the transformed reactive compounds. The Michael acceptor model allows for a theoretical prediction of RC50 values, even where experiments are limited by the solubility of particular chemical compounds. The Michael addition reactivity model can be used in the reactivity or toxicity profiling of electrophilic compounds. It can support mechanistic grouping to order compounds according to reactivity and also confirm that a compound is reactive. Further work will focus on the applications of the model in qualitative and quantitative toxicity prediction, for example, the acute toxicity and skin sensitization of potential Michael-type acceptors. Acknowledgment. This study was financially supported by the European Union through the InSilicoTox Marie Curie Project (Contract MTKD-CT-2006-42328) and OSIRIS Integrated Project (Contract GOCE-CT-2007-037017).

References (1) Aptula, A. O., and Roberts, D. W. (2006) Mechanistic applicability domains for nonanimal-based prediction of toxicological end points: General principles and application to reactive toxicity. Chem. Res. Toxicol. 19, 1097–1105. (2) Schultz, T. W., Carlson, R. E., Cronin, M. T. D., Hermens, J. L. M., Johnson, R., O’Brien, P. J., Roberts, D. W., Siraki, A., Wallace, K. D., and Veith, G. D. (2006) A conceptual framework for predicting toxicity of reactive chemicals: Models for soft electrophilicity. SAR QSAR EnViron. Res. 17, 413–428. (3) Roberts, D. W., Aptula, A. O., and Patlewicz, G. (2007) Electrophilic chemistry related to skin sensitization. Reaction mechanistic applicability domain classification for a published data set of 106 chemicals tested in the mouse local lymph node assay. Chem. Res. Toxicol. 20, 40–60. (4) Roberts, D. W., Patlewicz, G., Kern, P. S., Gerberick, F., Kimber, I., Dearman, R. J., Ryan, C. A., Basketter, D. A., and Aptula, A. O. (2007) Mechanistic applicability domain classification of a local lymph node assay dataset for skin sensitization. Chem. Res. Toxicol. 20, 1019– 1030. (5) Von der Ohe, P. C., Ku¨hne, R., Ebert, R.-U., Altenburger, R., Liess, M., and Schu¨u¨rmann, G. (2005) Structural alerts: A new classification model to discriminate excess toxicity from narcotic effect levels of organic compounds in the acute daphnid assay. Chem. Res. Toxicol. 18, 536–555. (6) Cronin, M. T. D., Dearden, J. C., Duffy, J. C., Edwards, R., Manga, N., Worth, A. P., and Worgan, A. D. P. (2002) The importance of hydrophobicity and electrophilicity descriptors in mechanistically-based QSARs for toxicological endpoints. SAR QSAR EnViron. Res. 13, 167– 176. (7) Schultz, T. W., Yarbrough, J. W., Hunter, R. S., and Aptula, A. O. (2007) Verification of the structural alerts for Michael acceptors. Chem. Res. Toxicol. 20, 1359–1363. (8) Koleva, Y. K., Madden, J. C., and Cronin, M. T. D. (2008) Formation of categories from structure-activity relationships to allow read-across for risk assessment: Toxicity of R,β-unsaturated carbonyl compounds. Chem. Res. Toxicol. 21, 2300–2313. (9) Tanii, H., and Hashimoto, K. (1982) Structure-toxicity relationship of acrylates and methacrylates. Toxicol. Lett. 11, 125–129. (10) Potter, D. W., and Tran, T.-B. (1992) Rates of ethyl acrylate binding to glutathione and protein. Toxicol. Lett. 62, 275–285.

Schwo¨bel et al. (11) Schultz, T. W., Yarbrough, J. W., and Johnson, E. L. (2005) Structureactivity relationships for reactivity of carbonyl-containing compounds with glutathione. SAR QSAR EnViron. Res. 16, 313–322. (12) Aptula, A. O., Patlewicz, G., Roberts, D. W., and Schultz, T. W. (2006) Non-enzymatic glutathione reactivity and in vitro toxicity: A nonanimal approach to skin sensitization. Toxicol. in Vitro 20, 239–247. (13) Freidig, A. P., Verhaar, H. J. M., and Hermens, J. L. M. (1999) Comparing the potency of chemicals with multiple modes of action in aquatic toxicology: Acute toxicity due to narcosis versus reactive toxicity of acrylic compounds. EnViron. Sci. Technol. 33, 3038–3043. (14) Cronin, M. T. D., Bajot, F., Enoch, S. J., Madden, J. C., Roberts, D. W., and Schwo¨bel, J. (2009) The in chemico-in silico interface: Challenges for integrating experimental and computational chemistry to identify toxicity. ATLA, Altern. Lab. Anim. 37, 513–521. (15) Schaafsma, G., Kroese, E. D., Tielemans, E. L. J. P., van de Sandt, J. J. M., and van Leeuwen, C. J. (2009) REACH, non-testing approaches and the urgent need for a change in mind set. Regul. Toxicol. Pharmacol. 53, 70–80. (16) Bo¨hme, A., Thaens, D., Paschke, A., and Schu¨u¨rmann, G. (2009) Kinetic glutathione chemoassay to quantify thiol reactivity of organic electrophiles: Application to R,β-unsaturated ketones, acrylates, and propiolates. Chem. Res. Toxicol. 22, 742–750. (17) Yarbrough, J. W., and Schultz, T. W. (2007) Abiotic sulfhydryl reactivity: A predictor of aquatic toxicity for carbonyl-containing R,βunsaturated compounds. Chem. Res. Toxicol. 20, 558–562. (18) Esterbauer, H., Zollner, H., and Scholz, N. (1975) Reaction of glutathione with conjugated carbonyls. Z. Naturforsch. 30c, 466–473. (19) McCarthy, T. J., and Witz, G. (1997) Structure-activity relationships in the hydrolysis of acrylate and methacrylate esters by carboxylesterase in vitro. Toxicology 116, 153–158. (20) Freidig, A. P., Verhaar, H. J. M., and Hermens, J. L. M. (1999) Quantitative structure-property relationships for the chemical reactivity of acrylates and methacrylates. EnViron. Toxicol. Chem. 18, 1133– 1139. (21) Schultz, T. W., Rogers, K., and Aptula, A. O. (2009) Read-across to rank skin sensitization potential: Subcategories for the Michael acceptor domain. Contact Dermatitis 60, 21–31. (22) Clarke, E. D., Greenhow, D. T., and Adams, D. (1998) Metabolismrelated assays and their application to agrochemical research: Reactivity of pesticides with glutathione and glutathione transferases. Pestic. Sci. 54, 385–393. (23) Blandamer, M. J., Burgess, J., Robertson, R. E., and Scott, J. M. W. (1982) Dependence of equilibrium and rate constants on temperature and pressure. Chem. ReV. 82, 259–286. (24) Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., Montgomery, J. A., Vreven, T., Kudin, K. N., Burant, J. C., Millam, J. M., Iyengar, S. S., Tomasi, J., Barone, V., Mennucci, B., Cossi, M., Scalmani, G., Rega, N., Petersson, G. A., Nakatsuji, H., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Klene, M., Li, X., Knox, J. E., Hratchian, H. P., Cross, J. B., Adamo, C., Jaramillo, J., Gomperts, R., Stratmann, R. E., Yazyev, O., Austin, A. J., Cammi, R., Pompelli, C., Ochterski, J. W., Ayala, P. Y., Morokuma, K., Voth, G. A., Salvador, P., Dannenberg, J. J., Zakrzewski, V. G., Dapprich, S., Daniels, A. D., Strain, M. C., Farkas, O., Malick, D. K., Rabuck, A. D., Raghavavachari, K., Foresman, J. B., Ortiz, J. V., Cui, Q., Baboul, A. G., Clifford, S., Cioslowski, J., Stefanov, B. B., Liu, G., Liashenko, A., Piskorz, P., Komaromi, I., Martin, R. L., Fox, D. J., Keith, T., Al-Laham, M. A., Peng, C. Y., Nanayakkara, A., Challacombe, M., Gill, P. M. W., Johnson, B., Chen, W., Wong, M. W., Gonzalez, C., and Pople, J. A. (2003) Gaussian 03, revision C.02, Gaussian, Inc., Pittsburgh, PA. (25) Reed, A. E., Weinstock, R. B., and Weinhold, F. (1985) Natural population analysis. J. Chem. Phys. 83, 735–746. (26) Foster, J. P., and Weinhold, F. (1980) Natural hybrid orbitals. J. Am. Chem. Soc. 102, 7211–7218. (27) Klamt, A. (1993) Estimation of gas-phase hydroxyl radical rate constants of organic compounds from molecular orbital calculations. Chemosphere 26, 1273–1289. (28) Schwo¨bel, J., Ebert, R.-U., Ku¨hne, R., and Schu¨u¨rmann, G. (2009) Modeling the H bond donor strength of -OH, -NH, and -CH sites by local molecular parameters. J. Comput. Chem. 30, 1454–1464. (29) Bo¨hnhardt, A., Ku¨hne, R., Ebert, R.-U., and Schu¨u¨rmann, G. (2008) Indirect photolysis of organic compounds: Prediction of OH reaction rate constants through molecular orbital calculations. J. Phys. Chem. A 112, 11391–11399. (30) Schwo¨bel, J., Ebert, R.-U., Ku¨hne, R., and Schu¨u¨rmann, G. (2009) Prediction of the intrinsic hydrogen bond acceptor strength of chemical substances from molecular structure. J. Phys. Chem. A 113, 10104– 10112. (31) Sanner, M. (1996) MSMS, Molecular Graphics Laboratory, La Jolla, CA.

Michael Addition ReactiVity (32) Connolly, M. L. (1992) The molecular surface package. J. Mol. Graphics 11, 139–141. (33) Parr, R. G., Szentpa´ly, L. v., and Liu, S. (1999) Electrophilicity index. J. Am. Chem. Soc. 121, 1922–1924. (34) Chattaraj, P. K., and Roy, D. R. (2007) Update 1 of Electrophilicity Index. Chem. ReV. 107, PR46–PR74. (35) Domingo, L. R., Aurell, M. J., Pe´rez, P., and Contreras, R. (2002) Quantitative characterization of the global electrophilicity power of common diene/dienophile pairs in Diels-Alder reactions. Tetrahedron 58, 4417–4423. (36) Domingo, L. R., Pe´rez, P., and Contreras, R. (2003) Electronic contributions to the σp parameter of the Hammett equation. J. Org. Chem. 68, 6060–6062. (37) Enoch, S. J., Cronin, M. T. D., Schultz, T. W., and Madden, J. C. (2008) Quantitative and mechanistic read across for predicting the skin sensitization potential of alkenes acting via Michael addition. Chem. Res. Toxicol. 21, 513–520. (38) Enoch, S. J., Roberts, D. W., and Cronin, M. T. D. (2009) Electrophilic reaction chemistry of low molecular weight respiratory sensitizers. Chem. Res. Toxicol. 22, 1447–1453. (39) Campodo´nico, P. R., and Contreras, R. (2008) Structure-reactivity relationships for electrophilic sugars in interaction with nucleophilic biological targets. Bioorg. Med. Chem. 16, 3184–3190. (40) Campodo´nico, P. R., Fuentealba, P., Castro, E. A., Santos, J. G., and Contreras, R. (2005) Relationships between the electrophilicity index and experimental rate coefficients for the aminolysis of thiolcarbonates and dithiocarbonates. J. Org. Chem. 70, 1754–1760. (41) Wondrousch, D., Bo¨hme, A., Thaens, D., Ost, N., and Schu¨u¨rmann, G. (2010) Local electrophilicity predicts the toxicity-relevant reactivity of Michael acceptors. J. Phys. Chem. Lett. 1, 1605–1610. (42) Bo¨hnhardt, A., Ku¨hne, R., Ebert, R.-U., and Schu¨u¨rmann, G. (2010) Predicting rate constants of OH radical reactions with organic substances: advances for oxygenated organics through a molecular

Chem. Res. Toxicol., Vol. 23, No. 10, 2010 1585

(43)

(44) (45)

(46) (47) (48)

(49) (50) (51)

orbital HF/6-31G** approach. Theor. Chem. Acc. DOI: 10.1007/ s00214-009-0724-8. Schwo¨bel, J., Ebert, R.-U., Ku¨hne, R., and Schu¨u¨rmann, G. (2009) Prediction of the intrinsic hydrogen bond acceptor strength of organic compounds by local molecular parameters. J. Chem. Inf. Model. 49, 956–962. Schwo¨bel, J., Madden, J. C., and Cronin, M. T. D. (2010) Examination of Michael addition reactivity towards glutathione by transition state calculations. SAR QSAR EnViron. Res. (in press). Schu¨u¨rmann, G., Ebert, R.-U., Chen, J., Wang, B., and Ku¨hne, R. (2008) External validation and prediction employing the predictive squared correlation coefficient - Test set activity mean vs training set activity mean. J. Chem. Inf. Model. 48, 2140–2145. Aptula, N., Roberts, D. W., Schultz, T. W., and Pease, C. (2007) Reactivity assays for non-animal based prediction of skin sensitisation potential. Toxicology 231, 104–119. Chan, K., and O’Brien, P. J. (2008) Structure-activity relationships for hepatocyte toxicity and electrophilic reactivity of R,β-unsaturated esters, acrylates and methacrylates. J. Appl. Toxicol. 28, 1004–1015. Chan, K., Poon, R., and Brien, P. J. (2008) Application of structureactivity relationships to investigate the molecular mechanisms of hepatocyte toxicity and electrophilic reactivity of R,β-unsaturated aldehydes. J. Appl. Toxicol. 28, 1027–1039. Harder, A., Escher, B. I., and Schwarzenbach, R. P. (2003) Applicability and limitation of QSARs for the toxicity of electrophilic chemicals. EnViron. Sci. Technol. 37, 4955–4961. ChemIDplus NLM (2009) U.S. National Library of Medicine, Bethesda, MD. EPISUITE WSKOW, version 1.40 (2000) U.S. Environmental Protection Agency, Washington, DC.

TX100172X