A Computational Study of Detoxification of Lewisite Warfare Agents by

Article pubs.acs.org/JPCA

A Computational Study of Detoxification of Lewisite Warfare Agents by British Anti-lewisite: Catalytic Effects of Water and Ammonia on Reaction Mechanism and Kinetics Chandan Sahu, Srimanta Pakhira, Kaushik Sen, and Abhijit K. Das* Department of Spectroscopy, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700032, India S Supporting Information *

ABSTRACT: trans-2-Chlorovinyldichloroarsine (lewisite, L agent, Lew-I) acts as a blistering agents. British anti-lewisite (BAL, 2,3-dimercaptopropanol) has long been used as an L-agent antidote. The main reaction channels for the detoxification proceed via breaking of As−Cl bonds and formation of As−S bonds, producing stable, nontoxic ring product [(2-methyl-1,3,2-dithiarsolan-4-yl)methanol]. M06-2X/GENECP calculations have been carried out to establish the enhanced rate of detoxification mechanism in the presence of NH3 and H2O catalysts in both gas and solvent phases, which has been modeled by use of the polarized continuum model (PCM). In addition, natural bond orbital (NBO) and atoms in molecules (AIM) analysis have been performed to characterize the intermolecular hydrogen bonding in the transition states. Transitionstate theory (TST) calculation establishes that the rates of NH3-catalyzed (2.88 × 10−11 s−1) and H2O-catalyzed (2.42 × 10−11 s−1) reactions are reasonably faster than the uncatalyzed detoxification (5.44 × 10−13 s−1). The results obtained by these techniques give new insight into the mechanism of the detoxification process, identification and thermodynamic characterization of the relevant stationary species, the proposal of alternative paths on modeled potential energy surfaces for uncatalyzed reaction, and the rationalization of the mechanistic role played by catalysts and solvents.

■

INTRODUCTION Contamination of environment by chemical warfare agents (CWA) can cause serious health hazards. Among CWAs, trans2-chlorovinyldichloroarsine, also known as lewisite (Lew-I), is a powerful blistering chemical warfare agent, first produced near the end of World War I, capable of inflicting severe chemical burns of the eyes, skin, and lungs.1 Lewisite is prepared by the reaction of acetylene with arsenic trichloride.2 Trivalent arsenic is toxic because of its reactivity with important biological sulfhydryls.3 Peters et al.4 analyzed lewisite compounds of a product made from the protein in hair (or skin), known as keratin, and found that the arsenical compound combined with the sulfhydryl groups in the hair keratin. Moreover, one arsenic atom combined with two sulfur groups, forming a chemical ring. The idea arose that if a compound containing two SH groups could be made capable of forming a more tightly fixed compound with the arsenic in the lewisite than that in the tissue, the poisonous arsenic could be removed. The only known antidotes are the dithiol compounds that compete and extract the arsenical from tissue sulfhydryl. The efficiency of dithiols depends on their ability to form cyclic adducts with arsenic, which are more stable than the adducts formed with monothiols.5 A requisite feature of antidotes for heavy metals such as arsenic is the ability of antidotes to sequester the metal and eventually excrete the adduct. Thus, the pharmacokinetics of the antidote and the stability of adduct formed are most important in the first step of the detoxification process. Webb et al.3 and Waters and Stock5 synthesized the dithiol-containing © 2013 American Chemical Society

compound British anti-lewisite (BAL, 2,3-dimercaptopropanol) in the search for an relevant antidote to Lew-I.4,5 BAL is a mobile oil, which passes readily through skin, where two OH groups of glycerol are replaced by SH. The BAL−arsenic compound is removed from the skin and excreted in the urine. Although a number of experimental investigations of detoxification procedures for lewisite-type compounds were carried out carefully, no theoretical study is devoted to establish the catalytic efficiency in detoxification mechanisms of Lew-I warfare agent. The objective of the present article is to provide mechanistic details of the detoxification of lewisite in presence of catalyst, and in a forthcoming paper we will focus on the modification of BAL to decrease its toxicity without compromising its efficiency. The nucleophilic reaction of BAL may proceed via two different routes having two different transition states (Scheme 1): (1) nucleophile attack by a 3mercapto group (SA) followed by a 2-mercapto group (SB), and (2) nucleophile attact by a 2-mercapto group (SB) followed by a 3-mercapto group (SA). The addition process to arsenic proceeds via cyclic transition states through -H transfer from S to Cl, and elimination preferentially takes place as HCl. Therefore, determination of the preferred route (sketched in Scheme 1) is quite challenging. Received: December 13, 2012 Revised: March 28, 2013 Published: March 29, 2013 3496

dx.doi.org/10.1021/jp312254z | J. Phys. Chem. A 2013, 117, 3496−3506

The Journal of Physical Chemistry A

Article

Scheme 1. (A) Possible Competitive Pathways for Uncatalyzed Detoxification of Lewisite and (B) Possible Mechanisms for Catalytic Detoxification Pathways of Lewisite

■

COMPUTATIONAL DETAILS The biological activity of molecules is known to be crucially dependent on electronic structures of the active part of the given compound and its conformation. Therefore, we first performed a conformational analysis of BAL using molecular dynamics (MD) and molecular mechanics (MM) approaches. For computational investigation of conformational minima of BAL, it is subjected to a MD conformational search with an unconstrained MD trajectory via the Verlet velocity algorithm and NVE thermostat with other default parameters in Gabedit v2.3.8.6 The minimum conformational geometries are obtained via the PM6 semiempirical method as implemented in MOPAC 2009.7 Additional molecular mechanics conformational studies are carried out with the MMFF94 force field and a systematic

