Molecular structure, dynamics and stabilities of di-gadolinium aminopolycarboxylate complexes

Eur J Med Chem (1991) 26,299-304 0 Elsevier, Paris

299

Molecular structure, dynamics and stabilities of di-gadolinium aminopolycarboxylate complexes R FossheimQ*, SG Dahll, H Dugstad2 ‘Institute of Medical Biology, University of Trams@, N-9001 Tromso; 2Nycomed AS, PO Box 4220 Torshov, N-0401 Oslo 4, Norway (Received

3 1 July 1990; accepted 5 November

1990)

Summary - The molecular structures, energies and dynamics of 2,3,5,6-tetrakis[N~-bis(carboxymethyl)aminomethyl]piperazineN,N’-diacetic acid (DIMER A), l,%bis[ 1 ,l,4-tris(carboxymethyl)-1,4-diazabutyl]-2,4-bis[[N~-bis(carboxymethyl)]amino]cyclohexane (DIMER B) and 1,4-bis[ 1,I ,4-tris(carboxymethyl)-1 ,4-diazabutyl]-2,5-bis[[N~-bis(carboxymethyl)]amino]cyclohexane (DIMER C) and their diGd-complexes were examined by molecular mechanical calculations and molecular dynamics simulations using the AMBER force field. The results indicate that DIMER A would form the most stable binuclear complex among the 3 ligands, mainly because the energy needed for the ligand to obtain the conformation in the complex is less in DIMER A than in DIMER B and DIMER C. The (RSRS), (RRSS), (ERRS), (SIRS) and (SSSS) isomers of DIMER A had similar complex stabilities in these calculations. The (SSRR), (RSSR) and (SW?) isomers of DIMER B and of DIMER C also had similar stabilities. Steric and electrostatic interactions with neighbouring groups determined the conformations of N-substituents in the ligands. RbumC - Structure, dynamique molCculaire et stabilit6 de complexes di-gadolinium aminopolycarboxylate. Les structures moleculaires, les energies et la dynamique mole’culaire de l’acide te’trakis[N,N-bis(carboxyme’thyl)aminomezine-N,N~diace’tique (DIMER A), du I ,5-bis[l ,I ,4-tris(carboxyme’thyl)-I ,I-diazabutyl]-2,4-bis[[N,N-bis(carboxymtthyl)]amino]cyclohexane (DIMER B) et du I ,4-bis[l ,I ,4-tris(carboxyme’thyl)-I ,4-diazabutyl]-2,5-bis[[N,N-bis(carboxyme’thyl)]amino]cyclohexane (DIMER C) et de leurs complexes di-Gd ont e’te’ e’tudies par calcul de me’canique mole’culaire et simulation de dynamique mole’culaire en utilisant le champ de forces AMBER. Les re’sultats montrent que le dim&e A repre’sente le complexe binuclt!aire le plus stable parmi les trois ligands, essentiellement du fait que l’e’nergie requise par le l&and pour adopter sa conformation complexe’e est plus faible dans le dim&e A que dans les dim&-es B et C. Les isomeres (RSRS), (RRSS), (RRRS), (SRRS) et (SSSS) du dim&e A ont des stabilite’s du msme ordre de grandeur dans nos calculs. Les isomeres (SSRR), (RSSR) et (SSSR) des dim&es B et C ont egalement des stabilite’s similaires. Des interactions ste’riques et e’lectrostatiques avec des groupements adjacents de’terminent la conformation des substituants des atomes d’azote des ligands. magnetic resonance namics simulations

imaging

/ gadolinium

complexes

/ complex

Introduction DTPA (1,4,7-triazaheptane- 1,1,4,7,7-pentaacetic acid) and other aminopolycarboxylates are used as chelating agents for paramagnetic ions in magnetic resonance imaging techniques [ 1, 21. Complexes of gadolinium (III) ions are of particular interest, due to their high magnetic moment. To be used as diagnostic agents, such paramagnetic complexes must have high

*Correspondence and reprints; Nycomed Torshov, N-0401 Oslo 4, Norway

AS, PO Box 4220

stability

/ molecular

mechanics

calculations

/ molecular

dy-

thermodynamic and kinetic stability, since lanthanide ions are highly toxic in man [2]. It is also desirable that the complexes have high water relaxivity in order to obtain a high image contrast [2]. We have previously examined various mononuclear aminopolycarboxylate ligands and their lanthanide ion complexes by molecular mechanical calculations and molecular dynamics simulations with the aim of gaining insight into factors that determine the stability of such complexes [3, 41. Since relaxation rates are linearly dependent on the concentration of paramagnetic ions in mononucleating complexes [2], it might be possible to obtain increased relaxivity without altering the osmomolarity by incorporating two paramagnetic ions in one ligand. Therefore, there is

300 considerable interest in bi- or higher nucleating ligands and it has been demonstrated that dimers of the chelating agent D03A (1,4,7,10-tetra-azacyclododecane-N,N’,N”&acetic acid) even have higher relaxivity per Gd-atom than the monomeric D03A complex [5]. We report on molecular mechanical calculations and molecular dynamics simulations on the three binucleating ligands and the isomers shown in scheme 1 and their di-Gd complexes, in order to evaluate their relative complex stabilities.

simulations were performed with step length 0.001 ps, and with constraining of all bonds involving hydrogen atoms. The coordinates were saved on a disk file 100 times during each simulation. Molecular graphics were performed with the MIDAS programs [7] on an Evans and Sutherland PS390 workstation with a MicroVax II/Ultrix system as host computer, and with the ALCHEMY [S] program on an MS-DOS personal computer. Atomic point charges

Methods

Atomic point charges, which were used in the molecular mechanical calculations and molecular dynamics simulations, were obtained from previous ab initio calculations on mononuclear ligands [3, 41, which were obtained with the QUEST 1.0 program [9]. QUEST calculates quantum mechanical electrostatic potentials in a number of layers around the molecule, and net atomic charges by least squares fitting to the electrostatic potentials. Electrostatic potentials in 4 layers 0.2 8, apart, the innermost layer corresponding to 1.4 times the van der Waals radii, and atomic point charges were calculated on a CRAY X-MP computer using an STO-3G basis set. Charges for the monomeric unit of DIMER A (2,3,5,6-tetrakis [N,N-bis(carboxy-methyl)aminomethyllpiperazine-NJ’-diacetic acid) were obtained from charge calculations on a DTPA analog having the central nitrogen atom constrained in a piperidine ring [4]. Equal charges were assigned to atoms related by mirror symmetry. The C-atom charges of the pipe-

Molecular mechanical calculations and molecular dynamics simulations were performed with the AMBER 3.0 software package [6] on a VAX 8600/VMS system. All energy minimizations and molecular dynamics simulations were performed with a distance dependent dielectric, E = rv, where rq is the interatomic distance. A non-bonded cut-off distance of 99 A was used in the energy minimizations and molecular dynamics simulations. The energy minimizations were performed by the steepest descent method for the initial 10 cycles, with an initial step length of 0.05, followed by conjugate gradient minimization until convergence. Convergence was assumed when the difference in energy or gradient was less than 0.05 between subsequent steps. The molecular dynamics simulations were run for 10 ps at constant temperature (300 K), after an initial 1 ps equilibration period starting at 0.1 K. The dynamics

DIMER B

I RSRS

Scheme 1.

II RRSS

Ill RRRS

IV SAW

v ssss

Y

X

0 X

Y

DIMER C

II RSSR

III SSSR

301 razine ring were adjusted in order to give the correct total charge of the ligand. The same charge model was used for all isomers of DIMER A. Charges for DIMER B (1,5-bis[ 1,1,4-tris(carboxymethyl)-l,4-diazabutyl]-2,4-bis[[N,N-bis(carboxymethyl)]amino]cyclohexane) and DIMER C (1,4bis[ 1,1,4-tris(carboxymethyl)-1,4-diazabutyll-2,5bis[[N,N-bis(carboxy-methyl)]amino]cyclohexane) were obtained by assigning the charges calculated for corresponding atoms in DTPA[3] to the two DTPAparts of the molecules. Charges for the two -CH,groups in the cyclohexane ring were obtained from a charge calculation on a DTPA derivative having the central nitrogen constrained in a piperidine ring [4], which contains -CH2- groups in the same bonding environment. The -CH,- atomic charges were adjusted to give each ligand a correct total charge. The same set of charges were used in calculations on the (SSRR), (RSSZ?) and (SSSR) isomers.

complexes of binucleating aminopolycarboxylate ligands would be expected to coordinate two water molecules. Models of dihydrated complexes were built with the graphics system by docking of water molecules into the vacant coordination sites, using the van der Waals surfaces of the complexes. The Gd,L*2H20 complexes (L: ligand) were then energy minimized. As a measure of complex stability of the binuclear complexes, the reaction energy of complexation (&) between a ligand and two cations was calculated by equation (1). ER=EML-EL

(1)

EML is the energy of the ligand-two cation complex. EL is the energy of the free ligand in its most stable conformation, which was found from molecular dynamics simulations in vacua and molecular mechanical calculations as described above.

Free ligands Starting geometries for all ligands were obtained from model building with the ALCHEMY program. The following calculations were then performed on each of the ligands in their unprotonated forms: - the geometries were initially optimized by energy minimization without electrostatic interactions, using the AMBER programs; - a further geometry optimization was performed including electrostatic interactions; - these refined ligand structures were used as starting points for the molecular dynamics simulations; - the 3 structures which had the lowest potential energy among those observed in each molecular dynamics simulation, were selected and refined by energy minimization. Gd(III) ion complexes Previously determined values for van der Waals well depth e* = 0.2 and van der Waals radius R* = 1.8 [3] were used as non-bonded Gd(II1) parameters in the calculations. Amine nitrogen charges of - 0.4 were used, as it has been shown that this value gives calculated structures with good agreement to Gd-N and Gd0 distances found in crystal structures [3]. Nitrogen charges calculated for the free ligands were usually in the range of - 0.1 to - 0.3. Starting coordinates for the diGd complexes were obtained with the graphics system using the MIDAS programs. The torsion angles of the free ligands were adjusted in order to give a reasonable coordination geometry of the centrally placed metal ions, and the diGd complexes were refined by energy minimization. Since DTPA and related octa-dentate ligands usually coordinate one water molecule in their Gd complexes,

Results Table I gives molecular mechanical potential energies of the lowest-energy conformer of each of the nine free ligands which are shown in scheme 1, together with the total energies of the dihydrated complexes. Table I also shows the energies of the ligand and two cations, interaction energies between one cation and the ligand, cation-cation interaction energies, the total energy of the ligand in the complex conformation and the calculated reaction energies of complexation. The five isomers of DIMER A shown in scheme 1 are the unique isomers of this compound. There are additional possible isomers for DIMER B and DIMER C than those shown in scheme 1. However, initial model building indicated that the (SRSR) and (SSSS) isomers are less likely to form stable complexes, either due to steric crowding or to poor coordination. Figure 1 shows stereoscopic drawings the DIMER A IV, DIMER B I, DIMER B II and DIMER C I complexes. Only one isomer of DIMER A is shown in figure 1 since there were no substantially different features in the structures of the various DIMER A isomers. Two isomers of DIMER B are shown in figure 1 in order to illustrate the differences in steric crowding in the complexes. Only one isomer of DIMER C is shown in figure 1 since the structural differences between these isomers were analogous to the differences between the DIMER B isomers. The piperazine ring in the free ligands had a chair conformation in all energy minimized conformers of isomers I-V of DIMER A. The molecular dynamics simulations indicated that the conformational flexibility of the piperazine ring is restricted by the 6 ring substituents. This is illustrated in the upper left part of

302

Table I. Total ligand energy (EL), total dihydrated complex energy (ET), complex energy (E&, ligand-cation interaction energies (Eci and EC*), cation-cation repulsion energy (.!&), energy of ligand in the complex conformation (&-) energy (&a) for complex formation (kcal/mol) of aminopolycarboxylate diGd-complexes. Isomer

Ligand

Complex

EL

E ML

ET

-

821.8 819.8 824.7 826.7 824.9

Reaction

E Cl -

E LC

E R-g

685.0 684.7 685.3 685.1 684.3

42.8 42.4 43.7 43.9 42.9

501.6 503.0 503.0 501.5 495.8

-

-

BI B II B III

248.0 261.1 257.0

- 882.9 - 867.8 - 872.8

- 846.5 - 830.3 - 836.6

- 682.4 - 691.2 - 685.8

- 687.5 - 686.2 - 688.8

33.8 37.1 34.8

489.5 510.0 503.1

- 1094.5 - 1091.4 - 1093.6

CI c II c III

247.5 261.7 260.6

- 882.4 - 871.3 - 874.4

- 846.8 - 834.4 - 838.5

- 684.1 - 686.7 - 691.1

- 687.8 - 685.4 - 695.0

34.0 35.2 38.8

491.1 502.4 508.8

- 1094.3 - 1096.1 - 1099.1

g Fig 1. Stereoscopic drawings of diGd complexes of DIMER A isomer IV, DIMER B isomer I and II, and DIMER C isomer I.

-

E Cl2

288.8 288.5 287.9 284.8 282.7

*

681.0 680.5 686.1 687.0 679.4

E c2

AI A II A III A IV AV

+!t$$$t

859.8 863.4 860.7 864.1 859.9

and reaction

1110.6 1108.3 1112.6 1111.5 1107.6

figure 2 which shows the structures of DIMER A IV observed at 0.5 ps intervals during 10 ps of molecular dynamics simulation. Inversion of the piperazine ring was not observed during the molecular dynamics simulation. The axial conformation, which was seen in all conformers of the complexes except one (DIMER A I), appears to be preferred for the Nacetate groups. In DIMER A I one N-substituent was equatorial in one of the 4 energy minimized conformers. This conformer had 1.4 kcal/mol higher energy than the lowest-energy conformer of this isomer. DIMER A I is the only isomer that had both the neighbouring ring 2,6-substituents in axial positions at one end of the piperazine ring. This indicates that steric and electrostatic interactions with the neighbouring substituents is important for the conformation of the N-substituents. The 1,3-diaxial interaction found at one end of the piperazine ring in DIMER A I does not appear to be energetically unfavourable, judged from the total ligand energies of the 5 DIMER A isomers (table I) The cyclohexane ring in DIMER B I and DIMER C I also has restricted conformational flexibility, as indicated by the molecular dynamics simulations. This is illustrated by the superimposed structures, observed at 0.5 ps intervals, which are shown in figure 2. All

303

Fig 2. Superimposed drawings of structures observed at 0.5 ps intervals during 10 ps of molecular dynamics simulations of DIMER A IV, DIMER B I, DIMER B II and DIMER C I.

eight energy minimized conformers of DIMER B I and DIMER C I had cyclohexane chair conformations with all 4 ring-substituents in equatorial positions. On the other hand, the cyclohexane ring in DIMER B II showed larger conformational flexibility, and the cyclohexane ring adopted chair, boat and twist-boat conformations during the 10 ps molecular dynamics simulations, as illustrated in figure 2. Only chair conformations were observed during the molecular dynamics simulation of DIMER C II and DIMER C III. DIMER B III on the other hand adopted boat and twist-boat conformations during the molecular dynamics simulation. In this isomer the cyclohexane ring also adopted a boat conformation in the lowestenergy conformer found, whereas the lowest-energy conformers of the other DIMER B and DIMER C isomers adopted chair conformations. All isomers of DIMER A have higher complex stabilities than the isomers of DIMER B and DIMER C. There were only small differences between the complex stabilities of the DIMER A isomers (table I). The total energies of the complexes (ET) and the ligands (EL) were similar. As indicated by the ligandcation interaction energies (Ec, and E,,) given in

table I, all isomers of DIMER A may obtain a favourable coordination geometry at both ends of the piperazine ring. This geometry may be attained without sterical hinderance, as indicated by the similar energies of the ligands in their complex conformations (EL,-). Furthermore, the cation-cation distances were similar in the five isomers, ranging from 8.28 to 8.53 A. The cation-cation repulsion energy (EC,J does not contribute therefore to distinguish the five isomers of DIMER A regarding their complex stabilities. The total energies (/ST) of the DIMER B I and DIMER C I complexes were lower than those of the corresponding isomer II and III complexes. As can be seen from table I, the main reason for this difference was that the energy of the ligand in the complex (&c) was higher for isomer I than for isomer II of DIMER B and DIMER C (table I). As illustrated in figure 1 for DIMER B I and II, there were generally more short H...H contacts between the ring H-atoms and sidechain H-atoms in isomers II and III than in isomers I (fig 1). DIMER B I has five such H...H contacts with distances between 2.20 and 2.10 A. DIMER B 11 and III have six such H...H contacts, one of which is 2.02 A in B II. Four H...H contacts in isomer I of

304 DIMER C had distances ranging from 2.10 to 2.20 A. In isomer II of DIMER C there were six short H...H contacts, two of which were 2.08 Ad whereas C III had eight H...H contact less than 2.20 A. The larger steric crowding in isomers II and III of DIMER B and DIMER C resulted in deformations of bond-angles, in particular those pertaining to the cyclohexane ring. The total angle strain energy was 16.4 kcal/mol for both DIMER B I and DIMER C I, 29.1 kcal/mol for DIMER B II, 26.4 kcal/mol for DIMER C II, 24.2 kcal/mol for DIMER B III and 23.2 kcal/mol for DIMER C III. The ligand-cation interaction energies (E,, and Ec2) and the cation-cation repulsion energies (E,,,) were similar for the isomers (table I), the largest ligand-cation interaction energies being found in DIMER C III. Since the free ligands of DIMER B I and DIMER C I had lower energies than the corresponding isomer II and isomer III ligands, the complex stabilities calculated according to equation (1) were similar for the six compounds within the accuracy of the calculations.

Discussion The calculated complex stabilities of 2,3,5,6tetrakis[NjVbis(carboxymethyl)aminomethyl]piperazine-N,N-diacetic acid (DIMER A), 1,5-bis[ 1,1,4tris(carboxymethyl)1,4-diazabutyll-2,4-bis[ [N,Nbis(carboxymethyl)]amino]cyclohexane (DIMER B) and 1,4-bis[ 1,1,4-tris(carboxymethyl)1,4-diazabutyll2,5-bis[[N,N-bis(carboxymethyl)]amino]cyclohexane (DIMER C) indicate that the piperazine derivative (DIMER A) forms the most stable di-Gd complex among the three, and that the two cyclohexane derivatives (DIMER B and DIMER C) form diGdcomplexes with similar stabilities. This difference was mainly caused by the energies required to bring the free ligands into their respective complex conformations, which generally were lower for DIMER A than for DIMER B and DIMER C. This may be seen from the difference in energy between the ligands in the complex conformation and the energy of the free ligand (&-EL), which was in the range of 212.8-216.7 kcal/mol for the DIMER A isomers and in the range of 230.1-248.2 kcal/mol for the DIMER B and DIMER C isomers. A qualitative explanation for this difference may be found by comparing the structures of the ligands. In DIMER A all the piperazine ring atoms are substituted with relatively short pendant carboxylate arms. In DIMER B and C the cyclohexane ring is 1,2,4,5-substituted with two relatively long pendant carboxylate arms. This forces the carboxylate groups closer together in DIMER A than in DIMER B and C. In this respect, the DIMER A free ligands are closer to

the complex conformation, where the carboxyl groups are brought closer by coordination, than the DIMER B and C free ligands. When the negatively charged carboxylate groups are brought closer together, the potential energies of the free ligands become larger, which leads to higher complex stability according to equation (1). These calculations do not consider the contributions of solvent effects or entropy effects to complex stabilities. However, from the structures of these compounds it seems likely that ligands and complexes will be equally well solvated. Both contain the same number of carboxylate and amine groups of which all are involved in cation-coordination in the complexes. Neither of the compounds may form intramolecular hydrogen bonds. For these reasons, effects of solvation on complex formation would not be expected to be significantly different in these compounds. Differences in reaction entropies are also deemed small. All compounds are likely to coordinate the same number of water molecules, which would produce a similar increase in the entropy by loss of coordinated water molecules from the Gd(II1) ions. Since all ligands contain a six-membered ring with side-chains, the entropy decrease by loss of conformational mobility of the ligands on complex formation would not be expected to be significantly different. The potential use of these binucleating ligands as chelating agents in magnetic resonance imaging techniques depends on several factors which may be assessed by synthesis and biological testing of the compounds. Since thermodynamic stability is one important factor for the use of such complexes in magnetic resonance imaging techniques, the performed calculations indicate that the piperazine derivative DIMER A should be considered as the first choice among the three as a synthetic target.

References Desreux JF, Barthelemy PP (1988) Nucl Med Biol 15, 9 Lauffer RB (1987) Chem Rev 87.901 Fossheim R,‘DahlSG (1990) Acta Chem Stand 44,698 Fossheim R, Dugstad H, Dahl SG (1991) J Med Chem 34, 819 Gries H, Platzek J, Raduchel B, Renneke FJ, Weinmann HJ (1990) Poster abstract, European Congress of NMR in Medicine and Biology, Strasbourg 2-5 Miy 1990 Sinnh UC. Weiner PK. Caldwell JW. Kollman PA (1986) AMBER (UCSF, version 3.0). Department of P&maceutical Chemistry, University of California, San Francisco Ferrin TE, Huang CC, Jarvis LE, Langridge R (1988) J Mol Graphics 6, 13 ALCHEMY: Distributed by Tripos Associates Inc, 1699 South Hanley Road, Suite 303, St Louis, MO 63144, USA Singh UC, Kollman PA ( 1984) J Comp Chem 5, 129