Charmm force field parameterization of rosiglitazone

CHARMM Force Field Parameterizationof Rosiglitazone ANDERS HANSSON, PAULO C. T. SOUZA, RODRIGO L. SILVEIRA,LEANDRO MARTI´NEZ, MUNIR S. SKAFInstitute of Chemistry, State University of Campinas—UNICAMP, C.P. 6154,Campinas, SP 13084-862, Brazil Received 10 December 2009; accepted 1 February 2010Published online 10 June 2010 in Wiley Online Library (wileyonlinelibrary.com).
DOI 10.1002/qua.22638 ABSTRACT: We develop a CHARMM-based interaction potential for rosiglitazone, awell-known selective ligand to the c isoform of the peroxisome proliferator-activatedreceptor (PPARc) and widely marketed antidiabetic drug of the thiazolidinedione (TZD)class. We derive partial atomic charges and dihedral torsion potentials for sevenrotations in the molecule, for which there are no analogs available in CHARMM. Thepotential model is validated by performing a series of molecular dynamics simulationsof rosiglitazone in neat water and of a fully solvated rosiglitazone-PPARc complex. Thestructural and dynamical behavior of the complex is analyzed in comparison withavailable experimental data. The potential parameters derived here are readilytransferable to a variety of pharmaceutically important TZD compounds. V Periodicals, Inc. Int J Quantum Chem 111: 1346–1354, 2011 Key words: rosiglitazone; TZD; CHARMM parameterization; nuclear receptors;PPARc; molecular dynamics differentiation and proliferation [2]. PPARc is also involved in inflammatory and immune responses and is highly expressed in various types of cancer he nuclear receptor superfamily comprises a and associated with type II diabetes [3]. PPARc group of roughly 48 proteins responsible for forms heterodimers with the retinoid X receptor regulation of gene transcriptional activity by (RXR), which is the necessary partner for DNA means of hormone binding [1]. The c isoform of binding and transcription. RXR, activated by its natural ligand 9-cis retinoic acid, serves as hetero- (PPARc) is a key receptor in the regulation of cell dimeric partner for many nuclear receptors [1].
The more promiscuous PPARc accommodates Correspondence to: M. S. Skaf; e-mail: skaf@iqm.unicamp.br various types of ligands, mostly agonists.
Contract grant sponsors: Brazilian Agencies Fapesp, CNPq.
Additional Supporting Information may be found in the ligands and to determine the crystal structure of International Journal of Quantum Chemistry, Vol 111, 1346–1354 (2011) CHARMM FORCE FIELD PARAMETERIZATION OF ROSIGLITAZONE FIGURE 1. Chemical structure of rosiglitazone. The arc-shaped arrows indicate the torsion (T) of the eight bonds forwhich full revolution is possible.
receptors are quite advanced. However, the dynamic interactions of receptor, DNA response elements, and ligands remain poorly understood.
dynamics (MD), admit studies of association and dissociation processes [4–10], as well as the molecular mechanisms involved in receptor acti- vation which give insight in ligand selectivity [11, 12]. The force fields primarily derived for molecu- lar mechanics of macromolecules, such proteins, saccharides, lipids, and nucleic acids, do not pro- vide parameters for specific chemical compounds,such as the thiazolidinediones (TZD), a well- where the sums extend to all bond stretchings, known group of PPARc agonists, of which rosigli- Urey-Bradley terms, angle bendings, dihedral tor- tazone among others have been clinically studied sions, improper dihedrals, and nonbonded van and produced by the pharmaceutical industry der Waals and Coulombic interactions. There are potential parameters developed for a great variety The TZDs appear in two enantiomeric forms of biomolecules [17–20]. The transferability of pa- (R)-(þ) and (S)-(À), due to the stereogenic center rameters between molecules is the basic principle at atom C5 of the thiazolidine ring (Fig. 1). Higher of these force fields, so parameters from similar antidiabetic activity have been predicted for the molecules should be used whenever possible.
(S)-(À) enantiomer of rosiglitazone [14], which However, for markedly flexible molecules, the also is the observed form in the available crystal- sampling over different conformational states lized PPARc structures in the Protein Data Bank dependent critically on torsion (dihedral angle) (PDB) [15]. Therefore, a complete set of parame- potentials, so obtaining accurate torsional parame- ters for the two enantiomers of rosiglitazone are ters derived from quantum mechanical potential derived here for future studies, especially the energy surface (PES) scans of a given particular interactions the (S)-(À) enantiomer with nuclear receptors. The particular parameterization chosen In this work, we propose a complete set of for these simulations is based on the CHARMM CHARMM compatible force field parameters for force field for biomolecular systems [16]. The rosiglitazone. Our main goal is to generate a force potential energy function is expressed by Eq. (1) field for this molecule which is suitable for INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY interactions between the drug and its known tar- For the self-consistent field (SCF) procedure, get receptor and other proteins. Partial point the default combination of the EDIIS and the CDIIS algorithms, with no damping or Fermi atom types assigned, and similar bond lengths, broadening, is used. The convergence criteria (in atomic units) are 1.0 Â 10À6 for the energy, 1.0 Â improper dihedral parameters for atom combina- 10À6 for the maximum of density matrix, and tions not included in CHARMM are adopted.
1.0 Â 10À8 for the root mean square (RMS) of den- Finally, novel torsional parameters for seven dihe- sity matrix. The geometry is relaxed with the dral rotations are derived based on CHARMM implemented version of the Berny geometry opti- procedures. The newly parameterized model for mization algorithm and the convergence criteria rosiglitazone is tested by performing molecular (in atomic units) are 1.5 Â 10À5 for maximum dynamics (MD) simulations of this molecule in force, 1.0 Â 10À5 for RMS force, 6.0 Â 10À5 for aqueous solutions and bound to the ligand bind- maximum displacement, and 4.0 Â 10À5 for RMS ing domain (LBD) of PPARc. We examine the displacement. All four criteria were simultane- simulated structure and dynamics of the rosiglta- zone-LBD complex in comparison with available Stationary points for RHF/6-31G(d) as well as experimental data and find that the proposed with the more accurate 6-311G(d,p) basis set model is very-well behaved as far as these prop- using both RHF and density functional theory DFT/B3LYP methods are obtained to evaluate the Because all TZDs share the same molecular minimum of the first method. The partial point backbone, the torsion potentials derived here are charge calculation is, however, based on the readily transferred to a series of pharmaceutically RHF/6-31G(d) geometry for consistency with important compounds used in the treatment of CHARMM. Net atomic charges are derived to fit type II diabetes, such as pioglitazone, troglitazone, the RHF/6-31G(d) electrostatic potential of the rivoglitazone, and ciglitazone, to name few. There are two specific features in the molecular structure selected according to the employed Merz-Singh- of the TZDs for which adequate CHARMM param- Kollman approach [26, 27] in 10 layers and eters were so far unavailable. One of them is the 17 grid point per unit area resulting in 95,155 five-member heterocyclic aromatic ring containing points (the default is four layers and one point S1 and N3 atoms for which there are neither per unit area, giving 1,850 points for rosiglita- CHARMM partial charges nor parameters for the zone.) The charges are also constrained to repro- T1 torsion (Fig. 1). On the other end of the mole- duce the molecular dipole moment. The atoms cule, there is an unusual pyridine ring containing a are classified in atom types of the CHARMM22 all-atom force field for proteins [17, 18] (release CHARMM, which greatly affects the potential bar- c35b2), based on the derived charges and the local chemical environment. Unavailable parametersare adopted primarily from similar groups of theall-atom force field for proteins CHARMM22(release CHARMM32 for esters [28, 29] (of the samerelease).
The force field parameters are developed at the The flexibility of the ligand is an important fac- restricted Hartree-Fock (RHF) ab initio level of tor for their binding modes in the active site of theory with the 6-31G(d) basis set in consistency the LBD and influence the ligand dissociation with the CHARMM parameterization of the c35b2 mechanism [4–7, 10]. The structure of rosiglita- release. All ab initio calculations are performed zone admits full revolution around eight bonds, with the Gaussian03 package revision E.01 [22], as depicted in Fig. 1. The dihedral angle parame- whereas classical force field potential energies are ters required for describing the torsional barriers computed with the NAMD package [23], which is are not included in the current releases of also used for the evaluative MD simulations. The CHARMM. Therefore, all torsional rotations, acid dissociation constant (pKa) of rosiglitazone is excluding the rotation of the CH3 group, are para- meterized. To avoid interference with the existing INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY CHARMM FORCE FIELD PARAMETERIZATION OF ROSIGLITAZONE dihedral parameters and to obtain unique set of dihedral angle parameters for the seven rotating hydration shell around the protein is at least 15 A bonds, a few atom type aliases are introduced.
thick and the entire system is enclosed in ortho- To calculate the rotational energy profiles, the rhombic box with periodic boundary conditions.
quadruples C4-C5-C6-C7, C5-C6-C7-C8, C9-C10- Langevin dynamics is employed to simulate the isothermal-isobaric ensemble at 300 K and 1.0 C14-C15-N16-C17, and C15-N16-C17-N18 define atm. The velocity Verlet algorithm is used for the dihedral torsions T1–T7, respectively. The time integration with a time step of 2.0 fs. Full energy profiles are obtained from rotational scans Coulomb forces are computed with the particle at 20 degree steps with a RHF/6-31G(d) geometry mesh Ewald algorithm [32], whereas van der relaxation where only the regarded rotational angle is fixed. These energy profiles correspond switching cutoff. The systems are initially relaxed to the total quantum mechanical potential energy, with the default conjugate gradient and line which in some cases differ substantially from the search algorithm (CG) of NAMD2.7. We used the dihedral torsion potential energy, despite the relaxed geometry, mainly due to the constant parameters for the protein and counterions.
point charge model. The missing classical dihe- The following protocols are used for equilibrat- ing the two systems. Rosiglitazone in water: (1) 300 steps of CG, keeping all atoms of rosiglita- zone fixed; (2) 600 steps of CG, without any restraint; (3) 1,000 ps of MD without any restraint.
The ligand–receptor complex in solvent: (1) 2,000 is obtained for each dihedral torsion by least- steps of CG, with all atoms of the protein fixed; square fitting this expression to the difference (2) 200 ps of MD with all atoms of the protein between the quantum and the molecular mechan- fixed; (3) 500 steps of CG, with the Ca atoms ics energy profiles. Since dihedral angle parame- fixed; (4) 200 ps of MD with the Ca fixed; (5) 400 ters are unknown for seven rotational barriers, the ps of MD without any restraint. Finally, 6 ns MD fitting procedure was conveniently done repeat- simulations were performed starting from the edly, in a self-consistent manner, until conver- equilibrated structures. The trajectories of first gence of all torsion potentials. The potential 1 ns of these simulations were discarded. This energy of the molecular mechanics force field is protocol was repeated 20 times for the ligand– calculated with the NAMD2.7 package [23] via receptor complex; thus, 20 independent 5 ns MD simulations were obtained for data analysis.
terms with multiplicity n 6 are considered and,as far as possible, phase angles d are multiples of180, consistent with CHARMM. The obtainedparameters for the dihedral angles defined by the aforementioned quadruples describe solely thetotal dihedral angle potential for each bond rota- First the neutral state of (S)-(À)-rosiglitazone in tion. Hence, the other force constants of the neat water is considered. The structure is relaxed involved dihedral angles are kept as zero.
(in atomic units) to a maximum force of 1.0 Â Finally, MD simulations are performed to eval- 10À6, RMS force of 1.0 Â 10À7, maximum dis- uate the developed parameters. Both rosiglitazone placement of 1.9 Â 10À5, and RMS displacement solvated in a box containing 1,828 water mole- of 0.3 Â 10À5. This geometry has total energy of cules at ambient conditions and rosiglitazone À1478.1721 atomic units, which has converged to bound to the LBD of PPARc (fully hydrated) are 0.37 Â 10À8, and its dipole moment is 3.014 D. On studied. The structure of the ligand–receptor the basis of this configuration, we estimated the complex, based on the ligand binding domain acid dissociation constants for rosiglitazone to of the PDB structure 1FM6, chain D [31], is fully 6.9 (base part) and 6.5 (acid part), as described solvated by 18,000 water molecules, 28 Naþ above. The values are consistent with the litera- and 23 ClÀ ions, reaching a concentration of ture [34]. At physiological pH around 7.4 [35], the $0.15 mol/L to obtain a neutral system. Struc- unprotonated, neutral form of rosiglitazone (sin- tural (crystallographic) water molecules within a glet state) predominates. Therefore, parameters INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY TABLE IAtom names, CHARMM atom types, and partial atomic charges derived in this work.
a CHARMM22 all-atom force field for proteins.
b À0.16 is used for the force field.
are derived for neutral rosiglitazone. Table I lists PES (black dots), the molecular mechanics PES the calculated Merz-Singh-Kollman partial atomic excluding the dihedral potential for the T3 torsion charges for the ground state geometry. Despite (crosses), the quantum and classical PES differ- the differentiated hydrogen charges of the CH3 ence (plus signs) and the adjusted dihedral poten- groups, the average 0.09 is adopted for all methyl tial (solid line). For comparison, the full molecular hydrogens, yielding À0.16 for the methyl carbons.
mechanics PES, now including the parameterized The ground state geometry at RHF and B3LYP torsional potential, is shown by empty circles. The level using the 6-311G(d,p) basis set confirm agreement between the quantum and the full clas- above calculated molecular conformation.
sical force filed PESs is excellent. In the remainder The atoms are classified in atom types of the panels, only the quantum and classical PESs dif- CHARMM22 all-atom force field for proteins ferences (plus signs) and the adjusted potentials (release c35b2) based on the derived charges and (lines) are shown. As mentioned in ‘‘Parameter- the local chemical environment (Table I). Bond, ization and Computational Details,’’ we have angle, and dihedral parameters (excluding T1–T7 applied the CHARMM restriction of using the dihedral angles discussed below), as well as first six terms in the fitting function [cf., Eq. (2)], Urey-Bradley distances and Lennard-Jones energy which limits quality of the fit. Moreover, the and distance parameters, which could be obtained adjustments of the torsion parameters have been by group analogy from the CHARMM force field, carried out in a self-consistent way due to the are provided as Supporting Information.
mutual dependence between the distinct rota- The calculated dihedral angle PES for the tor- tions. The complete set of fitted parameters thus sions T1–T7 and the corresponding fitting func- tions [Eq. (2)] are depicted in Fig. 2. The top panel The first test for the force field is performed (T3) also shows the ab initio quantum mechanical by examining the behavior of rosiglitazone in an INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY CHARMM FORCE FIELD PARAMETERIZATION OF ROSIGLITAZONE FIGURE 2. Torsion energies of the bonds C5-C6 (T1), C6-C7 (T2), C10-O13 (T3), O13-C14 (T4), C14-C15 (T5), C15-N16 (T6), and N16-C17 (T7). The black dots represent the ab initio PES, whereas the crosses (x) indicate the forcefield PES without the considered dihedral angle. The difference between them is shown by plus signs (þ) and theadjusted potential energy of the dihedral angle is depicted by a solid line. For comparison, the circles show the totalforce field, including the parameterized dihedral angle.
aqueous environment. The average geometry of gives a root mean square deviation (RMSD) of the molecule in water closely replicates the ab initio relaxed geometry after the initial relaxation maintained during the simulation, so no addi- of the system. Structural aligning of the geome- tional improper dihedral angles are therefore tries, where all atoms are taken into account, INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY TABLE IIDihedral parameters for seven torsional rotations.
FIGURE 3. (a) Average structure of the PPARc- rosiglitazone complex from simulations (light gray) superimposed to the crystallographic structure (dark gray). (b) Average structures of rosiglitazone for individual 5 ns simulations (light gray) and the crystalstructure (black), with the LBDs structurally aligned.
The available experimental data of rosiglita- the behavior of the ligand and its target nuclear zone regarding its dynamical behavior is limited receptor protein PPARc in the ligand–receptor to the crystallographic Debye-Waller or tempera- complex. In all currently reported structures in ture B-factors. A comparison between the crystal- PDB containing rosiglitazone (PDB IDs: 1FM6, lographic and simulated B-factors is presented in Fig. 4, showing very good agreement, apart from bounded to the LBD pocket of PPARc. One of the an overall multiplicative factor, which is common highest resolution structures of the ligand-LBD in this type of comparison because the crystallo- complex is the 1FM6 structure, with a resolution graphic B-factors depend on the structure resolu- ˚ . Therefore, our evaluative simulations of tion. The mobility of the different structural rosiglitazone interacting with the PPARc LBD elements of the LBD derived from diffraction data are well-reproduced by the motions of rosiglita- The average structure of the rosiglitazone-LBD zone obtained from the MD simulations with the complex over 20 independent 5 ns simulations [Fig. 3(a)]. The RMSD between the average simu-lated and crystal structures of rosiglitazone itself isonly $0.35 A ˚ . The average structure of rosiglitazone obtained from each 5 ns simulation and the crystalstructure are visualized in Fig. 3(b), after alignmentof the LBDs, indicating that the PPARc LBDremains very stable and structurally well-correlatedwith the crystallographic structure in the presenceof the parameterized ligand. In addition, the aver-age conformation and position of rosiglitazoneinside the ligand binding pocket (i.e., the ligandbinding mode), preserve the crystal conformationand the crystalline ligand-LBD contacts. Theseresults indicate that the ligand–protein and ligand FIGURE 4. Temperature B-factors obtained from intramolecular interactions are well accounted for the present simulations (A) and available from by the proposed interaction potential.
crystallography experiments (PDB ID: 1FM6) (B).
INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY CHARMM FORCE FIELD PARAMETERIZATION OF ROSIGLITAZONE Webb, P.; Skaf, M. S.; Polikarpov, I. Proc Natl Acad Sci 13. Houseknetch, K. L.; Cole, B. M.; Steele, P. J. Domest Anim for molecular simulations of rosiglitazone, an im- 14. Parks, D. J.; Tomkinson, N. C. O.; Villeneuve, M. S.; Blan- chard, S. G.; Willson, T. M. Bioorg Med Chem Lett 1998, 8, portant nuclear receptor ligand, with relevant pharmaceutical applications in the treatment of 15. The Protein Data Bank (PDB) archive, [online] Available at: type II diabetes. The proposed force field enables MD studies of the interactions of rosiglitazone 16. Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nils- and other TZD compounds with the nuclear re- son, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bar- ceptor PPARc, as well as with other proteins and tels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, other biomolecular systems under CHARMM. We A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.;Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; have specially focused on the energy profiles of Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; the rotating bonds, which give the molecule its Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, characteristic flexibility and are very significant D. M.; Karplus, M. J Comput Chem 2009, 30, 1545.
factors for ligand association/dissociation mecha- 17. MacKerell, A. D., Jr.; Feig, M.; Brooks, C. L. J Comput nisms and other features that depend on ligand conformational adaptations. MD simulations are 18. MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, being carried out for the rosiglitazone-PPARc R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.;Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, complex using this potential aiming to investigate K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, the concerted motions of different structural ele- T. D.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenk- ments of the LBD intermediated by rosiglitazone rich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; in the binding pocket and the pathways of ligand Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J Phys Chem dissociation from the PPARc LBD core.
19. Foloppe, N.; MacKerell, A. D., Jr. J Comput Chem 2000, 21, 20. MacKerell, A. D., Jr.; Banavali, N. K. J Comput Chem 2000, 21. Guvench, O.; MacKerell, A. D., Jr. J Mol Model 2008, 14, 1. Ribeiro, R. C. J.; Kushner, P. J.; Baxter, J. D. Annu Rev 22. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; 2. Roberts-Thompson, S. J. Immunol Cell Biol 2000, 78, 436.
Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyen-gar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; 3. Parker, J. C. Adv Drug Deliv Rev 2002, 54, 1173.
Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; 4. Blondel, A.; Renaud, J. P.; Fischer, S.; Moras, D.; Karplus, Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; 5. Kosztin, D.; Israilev, K. S.; Schulten, K. Biophys J 1999, 76, Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Strat- 6. Martı´nez, L.; Sonoda, M. T.; Webb, P.; Skaf, M. S.; Polikar- mann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, 7. Martı´nez, L.; Webb, P.; Polikarpov, I.; Skaf, M. S. J Med G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.;Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, 8. Sonoda, M. T.; Martı´nez, L.; Webb, P.; Skaf, M. S.; Polikar- J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cio- pov, I. Mol Endocrinol 2008, 22, 1565.
slowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, 9. Carlsson, P.; Burendahl, S.; Nilsson, L. Biophys J 2006, 9, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al- Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, 10. Martı´nez, L.; Polikarpov, I.; Skaf, M. S. J Phys Chem B M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision E. 01; 11. Bleicher, L.; Aparicio, R.; Nunes, F. M.; Martı´nez, L.; Dias, Gaussian, Inc.: Wallingford, CT, 2004.
S. M. G.; Figueira, A. C. M.; Santos, M. A. M.; Venturelli, 23. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhor- W. H.; Silva, R.; Donate, P. M.; Neves, F. A. R.; Simeoni, L.
shid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schul- A.; Baxter, J. D.; Webb, P.; Skaf, M. S.; Polikarpov, I. BMC ten, K. J Comput Chem 2005, 26, 1781.
24. Tetko, I. V.; Gasteiger, J.; Todeschini, R.; Mauri, A.; Living- 12. Martı´nez, L.; Nascimento, A. S.; Nunes, F. M.; Phillips, K.; stone, D.; Ertl, P.; Palyulin, V. A.; Radchenko, E. V.; Aparicio, R.; Dias, S. M. G.; Figueira, A. C. M.; Lin, J. H.; Zefirov, N. S.; Makarenko, A. S.; Tanchuk, V. Y.; Proko- Nguyen, P.; Apriletti, J. W.; Neves, F. A. R.; Baxter, J. D.; penko, V. V. J Comput Aided Mol Des 2005, 19, 453.
INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY 25. VCCLAB’s online software Alogps 2.1 (2005), [online] 31. Gampe, R. T., Jr.; Montana, V. G.; Lambert, M. H.; Miller, Available at: http://www.vcclab.org.
A. B.; Bledsoe, R. K.; Milburn, M. V.; Kliewer, S. A.; Will- 26. Singh, U. C.; Kollman, P. A. J Comput Chem 1984, 5, 129.
son, T. M.; Xu, H. E. Mol Cell 2000, 5, 545.
27. Besler, B. H.; Merz, K. M., Jr.; Kollman, P. A. J Comput 32. Darden, T.; York, D.; Pedersen, L. J. J Chem Phys 2003, 98, 28. Vorobyov, I.; Anisimov, V. M.; Greene, S.; Venable, R. M.; 33. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J Chem Moser, A.; Pastor, R. W.; MacKerell, A. D., Jr. J Chem 34. Block, J. H.; Beale, J. M., Jr., Eds. Wilson and Gisvolds Text- 29. Lee, H.; Venable, R. M.; MacKerell, A. D., Jr.; Pastor, R. W.
book of Organic Medicinal and Pharmaceutical Chemistry, 11th ed.; Lippincott Williams & Wilkins: Philadelphia, 2004.
30. Humphrey, W.; Dalke, A.; Schulten, K. J Mol Graph 1996, 35. Voet, D.; Voet, J. G.; Pratt, C. Fundamentals of Biochemis- INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

Source: http://leandro.iqm.unicamp.br/leandro/shtml/publications/2011_int_j_quantum_chem.pdf