3-Aryl-2-[1H-benzotriazol-1-yl]acrylonitriles: A novel class of potent tubulin inhibitors

Article abstract of DOI:10.1016/j.ejmech.2011.06.018

During a screening for compounds that could act against Mycobacterium tuberculosis, a series of new cellular antiproliferative agents was identified. The most cytotoxic molecules were evaluated against a panel of human cell lines derived from hematologica

Full text of DOI:10.1016/j.ejmech.2011.06.018

European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Contents lists available at ScienceDirect
European Journal of Medicinal Chemistry
journal homepage: http://www.elsevier.com/locate/ejmech
Original article
3-Aryl-2-[1H-benzotriazol-1-yl]acrylonitriles: A novel class of potent tubulin
inhibitors
,
a
a
a
b
b
Antonio Carta a , Irene Briguglio , Sandra Piras , Giampiero Boatto , Paolo La Colla , Roberta Loddo ,
Manlio Tolomeo c, Stefania Grimaudo c, Antonietta Di Cristina c, Rosaria Maria Pipitone c, Erik Laurini d,
Maria Silvia Paneni d, Paola Posocco d, Maurizio Fermeglia d, Sabrina Pricl d
*
a Dipartimento Scienze del Farmaco, Università degli Studi di Sassari, Via Muroni 23/a, 07100 Sassari, Italy
b Dipartimento di Scienze e Tecnologie Biomediche, Sezione di Microbiologia e Virologia Generale e Biotecnologie Microbiche, Università degli Studi di Cagliari,
Cittadella Universitaria, S.S. 554, Km 4,500, 09042 Monserrato, CA, Italy
c Centro Interdipartimentale di Ricerca in Oncologia Clinica (C.I.R.O.C), Università di Palermo, via del Vespro 129, 90127 Palermo, Italy
d Molecular Simulation Engineering (MOSE) Laboratory, Department of Industrial Engineering and Information Technology (DI3), University of Trieste, Piazzale Europa 1,
34127 Trieste, Italy
a r t i c l e i n f o
a b s t r a c t
Article history:
During a screening for compounds that could act against Mycobacterium tuberculosis, a series of new
cellular antiproliferative agents was identified. The most cytotoxic molecules were evaluated against
a panel of human cell lines derived from hematological and solid human tumors. In particular, (E)-2-(1H-
Received 20 March 2011
Received in revised form
10 June 2011
Accepted 11 June 2011
Available online 8 July 2011
benzo[d] [1,2,3]triazol-1-yl)-3-(4-methoxyphenyl)acrylonitrile (1) was found to be of
a potency
comparable to etoposide and greater than 6-mercaptopurine in all cell lines tested. Accordingly,
a synthesis of a new series of (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)-3-(4-R-phenyl)acry-
lonitriles was conducted in order to extend the studies of structure-activity relationship (SAR) for this
class of molecules. With the aim to evaluate if 3-aryl-2-[1H-benzotriazol-1-yl]acrylonitriles were able to
act like tubulin binding agents, the effects on cell cycle distribution of the most active compounds (1, 2a,
3 and 4) were analyzed in K562 cells. A detailed molecular modeling study of the putative binding mode
of this series of compounds on tubulin is also reported.
Keywords:
3-Aryl-2-[1H-benzotriazol-1-yl]
acrylonitriles
Antiproliferative activity
Tubulin
Homology studies
Molecular modeling study
Cell cycle distribution
Podophyllotoxin
Ó 2011 Elsevier Masson SAS. All rights reserved.
Colchicine
1. Introduction
elongation and shortening) and heavily involved in the process of
cell mitosis. Indeed, at the onset of mitosis, cytoplasmic microtu-
Microtubules are cytoskeletal filaments that play a key role in the
regulation of processes such as cell shape maintenance, segregation
of chromosomes during mitosis, location of membrane-bound
organelles, and transport. Structurally they resemble hollow
cylinders in which the basic structural unit is constituted by an
bules disassemble, and the free tubulin molecules rearrange to form
the mitotic spindle, which bridges between the chromosomes in the
center and the centrosomes at opposite poles of the cell. The spindle
remains in dynamic equilibrium with the pool of free tubulin, and
therefore must constantly add tubulin subunits to function properly
[1]. Accordingly, any alteration of this dynamic equilibrium can
induce cycle arrest and cell death (apoptosis) soon after.
a,b-tubulin heterodimer. Both globular polypeptides (approxi-
mately 55 kDa each), are homologous but not identical, and combine
stoichiometrically. Along the microtubule axis, tubulin dimers are
joined to form protofilaments that associate laterally into a polar
microtubule structure. By nature, these are structures characterized
by highly dynamic behavior (e.g., fluctuations between phases of
Attacking the microtubule system is
a common strategy
employed in tumor cell proliferation inhibition; for this reason,
the binding sites for compounds able to interfere with the micro-
tubules/tubulin dynamics on tubulin are distinct and well charac-
terized [2,3]. Accordingly, those molecules able to bind at the
two natural compounds vinblastine and colchicine binding sites
could inhibit tubulin assembly and are classified as microtubule
destabilizing agents. In contrast, those compounds that bind to the
* Corresponding author. Tel.: þ39 079 228722; fax: þ39 079 228720.
E-mail address: acarta@uniss.it (A. Carta).
0223-5234/$ e see front matter Ó 2011 Elsevier Masson SAS. All rights reserved.
doi:10.1016/j.ejmech.2011.06.018

4152
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
paclitaxel binding site prevent microtubule disassembly and are
classified as microtubule stabilizing compounds [4].
In order to check whether the nature of the substituents on the
benzene ring could act favorably toward the biological activity of
these molecules, we then synthesized the new series of (E)-
In the effort of developing new agents against Mycobacterium
strains, our group found that some compounds belonging of the
class of 3-aryl-1-benzotriazol-1-yl-acrilonitrile derivatives were
active against the proliferation of a number of liquid and solid
human tumors. The benzotriazole moiety characterizing this series
of compounds has been amply studied by Katritzky et al. as
a synthetic auxiliary, as it constitutes a good leaving group after
reaction with a variety of carbonyl groups [5e7], whereas Sparatore
et al. studied the biological properties [8,9] and described the
interesting pharmacological activities of several benzotriazole
derivatives bearing substituents at positions 1 or 2 [10e18]
Furthermore, some of us have developed a separate chemistry
and investigated the biological potential of benzotriazoles [19e23].
In particular, we described several 3-aryl-2-benzotriazol-2-yl-
acrilonitriles which were found endowed with antitubercular
activity [24e27], partly retained by the corresponding amides and
carboxylic acids [25]. In addition, we described 2-benzotriazol-(1)-
yl-3-phenylacrilonitriles, which possess an interesting anti-
proliferative activity [26,27]. All these findings are summarized in
Table 1, from which the promising activity in the nanomolar range
of these compound series against both hematological and solid
human tumors is evident.
2-(5,6-dichloro-1H-benzo[d]
[1,2,3]triazol-1-yl)-3-(4-R-phenyl)
acrylonitriles 1ae3a depicted in Fig. 1. The substituents on the
phenyl-moiety (indicated with R in Fig. 1) have been selected
among those showing the best antiproliferative activity in the
previous series. All new derivatives have been evaluated for cyto-
toxicity against MT-4 cells and for their antiproliferative activity
against a panel of cell lines derived from hematological and solid
human tumors.
With the aim of understanding the mechanism of action of this
class of antiproliferative agents, at the beginning we performed an
extensive analysis of the scientific literature to detect possible
structural analogies with other molecules endowed with similar
antiproliferative activity. Partial similarities were found with many
molecules, but the most compelling were those with derivatives of
resveratrol and trimethoxy-aryl-thiophene, well-known antimi-
totic agents acting binding at the colchicine-binding site of tubulin
[28]. Recent studies conducted by other groups on these derivatives
have confirmed their ability to inhibit the polymerization of tubulin
thereby preventing the formation of spindle cells by blocking cell
replication in its metaphase [29e31]. Further, these experimental
pathways have been paralleled and assisted by molecular modeling
investigations, aiming both at optimizing the synthetic strategy and
efforts, and at exploring the eventual molecular interactions
between antimitotic inhibitors and tubulin. Therefore, we decided
to attempt to develop a structure-activity relationship (SAR) for all
these compound sets by means of molecular modeling ansatz based
on a hierarchical procedure consisting of i) three-dimensional
pharmacophore development for the compound series, ii) dock-
ing studies guided by the pharmacophore requirements, and
On the basis of these results, (E)-2-(1H-benzo[d] [1,2,3]triazol-
1-yl)-3-(4-methoxyphenyl)acrylonitrile 1 was identified as lead
compound [26]. Interestingly, the benzotriazol-(1)-yl-3-
phenylacrilonitrile derivative
1 bears no substituents at the
benzene ring of the condensed heterocyclic moiety. In order to
further potentiate its antiproliferative activity, several 5,6-
dimethylbenzotriazol-(1)-yl-3-phenylacrilonitriles were synthe-
sized in previous efforts, although with disappointing results [27].
Table 1
Antiproliferative activity against hematological and solid human tumor-derived cell lines of (E)-2-(1H-benzo[d] [1,2,3]triazol-1-yl)-3-(3-R1-4-R2-phenyl)acrylonitriles 1e4.
Compd
R1
R2
aCC50
MT-4
bCC50
fSKM EL28
cCCR FCCM
dWIL2 NS
eCCR FSB
gMCF7
hSKM ES1
iHepG2
jDU145
1
2
3
4
H
H
H
OCH3
CH3
Br
0.05
0.5
0.4
0.1
0.1
0.2
0.2
0.4
0.9
1.0
0.1
0.2
0.6
0.6
3.0
0.2
0.09
0.07
0.3
0.6
1.1
0.2
0.3
1.1
3.3
15
0.1
0.7
0.7
1.3
3.2
1.0
0.6
1.7
2.4
3.3
58.0
0.3
0.8
2.2
2.1
3.1.
8.0
0.7
0.6
0.5
1.7
3.4
2.0
0.4
O-CH2-O
6MP
Etoposide
0.09
0.09
0.1
1.2
a
Compound concentration (
m
M) required to reduce the viability of mock-infected MT-4 cells by 50%, as determined by the MTT method. Data represent mean values for two
separate experiments. Variation among duplicate samples was less than 15%.
b
Compound concentration (mM) required to reduce cell proliferation by 50%, as determined by the MTT method, under conditions allowing untreated controls to undergo at
least three consecutive rounds of multiplication. Data represent mean values for two separate experiments. Variation among duplicate samples was less than 15%.
c
CD4þ human acute T-lymphoblastic leukemia.
Human splenic B-lymphoblastoid cells.
Human acute B-lymphoblastic leukemia.
Human skin melanoma.
Human breast adenocarcinoma.
Human lung squamous carcinoma.
Human hepatocellular carcinoma.
Human prostate carcinoma.
d
e
f
g
h
i
j

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4153
3. Microbiology
All new compounds were evaluated for cytotoxicity toward
MT-4 cells and for antiproliferative activity against a panel of
human cell lines derived from hematological (CCRF-CEM, WIL-2NS
and CCRFSB) and solid (SKMEL28, MCF7, SKMES-1, HepG2, and
DU145) human tumors. The same antitumor agents previously
used, 6-mercapto-purine (6 MP) and etoposide (a semisynthetic
derivative of podophyllotoxin [33]) were tested as reference drugs.
4. Cell cycle analysis
Fig. 1. Structure of the (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)-3-(4-R-
phenyl)acrylonitrile derivatives considered in this work.
Because molecules exhibiting activity on tubulin should cause
the alteration of cell cycle parameters leading to a preferential
G2/M blockade, the effects of the most active compound of the new
series (2a) on cell cycle distribution were analyzed in K562 cells.
The relevant data were then compared to those obtained with the
non-halogenated compounds 1, 3 and 4. Cells were exposed to each
compound used at the TGI (total growth inhibition concentration),
and after 24 h cell cycle analysis was performed as described in
materials and methods. As expected from the results shown in
Fig. 2, compound 2a induced a marked increase of the G2-M peak
(67% vs. 18% of the non-treated control), consistent with its tubulin
binding activity. As shown in Fig. 2, compounds 1 and 4 also act on
the G2/M phase of the cell cycle (70.5% and 66.34%, respectively),
while compound 3 caused an increase of both G0-G1 and G2-M
peaks, suggesting that its cell growth inhibition activity could be
only partially correlated to its tubulin binding activity.
iii) estimation of the affinity of each compound toward tubulin by
means of atomistic molecular dynamics simulations. According to
this computational recipe, a putative binding mode of our series of
compounds onto tubulin was obtained, and the structural bases of
the biological results discussed above were investigated in details.
2. Chemistry
The intermediate compound 2-(5,6-dichloro-1H-benzo[d]
[1,2,3]triazol-1-yl)acetonitrile 6a was prepared as described in
Scheme 1 by condensation of 5,6-dichloro-1H-benzo[d] [1,2,3]tri-
azole 5 with chloroacetonitrile in dimethylformamide (DMF) in the
presence of triethylamine (Et3N). Compound
prepared following a procedure previously reported [32]. The
5 was in turn
5. [3H]Colchicine competition binding assay
synthesis of the new series of (E)-2-(5,6-dichloro-1H-benzo[d]
[1,2,3]triazol-1-yl)-3-(4-R-phenyl)acrylonitriles
1ae3a
was
To further confirm the tubulin binding activity of the title
compounds, we assessed their capability to compete with colchi-
cine or paclitaxel for binding to tubulin by resorting to a competi-
tive scintillation proximity assay [34] using compound 4 as
reference. As shown in Fig. 3, we found that compound 4
competitively inhibited [3H]colchicine binding to biotinylated
accomplished by a straightforward condensation of 6a with the
appropriate commercially available aldehydes, following the
procedure early described for 1 [26]. In this new series, among
the two possible geometric isomers (E/Z), E-isomers were obtained
as the sole product. Spectral (IR, UVevis, 1H NMR) and analytical
(elemental analyses, MS) data obtained for all new compounds are
in accordance with those of the previously described counterparts
[26], and support the assigned chemical structure.
tubulin with an IC50 value of 0.85
value for colchicine was 1.02 M. At the same time, compound 4 did
not compete with [3H]paclitaxel. No stabilization of the colchicine
mM whereas the corresponding
m
Scheme 1. Synthetic route to (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)-3-(4-R-phenyl)acrylonitrile derivatives 1ae3a.

