4486
J. Phys. Chem. A 2010, 114, 4486–4493
Aziridinium Ion Ring Formation from Nitrogen Mustards: Mechanistic Insights from Ab Initio Dynamics David J. Mann* Sterling Heights, Michigan, 48312 ReceiVed: August 17, 2009; ReVised Manuscript ReceiVed: February 22, 2010
Explicit solvent ab initio molecular dynamics simulations have been used to study the activation of the nitrogen mustards mechlorethamine and phosphoramide mustard into the aziridinium ion intermediates known to be involved in the alkylation of DNA. The simulations predict a concerted reaction occurring by way of neighboring-group participation with the nearby nucleophilic tertiary nitrogen, in agreement with previous theoretical studies. The calculated activation free energy of 20.4 kcal/mol for mechlorethamine agrees well with the value of 22.5 kcal/mol determined from available experimental data. Furthermore, these simulations indicate a dynamic transition state characterized by pronounced changes in the local water structure within the first hydration shell. Trajectories initiated from the transition-state region toward products reveal a reaction that proceeds directly to the solvent-separated ion pair. Failure to cross the transition state was observed in a small number of trajectories where the degree of coordination with the leaving group is less favorable for reaction. Direct vibrational excitation of mechlorethamine alone is not sufficient to induce reaction, even with energies far in excess of the activation energy. The simulations suggest that reactivity is strongly dependent on the degree of coordination and orientation of water molecules within the vicinity of the leaving group. The complete mechanism involving solvent reorganization, ionization of the C-Cl bond, and internal cyclization of the aziridinium ion ring was captured from elevated temperature simulations. Rate constants for aziridinium ion ring formation from both mechlorethamine and phosphoramide mustard (the alkylating moiety of cyclophosphamide) are calculated and compared with available phamacokinetic data. 1. Introduction Nucleophilic substitution reactions remain one of the most captivating and widely studied chemical reactions, due in large part to their representing a fundamental class of substitution reactions encountered in organic, inorganic, and biological chemistry. Despite being one of the first reaction types encountered in the study or organic chemistry, new discoveries on their rich mechanism and dynamics are being made through high-level experimental and theoretical studies.1-5 The body of literature on computational studies of nucleophilic substitution reactions has been devoted to the gas-phase bimolecular SN2 reaction and the solution-phase SN1 and SN2 reactions.6-10 The intramolecular nucleophilic displacement reactions have received much less attention despite their involvement in synthetic and medicinal chemistry. These reactions differ from the traditional SN1 mechanism in that the presence of a neighboring nucleophile may help facilitate ionization of the carbon-halogen bond. One reaction, in particular, involves the formation of cyclic quaternary amines formed from alkylating agents used in traditional chemotherapy regimens. The class of alkylating agents known as the nitrogen mustards are among the first and remain one of the most widely used antineoplastic drugs. Their clinical use evolved from observations during World War I of lymph node shrinkage and bone marrow suppression in soldiers exposed to sulfur mustard gas.11,12 Synthetic variants with less reactivity and toxicity were subsequently developed and introduced to clinical trials in the 1940s.11,12 The two classes of nitrogen mustards consist of the more reactive aliphatic (mechlorethamine, cyclophosphamide, * Corresponding author.
Figure 1. (A) Structures of the aliphatic nitrogen mustards. (B) Scheme illustrating the activation of mechlorethamine, followed by monoalkylation of the nucleophilic site of a nucleotide.
ifosphamide) and less reactive aromatic (chlorambucil, melphalan) mustards. (See Figure 1.) These agents exert their chemotherapeutic and cytotoxic effects through the alkylation of DNA.12-14 The seven nitrogen and six oxygen of guanine, the one and three nitrogens of adenine, and the three nitrogen of cytosine are all proposed nucleophilic sites susceptible to alkylation by nitrogen mustards and other alkylating agents.11
10.1021/jp9079553 2010 American Chemical Society Published on Web 03/11/2010
AZ+ Ring Formation from Nitrogen Mustards Several studies have demonstrated that the N7 position of guanine is the preferred alkylation site.15-17 One feature of the nitrogen mustards is that they are capable of both monoalkylation and bisalkylation of DNA, the latter resulting in either interstrand or intrastrand cross-linking. The alkylation of DNA has been shown to lead to cell-cycle arrest (through slowing down the G2 phase of the cell-cycle) as well as apoptosis (naturally programmed cell death).18 A key step in DNA alkylation via nitrogen mustards is their conversion into highly reactive aziridinium ion (AZ+) rings (Figure 1B). This prodrug activation culminating in chemically unstable substrates limits both their ease of administration and degree of toxicity.11 Mechlorethamine is very unstable in vivo, undergoing complete activation within minutes of administration, whereas chlorambucil and melphalan react much less readily, secondary to aromatic stabilization of the carbonium ions.12 Early kinetic studies by Williamson and Witten19 revealed a first-order mechanism for the activation of both classes of nitrogen mustards, with a reduced rate of formation among the aromatic agents. Since this original work, few experimental or theoretical studies on this or related reactions have been reported in the literature.20-24 The focus of the present study is on examining the mechanism and dynamics of AZ+ ring formation from nitrogen mustards through the use of explicit solvent ab initio molecular dynamics (AIMD) simulations. It is already known that the internal cyclization of both nitrogen and sulfur mustards involves intramolecular catalysis through neighboring group participation. However, less is known regarding whether the reaction proceeds through a concerted or stepwise mechanism, the energetic requirements for reaction, and the role of microscopic solvation on the key ionization step and on the formation and stability of post-transition state ionic complexes. Answering fundamental questions such as these is amenable to the accuracy and time scales affordable to AIMD simulations. Most of the studies on nucleophilic substitution reactions have been limited to gas-phase SN2 and aqueous-phase SN1 reactions. To the best of our knowledge, very few theoretical studies on intramolecular nucleophilic substitution reactions have been reported in the literature.17,20-22 Computational studies by Colvin,20 using the B3LYP/6-31G* level of theory and the polarizable continuum model (PCM) to include the effects of solvation, identified a transition state (TS) with an energy of 14.0 kcal/mol. A second pathway, involving a 1,2 shift of the leaving group, was also identified with a slightly lower barrier height of 13.3 kcal/mol. The AZ+ (1-methyl-1-(β-chloroethyl)ethylenimonium ion) formed from mechlorethamine is a cyclic quaternary amine formed by way of a cyclization-nucleophilic displacement reaction. The mechanistic and dynamical aspects of this transformation are the subject of this article. 2. Methods 2.1. Ab Initio Electronic Structure Methods. The AIMD simulations reported in the present study were carried out using density functional theory (DFT) with a plane-wave basis, as implemented in the VASP code.27 The generalized gradient approximation, as parametrized by Perdew and Wang (PW91), is used for the exchange-correlation energy. The one-electron wave functions are expanded in a plane-wave basis set using the Vanderbilt ultrasoft pseudopotentials of Kresse and Hafner.26 An energy cutoff of 396 eV was chosen with gamma point sampling of the Brillouin zone. Stationary point structures and relative energies of mechlorethamine and the ethylenimonium ion were compared with results from local orbital electronic structure calculations using the GAMESS program.25
J. Phys. Chem. A, Vol. 114, No. 13, 2010 4487 The simulation model consists of a single molecule of mechlorethamine solvated by 74 water molecules in a rectangular periodic box. The box dimensions of 12.0 × 12.0 × 18.0 Å were chosen such that a density of 1.0 g/cm3 is produced as well as to accommodate the shape of the solute molecule. Initial configurations for the AIMD simulations were obtained from equilibrium molecular dynamics simulations employing readily available parametrized force fields. The water molecules were modeled using the TIP3P water potential with the solute-solvent interaction parameters taken from the OPLS force field. (Note that mechlorethamine was constrained to its gas-phase equilibrium geometry throughout this equilibration run.) The structure obtained after 100 ps of molecular dynamics at a constant temperature of 310 K was used as the starting point for the AIMD simulations. To generate a configuration sampled from the ab initio potential energy surface, an additional equilibration run was performed with a single unconstrained AIMD simulation. The classical equations of motion were integrated using the velocity verlet algorithm for a total of 2 ps with an integration time step of 1.0 fs. A temperature of 310 K was maintained using the Nose´ algorithm28 with initial velocities sampled from a Boltzmann distribution. 2.2. Reaction Path Free Energy Simulations. A common approach for studying chemical reactions in solution involves calculation of the free energy profile along the reaction coordinate. Several methods have been developed and described extensively in the literature, including free energy perturbation (FEP)29-33 and thermodynamic integration (TI).34-37,33 For an excellent review of free energy calculations, see the article by Kollman.33 The TI approach is founded upon the master equation
∆G )
λ) dλ ∫01 〈 ∂H(x, ∂λ 〉λ
(1)
where x is the reaction coordinate and λ is a coupling parameter between states A and B. (In our case, A and B are the reactant and product states, respectively.) When the λ states refer to differing values of a fixed internal coordinate, the resulting free energy profile is termed a potential of mean force (PMF).37,38 Using the PMF bond contribution method,39 Pearlman has shown that while constraining the reaction coordinate along the reaction path the PMF due to these constraints can be calculated.37,40 In the resulting constraint force (CF) method, the derivative in the integrand of eq 1 can be written as
∂Hconstraint rij ∂rij ) (-Ficonstraint + Fjconstraint) ∂λ rij ∂λ
(2)
where Fconstraint is the net force on atom i, rij is the constraint i distance, and rij is the bond vector. Evaluation of the derivative in eq 2 amounts to calculating the net force between atoms i and j projected onto the bond axis, subject to the constraint rij. This algorithm is an alternative to umbrella sampling with the weighted histogram analysis method (WHAM).41 The formation of the AZ+ is accomplished through a nucleophilic displacement reaction with the neighboring tertiary nitrogen. The final products of the reaction are the AZ+ and chloride anion. As with previous studies of SN1 reactions in solution, the dissociating carbon-halogen bond has been chosen as the reaction coordinate. The construction of the PMF along this reaction coordinate amounts to running a series of nonequilibrium constrained AIMD simulations at various C-Cl distances, beginning at the minimum and extending the C-Cl bond up to and beyond
4488
J. Phys. Chem. A, Vol. 114, No. 13, 2010
Mann
TABLE 1: Comparison of Selected Internal Coordinates for Both Mechlorethamine and the Aziridinium Cation at the PW91/ Planewave-Pseudopotential and PW91/6-31+G(d,p) Levels of Theorya mechlorethamine coordinate C-Cl N-C(ethyl-Cl) N-C(methyl) C-C ∠Cl-C-C ∠C-C-N a
aziridinium cation
PW91/PW-PP(6-31+G(d,p)) 1.82(1.82) 1.46(1.47) 1.46(1.45) 1.52(1.53) 110.4(110.6) 108.8(110.5)
coordinate
PW91/PW-PP(6-31+G(d,p))
N-C(ring) C-C(ring) ∠N-C-C(ring) ∠C-N-C(ring)
1.50(1.50) 1.48(1.48) 60.6(60.3) 59.2(59.4)
Bond lengths and angles are in angstroms and degrees, respectively.
the dissociation threshold. The CF in eq 2 is averaged over the duration of the simulation and inserted into eq 1 for evaluation of the PMF. The reaction path PMF reported in this article was generated from simulations performed at a temperature of 310 K to mimic physiologic conditions. From the resulting PMF profile, the activation free energy, ∆G‡, can be determined and compared with experiment.19 3. Results and Discussion 3.1. Potential of Mean Force Simulations and Reaction Path Properties. A comparison of stationary point structures is made between both plane-wave and local orbital-based DFT predictions. Local minima from the explicit solvent model were generated by minimizing the final configuration obtained at the end of the equilibration run. The resulting minimum energy structure of mechlorethamine was compared with local orbital PW91 calculations carried out with a 6-31+G(d,p) basis set, with the PCM used to incorporate the effects of solvation. The internal coordinates of mechlorethamine, obtained from the explicit solvent plane-wave PW91 geometry optimization, were in excellent agreement with the PW91/6-31+G(d,p) minimum energy structures. A similar comparative analysis with the equilibrium geometry of the AZ+ ring was performed using the same methods described for mechlorethamine. Results for selected internal coordinates, shown in Table 1, were found to be very similar to the gas-phase values reported from a previous study.17 The constrained AIMD simulations were carried out starting from the equilibrium C-Cl bond length of 1.82 Å and extending up to a final distance of 3.35 Å, in varying intervals up to and beyond the TS. Each simulation was equilibrated for 500 fs, followed by a 2 ps production run, during which the average CF and other properties were accumulated. The location of the free energy barrier corresponds to the inflection point along the plot of the average CF versus reaction coordinate distance. In addition, the average CF should equal zero at the TS. After carefully analyzing the data, we chose to perform two additional constrained simulations at distances of 2.37 and 2.385 Å to identify the location of the TS more accurately. We found that, at a distance of 2.37 Å, the average CF was almost precisely equal to zero (-0.8 kcal/mol/Å). The PMF, as shown in Figure 2, was calculated using eq 1, and the averaged CFs were evaluated at each of the C-Cl bond distances. The calculated activation free energy ∆G‡ ) 20.4 kcal/mol, located at a distance of 2.37 Å along the reaction coordinate, can be compared with the experiments of Williamson and Witten.19 From the data provided in Chart 2 of their paper, we’ve calculated a first-order rate constant of 9.0 × 10-4 s-1. Using the transition state theory expression for the rate constant42
kTST )
( )
k BT ∆G‡ exp h kBT
(3)
Figure 2. Plot of the potential of mean force (PMF) along the reaction coordinate defined by the rupturing C-Cl bond distance. The dashed line corresponds to interpolated data.
the experimental activation free energy at a temperature of 310 K is determined to be 22.5 kcal/mol, which agrees within 10% of our calculated value. Note that in the experiments of Williamson and Witten, they were measuring the rate of disappearance of mechlorethamine. Given that AZ+ ring formation is expected to be rate-limiting, the measured overall rate constant should be very close, if not an upper limit, to the actual rate of AZ+ ring formation. Both the calculated activation free energy and C-Cl bond distance at the TS are very similar to those for the hydrolysis of tert-butyl chloride,43-45 perhaps the most widely studied SN1 reaction. For the hydrolysis of tertbutyl chloride under ambient conditions, Rossky et al. have calculated an activation free energy of 23 kcal/mol with corresponding locations of the C-Cl bond global minimum and TS of 1.8 and 2.3 Å, respectively.43 The fact that the calculated activation free energy for solvolysis of mechlorethamine (a primary alkyl halide) is so similar to the tertiary alkyl halide tert-butyl chloride illustrates the involvement of intramolecular catalysis by the neighboring nucleophile. Furthermore, the positive charge of the AZ+ cation is delocalized throughout the ring, giving this product carbocation-like character and contributing to its overall stability. Intramolecular nucleophilic substitution reactions involving neighboring-group participation are common in organic chemistry47-51 and have been observed to occur at both sp3 and sp2 carbon.52,53 These reactions exhibit first-order kinetics and may proceed through a two step process, analogous to the SN1 mechanism, or in a concerted fashion. Because primary carbocations (as would be produced by direct solvolysis of mechlorethamine) are often too unstable to form in solution (with energies of 30-35 kcal/mol higher than tertiary cations46), this would imply a concerted intramolecular nucleophilic displacement reaction, which is precisely the mechanism
AZ+ Ring Formation from Nitrogen Mustards
J. Phys. Chem. A, Vol. 114, No. 13, 2010 4489
Figure 3. Plots of the C-N bond distance versus time for simulations in which the C-Cl bond is constrained at distances of 2.02 (- - -), 2.37 (s), and 2.62 ( · · · ) angstroms.
observed from our simulations. The simulations reveal that dissociation of the C-Cl bond is accompanied by simultaneous nucleophilic attack by the neighboring nitrogen, culminating in formation of the heterocyclic ring. Animations of the simulation in the TS region reveal oscillations in the trajectory between an open and closed ring configuration, as illustrated from the plots shown in Figure 3. The back and forth oscillations between the reactant and product configurations occur with an approximate period of 400 fs, corresponding to a frequency of 80 cm-1. This behavior is characteristic of what one would expect of the dynamics near the TS and further illustrates the degree of structural complexity in the TS region. Although the reaction being studied exhibits first-order kinetics, the mechanism differs from that of the traditional SN1 reaction. Mechanistically, the reaction shares many of the features more commonly associated with SN2 reactions: simultaneous bond formation and dissociation, back-side nucleophilic attack, and inversion of configuration about the carbon bearing the leaving group (β carbon).49 An important problem in solution phase nucleophilic substitution reactions is the influence of solvation on bond heterolysis.54-59 As charge builds up on the leaving group within the TS region, hydrogen bonds begin to form with the leaving group that serve to stabilize the TS. The formation of these hydrogen bonds requires reorientation of water molecules and the breaking of water-water hydrogen bonds.59,60 There is a smooth evolution in the release of charge to the leaving group as the system passes through the TS. This was observed in previous studies of tertbutyl chloride heterolysis43 as well as in our present study. Furthermore, this shifting in the electron distribution along the reaction coordinate and subsequent release of charge occurs rather abruptly near the TS. Shown in Figure 4A is a plot of the fraction of charge released to the leaving group along the reaction coordinate. This plot was constructed from atomic charges obtained from constrained minimum energy structures along the reaction coordinate using our explicit water model. The atomic charges were generated from the charge density grid using the Bader charge analysis program.61-63 3.2. Elevated Temperature and Post Transition State Dynamics. To illustrate the effects of charge release on the solvent, we’ve calculated the average water-leaving group coordination number along the reaction coordinate. The results are shown in Figure 4B. The coordination number is defined as the number of nearest-neighbor water molecules with a Cl · · · H(water) distance