In the present work we have performed an exhaustive theoretical study of the title reaction by use of density functional theory (DFT) with the main aim of accurately defining the barrier height and stability of the complexes. We first locate and characterize the stationary points (reactants, products, saddle point, and complexes in the entry and exit channels) by using high-level electronic structure calculations. We then calculate the reaction path, which describes a chemical reaction, using energies only in the region of configuration space along the reaction path. Finally, we have performed transition-state theory (TST) calculations to obtain kinetics information. 3497



Article

rotor search with Avogadro v1.0.3.8 For final validation of the conformational analysis, 25 representative minimum structures are selected from the rigorous conformational search and are treated with density functional theory (DFT). The geometries of all the molecular species involved in this study are fully optimized by employing DFT with the hybrid functional (M062X) of Truhlar and Zhao.9 The all-electron Ahlrichs’ TZVP10 basis set is used for C, O, S, H, N, and Cl atoms and the triple-ζ cc-pVTZ-PP basis set11 with small-core relativistic effective core potential (RECP), which is found to have performed better for arsenic-containing systems,12 is used for As. RECP is a useful means of replacing the core electrons in a calculation with an effective potential, thereby eliminating the need for the core basis function, which usually require a large set of Gaussians to describe them. In addition to replacing the core, they are used to represent relativistic effects, which are largely confined to the core. In this context, both the scalar (spin-free) relativistic effects and spin−orbit (spin-dependent) relativistic effects are included in the effective potentials.13 For convenience, the theoretical level is denoted as M06-2X/GENECP. Harmonic vibrational frequencies are determined at the M06-2X/ GENECP level to confirm whether the optimized structures are local minima (no imaginary frequency) or transition states (one imaginary frequency) on the potential energy surfaces (PES) and to evaluate the zero-point vibration energy (ZPVE) and thermal corrections to the Gibbs free energy at T = 298 K. The connecting first-order saddle points, which are the transition states between the equilibrium geometries, are obtained by the synchronous transit-guided quasi-Newton (STQN) method. A parallel intrinsic reaction coordinate (IRC) calculation is performed with all transition states to confirm whether these transition states connect the right minima or not.14,15 The gas-phase optimization and frequency calculations are performed initially with each computational method. The integral equation formalism polarizable continuum model (IEFPCM)16 and simulation-based minimum distance (SMD) solvent model,17 as implemented in Gaussian 09,18 have been used to take into account the solvent effects. The Gibbs free energies in solvent (Gsolvent) are calculated by adding the total energy in the solvent (Esolvent) and the gasphase thermal correction to the Gibbs free energy (TCG gas) at 298.15 K. Three solvents, cyclohexane (ε = 2.02), dimethyl sulfoxide (DMSO) (ε = 46.70), and water (ε = 78. 39),19 available in the PCM model, are used to study the effects of solvent polarity on the reaction mechanism. The difference between nonprotic and protic solvents is due to the difference in the ALPHA scaling parameter by which SMD sphere radii, 1.4 for nonprotic solvents and 1.2 for protic ones,19 are multiplied in the PCM model. Potential energy surfaces for all routes have been constructed. An extensive investigation is conducted to detect the most favorable pathways as well as to quantify the reaction progress via the most common modern strategies: activation energy, natural bond orbital (NBO) analysis,20 Wiberg bond order calculation, and Bader’s atoms in molecules (AIM) approach.21−23 The problem of temperature dependence, crucial to rate theories, was first studied quantitatively by Wilhelmy,24 and after that a number of other equations were proposed but the matter was by no means settled.25,26 However, the Arrhenius equation

k = Ae−Ea / RT

has been accepted for calculation of temperature-dependent reaction rate constant. The Arrhenius equation gives the dependence of the rate constant k of a chemical reaction on the absolute temperature T, where A is the pre-exponential factor (or simply the prefactor), and Ea is the activation energy. The natural logarithm of the Arrhenius equation yields

ln(k) = ln(A) −

Ea RT

(2)

The derivative of eq 2 with respect to T gives

E d ln k = a2 dT RT

(3)

Or, integrating eq 3 we get ln

E ⎡1 k2 1⎤ = a⎢ − ⎥ k1 R ⎣ T1 T2 ⎦

(4)

where k1 and k2 are the rate constants at temperatures T1 and T2, respectively. In this article, the kinetic study of uncatalyzed and catalyzed detoxifications is conducted under collisionless conditions by use of transition-state theory (TST)27−32 where another Arrhenius-like expression appears for the rate constant. In conventional TST, the forward thermal rate constant can be calculated with information about only two configurations, the reactants and transition state, leaving the symmetry numbers out of the rotational partition function.33 The expression for TST rate constant is KTST(T ) =

⎛ ΔG⧧ ⎞ TkB ⎟ exp⎜ − h ⎝ RT ⎠

(5) ⧧