4154
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Fig. 2. Effects of compounds 1 (panel b), 2a (panel c), 3 (panel d), and 4 (panel e) on DNA content/cell following the treatment of K562 cells for 24 h. The cells were cultured without
any compounds (panel a) or with each compound used at TGI (total growth inhibition concentration). Cell cycle distribution was analyzed by the standard propidium iodide
procedure, as described in Material and Methods. Sub-G0-G1 (A), G0-G1, S, and G2-M cells are indicated in panel a.
binding was observed, as it is found for Vinca site binders [35,36].
We could then conclude that binding to tubulin at the colchicine
binding site for compound 4 and all other compounds of this series
is highly probable. Accordingly, the capacity to interact with the
mitotic spindle contributes to the antiproliferative activity of these
compounds.
a variety of parameters such as IC50, EC50 or CC50 values, binding/
inhibition constant, or even uptake rates may be used. Because our
interest was primarily the study of high-affinity inhibitors of
tubulin in order to prevent cancer cell proliferation, we decided to
use the CC50 values of our compounds against MT-4 cells as the
biological parameter to develop the pharmacophore model. To
obtain a realistic 3D quantitative-structure activity relationship
(QSAR) adequately describing the interactions between a ligand
series and their target with a high predictability, however, several
aspects have to be considered. Firstly, a collection of minimum
15 chemically diverse compounds is necessary to constitute the
so-called training set; also, the biological activity of the training set
compounds should cover at least 4 orders of magnitude to avoid
model biasing. Thus, to enrich the compound sets listed in Tables 1
and 2 of the present paper, for the generation of the 3D pharma-
cophore model we also considered three series of previously
synthesizes compounds tested for activity as tubulin inhibitors
using the same methodology employed in the present work
[24e27]. These compounds are shown in Table 3, together with
their biological data.
6. Molecular modeling
6.1. Three-dimensional pharmacophore model generation
A prerequisite for the development of a three-dimensional (3D)
pharmacophore model for tubulin inhibitors is the correlation of
a characteristic and reproducible biological activity to the structural
information of the respective compound. As biological activity,
A three-dimensional pharmacophore model captures the three-
dimensional arrangement of the structural features shared by all
active molecules that are presumably essential for the desired
pharmacological activity. The HypoGen algorithm of Discovery
Studio (DS) Catalyst [37] allows a maximum of five features to be
present in the pharmacophore generation process. Preliminary
runs including, among others, the negative charge (NC) site, the
positive ionizable (PI), the hydrogen bond acceptor lipid (HBAl),
and the hydrophobic aliphatic (HYAl), confirmed that these
features were never used in the generation of the pharmacophore
models. Thus, all these features were removed from the list. In
summary, the following five chemical features were taken into
account for hypothesis generation with HypoGen: hydrogen
Fig. 3. [3H]Colchicine competition binding assay of compound 4 and colchicine. Radio
labeled colchicine, unlabeled compound, and biotin-labeled tubulin were incubated
together at 37 ꢀC for 2 h. Symbols legend: B, colchicine: 6, 4.

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4155
Table 2
Cytotoxicity and antiproliferative activity against hematological and solid human tumor-derived cell lines of compounds 1ae3a.
Compd
aCC50
MT-4
bCC50
fSKM EL28
cCCRF CEM
dWIL2 NS
eCCR FSB
gMCF7
hSKM ES1
iHepG2
jDU145
1a
2a
14
1.7
>100
2.8
>100
7.0
>100
2.8
>100
20
>100
27
>100
27
>100
21
>100
17
3a
5.4
5.3
18
12
10
8.4
10
12
7
6MP
Etoposide
0.1
0.09
1.0
0.09
3.0
0.2
1.1
0.1
15
1.2
3.2
1.0
58.0
0.3
8.0
0.7
2.0
0.4
a
Compound concentration (mM) required to reduce the viability of mock-infected MT-4 cells by 50%, as determined by the MTT method. Data represent mean values for two
separate experiments. Variation among duplicate samples was less than 15%.
b
Compound concentration (mM) required to reduce cell proliferation by 50%, as determined by the MTT method, under conditions allowing untreated controls to undergo at
least three consecutive rounds of multiplication. Data represent mean values for two separate experiments. Variation among duplicate samples was less than 15%.
c
CD4þ human acute T-lymphoblastic leukemia.
Human splenic B-lymphoblastoid cells.
Human acute B-lymphoblastic leukemia.
Human skin melanoma.
Human breast adenocarcinoma.
Human lung.
Human hepatocellular carcinoma.
Human prostate carcinoma.
d
e
f
g
h
i
j
bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic
aromatic (HYAr), ring aromatic (RA) and more generic hydrophobic
feature (HY).
predicted to be quite less active (estimated CC50 ¼ 445
mM,
experimental CC50 ¼ 375 M). Finally, Fig. 4(C) depicts 33, the least
m
active compound, mapped onto Hypo1. As can be seen from this
panel, the molecule 33 fails to map one of the HBD features and the
HYAr function. The remaining HBD is mapped again by the cyano
group and the RA is overlaid by the phenyl ring of the aromatic ring
of the benzotriazole. This poor mapping scheme fully accounts
for the very low activity of this compound (estimated
A total of 10 hypotheses were generated by the HypoGen algo-
rithm, all characterized by the same three features: HBA, HYAr, and
Ar. The total hypothesis cost of these ten best models varies between
100.1 for the best ranked model (Hypo1) to 108.4 for the lowest
ranked one (Hypo10), as shown inTable 4. Such a confined difference
(8.6 bits) reflects both the homogeneity of the generated hypotheses
and the adequacy of the molecular training set. The difference
between the null and the fixed costs, which should be higher than 70
to guarantee a robust correlation, is 104.2 in our case. This corre-
sponds to a chance of true correlation in the data greater than 90%
[38]. Further, the evidence that for all hypotheses the total costs are
sensibly closer to the fixed cost (109.7) than to the null cost (212.1)
also constitutes an indication that meaningful models are obtained.
As reported in Table 4, for Hypo1 the configurational cost value is
equal to 11.6, which is well below the threshold value of 17. Finally,
the root-mean-square deviations (RMSD) and the correlation coef-
CC50 ¼ 42,800
mM, experimental CC50 ¼ 20,000
mM).
To check further on the usefulness of the generated 3D phar-
macophore, we validated the model by mapping a test set of 11
compounds (see Tables 3 and 6). Indeed, a good correlation coef-
ficient (0.84) was observed when a regression analysis was per-
formed by mapping the test set onto the features of the best
pharmacophore hypothesis Hypo1. The predicted and the experi-
mental CC50 values for the test set along with the respective errors
are given in Table 6. The average error in predicting the affinity of
the test set molecules is 1.4; therefore, given the inherent simplicity
of the pharmacophoric approach, and considering the intrinsic
variability of the biological responses, we can conclude that the
ability of the present 3D pharmacophore model in predicting the
activity of these molecular series as tubulin inhibitor is satisfactory.
Fig. 5 shows the mapping of one test set molecule 14 to Hypo1.
A second test was performed to verify the statistical significance
of 3D pharmacophore model Hypo1, based on a randomization
procedure derived from the Fisher method [39] available in the DS
Catalyst suite of programs. According to the validation procedure,
the experimental affinities of the compounds in the training set
were scrambled randomly, and the resulting new training set was
used for a new HypoGen run. The parameters used in running these
calculations were the same employed in the initial HypoGen
calculations and, since a 99% confidence level was selected, 99
random hypothesis runs were performed. The resulting data clearly
indicate that all values generated after randomization produced
hypotheses with no predictive values. Indeed, none of the outcome
hypotheses had lower cost score, better correlation or smaller root-
mean-square deviation than the initial one. Table 7 lists the first 10
lowest total score values of the resulting 99 hypotheses for our test
set molecules. In conclusion, there is a 99% chance for the best
hypothesis to represent a true correlation in the training set affinity
data for the present classes of compounds.
ficients (r) between estimated and experimental affinities range
from 1.6 to 1.9, and from 0.83 to 0.88, respectively. Considering that
all the generated pharmacophores mapthe molecules of the training
set in a similar way, the first model (Hypo1), characterized by the
highest cost difference, the lowest RMDS, and the best
selected for further analysis.
r values was
The affinities of 25 compounds in the training set estimated
using Hypo1 are reported in Table 5, along with the experimental
values and the relevant errors (expressed as the ratio between
estimated and experimental values). This Table clearly shows that
all activities are well predicted, all errors being mostly below 3.6.
Fig. 4 shows the mapping of compounds 1, 33, and 36 onto Hypo1,
respectively. As can be seen from Fig. 4(A), in 1 the aromatic ring of
the benzotriazole moiety matches the HYAr feature, whilst the
cyano group and one nitrogen atom of the triazole ring nicely map
the two HBD functions, respectively. The phenyl group bearing the
eOCH3 substituent finally overlaps the RA feature of the pharma-
cophore. The estimated activity for 1 is 0.13
mM, while the corre-
sponding experimental CC50 is 0.05 M. Fig. 4(B) is an example of
m
a pharmacophore mapping of a compound that is less active than
the former one. Compound 36 does not map all the features
encoded in Hypo1. In fact, 36 maps the HYAr function, again by
means of the benzomoiety of the heterocycle; however, it does not
map one of the HBD feature, and the mapping of the remaining RA
is not perfect. According to this partial mapping, this compound is
As a last statistical test, the leave-one-out method, which consists
of re-computing the hypothesis by excluding from the training set
one molecule at a time, was performed. Basically, this test aims at

