Tautomeric transformation of temozolomide, their proton affinities and chemical reactivities: A theoretical approach

Journal of Molecular Graphics and Modelling 66 (2016) 76–84

Contents lists available at ScienceDirect

Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM

Tautomeric transformation of temozolomide, their proton affinities and chemical reactivities: A theoretical approach Wichien Sang-aroon a,b,∗ , Vithaya Ruangpornvisuti c , Vittaya Amornkitbamrung b,d,∗∗ a

Department of Chemistry, Faculty of Engineering, Rajamangala University of Technology Isan, KhonKaen Campus, KhonKaen 40000, Thailand Integrated Nanotechnology Research Center, Department of Physics, Faculty of Science, KhonKaen University, KhonKaen 40002, Thailand Supramolecular Chemistry Research Unit, Department of Chemistry, Faculty of Science, Chulalongkorn University, Pathumwan, Bangkok 10320, Thailand d Thailand Center of Excellence in Physics, CHE, Ministry of Education, Bangkok 10400, Thailand b c

a r t i c l e

i n f o

Article history: Received 16 November 2015 Received in revised form 21 March 2016 Accepted 24 March 2016 Available online 26 March 2016 Keywords: Temozolomide isomers Reaction energy Thermodynamic properties Rate constants Proton affinity Chemical reactivity

a b s t r a c t The gas-phase geometry optimizations of bare, mono- and dihydrated complexes of temozolomide isomers were carried out using density functional calculation at the M06 − 2X/6 − 31 + G(d,p) level of the theory. The structures and protonation energies of protonated species of temozolomide are reported. Chemical indices of all isomers and protonated species are also reported. Energies, thermodynamic quantities, rate constants and equilibrium constants of tautomeric and rotameric transformations of all isomers I1 ↔ TZM ↔ HIa ↔ HIb ↔ I2 ↔ I3 in bare and hydrated systems were obtained. © 2016 Elsevier Inc. All rights reserved.

1. Introduction Glioblastoma is one of the most severe forms of human cancer, which is a primary brain tumor [1]. Malignant gliomas including glioblastoma multiforme and anaplastic astrocytoma occur more frequently than other types of primary central nervous system tumours. Even if these cancers are treated with surgery, chemotherapy and radiation, survival is less than 1 year [2]. Additionally, the average lifespan for a MG patient postdiagnosis is 14.6 months with most patients experiencing tumor relapse and outgrowth within 7 months of initial radiation therapy [3–5]. Nowadays, temozolomide is one of the most widely used and effective chemotherapeutic drug for glioblastoma patients. Temozolomide is a prodrug of the alkylating agent 5-(3 methyltriazen-1-yl) imidazole-4-carboximide (MTIC) [6]. Temozolomide is an oral chemotherapy drug. Temozolomide readily crosses the blood-brain barrier and has a broad

∗ Corresponding author at: Department of Chemistry, Faculty of Engineering, Rajamangala University of Technology Isan, KhonKaen Campus, KhonKaen 40000, Thailand. ∗∗ Corresponding author at: Integrated Nanotechnology Research Center, Department of Physics, Faculty of Science, KhonKaen University, KhonKaen 40002, Thailand. E-mail addresses: [email protected] (W. Sang-aroon), [email protected] (V. Amornkitbamrung). http://dx.doi.org/10.1016/j.jmgm.2016.03.016 1093-3263/© 2016 Elsevier Inc. All rights reserved.

spectrum of antineoplastic activity [7–9]. However, although radiotherapy followed by accompanying and auxiliary temozolomide therapy showed encouraging results in patient survival [3], recent studies indicate a rather discomfiting result. Over 40% of patients undergoing chemotherapy and 55% of newly diagnosed cases do not benefit at all from temozolomide therapy [10,11]. Upon administration, temozolomide hydrolyzes spontaneously at physiological pH to the active metabolite (MTIC). MTIC is therefore hydrolyzed to methylhydrazine [12,13]. There are still a little amount of work has been done on temozolomide. The anticancer, antitumor activity and skin delivery potency of temozolomide ester prodrugs have been reported [14,15]. TiO2 nanostructured reservoir with temozolomide has been synthesized and structural evaluated [16]. The quantitative determination of temozolomide in bulk and capsule has been determined by UV spectrophotometric method [17]. XRD, FTIR and FT-Raman analysis have been used to characterize some cocrystals of temozolomide [18]. The possibilities of encapsulating temozolomide into polybutylcyanoacrylate nanoparticles by polymerisation was studied [19]. The physiochemical properties of temozolomide process related impurities and their structure were investigated [20]. Interaction of temozolomide in water has been investigated by means of theoretical approach [21]. The structure, vibrational and electronic spectra of temozolomide molecule has been elucidated by combined experimental and theoretical methods [22].

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

According to the vast pharmaceutical importance, an investigation has been made in this study to investigate transformation of temozolomide. The compound should be extensively investigated for their isomers, energetic, thermodynamic properties for their equilibria and kinetics. In the present work, the optimized molecular geometries of temozolomide isomers have been calculated using DFT method with M06 − 2X/6 − 31 + G(d,p) basis function. The tautomerization due to direct (bare system) and water − assisted proton transfers (hydrated system) were investigated. Protonation energies and stabilities of various protonated temozolomide species were also investigated. Energies of all temozolomide isomers, thermodynamic quantities, rate constants and equilibrium constants of their tautomeric reactions were reported. Relative stabilities and reactivities based on frontier molecular orbital energies (EHOMO and ELUMO ) of all isomers were also computed and discussed.