where kB is Boltzmann’s constant, h is Planck’s constant, ΔG is the Gibbs free energy of activation, and T is the temperature. Comparing eq 1 with eq 5, overall rate constants and Arrhenius parameters are calculated at 300 K. This is also supported by an Arrhenius plot (ln k vs T−1) obtained by numerical integration of eq 2 where slope and intercept can be used to determine Ea and A in the temperature range 150−1000 K. Rate constants are minimized as a function of position along the minimum energy path (MEP) to get the rate constant at each temperature. All transition state structures on the MEP contain a single imaginary frequency with the mode vibration corresponding to the motion along the bond breaking and forming coordinate. An AIM analysis is performed to evaluate the strength of the intermolecular hydrogen bond in the transition-state structure. For determination of the strength of the H-bond, the important topological parameters electron density (ρBCP) and its Laplacian (∇2ρBCP) at the bond critical point (BCP) are calculated with the AIM 2000 program.34 An NBO analysis is carried out for the evaluation of the bond order in transitionstate geometry. The NBO population analysis and Wiberg bond index calculation are performed with the NBO 3.1 program35 as implemented in Gaussian 09. All electronic structure calculations are performed using the Gaussian 09 suite of quantum chemistry programs.18

■

RESULTS AND DISCUSSION Conformational Search of Reactants and Transition States. Conformational search is a process of finding lowenergy conformations of molecular systems by varying

(1) 3498



Article

optimized geometries are presented in Figure 1, and the relative energies are given in parentheses below the structures. The conformer of lowest energy found here is C1. Another reactant, chlorovinyldichloroarsine, also has three isomeric forms: germinal, cis, and trans. Urban and von Tersch36 reported that the trans isomer is most favored compared to cis isomer (∼0.6 kcal·mol−1) and geminal isomer (∼1 kcal·mol−1) at the MP2/6-311+G(2d,p)//MP2/6-31+G(2d,p) level. For this reason, all calculations are carried out on the trans isomer of lewisite. In the case of detoxification of Lew-I, there are two possibilities: nucleophilic attack by 3-mercapto group (SA) and by 2-mercapto group (SB) to the As metal center. In an effort to understand the patterns of bond cleavage and to make comparisons between the two pathways, it is first necessary to understand the mechanism and transition-state structures of entire reaction channels. We have also performed conformational searches for transition states for the catalytic detoxification pathways (see Figures S3 and S4 in the Supporting Information for conformations and corresponding energies). To determine the most favorable one, both routes are also investigated here in detail. Moreover, the nucleophilic reactions have been studied in the presence of one additional H2O and NH3 as catalyst on favorable pathway. We first analyze the uncatalyzed detoxification pathways. Uncatalyzed Detoxification. In the present theoretical study, all reactions follow the SN2 pattern and proceed through direct displacement mechanism involving only one step, in which the expulsion of a leaving group occurs at the same time when the substituting nucleophiles enter; that is, in a concerted and asynchronous way, which is in good agreement with the

Figure 1. Gas-phase optimized geometries at M06-2X/TZVP level (relative energies are given in parentheses) for the most stable conformers of BAL with important geometrical parameters.

geometric parameters. The minima for 25 selected conformers of BAL are finally optimized at the M06-2X/TZVP level of theory. After DFT optimization, we get four lower-energy structures having large dissimilar spatial arrangements. The

Figure 2. M06-2X/GENECP potential energy profile (ZPE-corrected) for the uncatalyzed detoxification of lewisite by BAL in gas phase. Free energy values (ΔG298K) are given in parentheses. 3499



Article

Figure 3. M06-2X/GENECP potential energy profile (ZPE-corrected) for H2O-catalyzed detoxification of lewisite by BAL in gas phase. Free energy values (ΔG298K) are given in parentheses.

comparable single-step nucleophilic reaction mechanism. Scheme 1 shows a stepwise mechanism through neutral intermediate along with a concerted transition state. So, we now explore the pathways in the following discussion. For the concerted step, the reaction consists of one step in which all bond-forming and bond-breaking processes occur in concert: reactant (R) → transition state (TS) → product (P). Here, a proton from the −SH group is shifted to the Cl atom attached to the metal center; simultaneously there is a nucleophilic attack by lone-pair electrons of S in the SH group on the electrophilic As metal center, and breaking of the As−Cl and formation of the As−S bond occur. Optimized geometries of the related species for both pathways are presented in the potential energy surface (PES). In this report all analyses are carried out with respect to the M06-2X/ GENECP energy. From the reaction scheme shown in Scheme 1, it is observed that two different pathways may be possible for detoxification of Lew-I. We first analyze nucleophilic attack by the first pathway. First Pathway: Nucleophilic SA Group Substitution Followed by SB Group. The changes of energy (ΔE) and free energy (ΔG) relative to the isolated reactants for the stationary points optimized along the reaction pathway are included in the potential energy profile (Figure 2). Coordinates of the gas-phase optimized geometries are given in the Supporting Information. For the neutral mechanism, two

