- Email: [email protected]
The molecular structure of disulfiram and its complexation with silica. A quantum chemical study
Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
Contents lists available at ScienceDirect
Journal of Molecular Structure: THEOCHEM journal homepage: www.elsevier.com/locate/theochem
The molecular structure of disulfiram and its complexation with silica. A quantum chemical study Thomas Sandberg, Jessica Rosenholm, Matti Hotokka * Department of Physical Chemistry, Åbo Akademi University, Porthansgatan 3-5, 20500 Åbo (Turku), Finland
a r t i c l e
i n f o
Article history: Received 30 January 2008 Received in revised form 8 April 2008 Accepted 8 April 2008 Available online 16 April 2008 Keywords: Disulfiram Conformation Ab initio calculations DFT Proton affinity Silica complex
a b s t r a c t The complexation between disulfiram, more commonly known by trade names such as Antabus or Antabuse, and amorphous silica was studied by quantum chemical calculations due to extensive adsorption observed experimentally. The conformational behavior of the drug molecule was studied computationally at molecular mechanics (MM), Hartree-Fock (HF)/6-31G* and DFT/B3LYP/TZVP level. Two almost equally stable conformers were observed with the energy difference of only 0.8 kJ/mol when comparing the Gibbs’ free energies (DG) at 25 °C. The Mulliken and Löwdin charges and the charges based on an electrostatic potential fit (ESP) were studied at HF/6-31G* level. Protonation of the optimized structures and complexation with geminal silanols were studied at DFT/ B3LYP/TZVP level of theory in combination with the MARI-J approximation. Three groups of protonated structures were obtained and the proton affinity of the most stable protonated structure was calculated to 957.6 kJ/mol. In the complexation study four different kinds of complexes did converge. The BSSE-corrected electronic complexation energy of the most stable complex was calculated to 13.6 kJ/mol. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Disulfiram, full IUPAC name 1-(diethylthiocarbamoyldisulfanyl)-N,N-diethyl-methanethioamide, more commonly known as the alcohol deterrent by trade names such as Antabus or Antabuse, is a widely applied aldehyde dehydrogenase inhibitor to support treatment of chronic alcoholism. Disulfiram has also been shown to induce apoptosis in a number of tumor cell lines [1–5] and in the reversal of fungal drug-resistance, as recently reviewed [6]. As many other anticancer agents, disulfiram is practically insoluble in water. Whereas low water solubility of drugs impedes the ability of drugs to be administered through the intravenous route, novel delivery systems for hydrophobic therapeutic anticancer drugs has received significant attention [7]. Ceramic carrier materials are attracting increased attention due to easy preparation routes via the sol–gel process which allows them to be prepared with desired shape, size and porosity. The mostly employed sol–gel derived material in such applications is SiO2, silica, which is nontoxic, biocompatible and bioerodible and moreover an essential component of cells throughout the human body. These properties provide sol–gel processed porous silica materials considerable potential as carriers for drug delivery. Due to their great adsorption properties, mesoporous silicas synthesized using supramolecular surfactant aggregates as structure directing agents [8,9] have at* Corresponding author. E-mail address: mhotokka@abo.fi (M. Hotokka). 0166-1280/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.theochem.2008.04.007
tracted lots of recent interest in applications related to immobilization and release of active substances, especially drugs [10]. Preliminary adsorption tests carried out in our laboratory showed that as much as C = 0.93 lmol m2 could be adsorbed from cyclohexane to a model mesoporous silica with 3 nm pores, corresponding to a loading degree of about 30 wt% [11]. Cyclohexane provides a non-aquatic and non-polar environment, which would best correspond to theoretical calculations. While most drugs are weak acids or bases, the disulfiram molecule does not contain any typical functional groups, and consequently the interactions with solid surfaces are not easily foreseen. Hence, we have used quantum chemical calculations in order to investigate the molecular reasons for the strong adsorption of disulfiram to silica. To the best of our knowledge, neither theoretical nor experimental structural analysis has been published to date. It is important to determine the minimum energy structure in order to understand the behavior of disulfiram, e.g., in adsorption on silica surfaces. An analogous study of organic aromatic probes on silica has recently been published [12]. 2. Computational methods 2.1. The minimum energy structures In order to find the minimum energy structure of disulfiram, an initial potential energy surface scan about the three middle bonds in the NACASASACAN-skeleton was performed using a step size
58
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
of 30 degrees at molecular mechanics (MM) level using the Tripos Force Field [13] as implemented in SYBYL 7.3 [14]. The 11 energetically most favourable structures (MM energies below 42.3 kJ/mol) were optimized with the program GAMESS version 22 Feb 2006 [15] by using Hartree-Fock (HF) theory [16] with the basis set 631G* [17–19]. Two almost equally stable conformers of disulfiram, structures 1 and 2, were found. Also Mulliken and Löwdin charges, electrostatic potential fit (ESP) charges and bond orders were calculated with GAMESS at HF/6-31G* level of theory. The conformers were then reoptimized for vibrational analysis using the TURBOMOLE program package version 5.9 [20–22] and density functional theory (DFT) with the B3LYP hybrid exchangecorrelation functional [23–25] in combination with the multipole accelerated resolution of identity (MARI-J) approximation [26– 28] and the TZVP basis set [29] for all atoms, as implemented in the TURBOMOLE program package. Vibrational analysis was carried out for the B3LYP/TZVP optimized structures in order to prove that all optimized structures are real minima on the potential energy surface and to obtain the thermodynamic contributions. Frequencies were scaled with a factor of 0.9614 [30] when calculating the thermodynamic contributions at 25 °C. 2.2. Protonation Ten different protonated structures were constructed by placing the proton close to the nitrogen or sulfur atoms, which are expected to be the reaction sites. The protonated structures were optimized using TURBOMOLE 5.9 at DFT/B3-LYP/TZVP level of theory in combination with the MARI-J approximation. Three groups of protonated structures were obtained and the most stable one was studied in detail at the B3LYP/TZVP level taking into account thermodynamic contributions. As the third type of protonated structures resulted in a splitting into two molecular units, a vibrational analysis was performed also on that system. 2.3. Complexation As the aim was to study disulfiram interaction with amorphous silica surfaces, also the complexation between the energetically most favourable conformation of disulfiram and a silica model was studied in detail at the B3LYP/TZVP level taking the thermodynamic contributions into account. As the silica model we used a metoxy terminated dihydroxy substituted silanol Me O ½SiðOHÞ2 O n Me;
n¼4
which was proved to very well represent the powder surface in an earlier study [12]. Eight different complexes were constructed and the whole complexes were left fully flexible for optimization, so that it would adopt the most favourable conformation. A vibrational analysis was performed the four converged complexes. The basis set superposition error (BSSE) was estimated by counterpoise method [31] at the B3LYP/TZVP level using the MARI-J approximation. The BSSE varied between 2.3 and 5.8 kJ mol1 for the complexes between disulfiram and the silica surface, and this value is the most accurate one in discriminating the binding mode of disulfiram on silica. 3. Results and discussion 3.1. The minimum energy structures The two almost equally stable conformers of disulfiram (1 and 2) are presented in Fig. 1.
Fig. 1. Two conformers of disulfiram, 1 (left) and 2 (right). The distance is given in Ångström.
The electronic energy of structure 1 was only 0.3 kJ/mol lower than that of structure 2 but when comparing the Gibbs’ free energies (DG, 25 °C) structure 2 was found to be 0.8 kJ/mol more stable. The difference in the structures is the torsion angle of the SASbond. The dipole moments of structure 1 and 2 are 0.49 and 0.59 debye, respectively. The Mulliken and Löwdin charges and the charges based on an electrostatic potential fit (ESP) are presented in Table 1 for the 16 heavy atoms as the hydrogens are not found to be relevant in this study. The numbering of the atoms is illustrated in Fig. 2. The thioketone sulphur atoms (S7 and S11) are the ones that are most likely to interact with the hydroxyls of silica. In all three representations of the charge, the negative charge on these atoms is obvious. For the other two sulphur atoms (S8 and S9) the negative charge that makes the hydrogen bond in complex 2 (Section 3.3) possible is seen only in the ESP-fitting values. The bond orders and bond distances in disulfiram for the most relevant atom pairs in this study are presented in Table 2. They do not differ between the two conformers. All the other bonds than the presented ones are pure single bonds with bond orders in the range [0.93, 0.98]. The results show quite clearly that the CANACAS-system is almost conjugated and almost planar. The bonding electrons are moved towards the thioketone groups, affecting also the bond between the nitrogen atoms and the ethyl groups, lowering these bond orders below 0.85. Table 1 Mulliken and Löwdin charges and electrostatic potential fit (ESP) charges for the 16 heavy atoms Conformer
Atom C1 C2 N3 C4 C5 C6 S7 S8 S9 C10 S11 N12 C13 C14 C15 C16
Mulliken
Löwdin
ESP-fit
1
2
1
2
1
2
0.52 0.14 0.57 0.15 0.49 0.06 0.26 0.13 0.13 0.06 0.26 0.57 0.15 0.49 0.14 0.52
0.51 0.14 0.57 0.15 0.49 0.06 0.26 0.13 0.13 0.06 0.26 0.57 0.15 0.49 0.14 0.51
0.45 0.18 0.13 0.18 0.45 0.08 0.27 0.11 0.11 0.08 0.27 0.13 0.18 0.45 0.19 0.45
0.44 0.18 0.13 0.18 0.45 0.08 0.27 0.12 0.12 0.08 0.28 0.13 0.18 0.45 0.19 0.45
0.34 0.24 0.25 0.30 0.24 0.24 0.38 0.10 0.10 0.22 0.37 0.25 0.29 0.23 0.26 0.36
0.36 0.30 0.24 0.29 0.25 0.23 0.38 0.10 0.09 0.21 0.37 0.22 0.27 0.23 0.29 0.37
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
59
Fig. 2. Numbering of the heavy atoms in the disulfiram molecule.
Table 2 Bond orders and bond distances in disulfiram for the most relevant atom pairs Atoms S8 C6 S9 C6 C10 N3 C10 C2 N3 N12 N12
S9 S8 C10 S7 S11 C6 N12 N3 C4 C13 C15
Distance [Å]
Bond order
2.02 1.81 1.81 1.66 1.66 1.33 1.33 1.48 1.47 1.47 1.48
0.99 0.94 0.94 1.54 1.54 1.26 1.26 0.84 0.85 0.85 0.84
Fig. 3. The HOMO and LUMO of structure 1 of disulfiram.
However, the bond distances are quite typical for analogous molecules. Experimentally the SAS bond length in dimethyl disulfide H3CASASACH3 has been determined to be 2.04 Å and the CAS bond length 1.81 Å [32], which is identical with our calculated values for disulfiram. The CASASAC torsional angle was optimized to 84.8 degrees, which compares quite well with our value for disulfiram. In dithio-oxamide the C@S bond length has been tabulated as 1.66 Å [33], which is exactly the value we also obtained for disulfiram. The CAN bond lengths that we obtained of 1.47 and 1.48 Å are pure paraffinic bonds, but ditto of 1.33 Å are almost conjugated such as the CAN bond lengths in peptides. To ensure that we have really found the minimum energy structure an additional potential energy surface scan about the SAS torsion angle was performed using a step size of 30 degrees at DFT/ B3LYP/TZVP level. A torsional barrier of 28.8 kJ/mol was obtained at the angle of 180 degrees, where also a 12.8 kJ/mol less stable conformer than structure 1 converged. The torsion angle of this less stable conformer was 177 degrees, i.e., a difference of only 3 degrees. However, the difference in bond distances is about 5 pm for the CAS bonds and 8 pm for the SAS-bond. Because of the strong p–p-conjugation between the lone pairs of the sulphur atoms and the p-orbitals in the C@S-bonds, even a small difference in the torsion angle makes a great impact on the energy. When considering the shape of the molecular orbitals on the single bonded sulphur atoms, the minimum energy structure is quite logical, because with the torsion angle of 0 or 180 degrees the repelling orbitals would overlap the most. The HOMO and LUMO of structure 1 are displayed in Fig. 3. 3.2. Protonation One representative of each of the three groups of protonated structures are presented in Fig. 4. The protonation energies without thermodynamic contributions varied between 892 and 977 kJ/mol. The proton affinity of the most stable protonated structure was calculated to be
Fig. 4. One representative of each of the three groups of protonated structures.
957.6 kJ/mol. The proton affinity of a compound A is defined as the negative of the enthalpy change PA = DHprot for the following reaction A+H+ ? AH+.
60
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
In structure 1, the proton is attached to one of the thioketone sulphur atoms. One covalent bond is formed as well as a hydrogen bond, and additionally the aromaticity of the SACANAC-system is preserved on both sides of the molecule. In structure 2, the reactive p-orbitals of the planar nitrogen atom are protonated, and the hybridization of the nitrogen atom is changed. In this case, the protonation energies varied between 913 and 927 kJ/mol. The third and final group of optimized geometries was unexpected. In this case, the proton is attached to one of the single bonded sulphur atoms in the CASASAC-system and, consequently, the CAS-bond is broken. One of the thioketone sulphur atoms has become triple bonded, and the two molecular units are coupled by a hydrogen bond with the length of 2.86 Å as can be seen in Fig. 4. These two substances were studied further taking into account thermodynamic contributions. No negative frequencies were observed and the complex was 34.6 kJ/mol (DH) more stable than the units apart. 3.3. Complexation Out of the eight starting structures for the complexes between disulfiram and silica, only four did converge. The converged structure of the complexes as well as some distances (in Ångström) and the BSSE-corrected electronic complexation energies are presented in Fig. 5. In complex 1 and 2, the energies correspond to weak hydrogen bonds [34] but in complex 3 and 4 conventional hydrogen bonds are formed, which is also reflected in the energies of these two complexes. In complex 4 actually two new hydrogen bonds are formed,
but at the same time an internal hydrogen bond in silica is lost. The internal hydrogen bond distances vary between 1.8 and 2.2 Å. It needs to be pointed out that the outcome of this study is mainly the finding of adsorption sites between disulfiram and silica, as the disulfiram molecule does not contain any typical functional groups, and consequently the interactions with solid surfaces are not easily foreseen. It could be argued that more comparable results might be obtained by freezing the silica part. Considering amorphous silica, the solid surrounding might force the silica on the surface to keep its structure more or less unaltered on adsorption. In fully optimized silica internal hydrogen bonds may therefore be lost or formed contrary to the real life. Deciding how much energy is gained by forming new hydrogen bonds and how much is lost by breaking old internal ones is difficult. A vibrational analysis was also carried out for these complexes to obtain the thermodynamic contributions, but unfortunately the complexes had vibrations of such low frequency that the harmonic approximation could not describe the vibrations properly. In such cases large errors easily appear in the entropy term, which results in large errors in the differences in Gibb’s free energy (DG). For these complexes DG varied between 42.3 and 51.4 kJ/mol, but the differences in complexation enthalpy (DH) were 0.6 (1), 2.9 (2), 12.4 (3) and 11.5 kJ/mol (4), respectively, which corresponds quite well to the electronic complexation energies. 4. Conclusions The strong complexation between disulfiram and silica observed experimentally was verified with quantum chemical calculations. The conformational equilibria of the disulfiram molecule
Fig. 5. The complexes between disulfiram and silica. The distances are given in Ångström.
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
were studied and two conformers with a difference of 0.8 kJ/mol in Gibb’s free energy at 25 °C were observed. Protonation of the optimized structures was also studied. Three groups of protonated structures were obtained and the proton affinity of the most stable protonated structure was calculated to 957.6 kJ/mol. In the complexation study with geminal silanols four different kinds of complexes did converge. The BSSE-corrected electronic complexation energy of the most stable complex was calculated to 13.6 kJ/mol. Acknowledgements Computing resources provided by CSC, the Finnish IT Center for Science, are kindly acknowledged. References [1] D. Cen, R.I. Gonzalez, J.A. Buckmeier, R.S. Kahlon, N.B. Tohidian, F.L. Meyskens Jr., Mol. Cancer Ther. 1 (2002) 197–204. [2] S.-G. Shiah, Y.-R. Kao, F.Y.-H. Wu, C.-W. Wu, Mol. Pharmacol. 64 (2003) 1076–1084. [3] W. Wang, H.L. McLeod, J. Cassidy, Int. J. Cancer 104 (2004) 504–511. [4] M. Wickström, K. Danielsson, L. Rickardson, J. Gullbo, P. Nygren, A. Isaksson, R. Larsson, H. Lövborg, Biochem. Pharmacol. 73 (1) (2007) 25–33. [5] D. Chen, Q.C. Cui, H. Yang, Q.P. Dou, Cancer Res. 66 (21) (2006) 10425–10433. [6] Z.E. Sauna, S. Shukla, S.V. Ambudkar, Mol. Biosyst. 1 (2) (2005) 127–134. [7] J. Lu, M. Liong, J.I. Zink, F. Tamanoi, Small 3 (2007) 1341–1346. [8] J.S. Beck, J.C. Vartuli, W.J. Roth, M.E. Leonowicz, C.T. Kresge, K.D. Schmitt, C.TW. Chu, D.H. Olson, E.W. Sheppard, S.B. McCullen, J.B. Higgins, J.L. Schlenker, J. Am. Chem. Soc. 114 (1992) 10834–10843. [9] C.T. Kresge, M.E. Leonowicz, W.J. Roth, J.C. Vartuli, J.S. Beck, Nature 359 (1992) 710–712. [10] (a) M. Hartmann, Chem. Mater. 17 (2005) 4577–4593; (b) I.I. Slowing, B.G. Trewyn, S. Giri, V.S.-Y. Lin, Adv. Funct. Mater. 17 (2007) 1225–1236.
61
[11] J.M. Rosenholm, T. Sandberg, M. Hotokka, M. Lindén, in preparation. [12] B. Granqvist, T. Sandberg, M. Hotokka, J. Coll. Int. Sci. 310 (2007) 369–376. [13] M. Clark, R.D. Cramer III, N. van Opdenbosch, J. Comput. Chem. 10 (1989) 982. [14] SYBYL 7.3, Tripos Inc., 1699 South Hanley Road, St. Louis, Missouri, 63144, USA. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S.J. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347–1363. [16] Warren J. Hehre et al., Ab Initio Molecular Orbital Theory, John Wiley & Sons, New York, 1986. [17] W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [18] P.C. Hariharan, J.A. Pople, Theor. Chim. Acta 28 (1973) 213. [19] M.M. Francl, W.J. Pietro, W.J. Hehre, J.S. Binkley, D.J. DeFrees, J.A. Pople, M.S. Gordon, J. Chem. Phys. 77 (1982) 3654. [20] R. Ahlrichs, M. Bär, M. Häser, H. Horn, C. Kölmel, Chem. Phys. Lett. 162 (1989) 165. [21] M. Häser, R. Ahlrichs, J. Comput. Chem. 10 (1989) 104. [22] M. von Arnim, R. Ahlrichs, J. Comput. Chem. 19 (1998) 1746. [23] A.D. Becke, Phys. Rev. A 38 (1988) 3098–3100. [24] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785–789. [25] A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. [26] K. Eichkorn, O. Treutler, H. Öhm, M. Häser, R. Ahlrichs, Chem. Phys. Lett. 242 (1995) 652. [27] K. Eichkorn, F. Weigend, O. Treutler, R. Ahlrichs, Theor. Chem. Acc. 97 (1997) 119. [28] M. Sierka, A. Hogekamp, R. Ahlrichs, J. Chem. Phys. 118 (2003) 9136. [29] A. Schäfer, C. Huber, R. Ahlrichs, J. Chem. Phys. 100 (1994) 5829. [30] A. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502–16513. [31] C.J. Cramer, Essentials of Computational Chemistry – Theories and Models, Wiley, Chichester, 2002. [32] W.N. Hubbard, D.R. Douslin, C. Lee, J.P. Mccullaugh, D.W. Scott, S.S. Todd, J.F. Messerly, I.A. Hossenlopp, A. George, G. Waddington, J. Am. Chem. Soc. 80 (1958) 3547. [33] B. Long, P. Markey, P.J. Wheatley, Acta Cryst. 7 (1954) 140. [34] G.R. Desiraju, T. Steiner, The Weak Hydrogen Bond in Structural Chemistry and Biology, Oxford University Press, New York, 1999.
Contents lists available at ScienceDirect
Journal of Molecular Structure: THEOCHEM journal homepage: www.elsevier.com/locate/theochem
The molecular structure of disulfiram and its complexation with silica. A quantum chemical study Thomas Sandberg, Jessica Rosenholm, Matti Hotokka * Department of Physical Chemistry, Åbo Akademi University, Porthansgatan 3-5, 20500 Åbo (Turku), Finland
a r t i c l e
i n f o
Article history: Received 30 January 2008 Received in revised form 8 April 2008 Accepted 8 April 2008 Available online 16 April 2008 Keywords: Disulfiram Conformation Ab initio calculations DFT Proton affinity Silica complex
a b s t r a c t The complexation between disulfiram, more commonly known by trade names such as Antabus or Antabuse, and amorphous silica was studied by quantum chemical calculations due to extensive adsorption observed experimentally. The conformational behavior of the drug molecule was studied computationally at molecular mechanics (MM), Hartree-Fock (HF)/6-31G* and DFT/B3LYP/TZVP level. Two almost equally stable conformers were observed with the energy difference of only 0.8 kJ/mol when comparing the Gibbs’ free energies (DG) at 25 °C. The Mulliken and Löwdin charges and the charges based on an electrostatic potential fit (ESP) were studied at HF/6-31G* level. Protonation of the optimized structures and complexation with geminal silanols were studied at DFT/ B3LYP/TZVP level of theory in combination with the MARI-J approximation. Three groups of protonated structures were obtained and the proton affinity of the most stable protonated structure was calculated to 957.6 kJ/mol. In the complexation study four different kinds of complexes did converge. The BSSE-corrected electronic complexation energy of the most stable complex was calculated to 13.6 kJ/mol. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Disulfiram, full IUPAC name 1-(diethylthiocarbamoyldisulfanyl)-N,N-diethyl-methanethioamide, more commonly known as the alcohol deterrent by trade names such as Antabus or Antabuse, is a widely applied aldehyde dehydrogenase inhibitor to support treatment of chronic alcoholism. Disulfiram has also been shown to induce apoptosis in a number of tumor cell lines [1–5] and in the reversal of fungal drug-resistance, as recently reviewed [6]. As many other anticancer agents, disulfiram is practically insoluble in water. Whereas low water solubility of drugs impedes the ability of drugs to be administered through the intravenous route, novel delivery systems for hydrophobic therapeutic anticancer drugs has received significant attention [7]. Ceramic carrier materials are attracting increased attention due to easy preparation routes via the sol–gel process which allows them to be prepared with desired shape, size and porosity. The mostly employed sol–gel derived material in such applications is SiO2, silica, which is nontoxic, biocompatible and bioerodible and moreover an essential component of cells throughout the human body. These properties provide sol–gel processed porous silica materials considerable potential as carriers for drug delivery. Due to their great adsorption properties, mesoporous silicas synthesized using supramolecular surfactant aggregates as structure directing agents [8,9] have at* Corresponding author. E-mail address: mhotokka@abo.fi (M. Hotokka). 0166-1280/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.theochem.2008.04.007
tracted lots of recent interest in applications related to immobilization and release of active substances, especially drugs [10]. Preliminary adsorption tests carried out in our laboratory showed that as much as C = 0.93 lmol m2 could be adsorbed from cyclohexane to a model mesoporous silica with 3 nm pores, corresponding to a loading degree of about 30 wt% [11]. Cyclohexane provides a non-aquatic and non-polar environment, which would best correspond to theoretical calculations. While most drugs are weak acids or bases, the disulfiram molecule does not contain any typical functional groups, and consequently the interactions with solid surfaces are not easily foreseen. Hence, we have used quantum chemical calculations in order to investigate the molecular reasons for the strong adsorption of disulfiram to silica. To the best of our knowledge, neither theoretical nor experimental structural analysis has been published to date. It is important to determine the minimum energy structure in order to understand the behavior of disulfiram, e.g., in adsorption on silica surfaces. An analogous study of organic aromatic probes on silica has recently been published [12]. 2. Computational methods 2.1. The minimum energy structures In order to find the minimum energy structure of disulfiram, an initial potential energy surface scan about the three middle bonds in the NACASASACAN-skeleton was performed using a step size
58
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
of 30 degrees at molecular mechanics (MM) level using the Tripos Force Field [13] as implemented in SYBYL 7.3 [14]. The 11 energetically most favourable structures (MM energies below 42.3 kJ/mol) were optimized with the program GAMESS version 22 Feb 2006 [15] by using Hartree-Fock (HF) theory [16] with the basis set 631G* [17–19]. Two almost equally stable conformers of disulfiram, structures 1 and 2, were found. Also Mulliken and Löwdin charges, electrostatic potential fit (ESP) charges and bond orders were calculated with GAMESS at HF/6-31G* level of theory. The conformers were then reoptimized for vibrational analysis using the TURBOMOLE program package version 5.9 [20–22] and density functional theory (DFT) with the B3LYP hybrid exchangecorrelation functional [23–25] in combination with the multipole accelerated resolution of identity (MARI-J) approximation [26– 28] and the TZVP basis set [29] for all atoms, as implemented in the TURBOMOLE program package. Vibrational analysis was carried out for the B3LYP/TZVP optimized structures in order to prove that all optimized structures are real minima on the potential energy surface and to obtain the thermodynamic contributions. Frequencies were scaled with a factor of 0.9614 [30] when calculating the thermodynamic contributions at 25 °C. 2.2. Protonation Ten different protonated structures were constructed by placing the proton close to the nitrogen or sulfur atoms, which are expected to be the reaction sites. The protonated structures were optimized using TURBOMOLE 5.9 at DFT/B3-LYP/TZVP level of theory in combination with the MARI-J approximation. Three groups of protonated structures were obtained and the most stable one was studied in detail at the B3LYP/TZVP level taking into account thermodynamic contributions. As the third type of protonated structures resulted in a splitting into two molecular units, a vibrational analysis was performed also on that system. 2.3. Complexation As the aim was to study disulfiram interaction with amorphous silica surfaces, also the complexation between the energetically most favourable conformation of disulfiram and a silica model was studied in detail at the B3LYP/TZVP level taking the thermodynamic contributions into account. As the silica model we used a metoxy terminated dihydroxy substituted silanol Me O ½SiðOHÞ2 O n Me;
n¼4
which was proved to very well represent the powder surface in an earlier study [12]. Eight different complexes were constructed and the whole complexes were left fully flexible for optimization, so that it would adopt the most favourable conformation. A vibrational analysis was performed the four converged complexes. The basis set superposition error (BSSE) was estimated by counterpoise method [31] at the B3LYP/TZVP level using the MARI-J approximation. The BSSE varied between 2.3 and 5.8 kJ mol1 for the complexes between disulfiram and the silica surface, and this value is the most accurate one in discriminating the binding mode of disulfiram on silica. 3. Results and discussion 3.1. The minimum energy structures The two almost equally stable conformers of disulfiram (1 and 2) are presented in Fig. 1.
Fig. 1. Two conformers of disulfiram, 1 (left) and 2 (right). The distance is given in Ångström.
The electronic energy of structure 1 was only 0.3 kJ/mol lower than that of structure 2 but when comparing the Gibbs’ free energies (DG, 25 °C) structure 2 was found to be 0.8 kJ/mol more stable. The difference in the structures is the torsion angle of the SASbond. The dipole moments of structure 1 and 2 are 0.49 and 0.59 debye, respectively. The Mulliken and Löwdin charges and the charges based on an electrostatic potential fit (ESP) are presented in Table 1 for the 16 heavy atoms as the hydrogens are not found to be relevant in this study. The numbering of the atoms is illustrated in Fig. 2. The thioketone sulphur atoms (S7 and S11) are the ones that are most likely to interact with the hydroxyls of silica. In all three representations of the charge, the negative charge on these atoms is obvious. For the other two sulphur atoms (S8 and S9) the negative charge that makes the hydrogen bond in complex 2 (Section 3.3) possible is seen only in the ESP-fitting values. The bond orders and bond distances in disulfiram for the most relevant atom pairs in this study are presented in Table 2. They do not differ between the two conformers. All the other bonds than the presented ones are pure single bonds with bond orders in the range [0.93, 0.98]. The results show quite clearly that the CANACAS-system is almost conjugated and almost planar. The bonding electrons are moved towards the thioketone groups, affecting also the bond between the nitrogen atoms and the ethyl groups, lowering these bond orders below 0.85. Table 1 Mulliken and Löwdin charges and electrostatic potential fit (ESP) charges for the 16 heavy atoms Conformer
Atom C1 C2 N3 C4 C5 C6 S7 S8 S9 C10 S11 N12 C13 C14 C15 C16
Mulliken
Löwdin
ESP-fit
1
2
1
2
1
2
0.52 0.14 0.57 0.15 0.49 0.06 0.26 0.13 0.13 0.06 0.26 0.57 0.15 0.49 0.14 0.52
0.51 0.14 0.57 0.15 0.49 0.06 0.26 0.13 0.13 0.06 0.26 0.57 0.15 0.49 0.14 0.51
0.45 0.18 0.13 0.18 0.45 0.08 0.27 0.11 0.11 0.08 0.27 0.13 0.18 0.45 0.19 0.45
0.44 0.18 0.13 0.18 0.45 0.08 0.27 0.12 0.12 0.08 0.28 0.13 0.18 0.45 0.19 0.45
0.34 0.24 0.25 0.30 0.24 0.24 0.38 0.10 0.10 0.22 0.37 0.25 0.29 0.23 0.26 0.36
0.36 0.30 0.24 0.29 0.25 0.23 0.38 0.10 0.09 0.21 0.37 0.22 0.27 0.23 0.29 0.37
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
59
Fig. 2. Numbering of the heavy atoms in the disulfiram molecule.
Table 2 Bond orders and bond distances in disulfiram for the most relevant atom pairs Atoms S8 C6 S9 C6 C10 N3 C10 C2 N3 N12 N12
S9 S8 C10 S7 S11 C6 N12 N3 C4 C13 C15
Distance [Å]
Bond order
2.02 1.81 1.81 1.66 1.66 1.33 1.33 1.48 1.47 1.47 1.48
0.99 0.94 0.94 1.54 1.54 1.26 1.26 0.84 0.85 0.85 0.84
Fig. 3. The HOMO and LUMO of structure 1 of disulfiram.
However, the bond distances are quite typical for analogous molecules. Experimentally the SAS bond length in dimethyl disulfide H3CASASACH3 has been determined to be 2.04 Å and the CAS bond length 1.81 Å [32], which is identical with our calculated values for disulfiram. The CASASAC torsional angle was optimized to 84.8 degrees, which compares quite well with our value for disulfiram. In dithio-oxamide the C@S bond length has been tabulated as 1.66 Å [33], which is exactly the value we also obtained for disulfiram. The CAN bond lengths that we obtained of 1.47 and 1.48 Å are pure paraffinic bonds, but ditto of 1.33 Å are almost conjugated such as the CAN bond lengths in peptides. To ensure that we have really found the minimum energy structure an additional potential energy surface scan about the SAS torsion angle was performed using a step size of 30 degrees at DFT/ B3LYP/TZVP level. A torsional barrier of 28.8 kJ/mol was obtained at the angle of 180 degrees, where also a 12.8 kJ/mol less stable conformer than structure 1 converged. The torsion angle of this less stable conformer was 177 degrees, i.e., a difference of only 3 degrees. However, the difference in bond distances is about 5 pm for the CAS bonds and 8 pm for the SAS-bond. Because of the strong p–p-conjugation between the lone pairs of the sulphur atoms and the p-orbitals in the C@S-bonds, even a small difference in the torsion angle makes a great impact on the energy. When considering the shape of the molecular orbitals on the single bonded sulphur atoms, the minimum energy structure is quite logical, because with the torsion angle of 0 or 180 degrees the repelling orbitals would overlap the most. The HOMO and LUMO of structure 1 are displayed in Fig. 3. 3.2. Protonation One representative of each of the three groups of protonated structures are presented in Fig. 4. The protonation energies without thermodynamic contributions varied between 892 and 977 kJ/mol. The proton affinity of the most stable protonated structure was calculated to be
Fig. 4. One representative of each of the three groups of protonated structures.
957.6 kJ/mol. The proton affinity of a compound A is defined as the negative of the enthalpy change PA = DHprot for the following reaction A+H+ ? AH+.
60
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
In structure 1, the proton is attached to one of the thioketone sulphur atoms. One covalent bond is formed as well as a hydrogen bond, and additionally the aromaticity of the SACANAC-system is preserved on both sides of the molecule. In structure 2, the reactive p-orbitals of the planar nitrogen atom are protonated, and the hybridization of the nitrogen atom is changed. In this case, the protonation energies varied between 913 and 927 kJ/mol. The third and final group of optimized geometries was unexpected. In this case, the proton is attached to one of the single bonded sulphur atoms in the CASASAC-system and, consequently, the CAS-bond is broken. One of the thioketone sulphur atoms has become triple bonded, and the two molecular units are coupled by a hydrogen bond with the length of 2.86 Å as can be seen in Fig. 4. These two substances were studied further taking into account thermodynamic contributions. No negative frequencies were observed and the complex was 34.6 kJ/mol (DH) more stable than the units apart. 3.3. Complexation Out of the eight starting structures for the complexes between disulfiram and silica, only four did converge. The converged structure of the complexes as well as some distances (in Ångström) and the BSSE-corrected electronic complexation energies are presented in Fig. 5. In complex 1 and 2, the energies correspond to weak hydrogen bonds [34] but in complex 3 and 4 conventional hydrogen bonds are formed, which is also reflected in the energies of these two complexes. In complex 4 actually two new hydrogen bonds are formed,
but at the same time an internal hydrogen bond in silica is lost. The internal hydrogen bond distances vary between 1.8 and 2.2 Å. It needs to be pointed out that the outcome of this study is mainly the finding of adsorption sites between disulfiram and silica, as the disulfiram molecule does not contain any typical functional groups, and consequently the interactions with solid surfaces are not easily foreseen. It could be argued that more comparable results might be obtained by freezing the silica part. Considering amorphous silica, the solid surrounding might force the silica on the surface to keep its structure more or less unaltered on adsorption. In fully optimized silica internal hydrogen bonds may therefore be lost or formed contrary to the real life. Deciding how much energy is gained by forming new hydrogen bonds and how much is lost by breaking old internal ones is difficult. A vibrational analysis was also carried out for these complexes to obtain the thermodynamic contributions, but unfortunately the complexes had vibrations of such low frequency that the harmonic approximation could not describe the vibrations properly. In such cases large errors easily appear in the entropy term, which results in large errors in the differences in Gibb’s free energy (DG). For these complexes DG varied between 42.3 and 51.4 kJ/mol, but the differences in complexation enthalpy (DH) were 0.6 (1), 2.9 (2), 12.4 (3) and 11.5 kJ/mol (4), respectively, which corresponds quite well to the electronic complexation energies. 4. Conclusions The strong complexation between disulfiram and silica observed experimentally was verified with quantum chemical calculations. The conformational equilibria of the disulfiram molecule
Fig. 5. The complexes between disulfiram and silica. The distances are given in Ångström.
T. Sandberg et al. / Journal of Molecular Structure: THEOCHEM 861 (2008) 57–61
were studied and two conformers with a difference of 0.8 kJ/mol in Gibb’s free energy at 25 °C were observed. Protonation of the optimized structures was also studied. Three groups of protonated structures were obtained and the proton affinity of the most stable protonated structure was calculated to 957.6 kJ/mol. In the complexation study with geminal silanols four different kinds of complexes did converge. The BSSE-corrected electronic complexation energy of the most stable complex was calculated to 13.6 kJ/mol. Acknowledgements Computing resources provided by CSC, the Finnish IT Center for Science, are kindly acknowledged. References [1] D. Cen, R.I. Gonzalez, J.A. Buckmeier, R.S. Kahlon, N.B. Tohidian, F.L. Meyskens Jr., Mol. Cancer Ther. 1 (2002) 197–204. [2] S.-G. Shiah, Y.-R. Kao, F.Y.-H. Wu, C.-W. Wu, Mol. Pharmacol. 64 (2003) 1076–1084. [3] W. Wang, H.L. McLeod, J. Cassidy, Int. J. Cancer 104 (2004) 504–511. [4] M. Wickström, K. Danielsson, L. Rickardson, J. Gullbo, P. Nygren, A. Isaksson, R. Larsson, H. Lövborg, Biochem. Pharmacol. 73 (1) (2007) 25–33. [5] D. Chen, Q.C. Cui, H. Yang, Q.P. Dou, Cancer Res. 66 (21) (2006) 10425–10433. [6] Z.E. Sauna, S. Shukla, S.V. Ambudkar, Mol. Biosyst. 1 (2) (2005) 127–134. [7] J. Lu, M. Liong, J.I. Zink, F. Tamanoi, Small 3 (2007) 1341–1346. [8] J.S. Beck, J.C. Vartuli, W.J. Roth, M.E. Leonowicz, C.T. Kresge, K.D. Schmitt, C.TW. Chu, D.H. Olson, E.W. Sheppard, S.B. McCullen, J.B. Higgins, J.L. Schlenker, J. Am. Chem. Soc. 114 (1992) 10834–10843. [9] C.T. Kresge, M.E. Leonowicz, W.J. Roth, J.C. Vartuli, J.S. Beck, Nature 359 (1992) 710–712. [10] (a) M. Hartmann, Chem. Mater. 17 (2005) 4577–4593; (b) I.I. Slowing, B.G. Trewyn, S. Giri, V.S.-Y. Lin, Adv. Funct. Mater. 17 (2007) 1225–1236.
61
[11] J.M. Rosenholm, T. Sandberg, M. Hotokka, M. Lindén, in preparation. [12] B. Granqvist, T. Sandberg, M. Hotokka, J. Coll. Int. Sci. 310 (2007) 369–376. [13] M. Clark, R.D. Cramer III, N. van Opdenbosch, J. Comput. Chem. 10 (1989) 982. [14] SYBYL 7.3, Tripos Inc., 1699 South Hanley Road, St. Louis, Missouri, 63144, USA. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S.J. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347–1363. [16] Warren J. Hehre et al., Ab Initio Molecular Orbital Theory, John Wiley & Sons, New York, 1986. [17] W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [18] P.C. Hariharan, J.A. Pople, Theor. Chim. Acta 28 (1973) 213. [19] M.M. Francl, W.J. Pietro, W.J. Hehre, J.S. Binkley, D.J. DeFrees, J.A. Pople, M.S. Gordon, J. Chem. Phys. 77 (1982) 3654. [20] R. Ahlrichs, M. Bär, M. Häser, H. Horn, C. Kölmel, Chem. Phys. Lett. 162 (1989) 165. [21] M. Häser, R. Ahlrichs, J. Comput. Chem. 10 (1989) 104. [22] M. von Arnim, R. Ahlrichs, J. Comput. Chem. 19 (1998) 1746. [23] A.D. Becke, Phys. Rev. A 38 (1988) 3098–3100. [24] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785–789. [25] A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. [26] K. Eichkorn, O. Treutler, H. Öhm, M. Häser, R. Ahlrichs, Chem. Phys. Lett. 242 (1995) 652. [27] K. Eichkorn, F. Weigend, O. Treutler, R. Ahlrichs, Theor. Chem. Acc. 97 (1997) 119. [28] M. Sierka, A. Hogekamp, R. Ahlrichs, J. Chem. Phys. 118 (2003) 9136. [29] A. Schäfer, C. Huber, R. Ahlrichs, J. Chem. Phys. 100 (1994) 5829. [30] A. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502–16513. [31] C.J. Cramer, Essentials of Computational Chemistry – Theories and Models, Wiley, Chichester, 2002. [32] W.N. Hubbard, D.R. Douslin, C. Lee, J.P. Mccullaugh, D.W. Scott, S.S. Todd, J.F. Messerly, I.A. Hossenlopp, A. George, G. Waddington, J. Am. Chem. Soc. 80 (1958) 3547. [33] B. Long, P. Markey, P.J. Wheatley, Acta Cryst. 7 (1954) 140. [34] G.R. Desiraju, T. Steiner, The Weak Hydrogen Bond in Structural Chemistry and Biology, Oxford University Press, New York, 1999.