2. Computational methodology The dispersion-corrected functional M06–2X was used to perform calculation in this work. Even B3LYP is one of the most popular functional used in the calculation based on its reliable result and computational cost, however, it has been reported that B3LYP activation barriers overestimate experimental activation energy by 3 − 4 kcal/mol while M06-2X activation barrier underestimate the experimental value by 2.8–3.7 kcal/mol [23]. To understand the effect of basis set on the relative energy of temozolomide before selecting the suitable one, so, geometry optimization of the six isomers carried out at various basis sets were shown in Table 1. The zero − point corrected relative energies based on the M06 − 2X/6 − 31 + G(d,p), M06 − 2X/6 − 311 + +G(d,p) and M06 − 2X/aug − cc − pVQZ are in the same increasing order and the obtained values from the three basis sets are very close. The result is in agreement to the previous work [24] which the effect of basis function on relative energies of nitrosamine were computed based on B3LYP/6 − 31G(d), B3LYP/6 − 311 + +G(d,p) and B3LYP/aug-cc-pVQZ levels and the obtained values from these three levels were not much different (∼2 kcal/mol). Thus, in this work, the gas − phase structure of isomers of the temozolomide and the transition states of all transformation reactions were optimized using the hybrid functional of Truhlar and Zhao M06 − 2X [25] with 6 − 31 + G(d,p) basis function. The solvent effect was investigated by single–point calculations on the M06 − 2X/6 − 31 + G(d,p)-optimized gas-phase structures using the conductor–like polarizable continuum model (CPCM) and integrated-equation formalism polarizable continuum model (IEFPCM) [26–28]. All calculations were performed with the Gaussian 09 program [29]. The molecular graphics of all related species were generated with the MOLEKEL 4.3 program [30]. The reactivity indices including Mulliken electronegativity (␹), chemical hardness (), and electronic chemical potential () for all isomers of the temozolomide were computed using orbital energies of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) at the M06 − 2X/6 − 31 + G(d,p) level of theory. The chemical indices were derived from the first ionization potential (I) and electron affinity (A) of the N − electron molecular system with a total energy (E) and external potential (v(r៝ )) using the relations:  = −(∂E/∂N)vr៝ = 2 − ∼ = 1/2(I − A) and the first ion= 1/2(I + A),  = −(∂ E/∂N 2 ) ៝ ∼ vr

ization potential is I = E(N − 1) − E(N)and electron affinity is A = E(N) − E(N + 1)[31]. I and A were computed from the EHOMO and ELUMO using the relationship I = −EHOMO and A = −ELUMO according to the Koopmans theorem [32]. The rate constant at the standard temperature 298.15 K (k298 ), was computed based on transition state theory (TST). The rate constant derived from the Gibbs free

77

Fig. 1. Chemical structure and atomic labeling of temozolomide molecule.

energy of activation († G) computed based on Eq. (1) [33,34], and derived from activation energy († E) computed based on Eq. (2) [23]. k(T ) =

kB T exp(−† G/RT ) hc o

k(T ) = 

kB T QTS exp(−† E/RT ) h QREA

(1) (2)