similar types of concerted neutral transition state (TS) are calculated. The first TS (UTS1-SA), with 24.18 kcal·mol−1 activation energy with respect to the isolated reactants, has a four-member ring structure in which As−S (SA) bond-forming and As−Cl bond-breaking occur, along with the transfer of Hatom from S (SA) atom to Cl atom. UTS1-SA contains an imaginary frequency at 889.37 cm−1, and the transition vector for this mode corresponds mainly to the shuttle of proton from the nucleophile −SH group to the Cl atom. UTS1-SA finally yields the intermediate complex UINT-SA, which is located 20.12 kcal·mol−1 above the isolated reactants and 4.06 kcal·mol−1 below UTS1-SA in the PES. The last step of the pathway is an elimination reaction that occurs through the transition state UTS2-SB. UINT-SA transforms into the transition state UTS2-SB with 24.27 kcal·mol−1 activation barrier. It has an imaginary frequency at 574.81 cm−1 owing to the H-transfer transition vector from S (SB) to Cl atom, and at the same time the As−S (SB) bond is formed. Finally, the detoxified five-member ring product, PC, is formed with elimination of HCl. Second Pathway: Substitution of Nucleophilic SB Group Followed by SA Group. In this pathway, we focus mainly on elucidation of the mechanism of As−Cl bond cleavage and As−S bond formation during nucleophilic attack by the nucleophile SB group, followed by the SA group. The optimized geometry along with geometrical parameters and the 3500



Article

Figure 4. M06-2X/GENECP potential energy profile (ZPE-corrected) for NH3-catalyzed detoxification of lewisite by BAL in gas phase. Free energy values (ΔG298K) are given in parentheses.

frequency 468.71 cm−1 for UTS1-SB and 426.07 cm−1 for UTS2-SA, respectively, are due to the H shift from nucleophilic −SH to Cl atom. Several attempts have been made to check the formation of product by the attack of two nucleophiles on a single transition state, but all have been unsuccessful. From the above investigation, it is clear that the two pathways are roughly comparable. To get more information about different pathways, we have performed NBO and AIM analyses.

ZPVE-corrected gas-phase potential energy profile are presented in Figure 2. The IRC calculation shows that the reaction proceeds by forming prereactant (URC) and preproduct (UINT-SB) complexes with 7.56 kcal·mol−1 lower and 20.62 kcal·mol−1 higher energy than the separated reactants, respectively. In between URC and PC, there exist two transition states having four-member cyclic rings, with 25.84 kcal·mol−1 (UTS1-SB) and 22.63 kcal·mol−1 (UTS2-SA) energy of activation. The transition vectors for the imaginary 3501



Article

Figure 5. Structures of all transition states (TS) of uncatalyzed and catalyzed detoxification pathways of lewisite at M06-2X/GENECP level of theory.

Comparison of Different Pathways. Following the above discussion, the relative reactivity of two nucleophile groups needs special attention. We have evaluated some important conceptual parameters relative to nucleophilic reactivity to establish favorable pathways of detoxification. The most important factor that favors detoxification is the strength of intermolecular H-bonding in the corresponding transition states. To quantify the strength of intermolecular H-bonding, NBO analysis is performed on the transition states. The

occupation number (ON) and significant donor−acceptor interactions measured by second-order perturbation energies are collected in Table S4 in the Supporting Information. For the first pathway, second-order perturbation energies for interaction of the lone pair of donor chlorine (Cl 14) with acceptor σ* (H15−S7) for UTS1-SA and for donor Cl 15 with acceptor σ* (H14−S8) for UTS2-SB are 178.31 and 24.69 kcal·mol−1, respectively. Similarly, for the second pathway, energies for interaction of the lone pair of donor chlorine (Cl 3502



Article

in UTS2-SA and UTS2-SB is of the same nature. In both cases, the acceptor NBO (σ*S−H) is composed of about 66.3% hydrogen NHO, which is further made up of 1s orbital. The donor chlorine NBO is composed mainly of 3p orbital. These intermolecular H-bonds (Figure S1, Supporting Information) between chlorine (Cl 14 in UTS2-SA and Cl 15 in UTS2-SB) and hydrogen (H 15 in UTS2-SA and H 14 in UTS2-SB) are associated with 27.57 and 24.69 kcal·mol−1 second-order perturbation energies, respectively. The Wiberg bond order for the intermolecular H-bond in the respective TS shows that the H-bond in UTS1-SA is almost same as that in UTS1-SB (Table S4, Supporting Information), and also UTS2-SB is similar to UTS2-SA. As a second strategy, the theory of AIM has been applied to define and describe the H-bonding phenomena in the transition states formed by the attack of two nucleophiles on Lew-I in the gas phase. The necessary indicator for the existence of Cl···H interaction is the presence of a bond (3, −1)-type critical point (CP) (in which curvature of the density at the CP along the line connecting two nuclei is positive and the two curvatures along mutually perpendicular lines and perpendicular to the bond path are negative), on which other topological properties such as the electron density (ρBCP) gradient and Laplacian of the electron density (∇2ρBCP) are calculated for characterizing the nature and strength of the Cl···H interaction.37 The calculated AIM parameters are presented in Table S5 in the Supporting Information. According to Popelier and Bader,38 there exists a standard condition for H-bonding in AIM formalism. The extent of electron density, ρBCP, and its Laplacian, ∇2ρBCP, for a normal H-bond at the CP should occur at 0.002−0.035 and 0.024−0.139 au,39 respectively. From the values in Table S5 (Supporting Information), it can clearly be found that the ρBCP and ∇2ρBCP values are within the aforesaid range, and consequently the existence of H-bonding interaction is strongly confirmed. Topological analysis further suggests that the intermolecular H-bonding interaction in UTS1-SA is similar to that in UTS1-SB, like UTS2-SB and UTS2-SA, supporting the NBO outcomes as well. Finally, from the above analysis the reactivity trends for nucleophiles for the detoxification reactions with lewisite have undoubtedly been established. Therefore, the order of reactivity predicts that both pathways are similar. Owing to this, we try to find out the influence of solvent and characterize the reaction in the presence of one additional H2O or NH3 molecule explicitly as catalyst in the first reaction channel. We first explore solvent effects on the detoxification reaction via UTS1-SA and UTS2SB transition states, as they are involved in the first reaction path. Solvent Effects on Reaction Pathway. To take into account the influence of solvent (DMSO, C6H12, and H2O) on the reaction mechanism, an SMD model has been employed. The standard Gibbs free energies relative to the isolated reactants for the stationary points optimized along the detoxification pathway at the M06-2X/GENECP-PCM level of theory are reported in Tables S1−S3 in the Supporting Information. Three separate sets of data are obtained because of different SMD radii for protic polar, nonprotic polar, and nonpolar solvents. The study of solvent effects on the concerted nonionic detoxification reaction mechanism indicates that, regardless of which functional is used or which approach is followed to account for solvent effects, the Gibbs free energy barriers in all the processes considered increase with solvent polarity. Our results show that the most favorable reaction

Figure 6. Arrhenius plots of calculated rate constants: (a) combined plots for NH3-catalyzed and uncatalyzed and (b) combined plots for H2O-catalyzed and uncatalyzed detoxification of lewisite by BAL. Combination of two plots (a and b) is given in Supporting Information.

15) with acceptor σ* (H14−S8) for UTS1-SB and for donor Cl 14 with acceptor σ* (H15−S7) for UTS2-SA are 180.21 and 27.57 kcal·mol−1, respectively. These results suggest that the interaction is similar in nature for the two pathways. The comparison of interaction is shown in the contour plot of orbital overlap depicted in Figure S1 in the Supporting Information. In the contour plot, the acceptor NBO (σ*S−H) is composed mainly of hydrogen’s natural hybrid orbital (NHO) (74.5% in UTS1-SA and 73.7% in UTS1-SB), and these hydrogen NHOs are further composed mostly of 1s orbital. Analysis of donor NBOs (LpCl) reveals that it is composed mainly of 3p orbital of Cl. The above intermolecular σ-type H-bonding between chlorine (Cl 14 in UTS1-SA and Cl 15 in UTS1-SB) and hydrogen (H 15 in UTS1-SA and H 14 in UTS1-SB) is shown in the plot. The intermolecular H-bonding in UTS1-SA and UTS1-SB is manifested in NBO calculation as LpCl → σ*S−H transitions. Also, the intermolecular H-bonding 3503



