Abstract
The intermolecular interactions of metallothionein with nitrogen mustard drugs were studied by molecular dynamics simulations. Previous laboratory experiments have defined selective alkylation of two cysteine residues, and selective binding was proposed to precede alkylation. The present study provides information about accessibility to cysteines based on evaluating the intermolecular energies and distances in the first few ps of dynamics simulations. A series of dynamics simulations was performed with three drug molecules positioned at the eight most solvent accessible cysteine residues of the dimeric form of the protein. Sites proximal to the sulfhydryl groups of Cys-33 and Cys-48 were found to be the most favorable for complexing the aziridinium forms of chlorambucil, melphalan, and mechlorethamine. The sites for preferential binding are in qualitative agreement with the sites of selective alkylation defined experimentally.
Metallothioneins (MTs)1 are small proteins consisting of 61 to 62 amino acid residues, among them 20 cysteines. Owing to this unique amino acid composition, MTs have the ability to bind seven or more divalent metal ions. This is the basis of the hypothesis that MTs have an important role in regulation of intracellular concentration of metal ions, such as Zn2+ (Vallee, 1995; Klaassen et al., 1999;Suhy et al., 1999). MTs similarly bind Cd2+ and Hg2+ ions; thus these proteins may be important in detoxification of heavy metal ions. It is very unlikely that this is an evolutionary result, and it may be only a property of these proteins, not their function (Vallee, 1987). Genetic experiments show that MTs are not essential or have a redundant function, because animals can reproduce in their absence (Michalska and Choo, 1993;Masters et al., 1994). However, conservation of the amino acid sequence of MTs through a broad phylogenetic range indicates that MTs have an important function. To date it has not been identified unambiguously (Palmiter, 1998).
The chelating thiolate groups in MTs create an environment of specific redox properties, so the MTs can be oxidized by cellular oxidants, releasing zinc ions, and the antioxidant activity of MTs is under study (Quesada et al., 1996; Lazo et al., 1998; Maret and Vallee, 1998). Experimental results also indicate that MTs are involved in cellular defense mechanisms. Overexpression of MTs is readily induced by a variety of endogenous and exogenous chemicals (Lazo and Pitt, 1995) and elevated MT levels have been correlated with resistance to electrophilic molecules used in chemotherapy (Kelley et al., 1988). Radioactive labeling experiments were interpreted to support covalent sequestration of chlorambucil (Endresen et al., 1983). More recent investigations have confirmed this suggestion in the case of that same drug and three related therapeutic alkylating agents: melphalan, mechlorethamine, and phosphoramide mustard (Yu et al., 1995; Zaia et al., 1996; Antoine et al., 1998; Wei et al., 1999). The sequencing of proteolytic products of MT-drug adducts by mass spectrometry also revealed that the alkylation of MT by these therapeutic mustards occurs preferentially at two of the 20 cysteine residues of the protein. The results of simple docking calculations (melphalan) and molecular dynamics simulations (chlorambucil, mechlorethamine) indicated a similar preference for reversible binding, which was suggested to occur before the alkylation reaction (Yu et al., 1995; Zaia et al., 1996;Antoine et al., 1998).
The most detailed information about the structure of MTs has been obtained by X-ray crystallography (Robbins et al., 1991) and NMR spectroscopy (Messerle et al., 1990; Narula et al., 1993; Peterson et al., 1996; Schultze et al., 1988; Wüthrich, 1991; Zhu et al., 1994). A detailed comparison of the structures obtained with the two methods pointed out that the crystal and solution structures of MT have consistent molecular architecture (Braun et al., 1992). The structure of MTs consists of two domains: the N-terminal (β) domain contains four, the C-terminal (α) domain contains three zinc ions. The major difference between the solution (NMR) and crystal (X-ray) structures is that although NMR provided information about “individual” domains, X-ray crystallography showed the relative orientation of the domains as well. Moreover, by X-ray crystallography the crystal packing could also be evaluated, which reveals an intimate association of two MT molecules having an approximately 2-fold symmetry with head-to-tail orientation.
The existence of dimeric MT in solution has also been studied. Ultracentrifugation results were interpreted to reveal a dimer (Palumaa et al., 1992), as was the behavior of MT during isolation by gel filtration. It was pointed out that the monomer-to-dimer ratio in solution depends on the concentration of the protein and also the concentration of excess metal ions (Palumaa et al., 1992; Otvos et al., 1993). The NMR studies emphasized the dynamic aspects of the MT dimer in solution, and cross-peaks corresponding to intermolecular proton-proton interactions were identified only in methanol-water solution (Otvos et al., 1993).
The aim of the present simulation study is 2-fold. On the one hand we have performed a series of docking simulations with three nitrogen mustard drugs (mechlorethamine, chlorambucil, and melphalan) using a more rigorous description of the protein-drug complex (including explicit representation of the solvent molecules) than in previous work. On the other hand, both the monomer and the dimer forms of the protein were modeled and the effect of dimerization on the binding preferences of drug molecules near cysteine residues is discussed and compared with earlier experimental observations.
Materials and Methods
In the calculations, the crystal structure MT isolated from rat liver was used (Robbins et al., 1991). The structure obtained from Brookhaven Protein Data Bank contains the atomic coordinates of the dimeric form of MT and explicitly shows the orientation of the monomers to each other (Fig. 1). Protons were assigned to the system assuming pH 7. This results in deprotonated acidic side chains and C-terminus carboxylate, and protonated lysine side chains in the sequence Ac-MDPNCSCATD10GSCSCAGSCK20 CKQCKCTSCK30KSCCSCCPVG40 CAKCSQGCIC50KEASDKCSCC60 A61. Subsequently all cysteines were deprotonated “manually” and the partial charges of Sγ and Cβ atoms of the cysteine residues were changed to −0.8 and −0.4, respectively (MacKerell et al., 1998). The metal ion composition of the original crystal structure is Cd5Zn2; in our calculation, however, Zn7-MT was used. Experimental (Messerle et al., 1992) and molecular modeling (Fowle and Stillman, 1997) studies have shown that the metal-sulfur connectivities and the overall structures of the metal binding sites are very similar for Cd- and Zn-MT. Besides these modifications, the Lys-51 side chain missing in the crystallographic structure was reconstructed in both monomers. The Na+ ions (one per monomer) and 174 “crystal” water molecules were retained.
Periodic boundary conditions were used in the simulations using a 32 × 45 × 55 Å cell in which the protein was solvated by approximately 1750 water molecules. This explicit representation of the solvent molecule and the Ewald summation method provides the most straightforward evaluation of electrostatic forces, which are among the most important determinants of ligand-receptor interactions. The charge distribution of the modeled system has been defined without counter ions. This choice was made because the location of the counter ions is not straightforward and it would have been an additional uncertainty in the simulations. The neutrality of the cell was ensured by the homogenous background charge density correction of the modeling software.
The standard consistent valence force field file, which had been extended with parameters of Zn2+ ions, was used. In the case of metal ions this modification is fairly easy, because only the nonbond (Lennard-Jones) parameters and charges have to be specified. The nonbond parameters used in this study had been determined by modeling the active site of HCAII (human carbonic anhydrase II) using a molecular orbital method (Hoops et al., 1991). In the active site of this enzyme, the Zn2+ ion is coordinated by three histidine residues and a water molecule, which is a different environment than Zn2+ chelation by four cysteine residues in MT. However, the utility of this approach is argued for dynamics simulations, because not the internal energy of the protein, but the energies of intermolecular interactions between the protein and the drugs are calculated.
Before the dynamic simulation the energy of the drug molecules was minimized. In this study the aziridinium ring structures of mechlorethamine, chlorambucil, and melphalan were used because previous experiments have shown that aziridinium ions are the predominant form of nitrogen mustards in solution at pH 7, and they are the alkylating agents (Golumbic et al., 1946; Colvin et al., 1976). These structures are shown in Fig. 2. Each minimized drug molecule was positioned close to each of the solvent-accessible (Robbins et al., 1991) cysteine residues (Cys-5, 7, 13, 33, 37, 41, 48, 57). In the initial orientation, the distance between the carbon atoms of the aziridinium ring of the drug molecule and the sulfur atom of the thiolate group was 2 to 3 Å. After solvation of the protein dimer-drug system, subsets of movable atoms were defined (amino acid residues within a 6-Å radius environment of the given sulfhydryl group; solvent molecules within an 8-Å radius environment; and the atoms of the drug molecule). The other atoms of the system were fixed during the simulation. The system thus defined was subjected to energy minimization (200 steps, conjugate gradient method) and dynamics simulation (temperature: 298 K; simulation time: 1 ps; timestep: 1 fs).
For model building and evaluation of the results, various modules of InsightII (version 97.0, Molecular Simulation Inc., San Diego, CA) molecular modeling software package (Builder, DeCipher, Analysis, etc.) were used; the energy minimizations and dynamics simulations were performed using Discover 3.
Results
The first step of the calculations described in this paper was an energy minimization of the solvated MT dimer. In this case a different solvation method (5 Å thick layer of water molecules around the protein) was used, and the minimization itself was longer (6000 steps, final divergence was 1 kcal/mol) than those performed before the dynamics simulations (200 steps). The aim of this longer energy minimization was to check the applicability of the Zn2+ nonbond parameters optimized for a different metal binding site, so it was performed without any constraint.
The result of this minimization shows that the backbone structure of the dimer of the protein remains very similar to that determined by X-ray crystallography. Only the side chains show significant movement compared with the starting geometry. The metal-thiolate clusters also retained their overall structure, only one of the zinc-sulfur “bonds” (Cys-44) was stretched.
In our calculations, flexible binding sites were used. The binding sites are defined as a given environment of a selected thiolate group consisting of 11 to 16 amino acid residues. One exception was made for Cys-33 when the binding site was defined as a set of 26 residues (13 residues in both domains) because of the “central” location of this particular cysteine in the structure of protein dimer (Fig. 1). Because these residues are not sequential and the remaining atoms of the MT dimer are fixed to their positions (according to the crystal structure) a “natural constraint” is used for the atoms of the binding site. The root mean square deviation curve (Fig.3a) shows that the binding subset is stable, even if the simulation time is longer (10 ps).
Figure 4 summarizes the results of 1-ps dynamics simulations performed with mechlorethamine. The average distance between the sulfur atom of the thiolate group and the carbon atom of the aziridinium ring of the drug that is closer to the sulfur atom during the simulation, and the energy of the intermolecular interaction between these atoms are shown for each cysteine. The energies were calculated taking into account both the monomer and the dimer forms of the protein.
Two other nitrogen mustard drug molecules, chlorambucil and melphalan, were studied using the same approach. The “binding” energies and interatomic distances obtained by dynamics simulations for the aziridinium forms of chlorambucil and melphalan are summarized in Figs.5 and 6.
Discussion
During energy minimization the structure of the dimer of the protein remains very similar to its crystal structure. The metal-thiolate clusters also retain their overall structures, lending validity to the use of the crystal structure for docking studies.
The modeling study contains some uncertainties, owing to the unique properties of MT-2. As it was mentioned above, the nonbond parameters for Zn2+ ions had been developed for a different ligand field. However, we believe that the nonbond parameters are reasonably applied to MT, in view of the highly dynamic structure of the metal-thiolate clusters, as pointed out by NMR (Palumaa et al., 1992; Otvos et al., 1993). The second uncertainty occurs in the protonation of thiolate groups. Accurate mass measurements using a double focusing mass spectrometer indicate that 6 of the 20 cysteine residues of MT are protonated at a pH value near 7 (Fabris et al., 1996). The most probable distribution of the six protons is that three thiolate groups are protonated in each domain; however, it is very likely that their locations change in a dynamic way. This phenomenon cannot currently be taken into account in molecular simulations; thus in our calculations all cysteine side chains are deprotonated. Neutrality of the cell was ensured by the homogeneous background correction of the Discover 3 software.
In general, the docking of substrates to proteins has two major difficulties. One of them is the inaccuracy of the molecular force fields. This has been mentioned above in connection with the nonbond parameters, but this problem cannot be avoided, even if the system does not contain atoms as exotic as zinc (Zimmer, 1995). The other difficulty is the high dimensional space in which the ligand-protein orientation has to be found at the global energy minimum (Nakajima et al., 1997). The binding site at each of the eight thiolate groups was characterized by the evaluation of a series of 1 ps dynamics simulations. With this very short simulation (without any constraint) it is not possible to find a global minimum of the energy surface of binding, but it may reveal the shortest distance possible between the binding site and the drug. In other words, the first few picoseconds of the dynamics simulation may show the accessibility of the selected sulfhydryl group by the drug. This accessibility can be characterized by the energy of the intermolecular interaction and the distance between the drug and the protein. The applicability of this approach was tested by a 10-ps simulation from which distances and energies are plotted against simulation time, as can be seen in Fig. 3, b and c. Both the energy and the distance values indicate that alkylation may occur in the first two picoseconds.
The energies of intermolecular interactions between mechlorethamine and the protein were calculated, taking into account both the monomer and the dimer forms of the protein. The difference between the results of the two evaluations is significant. Figure 4 shows the “binding” for mechlorethamine is stronger (lower energy) in the dimer form. This is not surprising, because the charge of the drug is +1 and the net charge of the dimer is −6. Another interesting feature of the energy curves is that the intermolecular energies corresponding to the binding sites in the β-(Cys-5, −7, −13) and α-domains (Cys-33, −37, −41, −48, −57) fall into two separate energy regions. For the monomer these are −60 to −95 and −100 to −140 kcal/mol, and for the dimer −80 to −110 and −140 to −190 kcal/mol). This is due to the “distribution” of charges between the two domains. The protonation of the system, as described above, results in a neutral β-domain and an α-domain of −3 net charge, so the electrostatic attraction between the drug and the α-domain is larger than that between the drug and the β-domain. Energy curves obtained for both the monomer and the dimer form of the protein show this tendency, with differences between the two electrostatic energy values of 20 to 30 kcal/mol for most of the sulfhydryl groups. The difference between the two curves exceeds this range in case of two cysteines, Cys-33 and Cys-48 (90 and 45 kcal/mol, respectively). When the dimer form of the protein is considered, the strongest intermolecular attraction can be observed for the sulfhydryl groups of these two cysteine residues. When the intermolecular energy is calculated between the monomer and aziridinium form of mechlorethamine, the lowest energy values are obtained for binding sites near Cys-37 and Cys-48. However, Cys-37 can not be considered as one of the most favorable sites for alkylation because the interatomic distance exceeds 4 Å and is the largest of those obtained for the cysteines of the α-domain. Experimental results obtained for the covalent sequestration of mechlorethamine by MT (Antoine et al., 1998) suggest that the initial alkylation takes place on Cys-48 with high selectivity. (Subsequently cross-linking occurs primarily at Cys-57.) Although the experimental work was carried out using rabbit MT-2 and the theoretical work is based on a crystallographic structure of rat MT-2, the comparison is considered to be valid because NMR studies have shown that these two mammalian MTs are isostructural (Vasak et al., 1987).
The antitumor drug, chlorambucil, was also modeled as its aziridinium ring (Williamson and Witten, 1967). Because it also has a deprotonated carboxylic group, it is neutral, and the evaluation of the series of dynamics simulations shows smaller attractions between the drug and the protein (Fig. 5). For the same reason the differences between the intermolecular energies calculated using the monomer and the dimer forms of the protein are also smaller than those obtained with the positively charged drug, mechlorethamine. Moreover, the calculated energy values are not clearly separated in contrast with those calculated for the two domains of the protein in case of mechlorethamine. If the energy values alone are taken into account, the differences among them are small and the most favorable binding site(s) is unclear. However, when the distances between the aziridinium group and the protein sulfhydryl group are also considered, the results indicate that the most probable alkylation sites are Cys-33 and Cys-48. This is in a good agreement with previous laboratory experiments (Zaia et al., 1996), which demonstrated that the first alkylation by chlorambucil occurs predominantly at these two cysteines (distribution of the alkylated products estimated by HPLC as 26 and 66%, respectively).
The third therapeutic alkylating drug studied was melphalan. This molecule is also a homolog of mechlorethamine and was modeled as its aziridinium form (Dulik and Fenselau, 1986). Its structure also includes a carboxylate and a protonated amino group resulting in a net charge of +1. The energies of intermolecular interactions calculated for the monomer and the dimer forms of MT, and the average interatomic distances between the reactive groups of the protein and the drug are shown in Fig. 6. Similar to mechlorethamine, which also has a net charge of +1, the difference between the calculated intermolecular energies for the monomer and the dimer forms is 15 to 30 kcal/mol for most of the cysteine sites studied. The two exceptions are Cys-33 and Cys-48, where the differences are 120 and 45 kcal/mol, respectively. The energies calculated for the dimer form represent the strongest binding of the drug. The dynamics simulations performed with melphalan placed close to each of the eight sulfhydryl groups result in the shortest intermolecular distances at Cys-33 and Cys-48 as well. Thus the results of the present simulation study are in good agreement with previous experiments, which show that Cys-33 and Cys-48 account for 23 and 66% of the initial alkylation sites, respectively (Yu et al., 1995).
The large energy difference between binding at Cys-33 in monomeric versus dimeric protein forms is observed for all three drug molecules. This probably reflects the central location of Cys-33 in the dimeric protein studied. Owing to the fairly compact dimer structure, binding and alkylation of Cys-33 requires that the drug penetrates into the central cavity formed by the four domains of the two monomers, where the positively charged aziridinium ring is close to the partial negative charge on Cys-33 residues of both monomers, resulting in larger electrostatic attractions than those obtained at other sulfhydryl groups. It is important to note that the sites of alkylation determined experimentally agree qualitatively with binding sites favored by present simulation study, only if the dimer form of MT was taken into account. One possible exception is Cys-33, which has not yet been identified experimentally as an alkylation site for mechlorethamine. The preferences determined taking into account monomeric MT do not converge readily with the experimental observations. For example, the binding of melphalan near Cys-33 appears less probable if MT monomer is the binding species. These observations strongly suggest that the dimer form of MT plays an important role in binding nitrogen mustard drugs. The qualitative agreement between theoretical and experimental results is also consistent with the hypotheses that selective binding precedes selective alkylation of MT by therapeutic nitrogen mustards.
Although glutathione has been proposed to bind to the amino terminus of MT (Brouwer et al., 1993) it has generally been difficult to detect or demonstrate nonmetal binding partners for MT. The nitrogen mustards may be viewed as models for endogenous substrate(s) that bind to the carboxy terminus domain, leaving evidence of their binding through alkylation, analogous to the photoaffinity strategy. The theoretical study reported here suggests that MT dimers should be used when searching for endogenous binding substrates.
The qualitative correlation shown here between dynamics simulations of binding complexes and experimental alkylation reactions also affirms the value of theoretical approaches in support of drug development.
Acknowledgments
We thank Ákos Bencsura (Institute of Chemistry of the Hungarian Academy of Sciences, Budapest, Hungary) for useful discussions.
Footnotes
-
Send reprint requests to: Catherine Fenselau, University of Maryland, Department of Chemistry and Biochemistry, College Park, MD 20742-2021. E-mail: fenselau{at}wam.umd.edu
-
This work was supported by National Institutes of Health Grant GM 21248.
- Abbreviation used is::
- MT
- metallothionein
- Received August 12, 1999.
- Accepted October 27, 1999.
- The American Society for Pharmacology and Experimental Therapeutics