= A exp(−† E/RT ) where kB is the Boltzmann’s constant, h is Plank’s constant, T is the absolute temperature, c◦ is the concentration factor, and R is the gas constant. QTS and QREA are the total partition functions of the transition states and reactants which composed of electronic, translational, rotational, and vibrational partition functions. The standard thermodynamic quantities including reaction enthalpy o ) and Gibbs free energy changes (Go ) for all reaction (H298 298 pathways were derived from the frequency calculations at the M06 − 2X/6 − 31 + G(d,p) level of theory. The equilibrium constant at standard temperature and pressure (K298 ) is computed based on o ) using a thermodynamic Eq. (3): Gibbs free energy changes (G298 K298 = exp(−Go /RT )

(3)

3. Results and discussion 3.1. Energetics and thermodynamics of transformational reaction of temozolomide The molecular structure and atomic labeling of temozolomide is shown in Fig. 1. For tautomeric isomers, the original temozolomide is named TZM. Imide isomers corresponding to the amide proton transfer to the three aromatic nitrogen atoms N1, N2 and N3 are defined as I1, I2 and I3 while two rotameric hydroxy imide isomers are called HIa and HIb, respectively. The M06 − 2X/6 − 31 + G(d,p) − optimized gas − phase structures of temozolomide isomers, their transformation reactions in bare, monohydrated and dihydrated systems are depicted in Fig. 2. The zero − point corrected relative energies derived from the M06 − 2X/6 − 31 + G(d,p), M06 − 2X/6 − 311 + +G(d,p) and M06 − 2X/aug − cc − pVQZ are in the same decreasing order. The relative stabilities based on the total energies derived from M06 − 2X/6 − 31 + G(d,p) level of the six isomers are in decreasing order: TZM (0.00) > HIa (12.15) > HIb (15.33) > I1 (45.36) > I2 (59.19) > I3 (67.76) kcal/mol. The relative energies due to single − point calculation using CPCM and IEFPCM at M06 − 2X/6 − 31 + G(d,p) level are very close. However, the relative stabilities of gas-phase optimized structures and single-point calculations in both solvation models

78

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

Table 1 Relative energy of tautomeric isomers of the temozolomide computed at various levels of theory in gas phase and aqueous phase. Isomers

e

TZM I1 HIa HIb I2 I3 a b c d e

In vacuo

In aqueous

EZPE a,b

EZPE a,b

M06-2X/6-31+G(d,p)

M06-2X/6–311 ++G(d,p)

M06 − 2X/aug − cc − pVQZ

CPCMc

IEF-PCMd

0.00 45.36 12.15 15.33 59.19 67.76

0.00 46.12 13.45 14.95 58.65 66.60

0.00 44.20 11.83 13.50 57.06 66.13

0.00 28.67 10.62 13.52 42.52 39.33

0.00 28.74 10.65 13.54 42.61 39.57

Relative to the most stable structure (TZM). In kcal/mol. The single-point CPCM calculation (␧ = 78.4) at M06-2X/6 − 31 + G(d,p) level. The single-point IEFPCM calculation (␧ = 78.4) at M06-2X/6 − 31 + G(d,p) level. The most stable neutral isomer.

are a bit different. The relative stabilities based on the total energies derived from CPCM//M06 − 2X/6 − 31 + G(d,p) (in square bracket) and IEF-PCM//M06 − 2X/6 − 31 + G(d,p) (in bracket) levels of the six isomers in solvent phases are in decreasing order: TZM [0.00]{0.00} > HIa [10.62]{10.65} > HIb [13.52]{13.54} > I1 [28.67]{28.74} > I3 [39.33]{42.61} > I2 [42.52]{42.61} kcal/mol. In addition, the relative energies computed in solvent media are more close to the most stable isomer TZM. I3 isomer is more stable than I2 isomer in solvent aqueous media. The competitive proton transfer reaction of TZM resulted in existence of either I1 or HIa isomers. It is obviously shown that HIa and HIb isomers are more stable than those of I1 isomer. Thus, tautomerization within the same functional group (amide↔hydroxy imide) is more energetic stable than those of formation of I1 isomer. A rotameric isomer; HIb is a bit less stable than those of HIa. Transformations of all isomers are presented as two reaction pathways namely (i) TZM↔I1 and (ii) TZM↔HIa↔HIb↔I2↔I3. The tautomerization between TZM↔I1 pair has occurred via transition

state TSt TZM I1 while TZM↔HIa, HIa↔HIb, HIb↔I2 and I2↔I3 pairs have occurred via transition states TSt TZM HIa, TSr HIa HIb, TSt HIb I2 and TSt I2 I3, respectively. The letter “t” and “r” in the transition state names are indicated to the proton transfer and rotational reactions, respectively. The letters “w” and “w2” in the names of isomers and transition states are indicated to the species in their monohydrate and dehydrate systems, respectively. The transition states correspond to rotational isomerization TSr HIaw HIbw and TSr HIaw2 HIbw2 for hydrated systems are not existed in the calculation. The M06 − 2X/6 − 31 + G(d,p) − optimized gas − phase structures of transition states involved in transformation reactions in bare, monohydrated and dihydrated systems are depicted in Fig. 3. The zero-point corrected activation energies, enthalpies and free energies and rate constant (k298 ) derived from transition state theory based on activation free energy (eq. 1) and activation energy (eq.2) are reported in Table 2. The rate constants computed based on eq. 1 and eq. 2 of all steps are in the same trend and not much dif-

Fig. 2. The M06 − 2X/6 − 31 + G(d,p) − optimized structures temozolomide isomers and their transformation reactions in (a) bare, (b) monohydrated (c) dihydrated systems.

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

79

Fig. 3. The M06 − 2X/6 − 31 + G(d,p)-optimized structures of transition states involved in transformation reactions in (a) bare, (b) monohydrated and (c) dihydrated systems.

ferent. The activation barriers and corresponding rate constants in hydrated systems are obviously different to those of bare systems. In the bare system, the activation energies of the reactions are in decreasing order: TZM↔I1 (53.62) > HIb↔I2 (52.75) > TZM↔HIa (44.81) > I2↔I3 (25.50) > HIa↔HIb (8.68) kcal/mol. For the complete set of the reaction ii; TZM↔HIa↔HIb↔I2↔I3, the rate determining step (RDS) is assigned to the HIb↔I2 step due its highest activation barrier (52.75 kcal/mol) among the other steps.

In the hydrated systems, the activation energies of the reactions are in increasing order: I2w↔I3w (32.02) > HIbw↔I2w (30.21) > TZMw↔I1w (26.30) > I1w↔HIaw (16.02) and I2w2↔I3w2 (88.21) > HIbw2↔I2w2 (27.59) > TZMw2↔I1w2 (26.35) > I1w2↔HIaw2 (12.05) kcal/mol, for mono- and dihydrated systems, respectively. In hydrated systems, isomer TZM can form two stable complexes based on the corresponding proton transfer paths. The forms wa/w2a and wb/w2b indicate the water

Table 2 Rate constants of tautomerization of all involved temozolomide isomers computed M06 − 2X/6 − 31 + G(d,p) level of theory. Reactionsa /systems

† E b,c

† H

Bare system TZM→TSt TZM I1→I1 I1→TSt I1 HIa→HIa HIa→TSr HIa HIb→HIb HIb→TSt HIb I2→I2 I2→TSt I2 I3→I3

53.62 44.81 8.68 52.75 25.50

Monohydrated system TZMw→TSt TZMw I1w→I1w I1w→TSt I1w HIaw→HIaw HIaw→TSr HIaw HIbw→HIbw HIbw→TSt HIbw I2w→I2w I2w→TSt I2w I3w→I3w Dihydrated system TZMw2→TSt TZMw2 I1w2→I1w2 I1w2→TSt I1w2 HIaw2→HIaw2 HIaw2→TSr HIaw2 HIbw2→HIbw2 HIbw2→TSt HIbw2 I2w2→I2w2 I2w2→TSt I2w2 I3w2→I3w2 a b c d e f

† Gb,d

k298 e



QTS /QREA

A

k298 e

52.33 44.03 9.04 52.33 25.87

56.41 46.33 8.12 53.85 25.53

2.75 × 10−29 6.72 × 10−22 6.97 × 106 2.07 × 10−27 1.19 × 10−6

2.60 4.75 1.21 1.36 3.83

9.56 × 10−3 9.17 × 10−1 2.61 × 10◦ 1.60 × 10−1 1.02 × 10◦

5.94 × 1010 5.70 × 1012 1.62 × 1013 9.96 × 1011 6.35 × 1012

7.65 × 10−29 3.81 × 10−20 8.43 × 106 2.91 × 10−27 4.93 × 10−6

26.30 16.02 −f 30.21 32.02

24.69 14.98 −f 28.92 32.10

28.80 17.41 −f 31.63 31.40

4.78 × 10−9 1.07 × 10◦ −f 4.05 × 10−11 5.91 × 10−11

2.70 3.35 −f 1.70 1.32

1.37 × 10−2 2.81 × 10−2 −f 9.86 × 10−2 1.10 × 10◦

8.50 × 1010 1.75 × 1011 −f 6.12 × 1011 6.84 × 1012

1.21 × 10−8 1.05 × 10◦ −f 7.45 × 10−11 3.05 × 10−11

26.35 12.05 −f 27.59 88.21

24.98 10.58 −f 25.68 85.66

28.69 14.04 −f 29.73 89.34

5.75 × 10−9 3.18 × 102 −f 9.92 × 10−10 1.99 × 10−53

1.03 2.95 −f 2.31 4.07

1.79 × 10−2 5.03 × 10−2 −f 1.74 × 10−1 1.85 × 10−2

1.11 × 1011 3.13 × 1011 −f 1.08 × 1012 1.15 × 1011

5.50 × 10−9 1.36 × 103 −f 1.50 × 10−10 1.01 × 10−53

b,c

The ending letter “r” and “t” of the name of transition state specify for rotational and proton transfer reaction, respectively. In kcal/mol. Total energy includes zero-point energy correction. Due to frequency calculation at M06 − 2X/6 − 31 + G(d,p) level. In s−1 . Undeterminable.

80

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

Table 3 Thermodynamic quantities and equilibrium constants of tautomerization of all involved temozolomide isomers computed M06 − 2X/6 − 31 + G(d,p) level of theory. Reactionsa /systems

E b,c

H b,c

Gb,c

K298

Bare system TZM→TSt TZM I1→I1 I1→TSt I1 HIa→HIa HIa→TSr HIa HIb→HIb HIb→TSt HIb I2→I2 I2→TSt I2 I3→I3

45.36 12.15 3.18 43.86 8.57

44.66 11.41 3.19 43.17 10.15

46.82 13.59 3.37 44.67 5.91

4.7 × 10−35 1.1 × 10−10 3.4 × 10−3 1.8 × 10−33 4.7 × 10−5

Monohydrated system TZMw→TSt TZMw I1w→I1w I1w→TSt I1w HIaw→HIaw HIaw→TSr HIaw HIbw→HIbw HIbw→TSt HIbw I2w→I2w I2w→TSt I2w I3w→I3w

33.71 7.89 −d 36.81 23.59

33.16 7.38 −d 36.21 23.64

34.15 8.92 −d 37.38 23.48

9.2 × 10−26 2.9 × 10−7 −d 4.0 × 10−28 6.1 × 10−18

Dihydrated system TZMw2→TSt TZMw2 I1w2→I1w2 I1w2→TSt I1w2 HIaw2→HIaw2 HIaw2→TSr HIaw2 HIbw2→HIbw2 HIbw2→TSt HIbw2 I2w2→I2w2 I2w2→TSt I2w2 I3w2→I3w2

28.63 4.65 −d 42.28 13.71

28.12 4.05 −d 42.51 13.93

30.24 6.07 −d 41.02 14.43

6.8 × 10−23 3.6 × 10−5 −d 8.4 × 10−31 2.7 × 10−11

a The ending letter “r” and “t” of the name of transition state specify for rotational and proton transfer reactions, respectively. b In kcal/mol. c Total energy includes ZPE correction. d Undeterminable.

complexes of isomer TZM for TZM↔I1 and TZM↔HIa reaction paths, respectively. Interaction between temozolomide and water has been theoretically observed [21]. It was found that the water prefers to bind to temozolomide simultaneously. The most potent proton donor site on temozolomide is the terminal amino group which corresponding to 1wb and 1w2b features in this work. As the reaction path which states that isomer TZM has two competitive proton transfer reactions; TZM↔I1 and TZM↔HIa, however, the faster rate is found in the second path either bare and hydrated systems. Thus, the reaction path ii is more energetically favorable than path i. The reaction HIa↔HIb has the largest rate constant k298 = 6.97 × 106 s−1 in bare system. As listed Table 2, the tautomeric transformations via direct (bare system) and waterassisted (monohydrated and dihydrated system) proton transfers of their corresponding transition states are remarkable different in their energies. The positive effect due to water assistance is obviously found in the reaction TZM↔I1 which the barrier is deduced from 53.62 to 26.30 and 26.35 kcal/mol from bare to monohydrate and dihydrated systems, respectively. Similarly, TZM↔I1 reaction possesses consecutively lowering barrier from 44.81 in bare system to 16.02 and 12.05 kcal/mol in monohydrated and dihydrated systems, respectively. Conversely, the consecutive increasing of barrier from 25.50 (bare) to 32.02 (monohydrated) and 88.21 (dihydrated) kcal/mol is found in the I2↔I3 reaction. Thermodynamic quantities and equilibrium constants of transformation reactions in all systems are reported in Table 3. All reactions steps are exothermic process. The HIa↔HIb step possesses the largest equilibrium constant (K298 = 3.4 × 10−3 ) in bare system while TZM↔HIa step has the largest K298 in both monohydrated (K298 = 2.9 × 10−7 ) and dihydrated (K298 = 3.6 × 10−5 ) systems. Based on the optimized transition state structure as depicted in Fig. 3, TSt I2 I3 shows the proton transfer from N2 to N3 in out-of-plane feature. It can be said that transformation of temozolomide in monohydrated system are more kinetically preferential than those of dihydrated system due to their activation barrier in each reaction steps in both systems are not much different. However, the transition state TSt I2w2 I3w2 has highest barrier may be due to the repulsive interaction between two water molecules confirms that mono water molecule is mostly suitable for transformational reaction.

The relative energy profile of the transformational reaction of temozolomide in various systems is depicted in Fig. 4. According to activation barrier and rate constant listed in Table 2 and thermodynamic quantities listed in Table 3 of all transformational reaction steps in bare and hydrated systems, the reaction path ii (TZM↔HIa) is more kinetically and thermodynamically preferential those of path i (TZM↔I1). So, the tautomerization of temozolomide is mainly according to TZM↔HI↔I2↔I3. The tautomerization via proton transfer based on mono-water assistance (monohydrated system) is found to be mostly favorable. The highest relative barrier is found in I2↔I3 for the complete set of TZM↔I3 reaction pathway both bare and hydrated systems. The rate determining step for the reaction in hydrated system is HIbw↔I2w and HIbw2↔I2w2 steps. Based on the M06 − 2X/6 − 31 + G(d,p) computations of different structural models (the bare, monohydrated and dihydrated structures), the conformations for the isomers TZM (in hydrated form b), HIa and I3 are not changed as shown in Fig. 5. Most of the isomers show clearly superimposed of their planar structures in different models except isomer I2. The most deviation in side chain amide (−CONH2 ) group is found in I1 and HIb isomers. Due to the change of bare to monohydrated and to dihydrated system, Fig. 5 illustrates that the side chain amide group of I1 and HIb isomers is mainly due to water complex. On the other hand, the protonation at N2 atom in I2 isomer induces the change in molecular planarity but not in the change of side chain. Similar result is found for TZM in their hydrated wa and w2a forms. Thermodynamic quantities of the hydration reaction of the hydrated temozolomide isomers are listed in Table 4. The pre − organization energies, enthalpies and free energies of temozolomide isomers computed based on the structure changes from the isolated to the hydrated forms. The binding energies, enthalpies and free energies of hydration reaction of temozolomide isomers are also reported. The highest pre-organization energies found in I3→I3w (19.89) and I3→I3w2 (24.23) kcal/mol which states that I3 isomer has significantly changed its structure due to hydration reaction. The other isomers possess pre-organization energies in narrow range 0.31 − 5.81 kcal/mol and 0.80 − 4.53 kcal/mol, respectively for monohydrated and dihydrated systems. Similarly, the reactions I3→I3w and I3→I3w2 possess highest binding energies 19.09 and 22.13 kcal/mol, respectively. This indicates that temozolomide isomers (except I3) can stabilize their structural planarity during hydration process with two water molecules which a little bit value of pre − organization may be due to a little change of amide moiety. It is found that only I1→I1w is exothermic reaction but the remains are endothermic. As the binding free energies of all hydration reaction of temozolomide isomers, the hydration reaction of all isomers is non − spontaneous process. The negative values of pre − organization free energies of all isomers in dihydrated systems indicate spontaneity of the pre − organization of temozolomide isomers but the results are converse in monohydrated systems. 3.2. Protonated species and their protonation energies There are six protonated species found and their relative energies are shown in Table 5. The name of protonated species are therefore defined as their corresponding protonated atoms (see Fig. 1). The M06 − 2X/6 − 31 + G(d,p) − optimized gas-phase structures of the six protonated species are also depicted in Fig. 6. The relative stabilities based on the total energies in gas-phase are in decreasing order: O1 (0.00) > N2 (9.60) > N1 (15.25) > N4 (18.92) > N3 (26.93) > O2 (34.33) kcal/mol. The relative stabilities based on the total energies in derived from CPCM//M06 − 2X/6 − 31 + G(d,p) (in square bracket) and IEF − PCM//M06 − 2X/6 − 31 + G(d,p) (in bracket) levels of the

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

81

Fig. 4. Relative energy profile of transformation reactions of temozolomide in bare (black line), monohydrated (red line) and dihydrated systems (purple line). For reaction TZM↔I1 and TZM↔HIa pathways, the relative energies are related to TZM wa/TZM w2 form and TZM wb/TZM w2b form, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

six protonated species in solvent phases are in decreasing order: O1 [0.00]{0.00} > N1 [0.89]{1.04} > N2 [2.00]{2.08} > N3 [5.41]{5.58} > N4 [14.92]{14.94} > O2 [18.31]{18.50}. The most stable protonated species is species O1 which is protonated at amide oxygen O1 atom. The result is in agreement to the stability of tautomeric isomers which HIa and HIb are the second and the third most stable isomers due to proton transfer from amide nitrogen to amide oxygen. The protonation energies of the six protonated species are also listed in Table 5. Based on the protonation energies of all six species, the order of protonation preference at various N and O atoms in gas phase are in decreasing order: O1 (−217.02) > N2 (−207.42) > N1 (−201.77) > N4 (−198.10) > N3 (−190.09) > O2 (−182.69) kcal/mol. The protonation energies of

all six species derived from CPCM//M06 − 2X/6 − 31 + G(d,p) (in square bracket) and IEF-PCM//M06 − 2X/6 − 31 + G(d,p) (in bracket) levels of the six isomers in solvent phases are in decreasing order: O1 [−274.81]{−274.79} > N1 [−273·77]{ − 273.90} > N2[−272·73]{ − 272.79} > N3 [−269·23]{ − 269.38} > N4 [−259.88]{−259.87} > O2 [−256·31]{ − 256.48} kcal/mol, respectively. Protonation at carbonyl oxygen O2 is the least preferential protonation site of temozolomide in both media. To elucidate the preferential protonation site, the natural population analysis (NPA) of atomic charge of temozolomide molecule was performed. As the NPA charge of N and O atoms of temozolomide molecule listed in Fig. 6, the most negative charge is N4 (−0.899) atom followed by O1 (−0.612), O2

Fig. 5. Superimpose pictures of the geometrical changes of the tautomeric isomers in bare, monohydrated and dihydrated systems. For clarity, water molecules are not displayed.

82

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

Table 4 Thermodynamic quantities of the hydration reaction of the hydrated temozolomide isomers computed M06 − 2X/6 − 31 + G(d,p) level of theory. Reactions/systems

Binding

Pre–organization E298,preorg a,b

H◦ 298,preorg a

Monohydrated system TZM→TZMwa TZM→TZMwb I1→I1w HIa→HIaw HIb→HIbw I2→I2w I3→I3w

0.82 0.31 4.01 0.67 3.17 5.81 19.89

0.19 0.34 4.04 0.56 3.11 6.51 19.03

Dihydrated system TZM→TZMw2a TZM→TZMw2b I1→I1w2 HIa→HIaw2 HIb→HIbw2 I2→I2w2 I3→I3w2

3.19 0.80 4.24 2.13 4.50 4.34 24.23

11.09 8.18 11.66 9.35 12.21 11.76 22.26

a b c

G◦ 298,preorg a

EZPE,binding a,b,c

H◦ ZPE,binding a,b,c

G◦ ZPE,binding a,b,c

2.29 1.82 4.10 0.91 3.16 5.01 21.67

9.84 7.32 −1.81 5.58 11.12 4.07 19.09

9.25 6.51 −2.25 5.22 11.08 4.12 17.61

18.88 16.79 6.21 14.21 19.22 11.93 29.50

−29.02 −30.91 −26.79 −28.69 −26.72 −25.83 −39.98

18.71 14.11 1.99 11.20 18.57 17.00 22.13

17.47 12.69 0.93 10.11 17.46 16.79 20.57

36.05 32.11 19.46 28.52 36.03 32.38 40.91

In kcal/mol. The ZPE corrected energies. Defined as E[Iw] − EI, Iw = hydrated from, I = free form.

Table 5 Relative energies and protonation energies of all involved protonated temozolomide species computed at M06 − 2X/6 − 31 + G(d,p) level of theory. Protonated species

Relative energies In vacuo

N1 N2 N3 N4 O1d O2 a b c d

Protonation energies In aqueous

In vacuo

In aqueous

C-CPM

IEF-PCM

C-CPM

IEF-PCM

EZPE a,b

Ewater a,b,c

Ewater a,b,d

EZPE a

Ewater a,b

Ewater a,c

15.25 9.60 26.93 18.92 0.00 34.33

0.89 2.00 5.41 14.92 0.00 18.31

1.04 2.08 5.58 14.94 0.00 18.50

−201.77 −207.42 −190.09 −198.10 −217.02 −182.69

−273.77 −272.73 −269.23 −259.88 −274.81 −256.31

−273.90 −272.79 −269.38 −259.87 −274.79 −256.48

In kcal/mol. The single − point CPCM calculation (␧ = 78.4) at M06 − 2X/6 − 31 + G(d,p). The single − point IEFPCM calculation (␧ = 78.4) at M06 − 2X/6 − 31 + G(d,p). The most stable species.

(−0.594), N1 (−0.473), N2 (−0.163) and N3 (0.032) atoms, respectively. The most favorable site for protonation is O1 (species O1) which possesses the second most negative charge. Even though N4 has the most negative charge but however protonation at N4 led its protonated species is the third least stable species due to unfavorable NH3 + group. In addition, protonation at N3 which is the most positive charged atom, shows that species N3 is second least protonated species. It is suggested that charges of N and O atoms are not related to the stabilities of protonated species of temozolomide. However, the atomic charges of N and O atoms and tautomerization reaction seem to be related. The atomic charge of O1 (−0.612) which is more negative than N1 (−0.473) led to faster rate of TZM↔HIa than those of TZM↔I1 step.

3.3. Frontier molecular orbital energies and chemical indices

Table 6 The ELUMO and EHOMO energies, EHOMO−LUMO (frontier molecular orbital energy gap) and their chemical indices computed at M06 − 2X/6 − 31 + G(d,p) level of theory. EHOMO a

EHOMO−LUMO a

a,b

a,b,c

a,d

Neutral form −5.31 TZM −5.25 I1 −5.31 HIa −5.31 HIb −5.36 I2 −5.80 I3

−8.71 −8.71 −8.71 −8.68 −8.68 −8.54

3.40 3.46 3.40 3.37 3.32 2.75

1.70 1.73 1.70 1.69 1.66 1.37

−7.01 −6.98 −7.01 −6.99 −7.02 −7.17

7.01 6.98 7.01 6.99 7.02 7.17

Protonated form −5.28 N1 −5.31 N2 −5.69 N3 −5.58 N4 −5.44 O1 −5.41 O2

−8.73 −8.71 −8.63 −8.79 −8.71 −8.92

3.46 3.40 2.94 3.21 3.27 3.51

1.73 1.70 1.47 1.61 1.63 1.76

−7.01 −7.01 −7.16 −7.18 −7.07 −7.17

7.01 7.01 7.16 7.18 7.07 7.17

Isomers

a

The energies of the lowest unoccupied molecular orbital (ELUMO ) and highest occupied molecular orbital (EHOMO ) and frontier molecular orbital energy gap (EHOMO−LUMO ) at the M06 − 2X/6 − 31 + G(d,p) level of tautomeric and protonated species of temozolomide are listed in Table 6. Their chemical indices, chemical hardness, chemical potential and Mulliken electronegativity of all temozolomide isomers are also shown in Table 6. These chemical indices are important tools to study the relative stabilities of different isomers of a molecular system [32]. How-

b c d

ELUMO a

In eV. Chemical hardness, EHOMO−LUMO /2. Electronic chemical potential,  = (EHOMO + ELUMO )/2. The Mulliken electronegativity,  = (EHOMO + ELUMO )/2.

ever, even though Koopmans theorem is not valid for DFT method, but the previous paper has clarified why long-range corrected (LC) density functional theory gives orbital energies quantitatively [35]. It has been found that only LC exchange functionals

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

83

Fig. 6. The M06 − 2X/6 − 31 + G(d,p) − optimized structures of various protonated temozolomide species. NPA charges on N and O atoms of temozolomide (TZM) are displayed.

were found to give the orbital energies close to the minus ionization potentials (IPs) and electron affinities (EAs) while other functionals considerably underestimate them. Additionally, the other work reported that the LC density functionals; LC-BOP and LCgau-BOP reproduced HOMO, LUMO and energy gaps better than other density functionals [36]. The negative of HOMO and LUMO energies were compared with the IPs and EAs, respectively, using CCSD(T)/6 − 311 + +G(3df,3pd) for 113 molecules, and found LC functionals to satisfy Koopmans’ theorem. Utilization of DFT methods to compute chemical indices derive from HOMO, LUMO and energy gaps were found in various molecular systems [24,31,37,38]. According to the Pearson’s maximum hardness principle (MHP) which states that the minimum energy structure has the maximum chemical hardness, so, TZM should has the maximum chemical hardness but it is the second most maximum value. The relative stabilities and chemical hardness of all isomers of temozolomide based on their frontier molecular orbital energy gap are in the same order and, are in decreasing order: I1 > TZM ∼ HIa > HIb > I2 > I3. Similarly, the protonated species O1 should has the maximum chemical hardness among other protonated species, but it is the third minimum value. The relative stabilities and chemical hardness of all protonated species of temozolomide based on their frontier molecular orbital energy gap are in the same order and, are in decreasing order: O2 > N1 > N2 > O1 > N4 > N3. The Mulliken electronegativity of all temozolomide isomers are very high, higher up to 7.0, but however not much different. Isomer I3 has the highest Mulliken electronegativity which is higher than the others by 0.15 − 0.19 eV. Protonated species of temozolomide possess electronegativity close to those of tautomeric isomers. Isomers N4 followed by O2 and N3 have the most, second and the third highest electronegativity which is higher than the others by 0.11 − 0.17 eV. Fig. 7 illustrates the change of energy gap and the orbital energies at HOMO − 5 to LUMO + 2 of tautomeric reaction isomers (I1↔TZM↔HIa↔HIb↔I2↔I3) and protonated species (species O1,

Fig. 7. The orbital energies at HOMO − 5 to LUMO + 2 of (a) tautomeric isomers and (b) protonated species.

O2 and N1 − N4). The LUMO + 2, LUMO + 1, LUMO and HOMO energies of all isomers are mostly stable except for LUMO of I3 which is close to −6.0 eV. The HOMO − 1 of the neighboring isomers I1, TZM, HIa and HIb are almost the same values which close to −9.0 eV but of I2 and I3 are deduced to −10.0 eV. The HOMO − 1 and HOMO − 2 of all isomers are remarkably separated except for I2 and I3 when their HOMO − 1 to HOMO − 5 are almost the same values. It is shown that these four isomers I1, TZM, HIa and HIb whose their HOMO and HOMO − 1 are hardly ever different, are is obviously more stable than the rest of two isomers I2 and I3. The stabili-

84

W. Sang-aroon et al. / Journal of Molecular Graphics and Modelling 66 (2016) 76–84

ties due to frontier MO energies are in agreement to the results in Table 6 which indicates that I3 and I2 isomers are the least and second least stable isomers. Similarly, HOMO and HOMO − 1 of all protonated species are almost close but these two states are obviously separated which found in species N3 and N2. It is indicated that species N2 and N3 are lesser stable than the others. As listed in Table 6, species N3 is the least stable isomer but however, N2 is the most third stable species. 4. Conclusion The bare, mono− and dihydrated complexes of temozolomide isomers were carried out using density functional calculation at the M06 − 2X/6 − 31 + G(d,p) level of the theory. The relative stabilities based on the total energies of the six isomers are in decreasing order: TZM > HIa > HIb > I1 > I2 > I3. The relative stabilities and chemical hardness of all isomers of Temozolomide based on their frontier molecular orbital energy gap are in decreasing order: I1 > TZM ∼ HIa > HIb > I2 > I3. The transformation of temozolomide follows as I1↔TZM↔HIa↔HIb↔I2↔I3. The relative activation energies of the reactions in bare systems are in increasing order: TZM↔I1 > HIb↔I2 > TZM↔HIa > I2↔I3 > HIa↔HIb. The barriers of the transformation reactions in mono− and dihydrated systems are much lower than those obtained via the bare system except TSt I2w I3w2. The relative activation energies of the reactions in mono− and dihydrated systems are in increasing order: HIb↔I2 > TZM↔I1 > I2↔I3 > TZM↔HIa. The numbers of explicit water and the location of water molecules bound to various site of Temozolomide show significantly affect to the tautomerization of temozolomide. In the presence of water molecule or in solution media, temozolomide can convert preferentially from original form (TZM) to hydroxy imide form (HIa,HIb) led to decreasing amount of the original active form for use in pharmaceutical aspect. In addition, the protonation behavior of Temozolomide varies depend on the position of N and O atoms. The most stable protonated species is hydroxy amide (O1) while protonation at N3 and N2 causes the species N3 and N2 possess the least and the third most stable compare to the others. It suggests that temozolomide becomes lesser active when presented in protonated form especially as its species N3. Acknowledgements This work was financially supported by National Research University Project (Khon Kaen University) through the research grant no. PD.57001 to WS was acknowledged. Institute of Research and Development, Rajamangala University of Technology Isan and Department of Chemistry, Faculty of Engineering, Khonkaen campus was also acknowledged. References [1] D. Zhang, A. Tian, X. Xue, M. Wang, B. Qiu, A. Wu, Int. J. Mol. Sci. 13 (2012) 1109–1125. [2] H.S. Friedman, T. Kerby, H. Calvert, Clin. Cancer Res. 6 (2000) 2585–2597.

[3] R. Stupp, W.P. Mason, M.J. Van Den Bent, et al., New Engl. J. Med. 352 (10) (2005) 987–996. [4] L.M. DeAngelis, N. Engl, J. Med, (2) (2001) 114–123. [5] E.R. Laws, I.F. Parney, W. Huang, et al., J. Neurosurg. 99 (3) (2003) 467–473. [6] G.C. De Gast, D. Batchelor, M.J. Kersten, et al., Br. J. Cancer 882 (2003) 175–180. [7] E.S. Newlands, M.F. Stevens, S.R. Wedge, R.T. Wheelhouse, C. Brock, Cancer Treat. Rev. 23 (1997) 35–61. [8] M.J.M. Darkes, G.L. Plosker, B. Jarvis, Am. J. Cancer 1 (2002) 55–80. [9] L. Tentori, G. Graziani, Curr. Med. Chem. 16 (2009) 245–257. [10] M.M. Mrugala, M.C. Chamberlain, F. Hutchinson, Nat. Clin. Practice Oncol. 5 (8) (2008) 476–486. [11] R.O. Mirimanoff, T. Gorlia, W. Mason, et al., J. Clin. Oncol. 24 (16) (2006) 2563–2569. [12] M.M. Mrugala, M.C. Chamberlain, F. Hutchinson, Nat. Clin. Pract. Oncol. 5 (8) (2008) 476–486. [13] B.J. Denny, R.T. Wheelhouse, M.F.G. Stevens, L.L.H. Tsang, J.A. Slack, Biochemistry 33 (31) (1994) 9045–9051. [14] M.J.M. Darkes, G.L. Plosker, B. Jarvis, Am. J. Cancer 1 (2002) 55–80. [15] P. Suppasansatorn, G. Wang, B.R. Conway, W. Wang, Y. Wang, Cancer Lett. 244 (2006) 42–52. [16] T. Lopez, J. Sotelo, J. Navarrete, J.A. Ascencio, Opt. Mater. 29 (2006) 88–94. [17] A.A. Razak, Sk. Masthanamma, B. Omshanthi, V. Suresh, P. Obulamma, Int. J. Pharm. Sci. Res. 4 (2013) 1419–1423. [18] N.J. Babu, P. Sanphui, A. Nangia, Chem. Asian J. 7 (2012) 2274–2285. [19] Xin-Hua Tian, Xiao-Ning Lin, F. Wei, W. Feng, Zhi-Chun Huang, P.B. Wang, L. Ren, Y. Diao, Int. J. Nanomed. 6 (2011) 445–452. [20] M. Laszcz, M. Kubiszewski, L. Jedynak, M. Kaczmarska, L. Kaczmarek, W. Luniewski, K. Gabarski, A. Witkowska, K. Kuziak, M. Malinska, Molecules 18 (2013) 15344–15356. [21] O.E. Kasende, A. Matondo, M. Muzomwe, J.T. Muya, S. Scheiner, Comput. Theor. Chem. 1034 (2014) 26–31. [22] S.A. Bhat, S. Ahmad, J. Mol. Struct. 1099 (2015) 453–462. [23] H. Dacosta, M. Fan, Rate Constant Calculation of Thermal Reactions: Methods and Application, Wiley, Hoboken, NJ, 2012. [24] V. Ruangpornvisuti, Int. J. Quant. Chem. 109 (2009) 275–284. [25] Y. Zhao, D.G. Trulhar, Thero. Chem. Acc. 120 (2008) 215–241. [26] J. Tomasi, B. Mennucci, E. Cances, J. Mol Struc, THEOCHEM 464 (1999) 211–226. [27] E. Cancès, B. Mennucci, J. Tomasi. J. Chem. Phys. 107 (1997) 3032–3039. [28] B. Mennucci, J. Tomasi, J. Chem. Phys. 106 (1997) 5151–5158. [29] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr, J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, O. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian09 Revision A.02, Gaussian Inc., Wallingford, CT, 2009. [30] P. Flükiger, H.P. Lüthi, S. Portmann, J. Weber, MOLEKEL 4.3, Swiss Center for Scientific Computing, Manno, 2000. [31] B. Wanno, V. Ruangpornvisuti, Chem. Phys. Lett. 415 (2005) 176–182. [32] T. Koopmans, Physica 1 (1933) 104–113. [33] J.W. Ochterski, Thermochemistry in Gaussian, Gaussian, Pittsburgh, 2000. [34] G. Bravo-Perez, J.R. Alvarez-Idaboy, A. Cruz-Torres, M.E. Ruiz, Phys. Chem. A 106 (2002) 4645–4650. [35] T. Tsuneda, J.-W. Song, S. Suzuki, K. Hirao, J. Chem. Phys. 133 (2010), 174101–174109. [36] R. Kar, J.W. Song, K. Hirao, J. Comput. Chem. 34 (11) (2013) 958–964. [37] B. Wanno, V. Ruangpornvisuti, J. Mol. Struct (Theochem.) 775 (2006) 113–120. [38] J. Baldenebro-López, J. Castorena-González, N. Flores-Holguín, J. Almaral-Sánchez, D. Glossman-Mitnik, Int. J. Mol. Sci. 13 (2012) 16005–16019.