Article

Table 1. Calculation of Overall Reaction Rate Constants and Arrhenius Parametersa reaction

rate constant at 300 K (s−1)

Ea (kJ/ mol)

5.44 × 10−13

146.70

koverallb at 300 K (s−1)

Eoverall (kcal/mol)

A (cm3/ molecule·s)

fitted expression for rate constantc

35.04

1.9072 × 1013

1.9072 × 1013 exp(−17652/T)

32.67

1.9034 × 1013

1.9034 × 1013 exp(−16458/T)

1.9018 × 1013

1.9018 × 1013 exp(−16508/T)

BAL + Lew-I k1

URC → UTS1‐SA k2

7.28 × 102

59.83

k3

4.32 × 10

20.03

UINT‐SA → UTS1‐SA UINT‐SA → UTS2‐SB overall

13

5.44 × 10−13 BAL + Lew-I + NH3

2.87 × 10−11

136.80

k 2A

8.47 × 10−1

76.68

k3A

2.48 × 10

39.53

k1A

CRC‐A ⎯→ ⎯ CTS1‐A

CINT‐A ⎯⎯⎯→ CTS1‐A

CINT‐A ⎯→ ⎯ CTS2‐A overall

k1W

CRC‐W ⎯⎯⎯→ CTS1‐W

6

2.88 × 10−11

2.42 × 10−11

BAL + Lew-I + H2O 137.23

k 2W

1.38 × 101

69.71

k3W

2.46 × 108

28.08

CINT‐W ⎯⎯⎯→ CTS1‐W

CINT‐W ⎯⎯⎯→ CTS2‐W overall

2.42 × 10−11

a

32.77 b

Steady-state conditions for catalyzed and uncatalyzed detoxification of lewisite are assumed. koverall = k1k3/(k3 + k2). cFitted expression for rate constant at 150 K ≤ T ≤ 1000 K.

Transition state CTS1-W corresponds to the cleavage of As−Cl and H−S (SA) bonds and the formation of As−S and Cl−H bonds. The transitions state then converted to the stable intermediate CINT-W. In CINT-W, the As−S (SA) and Cl−H bonds are completely formed, and the catalytic H2O molecule is situated in a typical orientation stabilized through H-bond formation in the Cl−H−O−H extended moiety. In the final step, the intermediate CINT-W forms the product complex PC through the transition state CTS2-W, associated with an activation barrier of 6.20 kcal·mol−1 with respect to CINT-W. CTS2-W basically possesses the cleavage of As−Cl single bond and formation of As−S (SB) bond with simultaneous proton transfer from S−H to −Cl part via the catalytic H2O molecule. It is important to note from Figure 3 that the presence of two H2O molecules lowers the energy barrier. So, it can be concluded that inclusion of more explicit solvent molecules lowers the barrier further. From AIM analysis (numerical values given in Table S5 in the Supporting Information), it can be found that in both transition states there is H bonding interaction, which stabilizes the TS and enhances the reaction process. Also, the effect of solvent polarity on the catalytic step has been reported in Table S2 in the Supporting Information. It is concluded from Table S2 that the ΔG⧧ increases with increase of solvent polarity. NH3-Catalyzed Detoxification. The next step in this study is to investigate the catalytic effect of an additional NH3 molecule in the mechanism considered for the detoxification of lewisite. The reaction is sketched in Scheme 1. In comparison with the H2O-assisted process, the NH3 molecule can be designated as a catalyst. All the optimized geometries and activation parameters along with the reaction energy profile are shown in Figure 4. A detailed analysis of IRC for the concerted transition state along both sides shows the existence of prereactive and preproduct complexes; their structures and intermolecular hydrogen bonds are shown in

