10412
J. Phys. Chem. 1995, 99, 10412-10416
Monte Carlo Simulations of the Interaction between the Dodecanucleotide d(CAATCCGGATTG)z and Tris(ethylenediamine)cobalt(III) Cations B. Svensson, C. E. Woodward,* A. P. Arnold, J. G. Collins, and J. Lafitani Department of Chemistry, University College, University of New South Wales, Australian Defence Force Academy, Canberra, ACT 2600, Australia Received: October 25, 1994; In Final Form: April 4, 1995@
The interaction between a synthetic DNA dodecanucleotide and tris(ethylenediamine)cobalt(III) cations has been studied by Monte Carlo simulations. A simple two-state model with a single adjustable parameter is used to model the 'H-NMR chemical shift movements obtained in a parallel experimental study of the metalion complex-dodecamer interaction. The observed chemical shift movements, due to changes in the metalion complex and electrolyte concentrations, are quantitatively reproduced by this model. The simulations indicate that the binding of the A-enantiomer of tris(ethylenediamine)cobalt(III) cation can be mainly attributed to electrostatic forces. Furthermore, the results suggest that the enantiomeric selectivity observed in the experimental study can be explained in terms of a suitable complementarity of hydrogen-bond donors on the metal-ion complex with hydrogen-bond acceptors on the dodecamer.
Introduction The importance of metal ion interactions with nucleic acids has long been recognized. For example, metal ions are believed to play a crucial role in the replication, transcription, and translation of the genetic code.' Other studies have shown that metal ions are able to influence the flexibility and conformational properties of DNAs2 Recent work has also focused on the interaction of metal ions with smaller DNA fragments such as oligonucleotides and even mononucleotides. The study of smaller fragments facilitates a more detailed analysis which may help rationalize observations on natural DNA. For example, NMR studies have concluded that Mg2+ and Ca2+may display specific interactions with DNA3-5 and that Mn2+ and Zn2' possibly bind selectively to some oligon~cleotides.~~~ Recently, the interaction between nonintercalatingcobalt(III) cation complexes and the self-complimentarydodecanucleotide d(CAATCCGGATTG)2 has been monitored by proton NMR.s Modeling the dodecamer as a cylinder with a charge of -22 and the complex cation as a simple trivalent monopole, one would expect a significant electrostatic attraction. One interesting result of the experimental work, however, was that the binding of the chiral cobalt tris(ethylenediamine)cobalt(III) cation to the dodecamer was enantiomerically selective. The binding constant for the A enantiomer of the complex was estimated to be some 5-10 times higher than that obtained for the A form. Clearly the simple model described above would be unable to explain this result. Nevertheless, one would expect that the electrostatic interactions, as embodied in such a model, would play a large role in stabilizing the association between the complex and dodecamer. It is of interest, therefore, to initially establish the predictions of the simple electrostatic model for this system. Chirally discriminating interactions, irrespective of their origin, are expected to be relatively shortranged. Once the properties of the simple electrostatic model have been determined, the addition of some short-ranged potential may well be sufficient to describe the chiral discrimination. In this work, we present a theoretical investigation of the chemical shift movements of the dodecamer protons, brought @
Abstract published in Advance ACS Abstracts, May 15, 1995.
0022-365419512099-10412$09.0010
about by varying the concentration of the cobalt cation complex. Our results do indeed suggest a plausible explanation for stereoselective binding, which could be modeled by an additional short-ranged potential. Most previous theoretical analyses of ion-DNA interactions have been based on the Poisson-Boltzmann (PB) approximation or its linearized v e r s i ~ n . ~ .The ' ~ DNA molecule has been usually modeled as an infinite cylinder with a uniform surface charge density. However, the PB approximation has also been solved for DNA on a cubic lattice, wherein the molecular details have been explicitly accounted for." In this work we use the Monte Carlo (MC) simulation technique which has several advantages over the PB approximation. For instance, within statistical uncertainties, exact properties of the model are obtained. In addition, simulations may easily be adapted to complex systems without symmetry.
Theory In the experimental work reported here and elsewhere,8 selected protons on the dodecamer were monitored using NMR spectroscopy. One of the protons studied is found on the guanine of the eighth base pair of the sequence and is denoted here as the G8H8 proton. This proton lies close to the middle of the dodecamer chain. The chemical shift of the G12H8 proton, which lies on the guanine at one end of the dodecamer, was also studied. It is reasonable to assume that the chemical shifts of the nucleotide protons are altered by changes to the proton-shielding constant caused by the presence of the complex ions. The complex is strongly attracted to the oligonucleotide by electrostatic interactions; however, the chemical shift changes are most likely due to a conformational change of the dodecamer. One proposal is that this is effected by hydrogen bonding to N7 and 0 6 acceptors on guanine residues.'* In this work, we are not concerned with the specific mechanisms responsible for the chemical shift changes, but simply assume that they are proportional to some measure of the concentration of the complex ion close to the monitored protons. The concentration measure, denoted by clot, was obtained as the average concentration in a sphere of radius 10 8, centered about the N7 nitrogen. Thus, in our model 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 25, 1995 10413
d(CAATCCGGATTG)2 and Cobalt(II1) Cations
with 6 the observed proton chemical shift in the presence of the complex and d~ the proton chemical shift of the “free” (no cobalt complex present) dodecamer. Note that it has been verified that the complex ions are in fast exchange with the bulk solution on the NMR time scale.8 The parameter il is unknown and is empirically determined in this work. The value that best fits the experiments obviously depends on the radius of the sphere within which the complex concentration is averaged. We assumed that beyond 10 8, the complex ions would have a negligible effect on the proton chemical shift. This value is quite large when compared to similar theoretical treatments in other systems;I0 however, smaller radii gave a similar quality of fit to experiment. Although the center of the averaging sphere was chosen to be the N7 atom of the guanine residue, tests also showed that equally good fits to experiment were obtained by centering the sphere at any site close to the G8H8 proton. Our approach can also be interpreted in terms of a two-state binding model. In this model a complex ion is assumed to “bind” to an appropriate site close to the affected proton. When bound it causes the proton chemical shift to change. Given fast exchange, the observed chemical shift may be written as
where d~ is the chemical shift of the bound complex and the X, are the fractions of bound and free dodecamer. Equation 1 is recovered by assuming that XBis proportional to clK. Implicit in the two-state model is the assumption that it is the formation of a bound complex between the dodecamer and the complex ion which causes the chemical shift to change. An associated “binding constant” for the bound complex can then be fitted to the measured chemical shifts. This type of calculation was used to obtain “binding constants’’ for the different enantiomers of the cobalt complex.8 The model used in this work, eq 1, assumes no specific binding and, furthermore, allows any number of complex ions to associate with the fragment. The quantity, clot, arises as the appropriate ensemble average over all ion configurations.
Simulation Technique and Model System The simulations were performed in the canonical ensemble wherein the number of particles, temperature, and volume are constant. All interactions were described using the primitive model13wherein the ions were modeled as charged hard spheres and the solvent was replaced by a structureless continuum with an appropriatedielectric constant. The interaction energy, u(r), between two particles with valencies qi and % is thus given by
+ rij < (a,+ ay2
rij I(ai ajy2
(3)
where e is the absolute electronic charge, EO is the permittivity of free space, E~ is the dielectric constant of the solvent, and ai and Oj are the hard core diameters of the particles. The system was enclosed within a spherical cell containing an appropriate number and type of ions and the dodecamer at its center. A schematic picture of the model system is shown in Figure 1. The cell radius was determined by the concentration of the dodecamer. The dodecamer has been shown to adopt a B-type conformation in s01ution.I~During the simulation the dodecamer conformation was assumed fixed at a standard B-type conformation, built using the BIOSYM package Insight 11using
Figure 1. Shematic picture of the model system. The dodecamer atoms are fixed at the center of the cell and are shown as shaded spheres. The open spheres of different sizes represent the various mobile ions.
the CVFF force field.15 Different charge distributions on the dodecamer atoms were tried. In our most elaborate model, all the atoms were given partial charges according to the BIOSYM parameters.15 A more simple model was also tried wherein the phosphate oxygens were each given a negative charge and all other atoms were assumed neutral. It was found that both models gave the same ion density profiles, within statistical fluctuations. As the simple model was computationally cheaper, it was used to generate all the results reported below. A number of different types of ions had to be included in the simulation, in order to mimic the experimental system. These consisted of sodium and chloride ions as well as a 10 mM phosphate buffer. We assumed that the buffer dissociated to give equal amounts of H2PO4- and HP042- ions together with monovalent cations. The cobalt complex contributed one trivalent cation and three monovalent anions. The hard core radius for all atoms and ions (except the cobalt complex) was chosen to be 2 A. The cobalt complex was given a radius of 4 A. The choice of a 2 A radius for the smaller ions has been shown to agree with experimentally determined ionic activity coefficients for sodium ch10ride.I~ It has previously been shown, for a charged protein in solution, that simulation results of the type presented below are rather insensitive to the radii of the atomic species.17 The Monte Carlo simulations were performed using the algorithm of Metropolis et a1.I8 The displacement parameters for all the ions were adjusted throughout the simulation to give an acceptance ratio of 50%. The temperature was 298 K, and the dielectric constant was 78.3, correspondingto water at room temperature. Statistical uncertainty in the measured density of the complex ions close to the dodecamer surface can become large at low concentrations. A simple and computationally cheap way to improve the precision was to increase the frequency of complex ion displacements. For instance, when only one complex was present in the simulation cell, the probability of attempting to move this ion was increased by a factor of 20. This modification did not affect any of the calculated averages such as the partial energies or the ionic density profiles. The dodecamer model consisted of 486 atoms. As a simple means to decrease the computational demands of checking for hard core overlaps, the dodecamer was enclosed within a fictitious smallest possible sphere: an even more efficient way would be to consider the smallest possible cylinder. Whenever an attempted position of a displaced ion was further from the center of the cell than the sphere radius plus its own, it was not necessary to test for hard core overlaps with the dodecamer.
10414 J. Phys. Chem., Vol. 99, No. 25, 1995
Svensson et al.
7.8
a1
7.7
7.6 og
7.5
7.4
7.3
0.0
5.0
15.0
10.0
20.0
R
bl 4
7.90
i
7.88
I 7.86 0.0
I
10.0
20.0
R Figure 2. (a) Experimental and simulated chemical shift of the G8H8 proton at different cobalt concentrations. The system contains a 2.5 mM fragment, 100 mM sodium chloride, phosphate buffer, and EDTA. ( 0 )experiment; (0)simulations; i= -0.480. (b) Experimental and simulated chemical shift of the G12H8 proton at different cobalt concentrations. The same system was used in part a. (0)experiment; (0)simulations; 1 = -0.135.
Results and Discussion Given the differences in binding affinity between the different enantiomeric forms of the cobalt complex, one might expect closer agreement with the theoretical predictions for only one of the forms. Indeed it transpires that the more strongly bound A form fits best with the theoretical model; thus, all the experimental measurements reported below were for the A form of the cobalt complex. Chemical shift movements were recorded as the complex concentration was varied. The sodium chloride concentration was 100 mM, and the dodecamer concentration was kept close to 2.5 mM. The phosphate buffer, as well as EDTA, was present in the solution. The experimental shifts for the G8H8 and G12H8 protons, as well as the theoretical predictions, are shown in Figure 2a,b. We use R to denote the ratio of the concentrations of complex to dodecamer. The theoretical results contained a single adjustible parameter A (see eq l), which was chosen to fit the experimental observations for a given proton. The relevant values of l. are given in the figure captions, where the measure of cia: was the average number of complex ions in the sphere of radius 10 A about the site of interest. Overall there is excellent agreement between theory and experiment. It is particularly noteworthy that the simulations predict the quite different saturation rates for the G8H8 and G12H8 protons. Compared to G12H8, the G8H8 proton chemical shift plateaus at lower R, due to the larger electrostatic attraction of the cobalt complex to the central region of the dodecamer.
The sensitivity of our results to the size of the cobalt complex was investigated. A few simulations were performed with different complex ion radii. It was found that radii of 2 and 4 A gave similar results, while a radius of 6 A gave a significantly less rapid saturation for the chemical shifts. This is to be expected: as the distance of closest approach between the complex and dodecamer increases, the maximum attractive electrostatic interaction decreases. No attempt was made to locate the radius which gave the best agreement with experiment. That value would seem most likely to be somewhere between 4 and 6 A, which is in good accord with the estimated bare van der Waals radius of the complex. This result suggests one possible reason for the apparent poor binding of the A enantiomer to the dodecamer, compared to the A form. No explicit account is taken of water molecules or a dielectric discontinuity in our simulations. Thus, the good agreement with experiment for the A enantiomer suggests that this complex is able to readily lose its waters of hydration when it is close to the dodecamer. Apart from direct dipole interactions with the central ion, the hydration sphere is most probably further stabilized by hydrogen bonds (H-bonds) between water and the amine nitrogens on the ligands. It is then plausible that, for the A enantiomer, H-bonds with water are replaced with others to suitably sited nucleotide acceptors, but that the same situation is sterically impossible for the A enantiomer. One such H-bonding scheme is that suggested by Xu et a1.I2 for the interaction of C O ( N H ~ ) with ~ ~ + nucleotides containing adjacent guanine residues. They have proposed that c0(NH3)63+ can specifically H-bond to the N7 and 0 6 acceptors on both guanine residues. However, the possibility of an H-bonding scheme that incorporates the nucleotide phosphate groups must also be considered. Farmer et al.I9 have shown that metal amine complexes can form multiple H-bonds with guanosine 5'monophosphate, with H-bonds to the guanosine base and phosphate group (either directly or via a water molecule). Regardless of the identity of the H-bonds, we conjecture that the A complex can be modeled using a larger radius, accounting for a hydration sphere which remains intact. Other possible reasons for the difference in binding of the A enantiomer will be discussed in the Conclusions section. Of relevance to the discussion above is the possible error introduced by our assumption of a rigid dodecamer molecule. In reality we expect some conformational change in the dodecamer upon binding of the complex. Indeed Xu et aL20 found that d(CCCCGGGG)2 underwent a B -A transition upon c0(NH3)63+ binding. Although the chemical shift changes found in this study are also due to a B-A conformational shift, NMR experiments indicated that the extent of the conformational shift is relatively smaller for d(CAATCCGGATTG)2 than for d(CCCCGGGG)2. The greater B-A conformational shift observed by Xu et al. for d(CCCCGGGG)2, at a corresponding ratio of metal complex to oligonucleotide, is probably due to the larger number of adjacent guanine residues in d(CCCCGGGG)2. Further, implicit in our model is the assumption that any free energy changes induced by structural rearrangements are included in the short-ranged contribution to the free energy of binding, as they should be. From our discussion above, it is apparent that, for the A enantiomer, the short-ranged contributions are either unimportant or else canceled by the free energy of dehydration of the complex. This may not be the case for the A enantiomer. The predicted chemical shifts also depend upon the size of the sphere used to calculate clW. In other similar theoretical studies of NMR experiments,I0 it has been usual to use a radius of just a few angstroms. It is advantageous to use a large radius,
d(CAATCCGGATTG)z and Cobalt(II1) Cations
J. Phys. Chem., Vol. 99, No. 25, 1995 10415
7.6
[
7.5
I
0. 0
" 1
7.3
v
0.0
10.0
7.4
20.0
0
R
150
conc (mM)
Figure 3. Effect of the cutoff distance in the predicted shifts of the G8H8 proton. The system is the same as in Figure 2a. The parameter 1 has been chosen- to make the curve coincide at R = 10. Cutoff distance: (- - -) 3 A; (-) 6 A; (0) 10 A.
as this reduces the statistical noise. In Figure 3 the predicted shifts for the G8H8 proton are shown for different radii. The constant in eq 1 was chosen so that curves passed through the same point at R = 10. No dramatic differences are evident. Even a radius of 10 8, gives a slightly lower saturation rate than the other values, which gives marginally better agreement with experiment. Irrespective of the detailed mechanism effecting chemical shift changes, one would expect that the influence of the cobalt complex would be a fairly short-ranged function about some site near to the monitored proton. Although our choice of 10 8,appears quite large compared to other studies, results are essentially the same for smaller radii. Furthermore, the larger radius gives better numerical accuracy. The effect of sodium chloride on the chemical shifts was studied at a complex concentration of 3 mM and a dodecamer concentration of 1 mM. A large increase in the chemical shift of the G8H8 proton was observed as the sodium chloride concentration increased from 0 to 100 mM. The increase is easily rationalized on the basis of reduced electrostatic interactions between the cobalt complex and the DNA fragment, due to Debye screening. Similar behavior was also observed in the simulations. The experimental shift changes for the G12H8 proton were significantly smaller and in the opposite direction to what one would expect from simple screening arguments. A comparison between theory and experiment is shown in Figure 4a,b. The 1 parameter for a particular proton was not fitted to the experimental data but was the same as used in Figure 2a,b. Considering this, the agreement is very good. Nevertheless, there appears to be a small discrepency between theory and experiment for the rate of increase in the chemical shift for the G8H8 proton. A possible explanation may be that the ionic composition of the buffer may be affected by the salt concentration. This will have a larger influence at low ionic strength. The simulated numbers for the G12H8 proton carry significant statistical noise. This is because there is only a small probability of finding the complex at the ends of the dodecamer. As stated above, the experimental results seem contrary to simple electrostatic arguments. The most likely explanation is that at low salt concentration the experimental shifts are affected by increased fraying at the end of the fragment. This effect is not present in the simulations. At low salt concentration, simulations predict very little change in the chemical shift. This is due to two competing effects. At low salt, the difference in electrostatic potential between the center and ends of the dodecamer is large. This favors more complex accumulation
100
50
7.916
b] e
m 0
t-3
e
7.914
7.912
'
0
0
J 50
100
150
conc (mM) Figure 4. (a) Effect of sodium chloride additions on the observed shift of the G8H8 proton. The system contains 1 mM dodecamer, 3 mM cobalt complex, phosphate buffer, and EDTA. The same value of 1 was used as in Figure 2a. (0)experiment; (0)simulations. (b) Effect of sodium chloride additions on the observed shift of the G12H8 proton. The system is the same as in part a. The same value of l. was used as in Figure 2b. ( 0 )experiment; (0)simulations.
in the center at the expense of density at the ends: At higher salt, the complex is less strongly attracted to the end of the fragment, but more complex ions are available to bind, due to a reduced attraction to the center. We also investigated the effect of adding multivalent electrolyte on the chemical shifts. A series of simulations were performed for 1:-2 and 2:-1 electrolytes and compared with the results for the simple monovalent salt. In addition to the electrolyte, the system consisted of cobalt complex, 1 mM dodecamer, and 10 mM phosphate buffer (no EDTA was present in these solutions). In the fust series of simulations, the complex concentration was varied and the electrolyte concentration was kept constant. Either 10 mM 1:-2 or 2:-1 electrolyte or 30 mM 1:-1 electrolyte was used. These values corresponded to a fixed Debye screening length. Figure 5 shows the predicted shift on the G8H8 atom as R varies from 1 to 20. The 1:-1 and 1:-2 electrolytes induced similar changes, whereas the 2:-1 electrolyte showed markedly different behavior. Thus, the interaction between the complex and the fragment, at least for the 2: - 1 electrolyte, cannot be described by a single DebyeHuckel screening parameter, as is the case in simple approximate theories like the linearized PB equation. In Figure 6 the predicted shifts are shown as a function of the ionic strength for the different electrolytes; R is kept equal to 3. It is obvious that the difference in the predicted shifts depends markedly on the ionic strength of the solution. At low ionic strength addition of 1 mM 2:-1 salt, i.e., one divalent cation per dodecamer, has a significant effect on the shifts.
Svensson et al.
10416 J. Phys. Chem., Vol. 99, No. 25, 1995 7.8
7.7 A
7.6
e A
Lo
7.5
7.4
7.3
e
e
A A
e 5
9 15
10
20
R
Figure 5. Predicted shifts of the G8H8 proton induced by the addition of cobalt complex to the system containing various types of electrolytes at the same ionic strength. The system contains 1 mM dodecamer, 10 mM phosphate buffer, and added electrolyte. The concentrations of the 1:-2 and the 2: - 1 electrolytes are 10 mM, while the concentration of the 1:-1 electrolyte is 30 mM. (0)1:-1 electrolyte; (0)1:-2 electrolyte; (A) 2:- 1 electrolyte.
applies and that it provides an explanation for the differences in the binding strength of the A and A enantiomers. Simply put, the A enantiomer is unable to effectively replace its waters of hydration upon binding to the dodecamer, due to the lack of spatial complementarity of H-bond acceptors on the dodecamer. We suggest that this scenario can be modeled by different hard core radii for the enantiomeric forms. One cannot, however, dismiss the possibility that other mechanisms may be important in the binding of the A enantiomer. Our assumption has been that the chemical shift is proportional to the local concentration of cobalt complex, eq 1. If more complex binding mechanisms existed, e.g., cooperative binding to multiple sites, one could expect a more complicated functional dependence than eq 1, with possible implications on the apparent binding constant.21 Simulations also predicted that the chemical shift changes induced by different types of salt, at the same ionic strength, can be very different. Although there are no experimental results available for these systems, it is reasonable to believe that the predicted trends are correct as the simulations are able to predict the correct trends for the more difficult case of a 3:-1 electrolyte.
References and Notes (1) Eichhorn, G. L.; Marzilli, L. G. In Advances in Inorganic Biochemistry; ElsevierPJorth-Holland: New York, 1981; Vol. 3. (2) Record, M. T., Jr.; Mazur, S. J.; Melancon, P.; Roe, J.-H.; Shaner, S. L.; Unger, L. Ann. Rev. Biochem. 1980, 50, 997. (3) Rose, D. M.; Polnaszek, C. F.; Bryant, R. G. Biopolymers 1982, 21, 653. (4) Braunlin, W. H.; Drakenberg, T.; Nordenskioeld, L. Biopolymers 1987, 26, 1047. (5) Braunlin, W. H.; Drakenberg, T.; Nordenskioeld, L. Biopolymers 1989, 28, 1339. (6) Froeystein, N. Sletten, E. Acta Chem. Scand. 1991, 45, 219. (7) Froeystein, N. A,; Davis, J. T.; Reid, B. R.; Sletten, E. Acta Chem. Scand. 1993, 47, 649. (8) Watt, T. A.; Collins, J. G.; Arnold, A. P. Inorg. Chem. 1994, 33, 609. (9) Manning, G. S. Q. Rev. Biophys. 1978, 11, 179. (10) Gunnarsson, G.; Gustavsson, H. J . Chem. Soc., Faraday Trans. I 1982, 78, 2901. (11) Jayaram, B.; Sharp, K.; Honig, B. Biopolymers 1989, 28, 975. (12) Xu, Q.; Jampani, S. R. B.; Braunlin, W. H. Biochemistry 1993, 32, 11754. (13) Friedman, H. L. In Ionic Solution Theory; Interscience Publishers: New York, 1962. (14) Collins, J. G. Biochem. Int. 1988, 5, 819. (15) Biosym modelling package Insight IUDiscover, Biosym Technologies: San Diego, CA. (16) Card, D. N.; Valleau, J. P. J . Chem. fhys. 1970, 52, 6232. (17) Svensson, B.; Jonsson, B.; Woodward, C. E.; Linse, S. Biochemistry 1991, 30, 5209. (18) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. fhys. 1953, 21, 1087. (19) Farmer, P. J.; Cave, J. R.; Fletcher, T. M.; Rhubottom, J. A.; Walmsley, J. A. Inorg. Chem. 1991, 30, 3414. (20) Xu, Q.; Shoemaker, R. K.; Braunlin, W. H. Biophys. J . 1993, 65, 1039. (21) We thank the referee for pointing out this possibility.
A.;
Conclusions Despite the complexity of the system, the predictions of the theoretical model are in very good agreement with the experimentally observed chemical shifts induced by the A enantiomer of the cobalt complex. This suggests that the binding of the A form can be attributed to mainly electrostatic interactions. This requires that either H-bond and dispersion interactions are small compared to the electrostatic contributions included in the model (certainly not true at short separations) or else there is a fortuitous cancellation of short-ranged interactions upon binding of the complex. We have conjectured that the latter mechanism
JP942878G