4156
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Table 3
Table 4
Additional set of molecules used in the training set and test set for 3D pharmaco-
phore model generation.
Composition, ranking score, statistical parameters, and cost analysis (expressed in
bits) of the top 10 generated hypotheses.
Hypothesis
Total cost
D
Costa
RMDS
r
1
2
3
4
5
6
7
8
9
10
138.4
140.1
140.8
141.6
142.0
143.8
144.0
146.0
146.7
147.0
73.7
72.0
71.3
70.5
70.1
68.3
68.1
66.1
65.4
65.1
1.64
1.74
1.75
1.71
1.79
1.83
1.82
1.86
1.87
1.91
0.876
0.859
0.857
0.846
0.851
0.843
0.845
0.839
0.837
0.828
a
D
Cost ¼ Cost difference ¼ null costꢁtotal cost. Null cost: 212.1. Fixed cost:
107.9. Configurational cost: 11.6.
hypotheses generated according to this method we did not obtained
meaningful differences between Hypo1 and each hypothesis
resulting from the exclusion of one compound at a time, in terms of
correlation coefficients, feature composition of the pharmacophore,
and quality of the predicted affinity of the excluded molecule.
According to all evidences highlighted above, we believe that
our in silico 3D pharmacophore model accounts for tubulin inhib-
itory activity of our set of compounds and, despite its inherent
simplicity, its predictive power is quite robust and can then be
employed as a guide for discovering a possible binding site and for
the successive docking of the inhibitors focus of the present study
within on their putative binding site on human tubulin.
a
Compd
R1
R2
R3
R4
R5
R6
CC50
Series I
7
8
H
H
H
CH3
H
CH3
H
H
H
H
H
H
H
CH3
H
H
H
H
H
CH3
CF3
CH3
H
H
H
H
H
H
H
CH3
H
H
H
H
H
H
H
H
H
H
H
Cl
H
OCH3
H
H
I
NO2
H
H
OeCH2eO
H
NO2
NO2
NO2
OeCH2eO
OCH3
H
OCH3
H
OCH3
H
H
H
H
H
H
H
OCH3
H
H
H
H
H
H
H
H
H
H
H
0.1
0.4
0.4
0.5
0.7
1.3
2.5
4.0
4.7
6.0
7.5
8.3
17
Cl
H
H
H
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
OCH3
H
H
Cl
OCH3
H
COOH
Cl
H
H
H
H
H
6.2. Docking and calculations of binding free energy by MM/PBSA
method
24
25
31
52
H
CF3
CF3
H
H
OeCH2eO
H
Docking studies were carried out to obtain putative binding
modes of our series of compounds to tubulin, and to investigate the
structural basis of the biological results discussed above. In the
absence of any information about the docking modes for our
OCH3
200
Series II
25
26
27
28
29
30
31
32
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
H
NO2
H
Cl
COOH
CH3
F
OCH3
OCH3
H
H
H
H
H
H
H
OCH3
H
H
0.4
12
16
16
20
Table 5
Experimental and estimated activity values (CC50) of the training set compounds
used to develop pharmacophore hypothesis for tubulin inhibitors.
H
OCH3
OCH3
OeCH2eO
OCH3
32
360
7500
20000
12000
CC50 (mM)
Errora
33
34
Compd
Experimental
Estimated
Mapping
CH3
CH3
OCH3
1
7
8
0.05
0.1
0.4
0.7
1.3
2.5
4.7
6.0
7.5
8.3
24
31
200
12
0.13
0.38
0.58
1.1
1.6
2.3
2.1
9.9
2.1
13
19
56
258
8.0
33
920
9100
42800
99
445
1400
6700
18
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 * 1
1 1 1 1
1 1 1 1
1 1 * 1
1 * 1 1
1 1 1 1
1 1 1 1
1 * 1 1
1 * * 1
1 1 * 1
1 * * 1
1 * 1 1
1 * 1 1
1 * * 1
1 * * 1
1 1 1 1
1 * 1 1
1 1 1 1
2.6
3.8
1.5
1.6
1.2
ꢁ1.1
ꢁ2.2
1.7
a
Compd
R1
H
R2
H
ARYL
CC50
Series III
35
S
11
12
13
15
16
17
18
20
22
24
26
29
31
32
33
35
36
37
38
1a
2a
3a
132
375
O
36
37
H
H
H
H
ꢁ3.6
1.6
ꢁ1.3
1.8
ꢁ3.4
ꢁ1.5
1.7
2.6
1.2
2.1
7.5
5000
20
H
N
360
7500
20000
132
375
5000
10000
14
38
H
H
10000
a
Compound concentration (mM) required to reduce the viability of mock-infected
1.2
MT-4 cells by 50%, as determined by the MTT method.
ꢁ3.6
ꢁ1.5
1.3
1.6
ꢁ1.1
verifying whether the correlation is strongly dependent on one
particular compound in the training set. The test is positive if the
affinity of each excluded molecule is correctly predicted by the
corresponding one-missing hypothesis. For each of the 23 new
1.7
5.4
2.7
4.8
a
Values in the error column represent the ratio of the estimated to experimental
affinity, or its negative inverse if the ratio is less than one.

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4157
Fig. 4. Pharmacophore mapping of compound 1 (A), 36 (B), and 33 (C) in the training set. The hypothesis features are portrayed as mashed spheres, color-coded as follows: light
green, HBA; light blue, HYAr; light orange, RA. HBA is actually represented as a pair of spheres, the smaller sphere representing the location of the HBA atom on the ligand and the
larger one the location of an HB donor on the receptor. Compounds are depicted as atom-colored sticks-and-balls: carbon, dark gray; oxygen, red; nitrogen, blue, hydrogen, white.
(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
molecules, we decided to apply some blank tests to verify the reli-
ability of the overall computational procedure adopted. Accordingly,
we began by modeling, docking and calculating the free energy of
binding (and the corresponding IC50 values) for podophyllotoxin
In the light of these blank-test results, and of the fact that
a similar computational ansatz was already successfully applied by
our group to estimate the activity of other protein inhibitors for
which no structural information were available [44e47], the
conceived modeling/docking procedure was applied for predicting
the binding mode of the present series of compounds in the tubulin
binding site targeted by podophyllotoxin and colchicine.
Twenty nanoseconds molecular dynamics simulations were
then performed on colchicine and podophyllotoxin in complex
with tubulin. Overall, compared to the corresponding crystal
structures, the RMSDs of the ligands and the protein backbone are
about 0.5 and 0.8 Å, respectively. The same protocol was applied to
MDL 27048, CHAL101, and a subset of our title compounds. Fig. 8
shows an equilibrated snapshot of compound 1/tubulin complex
as an example. From this figure it is evident how compound 1
and colchicine to human b-tubulin. The best docked structures for
both podophyllotoxin and colchicine, i.e., the configurations with
lowest docking energy in each prevailing cluster, were then
compared with the corresponding crystal structures and subjected
to MM/PBSA calculations [40] to obtain the free energy of binding
and, hence, the corresponding IC50 value [41,42].
Fig. 6 shows a superposition between the co-crystallized confor-
mations of podophyllotoxin and colchicine into the allosteric binding
site of tubulin, and thecorrespondingdocked conformationsobtained
upon application of the computational strategy adopted in this work.
As it can be seen, the agreement between crystal and docked struc-
tures is excellent: the root-mean-square deviation (RMSD) between
the docked configurations and the relevant crystal structures of these
test molecules is equal to 0.5 Å and 0.3 Å, respectively.
To further test whether our computational recipe was feasible
to predict unknown ligand binding modes to tubulin, we decided
to apply it to MDL 27048 ([trans-1-(2,5-dimethoxypheny1)-
3-[4-(dimethylamino)phenyl]-2-methyl-2-propen-1-one]), and CH-
ALC101, two structurally different but potent destabilizers of
microtubule formation for which only biological data and/or
modeling speculations are available [43,44]. As depicted in Fig. 7,
MDL 27048 and CHALC101 are found to share the same binding site
of podophyllotoxin and colchicine, in harmony with previous
modeling results [44], the overlapping of the binding site accounting
for the experimentally verified competition among colchicine,
podophyllotoxin, MDL 27048 and CHALC101 in tubulin binding.
interacts extensively with
contacts (coherently with previous modeling evidences reported
for podophyllotoxin [48]), while only two residues (i.e., Ala 180
and Val 181) of -tubulin are engaged in hydrophobic interactions
b-tubulin mainly via van der Waals
a
a
a
with the inhibitor. More importantly, the predicted binding mode is
found to comply with the pharmacophore hypothesis requirements
(see Fig. 4A). In fact, the first hydrogen bond acceptor (HBA) feature
on the inhibitor, represented by the triazole N3, localizes hydrogen
bond acceptor structures that are in an ideal position for forming
a stabilizing hydrogen bond (average dynamic distances ꢂ 2.8 Å)
with the donor guanidinium group of Lysb350 on the receptor.
The benzotriazole moiety is embedded in a hydrophobic pocket
consisting of the side chains of Lys 252, Leu 253, Asn 256,
b
b
b
Table 6
Experimental and estimated activity values of the test set compounds used to
develop pharmacophore hypothesis for tubulin inhibitors.
CC50 (mM)
Compd
Experimental
Estimated
Errora
9
0.4
0.5
4
17
25
52
0.4
16
16
32
0.8
1.8
10
4
14
76
0.1
10
23
35
2.0
3.6
2.5
ꢁ4.3
ꢁ1.8
1.5
ꢁ4.0
ꢁ1.6
1.4
10
14
19
21
23
25
27
28
30
34
1.1
1.8
12000
22000
a
Values in the error column represent the ratio of the estimated to experimental
Fig. 5. Pharmacophore mapping of compound 14 in the test set. Hypothesis features
affinity, or its negative inverse if the ratio is less than one.
are color-coded as in Fig. 3.

4158
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Table 7
computational procedure. Table 8 shows the binding free energy
components averaged over the MD trajectory for each of these
complexes: all ligands are correctly ranked for affinity, and the
Output parameters of the 10 lowest cost hypotheses resulting from the statistical
evaluation according to the Fisher method validation procedure for the tubulin
inhibitors.
magnitudes of the calculated
DGbind values are in good agreement
Hypothesis
r
RMDS
Total cost
with the corresponding experimental evidences [43,44,48].
1
2
3
4
0.816
0.802
0.789
0.782
0.769
0.876
0.745
0.743
0.740
0.737
0.728
1.689
1.711
1.745
1.796
1.802
1.644
1.834
1.871
1.992
2.014
2.109
140.0
152.1
154.4
155.4
155.8
138.4
158.1
159.8
160.0
160.4
160.8
Furthermore, Table 8 shows the corresponding affinities of
compounds (1e4 and 1ae3a) to the protein. If we examine the
contributions to binding across our compounds we can see that,
although both the electrostatic and the van der Waals interac-
5
Hypo1
6
7
8
9
10
tion energy contributions are important for binding, the overall
D
Gnsopl) constitute the driving force to
vdW
non-polar interactions (
D
E
þ
MM
tubulin binding for these molecules. Interestingly, and somewhat
contrarily to what could be expected, the presence of the [1,3] dioxolo
moiety in compound 4 does not afford a substantial favorable contri-
bution to the dispersion energies compared to that of the other
molecules (see Table 8), but mainly contributes to a molecular
conformation amenable for fitting into the protein binding site.
Met
hydrophobic aromatic feature (HYAr) of the 3D pharmacophore.
The p-methoxyphenyl ring fits in a pocket lined by residues Ile 236,
Ser 239, Leu 240, Leu 246, Ala 248, Leu 253, the cyano N atom
being located within hydrogen bond distance (average dynamic
312, as required
b257, Thrb312, Alab314, Lysb350, which all concur to satisfy the
Further analysis of the DGbind components reveals that while the van
b
der Waals and non-polar solvation energies drive binding for all
ligands, sterics alone do not completelyexplain the affinity differential.
We also find that the overall electrostatic component to binding either
provides no thermodynamic benefit or even opposes association
despitethepresenceofsomepolargroupsintheactivesite. Amongthe
five modeled compounds, the two most effective ligands 1 and 4 have
b
b
b
b
b
length (ADL) ¼ 3.2 Å) of the hydroxyl group of Thr
b
by the second HBA feature of the pharmacophore. The C]C bond is
stabilized by favorable hydrophobic interactions due to the vicinity
of the side chains of Leu
bearing the eOCH3 substituent is also extensively engaged in
hydrophobic interactions with the side chains of residues Ile 376,
Leu 253, and Ser 239, thus satisfying the last (RA) pharmacophore
b253 andIleb236, whilst the phenyl group
overall unfavorable electrostatic component to binding
ele
D
( E
þ
D
Gesolel )equaltoþ2.9kcal/molandþ1.3kcal/mol,respectively.
b
MM
As mentioned above, it appears that 1 gains some internal electrostatic
energy (
b
b
ele
D
E
) by having the eCN and eOCH3 groups that interact
feature. Finally, along the 20 ns MD trajectory another albeit
intermittent hydrogen bond interaction is also detected between
the oxygen atom of the eOCH3 group and the eOH of Serb239
MM
favorably with the eOH group of Thr
whilst 4 counteracts the unfavorable desolvation term (
b
312 and Ser
b
239, respectively,
D
Gesolel ) by
interacting with water molecules in the binding site (data not shown).
(ADL ¼ 3.6 Å), which contributes to anchor the molecule in the
proper, productive orientation within the binding site, thus
contributing further to the high activity of this compound.
By resorting tothe MM/PBSA approach we subsequentlyestimated
6.3. Projection of the free energy to each residue of tubulin
the free energy of binding,
MDL 27048 and CHALC101 to human tubulin as a benchmark for our
DGbind, for podophyllotoxin, colchicine,
In order to gain extra insight into the possible mechanism of
tubulin-drug binding, the values of
DGbind calculated from the
Fig. 6. Comparison between the co-crystallized (brown) and docked (light blue) structure of podophyllotoxin (top) and colchicine (bottom) into the allosteric binding site of human
-tubulin. Hydrogen atoms, water molecules, and counterions are omitted for clarity. Note that human -tubulin and the crystallized Bos Taurus -tubulin exhibit 98% sequential
homology. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
b
b
b

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4159
Fig. 7. (A) Superposition of the average structures of 1 (magenta), podophyllotoxin (dark red), colchicine (dark green), MDL 27048 (dark slate blue), and CHALC101 (dark kaki), and
(B) detailed view of the superposed docked conformations of the same compounds in the allosteric binding site of tubulin. Hydrogen atoms, water molecules, and counterions are
omitted for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
MM/PBSA methodology were further decomposed into the
nonbonded contribution afforded by individual tubulin binding site
residues. As can be seen from Fig. 9, several residues belonging to
the binding site responded differently to the different compounds.
Further, the binding site residues not only impacted their own
binding free energy interactions with the inhibitors but also
influenced the interactions of other residues with the inhibitors by
inducing slight alterations in the geometry of the binding site.
Favorable electrostatic interactions opposed by the polar solva-
6.4. Van der Waals energy contribution from inhibitors
To further explore the mechanism of loss in free energy of
binding between the different inhibitors and tubulin, the van der
Waals contributions were calculated for each atom of each
compound 2e4 and 1ae3a, and compared to those relative to the
lead molecule 1. On average, each inhibitor has 35 atoms, 12 of
which are hydrogens which afford a limited contribution to the
dispersion energy term. Thus, we present and discuss data for the
heavy atoms of each inhibitor only.
From the structural standpoint, both molecular series (1e4 and
1ae3a) can be considered formed by three major moieties: A) the
benzotriazole group, B) the nitrile-substituted double bond, and C)
the substituted phenyl ring (Fig. 11). To compare the energy vari-
ation between these three moieties for the different inhibitors, we
defined a van der Waals energy variation index, VVI, as:
tion energy penalty also apply to individual residues. A change in
ele
electrostatic energy (
D
E
) of any residue is always associated with
MM
an almost equal but opposite compensation in solvation energy
Gsol), as illustrated in Fig.10(A) for compounds 1e4 as an example.
(
D
The correlation coefficients for all inhibitors range between 0.95 and
0.99. This high correlation between these two counteracting
contributions to
nonbonded energy term
D
Gbind made the change of the van der Waals
vdW
D
E
the largest factor in the loss of
ꢀ
ꢁ
MM
P
EvndW;i ꢁ Ev1dW;i
binding free energy between tubulin and all 3-Aryl-2-[1H-benzo-
triazol-1-yl]acrylonitrile derivatives. To highlight those residues
i
VVI ¼
ꢃ 100
(1)
P
Ev1dW;i
with a significant difference between the lead compound 1 and all
vdW
MM
i
remaining inhibitors 2e4 and 1ae3a, a cut-off of 1 kcal/mol of E
D
was considered [49]. Thus, only those residues with a loss of van der
Waalsenergygreater than thecut-offare plotted in Fig.10 (B) and(C).
As can be seen by comparing Fig. 10 (B) and (C), some residues
which are critical for the favorable binding in the case of the lead
compound 1 (as well as for the other member of the series 2e4),
become unfavorable in the analogs bearing two chlorine atoms as
substituents on the benzotriazole group (1ae3a). This information,
coupled with the corresponding quantification of the van der Waals
energy loss, as discussed below, can yield important information
both for rationalizing the activity of the present molecules and the
design of new microtubule destabilizing agents.
in which i is the i-th atom within a specific moiety. Interestingly, the
benzotriazole heterocycle has a relatively low and constant value of
VVI in the series 2e4, the same moiety being characterized by
somewhat higher VVI values in compounds 1ae3a (Table 9). On the
other hand, both the phenyl ring and the double bond have
significantly higher and more variable VVI values in compounds
1ae3a with respect to the corresponding counterparts in
compounds 2e4 (Table 9).
A deeper examination of the MD trajectories of each inhibitor in
complex with tubulin reveals that, particularly when unsub-
stituted, the benzotriazole group maintains a relatively stable
Fig. 8. Overall view of compound 1 in complex with tubulin (left) and detailed view of tubulin residues interacting with 1 (right). Tubulin residue color code: A
181, dark green; I 236, dark kaki; S 239, dark magenta; L 240, gold; L 246, sienna; A 248, spring green; K 252, orange red; L 253, dark red; N 256, dark slate blue; M
pink; T 312, white; A 314, navy blue; V 316, chartreuse; K 350, dodger blue; A 352, sky blue; I 376, dark olive green. (For interpretation of the references to colour in this figure
legend, the reader is referred to the web version of this article.)
a180, dark cyan;
V
a
b
b
b
b
b
b
b
b
b257, hot
b
b
b
b
b
b

4160
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Table 8
Free energy components, total binding energy
DGbind from molecular dynamics simulations on human tubuline in complex with compounds 1e4, and 1ae3a. The calculated
values for the test compounds podophyllotoxin, colchicine, MDL 27048 and CHALC101 are also reported for comparison and validation.
np
sol
ele
vdW
MM
ele
sol
Compd
D
E
D
E
D
EMM
D
G
D
G
D
Gsol
ꢁTD
S
DGbind
MM
1
2
3
4
1a
2a
3a
ꢁ22.8 ꢄ 0.3
ꢁ21.0 ꢄ 0.2
ꢁ20.9 ꢄ 0.2
ꢁ23.9 ꢄ 0.3
ꢁ21.1 ꢄ 0.1
ꢁ21.8 ꢄ 0.2
ꢁ21.9 ꢄ 0.3
ꢁ26.4 ꢄ 0.3
ꢁ27.4 ꢄ 0.4
ꢁ29.5 ꢄ 0.3
ꢁ26.2 ꢄ 0.2
ꢁ25.6 ꢄ 0.2
ꢁ23.5 ꢄ 0.1
ꢁ23.2 ꢄ 0.2
ꢁ21.9 ꢄ 0.2
ꢁ20.8 ꢄ 0.2
ꢁ22.0 ꢄ 0.1
ꢁ21.6 ꢄ 0.2
ꢁ15.8 ꢄ 0.2
ꢁ20.4 ꢄ 0.3
ꢁ23.8 ꢄ 0.3
ꢁ21.2 ꢄ 0.2
ꢁ48.4 ꢄ 0.2
ꢁ44.5 ꢄ 0.1
ꢁ44.1 ꢄ 0.2
ꢁ45.8 ꢄ 0.2
ꢁ41.9 ꢄ 0.2
ꢁ43.8 ꢄ 0.1
ꢁ43.5 ꢄ 0.2
ꢁ42.3 ꢄ 0.2
ꢁ47.8 ꢄ 0.2
ꢁ53.4 ꢄ 0.2
ꢁ47.4 ꢄ 0.1
25.7 ꢄ 0.3
23.3 ꢄ 0.3
23.7 ꢄ 0.3
25.2 ꢄ 0.3
24.5 ꢄ 0.3
23.8 ꢄ 0.1
23.5 ꢄ 0.4
27.2 ꢄ 0.4
32.8 ꢄ 0.6
30.5 ꢄ 0.3
26.6 ꢄ 0.4
ꢁ2.5 ꢄ 0.0
ꢁ2.8 ꢄ 0.0
ꢁ2.9 ꢄ 0.0
ꢁ2.9 ꢄ 0.0
ꢁ3.1 ꢄ 0.0
ꢁ2.9 ꢄ 0.0
ꢁ2.8 ꢄ 0.0
ꢁ4.1 ꢄ 0.0
ꢁ4.3 ꢄ 0.0
ꢁ3.0 ꢄ 0.0
ꢁ2.8 ꢄ 0.0
23.2 ꢄ 0.3
20.5 ꢄ 0.3
20.8 ꢄ 0.3
22.3 ꢄ 0.3
20.5 ꢄ 0.2
20.9 ꢄ 0.1
ꢁ22.8 ꢄ 0.3
23.1 ꢄ 0.4
28.5 ꢄ 0.6
27.5 ꢄ 0.3
23.8 ꢄ 0.4
16.9 ꢄ 0.5
16.2 ꢄ 0.4
16.3 ꢄ 0.5
15.8 ꢄ 0.4
16.0 ꢄ 0.4
16.6 ꢄ 0.5
16.7 ꢄ 0.4
10.2ꢄ
ꢁ8.3
ꢁ7.8
ꢁ7.0
ꢁ7.7
ꢁ4.5
ꢁ6.3
ꢁ6.1
ꢁ9.0
ꢁ10.2
ꢁ8.5
ꢁ6.2
Podophyllotoxin
Colchicine
MDL 27048
CHALC101
9.1ꢄ
17.4ꢄ
17.4ꢄ
structure compared to the other two inhibitor main structural
parts. This relative stability for inhibitors 1e4 with respect to
1ae3a confirms that an unsubstituted benzotriazole group is
essential to maintain important van der Waals interactions with the
tubulin residues in most of the conformations that the molecules
sample in their motions within the binding site. When two chlorine
atoms replace two hydrogens on this heterocyclic scaffold, the free-
energy decomposition coupled with VVI values clearly shows that
more tubulin residues present considerable loss of van der Waals
interactions with the inhibitors. This is particularly true when the
eOCH3 substituent on moiety C is no longer present to anchor the
confirmed by competition experiments, in which compound 4
strongly displaced radio labeled colchicine from its binding site on
the tubulin, showing an IC50 value lower than that of colchicine.
Accordingly, the putative binding modes of these compounds on
human b-tubulin were predicted using a combined 3D pharmaco-
phore development/docking procedure, which was tested against
the available crystal structures of tubulin in complex with known
inhibitors. Further, the binding affinities of the lead compounds
1e4 and new molecules 1ae3a were estimated using a computa-
tional ansatz that combines MM/PBSA and normal mode calcula-
tions. Importantly, the results stemming from the free energy
calculations highlighted the dominant role of the van der Waals
molecules via the important hydrogen bond with Serb239. A
comparison of the relevant MD structure snapshots of 1, 1a, 2 and
2a in complex with tubulin, taken as proof-of-concepts and illus-
trated in Fig. 12, shown that both moieties B and C in the a series of
compounds undergo significant geometry changes within the
binding site, that ultimately lead these groups to lose some of their
van der Waals contacts with the protein binding pocket residues.
energy component ( DGbind values, as the
D
EMvdMW ) to the overall
favorable electrostatic contribution to binding was almost coun-
terbalanced by an almost equal but opposite desolvation term.
Further, the decomposition of the free energy analysis on tubulin
binding site residues indicated that, while the benzotriazole group
in all compounds sustained their van der Waals interactions with
most of the protein binding pocket, the other two moieties (i.e., the
cyano-substituted double bond and the substituted phenyl ring)
were able to sample a larger conformational space within the
tubulin binding site, adopting several conformations that lose
favorable dispersive interactions with important, anchor residues.
The combination of experimental and computational results
presented in this paper support the hypothesis that this class of
compounds can act as antiproliferative agents by interacting with
tubulin in the colchicine binding sites. Consequently, they could
inhibit tubulin assembly and might accordingly be classified as
microtubule destabilizing agents. More, the modeling study
proposes a rationale for the lower antiproliferative activity of the
compound series 1ae3a with respect to the leading series 1e4, and
7. Conclusions
The results from antiproliferative screenings, reported in
Table 2, show that no (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]tri-
azol-1-yl)-3-(4-R-phenyl)acrylonitrile derivative 1ae3a resulted
more potent than its counterpart 1e3 as so as of 4, thus proving
that the introduction of chlorine atoms on the benzo-moiety
decreases the antiproliferative activity of this class of benzo-
triazole derivatives.
Cell cycle analysis revealed that these series of compounds act
on the G2/M phase of the cell cycle, an evidence consistent with
a possible tubulin binding activity. This evidence was further
vdW
ele
Fig. 9. Decomposition of energy from MM/PBSA per tubulin binding site residue. Total molecular mechanics (
D
E
þ
D
E
) energy component of each residue with (A)
MM
MM
compounds 1 (gray), 2 (blue), 3 (red), and 4 (green), and (B) compounds 1 (gray), 1a (violet), 2a (orange), and 3a (light green). (For interpretation of the references to colour in this
figure legend, the reader is referred to the web version of this article.)

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4161
ele
D
Fig. 10. (A) Correlation between the solvation ( Gsol) and electrostatic contribution (
D
E
) of each tubulin binding site residue to the binding of compounds 1e4. Tubulin residues
MM
with van der Waals energy component larger than 1.0 kcal/mol interacting with (B) compounds 1 (gray), 2 (blue), 3 (red), and 4 (green), and (C) compounds 1 (gray), 1a (violet), 2a
(orange), and 3a (light green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
suggests that the design of new tubulin-binding agents based on
similar scaffolds should consider a careful balance between the
choice of substituents on the benzotriazole moiety A and the re-
optimization of the remaining molecular moieties B and C, as
these two may not be longer able to maintain some critical inter-
actions with tubulin by virtue of small albeit important confor-
mational changes. Further biological and computational studies on
other compounds will be reported in due course.
petroleum refers to the fraction with bp 40e60 ꢀC. The analytical
results for C, H, N and chlorine were within 0.4% of the theoretical
values.
8.2. Preparation of 5,6-dichloro-1H-benzo[d] [1,2,3]triazole (5)
Compound 5 was prepared following the procedure previously
described [32].
8. Materials and methods
8.1. Melting points, analytical and spectroscopical data
Melting points were determined by a Kofler hot stage or Digital
Electrothermal apparatus, and are uncorrected. Infrared spectra are
for Nujol mulls and were recorded using a PerkineElmer 781
spectrophotometer. UV spectra are qualitative and were recorded
in nm for solutions in EtOH with a PerkineElmer Lambda 5 spec-
trophotometer. 1H NMR spectra were recorded on a Varian XL-200
instrument, using TMS as internal standard. The chemical shift
values are reported in ppm (d) and coupling constants (J) in Hertz
(Hz). Signal multiplicities are represented by: s (singlet) and
d (doublet). MS spectra were performed on combined HP 5790eHP
5970 GC/MS apparatus. Column chromatography was performed
using 230e400 mesh silica gel (Merck silica gel 60). Light
Fig. 11. Definition of the three main chemical moieties of compounds 1e4 and 1ae3a
illustrated for compound 1 as an example.

4162
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
Table 9
stirred at room temperature, a solution of 4-methoxybenzaldehyde
(2.31 g, 17 mmol) in the same solvent (3 mL) was added dropwise.
After addition was complete, the whole mixture was heated under
reflux for 24 h. The desired compound (only E isomer), was
obtained by chromatography on silica gel column, eluenting by
light petroleum-ethyl acetate 70:30 of the crude residue obtained
after evaporation of the toluene solution. Yield 40%; M.p.
van der Waals energy variation index (VVI, %) for the three main moieties A, B, and C
of compounds 2e4 and 1ae3a.
Compd
A
B
C
2
3
4
1a
2a
3a
3
3
5
13
9
10
7
8
8
29
21
23
4
7
5
22
18
15
193e194 ꢀC (acetone); IR (Nujol):
lmax 318, 353 nm; 1H NMR (CDCl3):
n
d
2219 (CN) cmꢁ1; UV (EtOH):
8.26 (s, 1H, H-4), 8.04 (s, 1H,
H-7), 7.93 (d, J ¼ 8.4, 2H, H-20 þ H-60), 7.81 (s, 1H, vinyl-H), 7.04 (d,
J ¼ 8.4, 2H, H-30 þ H-50), 3.92 (s, 3H, OCH3); MS: m/z 344, 346, 348
(Mþ); Anal. C16H10Cl2N4O(C, H, N).
8.3. Preparation of 2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)
acetonitrile (6a) and 2-(5,6-dichloro-2H-benzo[d] [1,2,3]triazol-2-yl)
acetonitrile (6b)
8.4.2. (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)-3-p-
tolylacrylonitrile (2a)
2 g (10 mmol) of 5,6-dichloro-1H-benzo[d] [1,2,3]triazole are
dissolved in 20 ml of DMFa. To this solution are added 1.18 g
(11.7 mmol) 1.62 ml of triethylamine and 1.2 g (16 mmol) 1.02 ml of
cloroacetonitrile. The whole mixture was heated under reflux for
12 h. After cooling, the triethylamine hydrochloride formed is
filtered, and the solution is evaporated to dryness. Residue, dis-
solved in ether, is subjected to a series of washes with water.
Evaporation of the extract is obtained a residue containing the two
isomers 6a and 6b, that are separated by flash cromatografy. Elution
with light petroleum-ethyl acetate 80:20 afforded desired
compounds:2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)aceto-
To a solution of 2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)
acetonitrile (6a) (1.32 mmol) and Et3N (4.0 mmol) in toluene (5 mL)
stirred at room temperature, a solution of 4-methylbenzaldehyde
(2.04 g, 17 mmol) in the same solvent (3 mL) was added drop-
wise. After addition was complete, the whole mixture was heated
under reflux for 24 h. The desired compound (only E isomer), was
obtained by chromatography on silica gel column, eluenting by
light petroleum-ethyl acetate 80:20 of the crude residue obtained
after evaporation of the toluene solution. Yield: 20%; mp:
105e107 ꢀC (acetone); IR (Nujol):
n
2214 (CN) cmꢁ1; UV (EtOH):
8.27 (s, 1H, H-4), 8.06 (s,
lmax 221, 297, 330 nm; 1H NMR (CDCl3):
d
nitrile (6a): yield: 23%; mp: 149e150 ꢀC; IR (Nujol):nmax2240 cmꢁ1
;
1H, H-7), 7.88 (s, 1H, vinyl-H), 7.84 (d, J ¼ 8.4, 2H, H-20 þ H-60) 7.36
(d, J ¼ 8.4, 2H, H-30 þ H-50), 2.47 (s, 3H, CH3); MS: m/z 328, 330, 332
(Mþ); Anal. C16H10Cl2N4 (C, H, N).
UV(EtOH): lmax 217, 271, 292 nm; 1H NMR (CDCl3):
d 8.27 (s, 1H, H-
4), 7.85 (s, 1H, H-7), 5.58 (s, 2H, CH2); MS: m/z 226, 228, 230 (Mþ).
Anal. C8H4Cl2N4 (C, H, N). 2-(5,6-dichloro-2H-benzo[d] [1,2,3]tri-
azol-2-yl)acetonitrile (6b):yield: 15%; mp: 122e123 ꢀC; IR (Nujol):
8.4.3. (E)-3-(4-bromophenyl)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]
triazol-1-yl)acrylonitrile (3a)
To a solution of 2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)
acetonitrile (6a) (0.3g, 1.32 mmol) and Et3N (4.0 mmol) in toluene
nmax 2240 cmꢁ1, UV(EtOH) 219, 290, 302 nm; 1H NMR (CDCl3):
d
8.06 (s, 2H, H-4 e H-7), 5.64 (s, 2H, CH2); MS: m/z 226, 228, 230
(Mþ). Anal. C8H4Cl2N4 (C, H, N).
(5 mL) stirred at room temperature,
a
solution of 4-
8.4. Preparation of (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-
1-yl)-3-(4-R-phenyl) acrylonitriles (1a, 2a, 3a)
bromobenzaldehyde (3.14 g, 17 mmol) in the same solvent (3 mL)
was added dropwise. After addition was complete, the whole
mixture was heated under reflux for 24 h. The desired compound
was obtained by filtration of the resulting precipitates as soon as
the reaction mixture reaches r.t. Yield 33%; m.p. 150e152 ꢀC
8.4.1. (E)-2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)-3-(4-
methoxyphenyl)acrylonitrile (1a)
To a solution of 2-(5,6-dichloro-1H-benzo[d] [1,2,3]triazol-1-yl)
acetonitrile (6a) (1.32 mmol) and Et3N (4.0 mmol) in toluene (5 mL)
(acetone); IR (Nujol):
n
2225 (CN) cmꢁ1; UV (EtOH): lmax 216,
8.36 (s,1H, H-4), 8.27 (s,1H, H-7), 8.14 (s,
297 nm; 1H NMR (CDCl3):
d
Fig. 12. Conformational space of (A) compound 1, (B) 1a, (C) 2, and (D) 2a sampled in their complex with tubulin by MD simulations. For each compound, only 10 representative
snapshots taken along the corresponding MD trajectories are shown in different colors for clarity.

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4163
1H, vinyl-H), 7.92 (d, J ¼ 8.2, 2H, H-20 þ H-60) 7.72 (d, J ¼ 8.6, 2H, H-
30 þ H-50), 3.27 (s, 3H, CH3); MS: m/z 392, 394, 396, 397, 399 (Mþ);
Anal. C15H7BrCl2N4(C, H, N).
Counter (PE-Wallac Boston, MA). IC50 values were calculated via
data non-linear regression using GraphPad Prism 5.
8.7. Molecular modeling
8.5. Cell-based assay
The extensive, parallel molecular dynamics analyses were per-
formed using 64 processors of the Tartaglia cluster at the University
of Trieste (Trieste, Italy), as well as the same number of processors
on the IBM/BCX cluster at the CINECA supercomputer center
(Bologna, Italy). The Catalyst module of the Discovery Studio plat-
form [37] was employed for 3D pharmacophore model develop-
ment. The software AutoDock [51] was used for ligand docking. The
parallel version of the Amber 9 suite of programs [52] was
employed in all molecular dynamics simulations. Molecular
graphics images were produced using the UCSF Chimera package
from the Resource for Biocomputing, Visualization, and Informatics
at the University of California, San Francisco (supported by NIH P41
RR-01081) [53].
Test compounds were dissolved in DMSO at 100 mM and then
diluted into culture medium.
8.5.1. Cells
Cell lines were purchased from American Type Culture Collec-
tion (ATCC). Hematological tumor-derived cell lines were grown in
RPMI-1640 medium supplemented with 10% FCS, 100 units/mL
penicillin G and 100mg/mLstreptomycin. Solid tumor-derived cell
lines were grown in their specific media supplemented with 10%
FCS and antibiotics. Cell cultures were incubated at 37 ꢀC in
a humidified, 5% CO2 atmosphere. The absence of mycoplasma
contamination was checked periodically by the Hoechst staining
method.
8.7.1. 3D pharmacophore model generation
8.5.2. Antiproliferative assays
For the automated pharmacophore generation with DS Catalyst,
a training set of 25 inhibitors of tubulin, with CC50 values between
0.05 and 20,000 mM, was defined (Table 5). The test set used for the
Exponentially growing cells derived from human hematological
tumors [CD4þ human T-cells containing an integrated HTLV-1
genome (MT-4); CD4þ human acute T-lymphoblastic leukemia
(CCRF-CEM); Human splenic B-lymphoblastoid cells (WIL-2NS);
Human acute B-lymphoblastic leukemia (CCRF-SB)] were seeded at
an initial density of 1 ꢃ105 cells/mL in 96 well plates in RPMI-1640
medium supplemented with 10% fetal calf serum (FCS), 100 units/
validation of the pharmacophore model consisted of further 11
derivatives, as reported in Table 6. The details of the pharmaco-
phore development procedure with Catalyst have been extensively
described in previous work [49,54e57]. Briefly, the model struc-
tures of all compounds were sketched and each molecular structure
was subjected to energy minimization using the generalized
CHARMM force field [58] until the gradient dropped below 0.05. A
conformational search was then carried out using the Poling algo-
rithm [59e61] and the CHARMM force field as implemented in the
Catalyst program. The “best quality” generation option was adopted
to select representative conformers over a 0e20 kcal/mol interval
above the computed global energy minimum in the conformational
space, and the number of conformers generated for each compound
was limited to a maximum of 250 as a good compromise between
speed and maximum coverage in the conformational space.
Based on the conformations for each compound, the HypoGen
module of DS Catalyst was used to generate three-dimensional
pharmacophore models. During hypotheses generation, the soft-
ware attempts to minimize a cost function of two main terms: the
first penalizes the deviation between the estimated affinities of the
training set molecules and their experimental values, whilst the
second penalizes the complexity of the hypothesis. The uncertain
factor for each compound represents the ratio range of uncertainty
in the affinity value based on the expected statistical straggling of
biological data collection. Uncertainty influences the first step e
also called the constructive phase e of the hypothesis generating
process. In this work, the default factor of 3.0 was selected, as the
experimental affinities of our compounds span more than the
required four orders of magnitude.
mL penicillin G and 100
mg/mL streptomycin. Human cell lines
derived from solid tumors [skin melanoma (SK-MEL-28); breast
adenocarcinoma (MCF-7); lung squamous carcinoma (SK-MES-1);
hepatocellular carcinoma (HepG-2) and prostate carcinoma (DU-
145)] were also seeded at 1x105 cells/mL in 96 well plates in specific
media supplemented with 10% FCS and antibiotics as above. Cell
cultures were then incubated at 37 ꢀC in a humidified, 5% CO2
atmosphere in the absence or presence of serial dilutions of test
compounds. Cell viability was determined after 96 h at 37 ꢀC by the
3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-tetrazolium
(MTT) method [50].
bromide
8.5.3. Analysis of cell cycle
Cells were washed once in ice-cold PBS and resuspended at
1 ꢃ 106 ml in a hypotonic fluorochrome solution containing pro-
pidium iodide (Sigma) 50 mg/ml in 0.1% sodium citrate plus 0.03%
(v/v) nonidet P-40 (Sigma). After 30 min of incubation the fluo-
rescence of each sample was analyzed as single-parameter
frequency histograms using a FACS can flow cytometer (Bec-
toneDickinson, San Jose, CA). Distribution of cells in cell cycle was
determined using the ModFit LT program (Verity Software House,
Inc.). Apoptosis was determined by evaluating the percentage of
hypoploid nuclei accumulated in the sub-G0eG1 peak after
labeling with propidium iodide.
As the DS Catalyst software can generate pharmacophore
hypotheses consisting of a maximum of five hypotheses, an initial
analysis revealed that chemical feature types such as hydrogen
bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic
aromatic (HYAr), ring aromatic (RA), and general hydrophobic (HY)
sites could effectively map all critical chemical features of all
molecules in the training and test sets. Accordingly, these five
feature types were used to generate 10 pharmacophores from the
training set.
The HypoGen module in DS Catalyst performs two important cost
calculations (represented in bit units) that determine the success of
any pharmacophore hypothesis. One is called the fixed cost, which
represents the simplest model that fits all data perfectly, while the
8.6. [3H]Colchicine competition-binding scintillation proximity
assay (SPA) [34]
[3H]Colchicine was diluted, and biotin-labeled tubulin (T333,
Cytoskeleton, Denver, CO) was reconstituted following the manu-
facturer protocol. The diluted compound 4 and [3H]colchicine were
transferred to a 96-well isoplate (PE-Wallac, Boston, MA), to which
the binding buffer and the reconstituted biotinylated tubulin were
added. The plates were incubated at 37 ꢀC for 2 h; then,
streptavidin-coated yttrium SPA beads (Amersham Pharmacia
Biotech, Piscataway, NJ) were added, and the bound radioactivity
was determined using a MicroBetaTriluxMicroplate Scintillation