mechanism in nonpolar solvents, which is nonionic and fully concerted, becomes more energetic when the solvent becomes more polar. After implicit treatment of solvent, we find that the solvent has little effect on reactions, even though sometimes it proceeds through a higher barrier. This is due to the fact that the separated reactants are more stabilized and possess greater solvation energy. Although this seems to be doubtful, it can be explained by considering the fact that in this reaction the observed transition states are neutral and the charge separation is negligible. Owing to this, we try to find out the influence of solvent on reaction in presence of one additional H2O and NH3 molecules as catalyst. We first explore the H2O-assisted catalytic detoxification reaction. H2O-Catalyzed Detoxification. Since water is a proton donor and acceptor, it is possible that a water molecule can act as a catalyst in the concerted transition states (UTS1-SA and UTS2-SB) proposed previously. This concerted pathway is also modeled in different solvents having different polarities. The activation energies of H2O-catalyzed detoxification for As−Cl bond-breaking pathways are represented in Scheme 1. The optimized geometry along with the geometrical parameters and the ZPVE-corrected potential energy profile for the reactions with and without catalyst are shown in Figure 3. According to the findings, the first stationary point observed is the prereactive complex, CRC-W, which is much more stable (6.01 kcal·mol−1) compared to the prereactive complex found in the uncatalyzed pathway. In CRC-W, the important As−S distance is 3.29 Å. CRC-W then transforms into the first transition state, CTS1-W, with 16.12 kcal·mol−1 barrier of activation. The activation energy is significantly lower (8.06 kcal·mol−1) in the catalyzed detoxification than in the uncatalyzed one. In CTS1-W, the formation of As−S bond is quite advanced and it possesses a six-member cyclic structure. The transition vector for the imaginary frequency (655.01 cm−1) mode clearly reveals that the proton transfer occurs. 3504



Article

Figure 4. The transition vector indicates that, in CTS1-A, proton transfer takes place through a six-member cyclic path associated with 16.97 kcal·mol−1 transfer barrier, which is 7.21 kcal·mol−1 lower than the uncatalyzed process. In CTS1-A, the catalytic NH3 molecule forms a typical Cl···H−N bond with the Lew-I chlorine atom and an N···H−S hydrogen bond with the catalyst NH3 molecule. CTS1-A then transforms into a stable intermediate, CINT-A, with the formation of As−S bond. At the final step, the process of six-member ring formation takes place via the transition state CTS2-A, with formation of the As−S bond and removal of HCl. The CTS2-A has a barrier of activation of 16.23 kcal·mol−1. As expected and also observed from the PES, these transition states are similar to those of H2O-catalyzed reactions, and inclusion of additional NH3 molecules lowers the energy barrier further. According to the results obtained for H2O- and NH3-catalyzed reactions, it is obvious that the detoxification mechanism as well as the ratedetermining step along the pathway involves a proton shift through six-member cyclic transition states. These structures are characterized by less ring strain than four-member rings in the case of uncatalyzed reaction, thus lowering the proton transfer barrier and facilitating the processes. From AIM analysis (Table S5 in the Supporting Information), it can be found that in both transition states there is H-bonding interaction, which stabilizes the TS and enhances the reaction process. Also, the effect of solvent polarity on the catalytic step has been reported in Table S3 in the Supporting Information. It is concluded that ΔG⧧ increases with increasing solvent polarity. All the transition-state structures for catalyzed and uncatalyzed reactions are shown in Figure 5 to compare geometrical changes during detoxification. TST Calculation of Rate Constants. In the DFT calculation we have shown that the key channels responsible for uncatalyzed and catalyzed detoxification of Lew-I are those where breaking of the As−Cl bond and formation of the As−S bond, along with elimination of HCl via hydrogen-bonded complexes, take place. These are the reactions that pass through many transition states. The rate constant (k) corresponding to all the studied reaction channels can be analyzed by TST:40,41

The Arrhenius plots (data of Arrhenius plot obtained by numerical integration is given in the Supporting Information) of overall rate constants for three different reactions are shown in Figure 6, and the rate constant expressions that resulted from plot are presented in Table 1. From Figure 6, Arrhenius behavior is observed and the formation of hydrogen bond in the TS possibly takes place at the origin of the graph, showing the normal temperature dependence of reaction rates. The lowest barrier is calculated for the catalyzed reaction, and the rate constants of NH3-catalyzed (kA) and H2O-catalyzed (kW) are higher than those of the uncatalyzed channels. For the 150−1000 K temperature range, an overall A factor of 1.9034 × 1013 cm3/(molecule·s) and apparent activation energy of 32.67 kcal·mol−1 are calculated for the NH3-catalyzed pathway. The similar reaction for H2O-assisted detoxification has an overall activation barrier of 32.77 kcal·mol−1 and A factor of 1.9018 × 1013 cm3/(molecule·s). The slowest reaction among the three channels is the uncatalyzed detoxification, having an A factor of 1.9072 × 1013 cm3/(molecule·s) for the rate constant (kuc) at 300 K and an apparent activation energy of 35.04 kcal/mol. It is evident from the results presented in Table 1 that NH3 and H2O significantly catalyze the detoxification reaction at 300 K.

■

The model detoxification reaction of Lew-I has been investigated in gas and solvent phases by the DFT (M06-2X) method. All the reaction mechanisms are involved in proton transfer. Two competing detoxification pathways are explored. From the outcomes, we reach the conclusion that both pathways are quite favorable. In this study we have treated the pathways separately, and from the analysis of barrier height we see that the probability of breaking for both bonds is almost the same. The reaction in the presence of catalyst passes through a lower energy barrier, which indicates the purely uncatalyzed detoxification of Lew-I is not competitive with the H2O- or NH3-catalyzed mechanisms. The barrier is further lower with inclusion of more catalyst molecule explicitly. The present computational investigation provides quantitative results for detoxification of an arsenic-containing warfare compound for the first time. The mechanism proposed in this work presents some novelty arising from the catalytic role of NH3 and H2O in the nucleophilic attack of 3- and 2-mercapto groups. The results obtained by PCM-SMD model allow us to propose a mechanism consistent with general features of a SN2 reaction. A topological analysis of electron density, carried out via the AIM formalism, reinforced some of the previous conclusions and provided further insight into the characteristics of the hydrogen bonds. From kinetic modeling of the gas-phase chemistry of this chemical system, Arrhenius behavior has been observed theoretically for the effect of temperature on the rate constants of the Lew-I detoxification mechanism. The accuracy of these calculations can be anticipated by the use of standard high-level methods and basis sets. Although there is no experimental support, our results and predicted reaction mechanisms will serve as useful information and basis for future experimental work to understand the detoxification by chelation of arsenicor heavy-metal-containing warfare agents in the presence of drugs containing basic amino acid, aliphatic amine, and hydroxyl groups as catalysts.

Step A k1

C3H 7S2 OH + AsC2 H 2Cl3 ⇌ [C3H 7S2 OH ··· AsC2 H 2Cl3] k2

Step B k3

[C3H 7S2 OH ··· AsC2 H 2Cl3] → AsS2 C4 H5OHCl + 2HCl

As shown in potential energy profiles, the mechanisms of BAL + Lew-I or BAL + Lew-I + NH3/H2O reactions are two-step mechanisms. The relatively high stabilization energies of the prereactant complex (PRC) suggest unambiguously that the mechanism is complex, with a first reversible step leading to PRC formation and a second irreversible step yielding the corresponding stable product and HCl. According to the reaction mechanism proposed above for Lew-I detoxification, k1 and k2 are the forward and reverse rate constants for the first step and k3 corresponds to the second step. Steady-state analysis leads to a rate coefficient for each overall reaction channel, which can be written as koverall =

k1k 3 k3 + k 2

CONCLUSION

(6) 3505



■

ASSOCIATED CONTENT

■

AUTHOR INFORMATION

Article