4164
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
second is known as the null cost, and represents the highest cost of
a pharmacophore with no features and which estimates activity to
be the average of the activity data of the training set molecules. A
meaningful pharmacophore hypothesis may result when the
difference between these two values is large; for instance, a value of
40e60 bits for a 3D pharmacophore hypothesis may indicate that it
has 75e90% probability of correlating the data. Also, the total cost of
any pharmacophore hypothesis should be close to the fixed cost to
provide any useful model. A further parameter that also determine
the quality of a given 3D pharmacophore model with possible
predicting value is the configurational cost, also known as the
entropy cost and depends on the complexity of the pharmacophore
hypothesis space. For a good hypothesis, this cost, that is the
magnitude of the hypothesis space for a given training set of
compounds, should be less than 17. If this last cost exceeds 17, there
are more degrees of freedom in the training set that the DS Catalyst
algorithm can properly handle and, consequently, the correspond-
ing pharmacophore is likely to be poorly meaningful. Finally, the
root-mean-square deviations (RMSDs) and the correlation coeffi-
of homology-based techniques [62]. Briefly, we used the Discovery
Studio software (v.2.1, Accelrys Inc., San Diego, CA, USA) to generate
a sequence alignment among model protein (human
SwissProt accession code Q9H4B7) and three reference proteins
(i.e., Bos Taurus, Sus Scrofa, and Ovis Aries tubulin -chains, PDB
b-tubulin,
b
entry codes 1SA1, 1IAO, and 3HKC, respectively). Subsequently, the
Modeler program [63,64] was employed to build a preliminary 3D
model of human tubulin b-chain based on the alignment obtained
in the previous step. The model generated was further refined using
statistical pair potentials for the loops. The geometry of protein was
refined for 200 steps (steepest descent) in vacuum, using the all-
atom the ff03 force field of Duan et al. [65,66] as defined in
parm99.dat and frcmod.ff03 parameter files [67]. Further protein
geometry refinement was carried out using the PMEMD module of
Amber 9 via a combined steepest descent e conjugate gradient
algorithm, using as a convergence criterion for the energy gradient
the root-mean-square of the Cartesian elements of the gradient
equal to 10ꢁ2 kcal/(mol Å). The quality of the model was assessed by
using different validation tools such as Procheck [68] and Whatif
[69]. Ramachandran plot statistics indicated that 98% of the main-
chain dihedral angles were found in the most favorable region,
thus confirming the exceptional quality of the 3D model of human
cients r represent de facto the quality of the correlation between the
estimated and the actual activity data.
Three validation procedures were followed to determine the
statistical relevance and the validity of the proposed 3D pharma-
cophore models: the test set prediction model, the randomization
method, and the leave-one-out procedure. In this work, the former
procedure consisted in the collection of further, different
compounds into a test set, and in performing a regression analysis
by mapping the test set molecules onto the best pharmacophore
hypothesis. The high correlation coefficients obtained using the test
set compounds revealed the good correlation between the actual
and estimated affinities and, hence, the predictive validity of the
corresponding 3D hypothesis. The randomization validation
procedure is based on Fisher’s randomization test [39]. The goal of
this type of validation is to check whether there is a strong corre-
lation between the chemical structures and the biological activity.
This is done by randomizing the affinity data associated with the
training set compounds, generating pharmacophore hypothesis
using the same features and parameters employed to develop the
original pharmacophore model. The statistical significance is
calculated according to the following formula:
b
-tubulin obtained. This optimized 3D model was then used as the
entry point for molecular dynamics (MD) simulations.
The model structures of the selected ligands were generated with
Discovery Studio. All molecules were subjected to an initial energy
minimization, withtheconvergencecriterionsetto10ꢁ4 kcal/(molÅ).
A conformational search was carried out using a well-validated, ad-
hoc developed combined molecular mechanics/molecular dynamics
simulated annealing (MDSA) protocol [45,46,62,70]. Accordingly, the
relaxed structures were subjected to 5 repeated temperature cycles
(from 310 K to 1000 K and back) using constant volume/constant
temperature (NVT) MD conditions. At the end of each annealing cycle,
the structures were again energy minimized to converge below
10ꢁ4 kcal/(mol Å), and only the structures corresponding to the
minimum energy were used for further modeling. The atomic partial
charges for the geometrically optimized compounds were obtained
using the RESP procedure [71], and the electrostatic potentials were
produced by single-point quantum mechanical calculations at the
Hartree-Fock level with a 6-31G* basis set, using the Merz-Singh-
Kollman van der Waals parameters [72,73]. Eventual missing force
field parameters for the inhibitor molecules were generated using the
Antichamber module of Amber 9.0.
h
i
~
Significance ¼ 100 ꢃ 1ð1 þ x=yÞ
where x is the total number of hypotheses having a total cost lower
than the original (best) hypothesis, and y is the total number of
HypoGen runs (initial þ random runs). Thus, for instance, 99
random spreadsheets (i.e.,99 HypoGen runs) have to be generated
to obtain a 99% confidence level. Should any randomized data set
result in the generation of a 3D pharmacophore with similar or
even better cost values, root-mean-square deviations, and corre-
lation coefficients, then it is likely that the original hypothesis does
reflect a chance correlation.
The optimized structures of the test compounds were docked
into the colchicin - podophyllotoxin binding site of tubulin by
applying a consolidated procedure [74], accordingly, it will be
described here only briefly. The software AutoDock4 [51] was
employed to estimate the possible binding orientations of all
compounds in the receptor. In order to encase a reasonable region
of the protein surface and interior volume, centered on the crys-
tallographic identified binding site (as determined by super posi-
tion with the highly homolog Bos Taurus b-tubulin in complex with
Finally, the leave-one-out test checks if the correlation between
experimental and computed affinities is heavily dependent on one
particular molecule of the training set by re-computing the phar-
macophore model with the exclusion of one molecule at a time.
Accordingly, 24 new training sets were built, each composed by 23
molecules, and 24 HypoGen calculations were launched under the
same conditions. For each run, the hypothesis characterized by the
lowest total cost was employed to predict the affinityof the excluded
compound and to estimate the new correlation coefficient.
colchicines, PDB entry 1SA0.pdb), the grids were 60 Å on each side.
Grid spacing (0.375 Å), and 120 grid points were applied in each
Cartesian direction so as to calculate mass-centered grid maps.
Amber 12e6 and 12e10 LennardeJones parameters were used in
modeling van der Waals interactions and hydrogen bonding (NeH,
OeH and SeH), respectively. In the generation of the electrostatic
grid maps, the distance-dependent relative permittivity of Mehler
and Solmajer was applied [75].
For the docking of each compound to the protein, three hundred
Monte Carlo/Simulated Annealing (MC/SA) runs were performed,
with 100 constant temperature cycles for simulated annealing. For
these calculations, the GB/SA implicit water model [76,77] was used
8.7.2. Molecular docking studies
The 3D model structure of human
b-tubulin, presently not
available in the Protein Data Bank (PDB), was built by a combination
to mimic the solvated environment. The rotation of the angles
4

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4165
and
f, and the angles of side chains were set free during the
D
Gbind
¼
D
D
EMM
þ
D
Gsolv ꢁ T
D
S
(2)
calculations. All other parameters of the MC/SA algorithm were
kept as default. Following the docking procedure, the structure of
all compounds were subjected to cluster analysis with a tolerance
of 1 Å for an all-atom root-mean-square (RMS) deviation from
a lower-energy structure representing each cluster family. In the
absence of any relevant crystallographic information, the structure
of each resulting complex characterized by the lowest interaction
energy in the prevailing cluster was selected for further evaluation.
Each best substrate/tubulin complex resulting from the auto-
mated docking procedure was further refined in Amber 9 using the
quenched molecular dynamics (QMD) method [70,78,79]. In this
case, 1 ns molecular dynamics (MD) simulation at 310 K were
employed to sample the conformational space of the sub-
strateetubulin complex in the GB/SA continuum solvation envi-
ronment [76,77]. The integration step was equal to 1 fs. After each
ps, the system was cooled to 0 K, the structure was extensively
minimized, and stored. To prevent global conformational changes
of the protein, the backbone of the protein binding site were con-
strained by a harmonic force constant of 100 kcal/Å, whereas the
amino acid side chains and the ligands were allowed moving
without any constraint.
where
Gbind is the binding free energy in water,
D
EMM is the
DGsol is the
interaction energy between the ligand and the protein,
solvation free energy, and ꢁT S is the conformational entropy
contribution to the binding. EMM is calculated from the molecular
mechanics (MM) interaction energies, according to:
D
D
ele
vdw
MM
D
EMM
¼
D
E
þ
D
E
(3)
MM
ele
ele
vdw
where
D
E
and
D
E
are electrostatic and van der Waals inter-
MM
action energies between the ligand and the receptor, which are
calculated using the MM/PBSA module in the Amber 9 software suite.
The solvation energy,
trostatic contributions,
DGsol, is divided into two parts, the elec-
np
sol
D
Gesolel, and all other contributions,
D
G
:
np
sol
ele
sol
D
Gsol
¼
D
G
þ
D
G
(4)
ele
sol
The electrostatic contribution to the solvation free energy,
D
G
,
is calculated using the DelPhi software package [85], which solves
the PoissoneBoltzmann equations numerically and calculates the
electrostatic energy according to the electrostatic potential. Interior
and exterior dielectric constant values e were set equal to 1 and 80,
The best energyconfiguration of each complex resulting from the
previous step was allowed to relax in an 80 Å ꢃ 80 Å ꢃ 80 Å box of
TIP3P water molecules [80]. The resulting system was minimized
with a gradual decrease in the position restraints of the protein
atoms. Finally, to achieve electroneutrality, a suitable number of
counterions were added, in the positions of largest electrostatic
potential, as determined by the module LEaP module within Amber
9. After energy minimization of the added ions for 1500 steps,
keeping the protein, the ligand, and the pre-existing waters rigid,
followed by an MD equilibration of the entire water/ion box with
fixed solute for 2 ns, further unfavorable interactions within the
structures were relieved by progressively smaller positional
restraints on the solute (from 25 to 0 kcal/(mol Å2) for a total of 4 ns.
Each hydrated complex system was gradually heated to 310 K in
three intervals, allowing a 2 ns interval per each 100 K, and then
equilibrated for 5 ns at 310 K, followed by 20 ns of data collection
runs, necessary for the estimation of the free energy of binding (vide
infra). The MD simulations were performed at constant T ¼ 310 K
using the Berendsen et al. coupling algorithm [81] with separate
coupling of the solute and solvent to the heat, an integration time
step of 2 fs, and the applications of the Shake algorithm [82] to
constrain all bonds to their equilibrium values, thus removing high
frequency vibrations. Long-range nonbonded van der Waals inter-
actions were truncated by using a dual cut-off of 9 and 13 Å,
respectively, where energies and forces due to interactions between
9 and 13 Å were updated every 20 time steps. The particle mesh
Ewald method [83] was used to treat the long-range electrostatics.
For the calculation of the binding free energy between tubulin and
each inhibitor in water, a total of 20,000 snapshots were saved
during the MD data collection period described above.
respectively. A grid spacing of 2/Å, extending 20% beyond the
dimensions of the solute, was employed. The non-polar component
np
D
D
b
G
was obtained using the following relationship [86]:
ꢃ SA þ in which
0.00542 kcal/(mol Å2),
sol
np
G
¼
g
b
,
g
¼
sol
¼ 0.92 kcal/mol, and the surface area SA was estimated by means
of the MSMS software [87]. The last parameter in Equation (1), i.e.
the change in solute entropy upon association eT S, was calculated
D
through normal mode analysis [88], using the N mode module of
Amber 9. In the first step of this calculation, an 8 Å sphere around the
ligand was cut out from an MD snapshot for each ligandeprotein
complex. This value was shown to be large enough to yield
converged mean changes in solute entropy. On the basis of the size-
reduced snapshots of the complex, we generated structures of the
uncomplexed reactants by removing the atoms of the protein and
ligand, respectively. Each of those structures was minimized, using
a distance-dependent dielectric constant
e
¼ 4r, to account for
solvent screening, and its entropy was calculated using classical
statistical formulas and normal mode-analysis. To minimize the
effects due to different conformations adopted by individual snap-
shots we averaged the estimation of entropy over 40 snapshots.
Finally, the IC50 values were calculated from the corresponding
binding free energies using the following relationship [41,42]:
D
Gbind ¼ RT lnKdiss ¼ RT lnðIC50 þ0:50CenzÞyRT lnIC50
(5)
The overall quality of the entire procedure described above e
i.e., protein and inhibitors modeling, conformational analysis,
docking and energetic calculations e was tested by carrying out the
same calculations on podophyllotoxin, colchicines, MDL 27048, and
CHALC101, for which the crystallographic structures of the relevant
tubulin complexes and/or the corresponding IC50 values were
available [40,41,45].
8.7.3. Free energy of binding
The binding free energy DGbind of each tubulin/drug complex in
water was calculated according to the procedure termed Molecular
Mechanic/PoissoneBoltzmann Surface Area (MM/PBSA), and orig-
inally proposed by Srinivasan et al. [40]. Since the theoretical
background of this methodology is described in details in the
original papers by Peter Kollman and his group [84], it will be only
briefly described below.
References
[1] B.G. Szczepankiewicz, G. Liu, H.S. Jae, A.S. Tasker, I.W. Gunawardana, T.W. von
Geldern, S.L. Gwaltney 2nd, J.R. Wu-Wong, L. Gehrke, W.J. Chiou, R.B. Credo,
J.D. Alder, M.A. Nukkala, N.A. Zielinski, K. Jarvis, K.W. Mollison, D.J. Frost,
J.L. Bauch, Y.H. Hui, A.K. Claiborne, Q. Li, S.H. Rosenberg, J. Med. Chem. 44
(2001) 4416e4430.
[2] O. Valiron, N. Caudron, D. Job, Cell. Mol. Life Sci. 58 (2001) 2069e2084.
[3] S. Sengupta, S.A. Thomas, Expert Rev. Anticancer Ther. 6 (2006) 1433e1447.
[4] Y.K. Chiang, C.C. Kuo, Y.S. Wu, C.T. Chen, M.S. Coumar, J.S. Wu, H.P. Hsieh,
C.Y. Chang, H.Y. Jseng, M.H. Wu, J.S. Leou, J.S. Song, J.Y. Chang, P.C. Lyu,
Y.S. Chao, S.Y. Wu, J. Med. Chem. 52 (2009) 4221e4233.
Basically, an MD simulation (typically in explicit solvent) is first
carried out which yields a representative ensemble of structures.
The average total free energy of binding between each drug and the
protein receptor can then be calculated as:

4166
A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
[5] A.R. Katritzky, Farmaco 43 (1988) 1175e1181.
[48] R.B. Ravelli, B. Gigant, P.A. Curmi, I. Jourdain, S. Lachkar, A. Sobel, M. Knossow,
Nature 428 (2004) 198e202.
[6] A.R. Katritzky, S. Rachwal, G.J. Hitchings, Tetraedron 47 (1991) 2683e2732.
[7] A.R. Katritzky, J. Jiang, J.V. Greenhill, P.J. Steel, J. Chem. Soc. Perkin Trans. 1
(1992) 3055e3060.
[8] F. Sparatore, F. Pagani, Farmaco Ed. Sci. 17 (1962) 414e429.
[9] F. Sparatore, F. Pagani, Farmaco Ed. Sci. 20 (1965) 248e258.
[10] F. Sparatore, M.I. La Rotonda, G. Paglietti, E. Ramundo, C. Silipo, A. Vittoria,
Farmaco Ed. Sci. 33 (1978) 901e922.
[49] G. Giliberti, C. Ibba, E. Marongiu, R. Loddo, M. Tonelli, V. Boido, E. Laurini,
P. Posocco, F. Fermeglia, S. Pricl, Bioorg. Med. Chem. 18 (2010) 6055e6068.
[50] R. Pauwels, J. Balzarini, M. Baba, R. Snoeck, D. Schols, P. Herdewijn,
J. Desmyster, E. De Clercq, J. Virol. Methods 20 (1998) 309e321.
[51] G.M. Morris, D.S. Goodsell, R.S. Halliday, R. Huey, W.E. Hart, R.K. Belew,
A.J. Olson, J. Comput. Chem. 19 (1998) 1639e1662.
[11] F. Sparatore, F. Pagani, G. Caliendo, E. Novellino, C. Silipo, A. Vittoria, Farmaco
Ed. Sci. 43 (1988) 141e151.
[12] A. Boido, I. Vazzana, F. Sparatore, B. Cappello, V. Diurno, C. Matera,
R. Marrazzo, S. Russo, M.L. Cenicola, E. Marmo, Farmaco 44 (1989) 257e277.
[13] M. Loriga, G. Paglietti, F. Sparatore, G. Pinna, A. Sisini, Farmaco 47 (1992)
439e448.
[52] D.A. Case, T.A. Darden, T.E. Cheatham III, C.L. Simmerling, J. Wang, R.E. Duke,
R. Luo, K.M. Merz, D.A. Pearlman, M. Crowley, R.C. Walker, W. Zhang, B. Wang,
S. Hayik, A. Roitberg, G. Seabra, K.F. Wong, F. Paesani, X. Wu, S. Brozell, V. Tsui,
H. Gohlke, L. Yang, C. Tan, J. Mongan, V. Hornak, G. Cui, P. Beroza,
D.H. Mathews, C. Schafmeister, W.S. Ross, P.A. Kollman, AMBER 9. University
of California, San Francisco, CA, USA, 2006.
[14] G. Paglietti, P. Sanna, A. Carta, F. Sparatore, I. Vazzana, A. Peana, M. Satta,
Farmaco 49 (1994) 693e702.
[53] (a) Chimera (v.1.4), Resource for Biocomputing, Visualization, and Informatics
at the University of California, San Francisco (CA), USA.
[15] I. Vazzana, A. Boido, F. Sparatore, R. Di Carlo, G. Mattace Raso, M. Pacilio,
Farmaco 52 (1997) 131e139.
[16] G. Paglietti, F. Sparatore, Farmaco Ed. Sci. 27 (1972) 380e396.
[17] G. Paglietti, F. Sparatore, Ann. Chim. 62 (1972) 128e143.
[18] A. Carta, P. Sanna, G. Paglietti, F. Sparatore, Heterocycles 55 (2001)
1133e1140.
[19] P. Sanna, A. Carta, G. Paglietti, S. Zanetti, G. Fadda, Farmaco 47 (1992)
1001e1019.
[20] P. Sanna, A. Carta, G. Paglietti, A. Bacchi, G. Pelizzi, J. Heterocyclic Chem. 34
(1997) 97e106.
[21] A. Carta, P. Sanna, H. Molinari, Heterocycles 45 (1997) 1391e1395.
[22] P. Sanna, A. Carta, G. Paglietti, Heterocycles 51 (1999) 2171e2181.
[23] A. Carta, P. Sanna, A. Bacchi, Heterocycles 57 (2002) 1079e1090.
[24] P. Sanna, A. Carta, M.E. Rahbar Nikookar, Eur. J. Med. Chem. 35 (2000) 535e543.
[25] P. Sanna, A. Carta, L. Gherardini, M.E. Rahbar Nikookar, Farmaco 57 (2002)
79e87.
[26] A. Carta, P. Sanna, M. Palomba, L. Vargiu, M. La Colla, R. Loddo, Eur. J. Med.
Chem. 37 (2002) 891e900.
[27] A. Carta, M. Palomba, G. Boatto, B. Busonera, M. Mureddu, R. Loddo, Farmaco
59 (2004) 637e644.
[28] Y. Schneider, P. Chabert, J. Stutzmann, D. Coelho, A. Fougerousse, F. Gossé,
J.F. Launay, R. Brouillard, F. Raul, Int. J. Cancer 107 (2003) 189e196.
[29] R. Romagnoli, P.G. Baraldi, M.G. Pavani, M.A. Tabrizi, D. Preti, F. Fruttarolo,
L. Piccagli, M.K. Jung, E. Hamel, M. Borgatti, R. Gambari, J. Med. Chem. 49
(2006) 3906e3915.
[30] R. Romagnoli, P.G. Baraldi, M.G. Pavani, M.A. Tabrizi, D. Preti, F. Fruttarolo,
L. Piccagli, M.K. Jung, E. Hamel, M. Borgatti, R. Gambari, J. Med. Chem. 49
(2006) 6425e6428.
[31] M. Belleri, D. Ribatti, S. Vicoli, F. Cotelli, L. Forti, V. Tannini, L.A. Stivala,
M. Presta, Mol. Pharmacol. 67 (2005) 1451e1459.
[32] A. Carta, M. Loriga, S. Piras, G. Paglietti, M. Ferrone, M. Fermeglia, S. Pricl, P. La
Colla, G. Collu, T. Sanna., R. Loddo, Med. Chem. 3 (2007) 520e532.
[33] Z. Xiao, K.F. Bastow, J.R. Vance, R.S. Sidwell, H.K. Wang, M.S. Chen, Q. Shi,
K.H. Lee, J. Med. Chem. 47 (2004) 5140e5148.
(b) E.F. Pettersen, T.D. Goddard, C.C. Huang, G.S. Couch, D.M. Greenblatt.,
E.C. Meng, T.E. Ferrin, J. Comput. Chem. 25 (2004) 1605e1612.
[54] E. Laurini, D. Zampieri, M.G. Mamolo, L. Vio, C. Zanette, C. Florio, P. Posocco,
M. Fermeglia, S. Pricl, Bioorg. Med. Chem. Lett. 20 (2010) 2954e2957.
[55] D. Zampieri, M.G. Mamolo, E. Laurini, C. Florio, C. Zanette, M. Fermeglia,
P. Posocco, M.S. Paneni, S. Pricl, L. Vio, J. Med. Chem. 52 (2009) 5380e5392.
[56] M. Tonelli, V. Boido, C. Canu, A. Sparatore, F. Sparatore, M.S. Paneni,
M. Fermeglia, S. Pricl, P. La Colla, L. Casula, C. Ibba, D. Collu, R. Loddo, Bioorg.
Med. Chem. 16 (2008) 8447e8466.
[57] R. Di Santo, M. Fermeglia, M. Ferrone, M.S. Paneni, R. Costi, M. Artico, A. Roux,
M. Gabriele, K.D. Tardif, A. Siddiqui, S. Pricl, J. Med. Chem. 48 (2005) 6304e6314.
[58] B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan,
M.J. Karplus, J. Comput. Chem. 4 (1983) 187e217.
[59] A. Smellie, S.L. Teig, P. Towbin, J. Comput. Chem. 16 (1994) 171e187.
[60] A. Smellie, S.D. Kahn, S.L. Teig, J. Chem. Inf. Comput. Sci. 35 (1995) 285e294.
[61] A. Smellie, S.D. Kahn, S.L. Teig, J. Chem. Inf. Comput. Sci. 35 (1995) 295e304.
[62] A. Carta, G. Loriga, S. Piras, G. Paglietti, M. Ferrone, M. Fermeglia, S. Pricl, P. La
Colla, B. Secci, G. Collu, R. Loddo, Med. Chem. 2 (2006) 577e589.
[63] A. Sali, T.L. Blundell, J. Mol. Biol. 234 (1993) 779e815.
[64] A. Sali, L. Potterton, F. Yuan, H. van Vlijmen, M. Karplus, Proteins 23 (1995)
318e326.
[65] Y. Duan, C. Wu, S. Chowdhury, M.C. Lee, G. Xiong, W. Zhang, R. Yang,
P. Cieplak, R. Luo, T. Lee, J. Comput. Chem. 24 (2003) 1999e2012.
[66] M.C. Lee, Y. Duan, Proteins 55 (2004) 620e634.
[67] J. Wang, P. Cieplak, P.A. Kollman, J. Comput. Chem. 21 (2000) 1049e1074.
[68] R.A. Laskowski, M.W. MacArthur, D.S. Moss, J.M. Thornton, J. Appl. Crystallogr.
26 (1993) 283e291.
[69] R. Rhodriguez, G. Chinea, N. Lopez,T. Pons, G. Vriend, CABIOS14(1998)523e528.
[70] For the application of the MDSA protocol to different bio(macro)molecular
systems, see for instance the selected list of recent publications and references
therein: (a) D. Zampieri, M.G. Mamolo, L. Vio, E. Banfi, G. Scialino,
M. Fermeglia, M. Ferrone, S. Pricl, Bioorg. Med. Chem. 15 (2007) 7444e7458;
(b) R. Toth, M. Ferrone, S. Miertus, E. Chiellini, M. Fermeglia, S. Pricl, Bio-
macromolecules 7 (2006) 1714e1719.
[34] S.K. Tahir, P. Kovar, S.H. Rosenberg, Ng Shi-Chung, BioTechniques 29 (2000)
156e160.
[35] R. Bai, M.C. Roach, S.K. Jayaram, J. Barkoczy, G.R. Pettit, R.F. Luduena, E. Hamel,
Biochem. Pharmacol. 45 (1993) 1503e1515.
[36] R. Bai, G.R. Pettit, E.J. Hamel, Biol. Chem. 265 (1990) 17141e17149.
[37] Discovery Studio v. 3.0, Accelrys, San Diego, CA, USA.
[38] J. Sutter, O.F. Guner, R.D. Hoffman, H. Li, M. Wadman, in: O.F. Guner (Ed.),
Pharmacophore Perception, Development, and Use in Drug Design, Interna-
tional University Line, La Jolla, CA, 1999, pp. 501e511.
[39] R. Fisher, The Design of Experiments. Hafner Publishing, New York, NY, USA,
1966.
[40] J. Srinivasan, T.E. Cheatham, P. Cieplak, P.A. Kollman, D.A. Case, J. Am. Chem.
Soc. 120 (1998) 9401e9409.
(c) E. Banfi, G. Scialino, D. Zampieri, M.G. Mamolo, L. Vio, M. Ferrone,
M. Fermeglia, M.S. Paneni, S.J. Pricl, Antimicrob. Agents Chemother. 58 (2006)
76e84.
[71] C.I. Bayly, P. Cieplak, W.D. Cornell, P.A. Kollman, C.I. Bayly, P. Cieplak,
W.D. Cornell, P.A. Kollman, J. Phys. Chem. 97 (1993) 10269e10280.
[72] U.C. Singh, P.A. Kollman, J. Comput. Chem. 5 (1984) 129e145.
[73] B.H. Besler, K.M. Merz, P.A. Kollman, J. Comput. Chem. 11 (1990) 431e439.
[74] For a further list of recent applications of the MM/PBSA methodology from our
group see, for instance: (a) M. Tonelli, I. Vazzana, B. Tasso, V. Boido,
F. Sparatore, M. Fermeglia, M.S. Paneni, P. Posocco, S. Pricl, P. La Colla, C. Ibba,
B. Secci, G. Collu, R. Loddo, Bioorg. Med. Chem. 17 (2009) 4425e4440;
(b) J.C. McAuliffe, W.L. Wang, G.M. Pavan, S. Pricl, D. Yang, S.S. Chen, A.J. Lazar,
R.E. Pollock, J.C. Trent, Mol. Oncol. 2 (2008) 161e163;
[41] M.B. Kroeger Smith, M.L. Lamb, J. Tirado-Rives, W.L. Jorgensen, C.J. Michejda,
S.K. Ruby, R.H. Smith Jr., Protein Eng. 13 (2000) 413e421.
[42] J. Wang, P. Morin, W. Wang, P.A. Kollman, J. Am. Chem. Soc. 123 (2001)
5221e5230.
[43] V. Peyrot, D. Leynardier, M. Sarrazin, C. Briand, M. Menendez, J. Laynez,
J.M. Andreu, Biochemistry 31 (1992) 11125e11132.
(c) T. Negri, G.M. Pavan, E. Virdis, A. Greco, M. Fermeglia, M. Sandri, S. Pricl,
M.A. Pierotti, S. Pilotti, E. Tamborini, J. Natl. Cancer Inst. 101 (2009) 194e204;
(d) M. Ferrone, F. Perrone, E. Tamborini, M.S. Paneni, M. Fermeglia, S. Suardi,
E. Pastore, D. Delia, M.A. Pierotti, S. Pricl, S. Pilotti, Mol. Cancer Ther. 5 (2006)
1467e1473;
(e) E. Tamborini, S. Pricl, T. Negri, M.S. Lagonigro, F. Miselli, A. Greco,
A. Gronchi, P.G. Casali, M. Ferrone, M. Fermeglia, A. Carbone, M.A. Pierotti,
S. Pilotti, Oncogene 25 (2006) 6140e6146;
[44] D.Y. Kim, K.-H. Kim, N.D. Kim, K.Y. Lee, C.K. Han, J.H. Yoon, S.K. Moon, S.S. Lee,
B.L. Seong, J. Med. Chem. 49 (2006) 5664e5670.
[45] For recent applications of our computational protocol see, for instance; (a)
D. Zampieri, M.G. Mamolo, E. Laurini, M. Fermeglia, P. Posocco, S. Pricl,
E. Banfi, G. Scialino, L. Vio, Bioorg. Med. Chem. 17 (2009) 4693e4707;
(b) M. Mazzei, E. Nieddu, M. Miele, A. Balbi, M. Ferrone, M. Fermeglia,
M.T. Mazzei, S. Pricl, P. La Colla, F. Marongiu, C. Ibba, R. Loddo, Bioorg. Med.
Chem. 16 (2008) 2591e2605;
(f) S. Pricl, M. Fermeglia, M. Ferrone, E. Tamborini, Mol. Cancer Ther. 4 (2005)
1167e1174;
(g) M.G. Mamolo, D. Zampieri, L. Vio, M. Fermeglia, M. Ferrone, S. Pricl,
G. Scialino, E. Banfi, Bioorg. Med. Chem. 13 (2005) 3797e3809.
[75] E.L. Mehler, T. Solmajer, Protein Eng. 4 (1991) 903e910.
[76] A. Onufriev, D. Bashford, D.A. Case, J. Phys. Chem. B. 104 (2000) 3712e3720.
[77] M. Feig, A. Onufriev, M. Lee, W. Im, D.A. Case, C.L. Brooks III, J. Comput. Chem.
25 (2004) 265e284.
[78] F. Felluga, G. Pitacco, E. Valentin, A. Coslanich, M. Fermeglia, M. Ferrone,
S. Pricl, Tetrahedron: Asymmetry 14 (2003) 3385e3399.
[79] V. Frecer, M. Kabelac, P. De Nardi, S. Pricl, S.J. Miertus, Mol. Graph. Model 22
(2004) 209e220.
(c) D. Zampieri, M.G. Mamolo, L. Vio, E. Banfi, G. Scialino, Fermeglia, M.,
M. Ferrone, S. Pricl, Bioorg. Med. Chem. 15 (2007) 7444e7458 (and references
therein).
[46] A. Carta, M. Loriga, G. Paglietti, M. Ferrone, M. Fermeglia, S. Pricl, T. Sanna,
C. Ibba, P. La Colla, R. Loddo, Bioorg. Med. Chem. 15 (2007) 1914e1927.
[47] S. Manfredini, N. Solaroli, A. Angusti, F. Nalin, E. Durini, S. Vertuani, S. Pricl,
M. Ferrone, S. Spadari, F. Focher, A. Verri, E. De Clercq, J. Balzarini, Antivir.
Chem. Chemother. 14 (2003) 183e194.
[80] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem.
Phys. 79 (1983) 926e935.

A. Carta et al. / European Journal of Medicinal Chemistry 46 (2011) 4151e4167
4167
[81] H.J.C. Berendsen, J.P.M. Postma, W.F. Van Gunsteren, A. DiNola, J.R. Haak,
J. Chem. Phys. 81 (1984) 3684e3690.
[85] M.K. Gilson, K.A. Sharp, B.H. Honig, Calculating the electrostatic potential of
molecules in solution e method and error assessment, J. Comput. Chem. 9
(1988) 327e335.
[86] D. Sitkoff, K.A. Sharp, B.H. Honig, J. Phys. Chem. 98 (1994) 1978e1988.
[87] M.F. Sanner, A.J. Olson, J.C. Spehner, Biopolymers 38 (1996) 305e320.
[88] E.B. Wilson, J.C. Decius, P.C. Cross, Molecular Vibrations. McGraw-Hill,
New York, NY, 1995.
[82] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327e341.
[83] A. Toukmaji, C. Sagui, J. Board, T. Darden, J. Chem. Phys. 113 (2000) 10913e10927.
[84] P.A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee, T. Lee,
Y. Duan, W. Wang, O. Donini, P. Cieplak, J. Srinivasan, D. Case,
T.E. Cheatham III, Acc. Chem. Res. 3 (2000) 889e897.