model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (18) Frisch, M. J.; et al. Gaussian 09, revision B.01; Gaussian, Inc., Wallingford, CT, 2009. (19) Frisch, Æ.; Frisch, M. J.; Trucks, G. W.; Clemente, F. R. Gaussian 03 User’s Reference, Gaussian Inc.: Pittsburgh, PA, 2003. (20) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chem. Rev. 1988, 88, 899−926. (21) Bader, R. F. W. Atoms in Molecules. A Quantum Theory; Oxford University Press: New York, 1994. (22) Cioslowski, J.; Nanayakkara, A.; Challacombe, M. Rapid evaluation of atomic properties with mixed analytical/numerical integration. Chem. Phys. Lett. 1993, 203, 137−142. (23) Cioslowski, J. A new robust algorithm for fully automated determination of attractor interaction lines in molecules. Chem. Phys. Lett. 1994, 219, 151−154. (24) Wilhelmy, L. Temperature dependence to rate theories. Pogg. Ann. 1850, 81, 422−499. (25) Logan, S. R. The origin and status of the Arrhenius equation. J. Chem. Educ. 1982, 59, 279−285. (26) Laidler, K. J. The development of the Arrhenius equation. J. Chem. Educ. 1984, 61, 494−504. (27) Eyring, H. The activated complex in chemical reactions. J. Chem. Phys. 1935, 3, 107−115. (28) Evans, M. G.; Polanyi, M. Some application of the transition state method to the calculation of reaction velocities, especially in solution. Trans. Faraday Soc. 1935, 31, 875−894. (29) Truhlar, D. G.; Isaacson, A. D., Garrett, B. C. Generalized transition state theory. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. 4, pp 65−137. (30) Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. J. Current status of transition-state theory. J. Phys. Chem. 1996, 100, 12771−12800. (31) Fernández-Ramos, A.; Ellingson, B. A.; Garrett, B. C.; Truhlar, D. G. Variational transition state theory with multidimensional tunneling. Rev. Comput. Chem. 2007, 23, 125−232. (32) Garrett, B. C., Truhlar, D. G. Transition state theory. In Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer, H. F., III, Eds.; Wiley: Chichester, U.K., 1998; Vol. 5, pp 3094−3104. (33) Fernández-Ramos, A.; Benjamin, A. E.; Meana-Pañeda, R.; Marques, M. C. J.; Truhlar, D. G. Symmetry numbers and chemical reaction rates. Theor. Chem. Acc. 2007, 118, 813−826. (34) Biegler-Konig, F.; Schonbohm, J.; Bayles, D. AIM2000, A program to analyze and visualize atoms in molecules. J. Comput. Chem. 2001, 22, 545−559. (35) Glendening, D. E.; Reed, A. E.; Carpenter, J. E.; Winhold, F. NBO, version 3.1; 1992. (36) Urban, J. J.; von Tersch, R. L. Conformational analysis of the isomers of lewisite. J. Phys. Org. Chem. 1999, 12, 95−102. (37) Barone, V.; Cossi, M.; Tomasi, J. Geometry optimization of molecular structures in solution by the polarizable continuum model. J. Comput. Chem. 1998, 19, 404−417. (38) Popelier, P. L. A.; Bader, R. F. W. The existence of an intramolecular C-H-O hydrogen bond in creatine and carbamoyl sarcosine. Chem. Phys. Lett. 1992, 189, 542−548. (39) Nowroozi, A.; Jalbout, A. F.; Roohi, H.; Khalilinia, E.; Sadeghi, M.; Leon, A. D.; Raissi, H. Hydrogen bonding in acetylacetaldehyde: Theoretical insights from the theory of atoms in molecules. Int. J. Quantum Chem. 2009, 109, 1505−1514. (40) Laidler, K. J. Chemical Kinetics, 3rd ed.; Harper and Row: New York, 1987. (41) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics; Prentice Hall: Englewood Cliffs, NJ, 1989.

S Supporting Information *

Seven tables and four figures showing conformational analysis, relative energy of conformers, Cartesian coordinates for optimized geometries in gas phase for all species, ZPVE, thermochemical parameters, and calculation of H-bond strength. This material is available free of charge via the Internet at http://pubs.acs.org. Corresponding Author

*E-mail [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS C.S. is grateful to the Council of Scientific and Industrial Research (CSIR), Government of India, for providing him research fellowships. C.S. is thankful to Dr. Bhaskar Mandal and Dr. Deepanwita Ghosh for their assistance and helpful discussions. A.K.D. is grateful to the Council of Scientific and Industrial Research (CSIR), Government of India, for a research grant under Scheme 03(1168)/10/EMR-II.

■

REFERENCES

(1) Marrs, T. C.; Maynard, R. L.; Sidell, F. R. Chemical Warfare Agents: Toxicology and Treatment; Wiley: New York, 1996. (2) Green, S. J.; Price, T. W. LVI.The Chlorovinylchloroarsines. J. Chem. Soc. Trans. 1921, 119, 448−453. (3) Webb, J. L. Arsenicals; Enzyme and Metabolic Inhibitors, Vol. 3; Academic Press: New York, 1966. (4) Peters, R. A.; Stocken, l. A.; Thompson, R. H. S. British AntiLewisite (BAL). Nature 1945, 156, 616−619. (5) Waters, L. L.; Stock, C. BAL (British Anti-Lewisite). Science 1945, 102, 601−606. (6) Allouche, A. R. Gabedit; http://gabedit.sourceforge.net/. (7) Stewart, J. J. P. MOPAC 2009, version 9.259W; Stewart Computational Chemistry; http://openmopac.net/. (8) Ali, S., et al. Avogadro: an open-source molecular builder and visualization tool, version 1.0.3; http://avogadro.openmolecules.net/ wiki/Main_Page/. (9) Zhao, Y.; Truhlar, D. G. Comparative DFT Study of van der Waals complexes: Rare-gas dimers, alkaline-earth dimers, zinc dimer, and zinc-rare-gas dimers. J. Phys. Chem. 2006, 110, 5121−29. (10) Schafer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (11) Peterson, K. A.; Puzzarini, C. Systematically convergent basis sets for transition metals. II. Pseudopotential-based correlation consistent basis sets for the group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements. Theor. Chem. Acc. 2005, 114, 283−296. (12) Xu, X.; Truhlar, D. G. Accuracy of effective core potentials and basis sets for density functional calculations, including relativistic effects, as illustrated by calculations on arsenic compounds. J. Chem. Theory Comput. 2011, 7, 2766−2779. (13) Bergner, A.; Dolg, M.; Kuechle, W.; Stoll, H.; Preuss, H. Ab initio energy-adjusted pseudopotential for elements of groups 13−17. Mol. Phys. 1993, 80, 1431−1441. (14) Gonzales, C.; Schlegel, H. B. An improved algorithm for reaction path following. J. Chem. Phys. 1989, 90, 2154−2162. (15) Gonzales, C.; Schlegel, H. B. Reaction path following in massweighted internal coordinates. J. Phys. Chem. 1990, 94, 5523−5527. (16) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3093. (17) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum 3506


A Computational Study of Detoxification of Lewisite Warfare Agents by

Recommend Documents