Brownian Dynamics Simulations of Polyalanine in Salt Solutions

Department of Physiology and Biophysics, Mount Sinai School of Medicine, New York, New York ... ReceiVed: January 10, 1996; In Final Form: June 25, 19...
0 downloads 0 Views 468KB Size
Download PDF
15280

J. Phys. Chem. 1996, 100, 15280-15289

Brownian Dynamics Simulations of Polyalanine in Salt Solutions Wenbin Yu† and Chung F. Wong* Department of Physiology and Biophysics, Mount Sinai School of Medicine, New York, New York 10029

John Zhang* Department of Chemistry, New York UniVersity, New York, New York 10003 ReceiVed: January 10, 1996; In Final Form: June 25, 1996X

The distribution and dynamics of the mobile monovalent ion atmosphere around polyalanine in aqueous solution at various ionic strengths have been simulated by the Brownian dynamics method. Initial simulations of pure sodium chloride solutions were performed to examine the influence of the choice of nonbonded cutoff distance on simulation results. Simulations were then performed to study how different conformations of polyalanine affect the distribution and dynamics of the mobile monovalent ion atmosphere. We found that it was necessary to use a nonbonded cutoff >30 Å in order to get reliable results for the ion pair radial distribution functions, the ionic polarizabilities, and the auto time-correlation functions of the collective dipole moments of the sodium chloride solutions. We also found that R-helical polyalanines preferentially bound sodium ions at 0.1 M sodium chloride concentration. The preferential binding of sodium ions was still noticeable but less pronounced when the polyalanines were in the extended conformation or when the concentration of the salt was increased. The polarizability of the ions was found to be insensitive to the presence of R-helical polyalanines but became smaller in the presence of extended polyalanines. The relaxation of the collective dipole moment of the ions was also found to be affected only by polyalanines in the extended conformation but not in the R-helical conformation. The biological significance of the simulation results is discussed.

I. Introduction It is well-known that the stability of proteins can be affected by the presence of salts. Timasheff and co-workers1,2 have suggested that the ability of a specific type of ions to stabilize or destabilize a protein can be correlated with the distribution of this type of ions around the protein. The study of the mobile ion atmosphere around a protein in different conformations may thus shed light on whether a specific ion is stabilizing the protein and which conformations of the protein it may preferentially stabilize. Theoretically, several methods have been used to study the distributions of ions around proteins and DNAs. A very popular and computationally-inexpensive approach to determining the ion distribution around a protein is via the solution of an appropriate Poisson-Boltzmann equation.3-5 Although the linearized approximation of the Poisson-Boltzmann equation has been used in many calculations, more accurate nonlinear versions are also becoming popular. Recent models for Poisson-Boltzmann calculations also account for the detailed shape of a biomolecule. Monte Carlo simulations have also been used to study ion distributions around biomolecules.6 In the Monte Carlo simulations, the ions are explicitly included in the simulations so that no linearized approximations of the Boltzmann distribution terms are used. The finite size of the ions are also accounted for explicitly. Neither the PoissonBoltzmann nor the Monte Carlo approach, however, provides information about the dynamics of ions arond a biomolecule. Although molecular dynamics simulations can be and have been carried out to study the equilibrium structures and dynamics of † Current address: Photonic Materials and Devices Laboratory, Department of Materials Science and Engineering, University of Southern California, Los Angeles, CA 90089. X Abstract published in AdVance ACS Abstracts, August 15, 1996.

S0022-3654(96)00124-4 CCC: $12.00

ions around biomolecules, these simulations are computationally very expensive because many explicit solvent molecules need to be included in the simulations. Since the number densities of ions for physiologically-relevant concentrations of salts are very small, it typically takes more than a few hundred nanoseconds of simulations to give statistically-significant descriptions of the distributions of ions around biomolecules. This is because the ions have a large amount of space to sample during the simulations. Since it is still difficult to carry out molecular dynamics simulations to such a time scale, we here explore the use of the Brownian dynamics simulation method to study the equilibrium distributions and dynamics of ions around biomolecules. The Brownian dynamics simulations carried out here also allow us to examine more easily the effects of the choice of nonbonded cutoff on the simulation results of ionic systems. In typical molecular dynamics simulations, it is usually difficult to use nonbonded cutoffs more than ∼15 Å due to the large computational costs involved. In this paper, we utilize the Brownian dynamics method to examine the extent to which the use of short nonbonded cutoff distances affects the simulated structural and dynamical properties of ionic systems. In addition to influencing the stability of proteins and DNAs, the altered behavior of salts around biomolecules may also reflect the conformational preference of biomolecules. Frequencydependent dielectric measurements of DNAs, for instance, have identified two dielectric dispersions in addition to the one arising from water at ∼17 GHz.7,8 The lowest frequency dispersion (j1 MHz), which is molecular weight sensitive, is thought to arise from the relaxation of bulk ions and the redistribution of bound ions over the whole length of a DNA molecule. The dielectric dispersion occurring between ∼1 MHz and ∼1 GHz is postulated to result from the fluctuation of counterions over short segments of DNAs. This fluctuation has been suggested © 1996 American Chemical Society

Brownian Dynamics Simulations of Polyalanine

J. Phys. Chem., Vol. 100, No. 37, 1996 15281

to be sensitive to the dynamic bending of DNAs and therefore the study of the dielectric dispersion in the 1 MHz-1 GHz region may yield useful information about the conformational preference of DNA molecules. Dielectric measurements can also probe the conformational properties of polypeptides and proteins. One of the goals of the Brownian dynamics study presented here is to examine whether the dielectric properties of polypeptides are sensitive to the conformation of the polypeptides. The study of the frequency-dependent dielectric properties of an ionic solution requires the calculation of the time-correlation function of the collective dipole moment of the solution. This information can only be obtained from models, such as the Brownian dynamics simulation method employed here, that also yield information about the dynamics of an ionic system. Monte Carlo simulations, for example, cannot provide dynamical information. The Brownian dynamics simulation models utilized here also remove some of the approximations employed in some earlier simplified models. For example, the interpretation of the experimental data on the dielectric properties of DNA solutions has so far been based on relatively simple models such as those introduced by Manning9,10 and Mandel.11,12 In these models, DNAs are treated as infinitely thin rigid rods with continuouslydistributed charges or charges placed symmetrically along the rods. In some cases, only the interactions between the DNAs and the counterions are included and the interactions among the counterions are ignored. The more realistic Brownian dynamics simulation models presented here can therefore provide a better microscopic understanding of the frequencydependent dielectric behavior of solutions of proteins and DNAs at various salt concentrations. In this paper, we choose to study a 13-residue polyalanine. Since the side chain of an alanine residue is uncharged, the polyalanine model allows us to first examine how the differences in the electrostatic fields arising from different arrangements of the backbone peptide groups in different conformations of the polypeptide can affect the structural and dynamical properties of the ion atmosphere around the polypeptide, without considering the further complications introduced by the presence of polar or charged side chains. II. Simulation Models and Method The systems we study are composed either of sodium and chloride ions only or with a 13-residue polyalanine immersed as well in a cubic or rectangular box. The screening effects of the solvent on the ion pair interactions are modeled by using a constant dielectric coefficient  of the solvent. This is the approximation used in most Monte Carlo simulations of proteins and DNAs carried out to date. Periodic boundary conditions are also used in all the simulations. The sodium and chloride ions interact through the Friedman-type solvent averaged potential:13,14

( )

r*i + r*j 9 zizje2 uij(r) ) Bij + r r

(1)

Here r is the distance between ions i and j, r* i is the radius of ion i, which is taken to be 0.97 Å for sodium and 1.81 Å for chloride in our calculations, zi is the valence of ion i, e is the electronic charge, and  is the solvent dielectric constant, which is taken to be 80 in our calculations. When r*i is expressed in Å, the coefficient Bij is given by βBij ) (5377.75/T)|zizj|/ -1 ) k T and T is the absolute temperature. (r*i + r* B j), where β The 13-residue polyalanine is put in the center of the simulation cell with the long axis roughly parallel to the y axis of the simulation cell. In the R-helical conformation, the long

axis of polyalanine is about 20 Å. The width of the R-helical polyalanine is about 8 Å. All the polar hydrogens of the polyalanine are explicitly included as is done in the GROMOS molecular modeling package.15 The atomic partial charges and the Lennard-Jones parameters for the polyalanine are also taken from GROMOS.15 The polyalanine interacts with the ions through a LennardJones plus Coulomb potential:

uib(r) ) 4ib

(( ) ( ) ) σib r

12

-

σib r

6

+

zizbe2 r

(2)

Here r is the distance between ion i and atom b of the polyalanine. ib and σib are the parameters of the Lennard-Jones potential between ion i and atom b of the polyalanine. zi is the valence of the ion and zb is the partial charge of the atom b in electrostatic unit. As for the Friedman-type potential,  ) 80 is used for the dielectric constant. In order to avoid overcountering interactions and save computational costs, the nonbonded interaction forces are usually truncated in molecular dynamics and Brownian dynamics simulations. In this paper, all the short-range interactions are truncated at 10 Å. Increasing the truncation length for the shortrange interactions to >10 Å has little effect on the simulation results. On the other hand, a large cutoff for the electrostatic interactions is required for giving reliable results. We have carried out simulations with a number of different cutoffs for the long-range electrostatic potential, and the results obtained by using two cutoffs, 20 and 40 Å, are presented in this paper. A propagation equation for following the Brownian motion of a system of N interacting particles has been developed by Ermak and McCammon16 and is used in this study. When hydrodynamic interactions are neglected, the Ermak-McCammon equation takes the following form:

ri(t + ∆t) ) ri(t) + βDiFi∆t + Ri

(3)

where ri(t) is the position of particle i at time t, β-1 ) kBT where T ) 300 K in this work, Di is the diffusion constant of particle i which is taken to be 1.483 × 10-5 cm2/s for both Na+ and Cl- (this is estimated from the experimental data of sodium chloride solutions at 300 K17), Fi is the force experienced by particle i at time t, ∆t is the time step, and Ri is a Gaussian random displacement having the characteristics 〈Ri〉 ) 0 and 〈RiRj〉 ) 2Di∆tδijI where I is a unit tensor. This random displacement is due to the constant bombardments of particle i by the solvent molecules. In the Ermak-McCammon algorithm, the time step ∆t should be larger than the momentum relaxation time of the particles but should be small enough so that the forces acting on the particles are almost constant during ∆t. ∆t ) 0.04 ps is used in this work. When using such a large time step, two atoms can occasionally get too close to each other and experience a large repulsive force due to the Lennard-Jones potential. One can remove this artifact by using approaches, such as the Langevin dynamics method, that model the inertial effects of ion movement as well. However, the use of the Langevin dynamics approach or the like significantly increases the computational time required for carrying out long simulations to study the long-time dynamics of the mobile ion atmosphere around biomolecules. On the other hand, it is reasonable to ignore the undesirable effects due to the occasional close encounters between two ions or between an ion and the polyalanine, since this event happens very infrequently because the number densities of ions are extremely low for physiological salt concentrations. We are using this approximation for the data presented in this paper. One can also simply model the averaged effects of the repulsive potentials in a Brownian

15282 J. Phys. Chem., Vol. 100, No. 37, 1996

Yu et al.

TABLE 1: Summary of Simulationsa type

C (M)

Rc (Å)

time (µs)

I I I H H H E E

0.1 0.1 0.5 0.1 0.1 0.5 0.1 0.5

20 40 40 20 40 40 40 40

1.0 1.0 0.1 1.0 1.0 0.1 1.0 0.1

1 2 3 4 5 6 7 8

a I: Pure salt simulations. H: Simulations with an R-helix. E: Simulations with an extended polypeptide. C: ion concentrations in M. Rc: Cutoff for Coulombic interactions. time: Total simulation length after equilibration.

dynamics simulation. The repulsive potentials mainly account for the fast librational motion of atoms in their solvent cages but are not responsible for the diffusive motion of the atoms. This fact is quite clear from molecular dynamics simulations of Lennard-Jones fluids.18 Therefore, to be consistent with the treatment of the averaged effects of the solvent, the repulsive potentials used in this work can actually be removed and replaced by some averaged effects, if one is not interested in the dynamics occurring faster than the ∼40 fs time scale. Anyhow, since the artifacts introduced by the repulsive potentials are so small in the applications presented here because of the low number densities of ions used in the simulations, these artifacts can be safely ignored in this study. In the simulations with the polyalanine, the polyalanine is placed in the center of a simulation box. In principle, one should also let the polyalanine undergo translational and rotational Brownian motion. However, the larger polyalanine diffuses much slower than the sodium and chloride ions so that it is a reasonable approximation to let the polyalanine fixed during the simulations. This makes it a bit easier to analyze the distribution of ions around the polyalanine. This approximation works even better for studying larger proteins and DNAs. If desired, it is straightforward to introduce translational and rotational motion of biomolecules in a Brownian dynamics simulation algorithm. In this paper, we present results from eight different simulations listed in Table 1. Simulations 1-6 are done by using a cubic simulation cell of dimension 100 × 100 × 100 Å3. For simulating a system with a 0.1 M sodium chloride concentration, 61 Na+ and 61 Cl- are put in the simulation cell. 301 of each type of ions are put in the cell when the sodium chloride concentration is 0.5 M. In simulations 7 and 8, a rectangular cell of dimension 83.3 × 128.0 × 82.5 Å3 is used. In simulation 7, 53 pairs of ions are put in the simulation cell to give a sodium chloride concentration of 0.1 M. In simulation 8, there are 265 pairs of ions in the cell to make up a sodium chloride concentration of 0.5 M. Periodic boundary conditions are used in all the simulations. In all the simulations with a 0.1 M sodium chloride concentration, the first 6 ns of simulation data are discarded for thermal equilibration. In all the simulations with a 0.5 M sodium chloride concentration, the first 0.6 ns data are discarded. Data are collected every 20 time steps (equivalent to 0.8 ps in real time) for analyses. III. Results A. Ion Pair Radial Distribution Functions. Figure 1a shows the Brownian dynamics simulation results of Na+-Na+, Cl--Cl-, and Na+-Cl- ion pair radial distribution functions of a 0.1 M NaCl solution using a long-ranged nonbonded cutoff of 20 Å (simulation 1). The mean-field Debye-Hu¨ckel (DH) results are also shown for comparison. The result are obtained by averaging over 20 ns of simulation data after equilibration.

Unlike the results of ion probabilities presented below, nanosecond simulations are sufficient to give satisfactory ion pair radial distribution functions. The statistical errors of the radial distribution functions presented here are less than 1% of the simulated values. Figure 1b is the same as Figure 1a except that a long-ranged nonbonded cutoff of 40 Å is used. The counterion and co-ion radial distribution function curves change less smoothly near the long-ranged nonbonded cutoff value. This is more severe when a 20 Å nonbonded cutoff is used instead of a 40 Å nonbonded cutoff. This is obviously an artifact resulting from the use of finite nonbonded cutoffs for the long-ranged interactions. We have tested the use of several nonbonded cutoffs ranging from 10 to 80 Å. It is necessary to use a nonbonded cutoff of g30 Å in order to give reasonably smooth ion pair radial distribution functions. If the potential of mean force between two ions is approximated by a sum of the repulsive part of the Friedmantype potential and a screened Coulomb potential of the DebyeHu¨ckel (DH) type

( )

9 zizje2 r* i + r* j U(r) ) Bij + exp(-br) r r

(4)

in which the constant b is determined from the equation

b2 )

4πe2 kBT

∑i Nizi2

where Ni is the average number of ions of kind i in unit volume in the solution, an ion pair radial distribution function can be obtained from this mean-field potential as

( )

g(r) ) exp -

U(r) kBT

(5)

The ion pair radial distribution functions calculated this way are also shown in Figure 1a and Figure 1b for comparison with the simulation results. It is clear that the results obtained by using the 40 Å long-ranged nonbonded cutoff agree better with the predictions from the Debye-Hu¨ckel theory than the results obtained by using the 20 Å long-ranged nonbonded cutoff. It is also apparent that even at this relatively low concentration of NaCl (0.1 M), the Debye-Hu¨ckel theory does not work as well at close distances. There seems to be an overscreening of ionpair interactions at close distances by the Debye-Hu¨ckel theory that the peak of the counterion radial distribution curve is underestimated whereas those of the coion radial distribution functions are overestimated at close distances. The Debye-Hu¨ckel theory badly fails for a 0.5 M NaCl solution as shown in Figure 1c in which the Brownian dynamics simulation results are obtained by using a 40 Å long-ranged nonbonded cutoff. The theory fails the most at short distances where it substantially underestimates the peak of the counterion radial distribution function curve. To further validate the Brownian dynamics simulation results, we also perform Metropolis Monte Carlo simulations for systems that have a 0.1 M sodium chloride concentration. The simulation cell is still 100 × 100 × 100 Å3. Each atom is randomly displaced according to a uniform distribution over the range (0, 2) Å in the x, y, and z directions. Since the concentration of sodium chloride is very dilute, the acceptance ratio of the trial moves is about 90%. The first 6.25 × 107 steps of the Monte Carlo simulations are discarded to allow for thermal equilibration and configurations taken from every 20 steps of the next 6.25 × 108 steps are used for analysis. The results in Figure 2, a and b, compare the ion pair radial

Brownian Dynamics Simulations of Polyalanine

J. Phys. Chem., Vol. 100, No. 37, 1996 15283

Figure 1. Na+-Na+, Cl--Cl-, and Na+-Cl- radial distribution functions from the Debye-Hu¨ckel theory (DH) and Brownian dynamics simulations of pure sodium chloride solutions (BD): (a, top) concentration ) 0.1 M, cutoff ) 20 Å; (b, middle) concentration ) 0.1 M, cutoff ) 40 Å; (c, bottom) concentration ) 0.5 M, cutoff ) 40 Å.

15284 J. Phys. Chem., Vol. 100, No. 37, 1996

Yu et al.

Figure 2. Na+-Na+, Cl--Cl-, and Na+-Cl- radial distribution functions from Monte Carlo simulations of a 0.1 M pure sodium chloride solution (a, top) with the Coulombic potential truncated at 40 Å, and (b, bottom) when the minimum image convention is used.

distribution functions obtained from the Monte Carlo simulations using two different schemes for truncating the Coulombic potential: (a) the Coulombic potential is truncated at 40 Å, and (b) the minimum image convention is used. From Figure 2a, it is easy to see the artifacts introduced by truncating the Coulombic potential at 40 Å. The ion pair radial distribution functions change drastically near the cutoff value for the Coulombic potential and these artifacts move to a region near 20 Å if the Coulombic potential is truncated at 20 Å. Clearly, these artifacts are associated with the choice of the long-range electrostatic cutoff. When the minimum image convention is used, the artifacts due to the use of a finite electrostatic cutoff are substantially reduced (see Figure 2b). The artifacts resulting from the use of finite spherical cutoffs for the long-ranged electrostatic interactions are, however, much less severe in the Brownian dynamics simulations. The more serious artifacts in

the Monte Carlo simulations probably result from not using information provided by the interaction forces. In the Metropolis Monte Carlo simulations, random coordinate displacements are made and the energies of the new configurations are subsequently checked to decide whether the configurations are kept or discarded. Moves that bring two oppositely-charged ions that are initially separated by a distance larger than the long-ranged nonbonded cutoff into a distance smaller than the cutoff will always be accepted as in the case when no cutoff is used. However, moves that separate two oppositely-charged ions from a distance smaller than the cutoff to one larger will have a smaller acceptance ratio than when no cutoff is used. This results in the overestimation of the counterion radial distribution function for distances that are slightly less than the cutoff but the underestimation of the same radial distribution function for distances that are slightly larger than the cutoff.

Brownian Dynamics Simulations of Polyalanine

J. Phys. Chem., Vol. 100, No. 37, 1996 15285

TABLE 2: Maximum Ion Probabilities in Different Simulations type

C (M)

Pmax(Na+)

Pmax(Cl-)

H H E E

0.1 0.5 0.1 0.5

33 ( 2 12 ( 2 10 ( 1 5(1

3.7 ( 0.8 6(2 4(1 5(1

Similar arguments apply to the explanation of the opposite effects observed for the co-ion radial distribution functions. The artifacts in the Brownian dynamics simulations are less severe probably because the information provided by the interaction forces are directly used in guiding the prediction of the new positions of the ions at each time step. B. Ion Probability Distributions around Polyalanine. In a pure salt solution, the probability of finding a specific type of ion in any location is the same. This is, however, no longer true when a polyalanine is introduced into the solution. The potential field around the polyalanine may also influence the probability distributions of positively and negatively-charged ions in different ways. We have examined the distributions of ions around a 13-residue polyalanine at two different conformations and at two different concentrations of NaCl. A normalized ion probability distribution function P(r) for finding an ion at position r is calculated as

P(r) )

n(r) V(r)‚(N/V)

(6)

Here V(r) is the volume of a tiny cube centered at r, n(r) is the average number of ions found in the cube, N is the total number of ions in the simulation cell, and V is the volume of the cell. P(r) ) 1 for r that is far away from the polyalanine. P(r) ) 3 means the probability of finding an ion at r is three times higher than the probability of finding the same ion in a region far away from the polypeptide. The results presented in this section are obtained from Brownian dynamics simulations of a cubic box of sodium and chloride ions with a polyalanine immersed in it. Periodic boundary conditions are applied as in the simulation of the pure NaCl solutions. Two ion concentrations are used: 0.1 and 0.5 M. In calculating the ion probability distributions, the tiny cubes have dimensions of 2 × 2 × 2 Å3. Since the number densities of ions are very small, very long simulations need to be carried out to get good statistics for calculating these ion probability distributions. For the 0.1 M solutions, simulations of 1.0 µs durations, in addition to 6 ns of equilibration, are carried out. The 0.5 M simulations do not need to be carried out as long to obtain comparable statistics as the 0.1 M simulations since there are 5 times more ions in the simulation cell than before. Only 0.1 µs simulations with 0.1 ns of equilibration are carried out for systems having a 0.5 M sodium chloride concentration. The results presented in this section are obtained by using a longranged nonbonded cutoff of 40 Å. The introduction of a 13-residue polyalanine significantly distorts the ion distributions near the surface of the polyalanine. The maximum values of the ion probabilities for the cation and the anion are listed in Table 2. Figure 3 shows the Na+ and Cl- probability distributions for a 0.1 M sodium chloride solution in which an R-helical polyalanine is immersed. The results for two slices cutting through the helix are shown. One slice cuts through the long axis of the helix. The other cuts through the center and is perpendicular to the long axis of the helix. In the figure, the darkness of each small square is related to the probability of finding a specific type of ion in a 2 × 2 × 2 Å3 cube associated with the square and is scaled between

Figure 3. Na+ and Cl- ion probability distributions P(r) from the simulation of a 0.1 M sodium chloride solution with R-helical polyalanines introduced. The xy slice is taken at z ) 0 with a thickness of 2 Å. The xz slice is taken at y ) 0 with a thickness of 2 Å. x and z range from -30 to 30 Å. y ranges from -40 to 40 Å. In the xy + slice, the maximum P values are Pmax ) 32 ( 2, and Pmax ) 2.1 ( + 0.5. In the xz slice, the maximum P values are Pmax ) 15 ( 2 and ) 2.2 ( 0.5. Pmax

[0, 1] according to the following rule:

darkness )

{

1.0 if P g 6.0 P/6.0 if P < 6.0

(7)

The darkest squares have darkness ) 1 and the highest P value for each kind of ion in each slice is shown in the figure caption. From Table 2 and Figure 3, one can see that the R-helix distorts the Na+ probability distribution much more than the Cl- probability distribution. The highest P value for the Na+ probability can be as high as 33, whereas the corresponding value for Cl- is only 3.7. Figure 3 also shows that, except at the N-terminal end of the helix, there is a higher probability of finding Na+ around the helical surface. The thickness of this ion “cloud” is about 4-6 Å. At the positively-charged N-terminus of the helix, the probability of finding Na+ is small, even smaller than 1.0. For the Cl-, there is also a slightly higher probability of finding the ions near the surface of the helix, but this preference is much less significant when compared to the case of the Na+. Analogous to the case of the cations, there is also a diminished probability of finding anions near the negatively-charged C-terminal end of the helix. The degree of inhomogeneity of the distributions of ions introduced by the presence of the polyalanine is strongest when the polyalanine is in the R-helical conformation and when the ion concentration is small. The effects are also more significant on the cations than on the anions. The preferential binding of cations by the R-helical polyalanine is probably due to the favorable interactions between the cations and the relativelyexposed carbonyl oxygens of the peptide groups. In the R-helical conformation, the amide hydrogens are buried and thus

15286 J. Phys. Chem., Vol. 100, No. 37, 1996

Figure 4. Same as Figure 3, but for a 0.1 M sodium chloride solution with the polyalanine in the extended form. x and z range from -40 to 40 Å. y ranges from -60 to 60 Å. In the xy slice, the maximum P + values are Pmax ) 7.1 ( 0.8 and Pmax ) 4 ( 1. In the xz slice, the + maximum P values are Pmax ) 6.2 ( 0.8, Pmax ) 2.4 ( 0.6.

may not interact as well with the anions. The degree of the heterogeneity of ion distributions is decreased when the polyalanine is in the extended conformation in which both the N-H and the C-O groups are exposed. Figure 4 presents the simulation results of the extended form of polyalanine at 0.1 M sodium chloride concentration (simulation 7 in Table 1). The corresponding results for a 0.5 M sodium chloride solution are also shown in Figure 5. At 0.1 M NaCl, there is still a relatively high probability of finding Na+ around the polypeptide surface, especially near the negatively-charged C-terminal end, but this preference is much less significant when compared to the case of the R helix. The Na+ “cloud” is also thinner than that in the helical case. At 0.5 M NaCl, although there is still a slightly higher probability of finding each type of ions near the polypeptide, the degree of heterogeneity of ion distributions introduced by the presence of the polypeptide is about the same for the cations and the anions. This is probably due to the larger screening effects of the interactions between the polypeptide and the ions at a higher salt concentration. This larger screening effect is also apparent in the R-helical case. The inhomogeneity of ion distributions around the R-helical polyalanine is also diminished when the concentration of the salt is increased to 0.5 M from 0.1 M (Figure 6). To show that the differences in the distributions of the two ions are not due to the difference in the size of the two ions, we have also carried out similar Brownian dynamics simulations and analyses by artificially changing the size of the chloride ion to that of the sodium ion and vice versa. The major differences in the distributions of the ions remain the same when these size changes are made, suggesting that the differences in the distributions of the cations and the anions introduced by the polyalanine do not result from the difference in the sizes of the two types of ions.

Yu et al.

Figure 5. Same as Figure 4, but the salt concentration is 0.5 M. In + the xy slice, the maximum P values are Pmax ) 3.9 ( 0.8, Pmax ) 4.2 + ( 0.9. In the xz slice, the maximum P values are: Pmax ) 4 ( 1, Pmax ) 3.6 ( 0.8.

Figure 6. Same as Figure 3, but the salt concentration is 0.5 M. In + the xy slice, the maximum P values are Pmax ) 11 ( 1, Pmax ) 3.4 ( + 0.8. In the xz slice, the maximum P values are Pmax ) 7 ( 1, Pmax ) 4 ( 1.

C. Fluctuations and Dynamical Properties. We have calculated the collective dipole moment of the ions and its fluctuation. The fluctuation of the collective dipole moment is related to the ionic polarizability of the system. The instanta-

Brownian Dynamics Simulations of Polyalanine

J. Phys. Chem., Vol. 100, No. 37, 1996 15287

Figure 7. Time autocorrelation functions 〈µi(0)µi(t)〉/〈µi(0)µi(0)〉 of the components of the collective dipole moment of ions in the x, y, and z directions for a pure 0.1 M sodium chloride solution with the Coulombic potential truncated at 20 and 40 Å.

TABLE 3: Collective Dipole Moment of Ions and Its Fluctuation in Different Simulations type C (M) Rc (Å) I I I H H H E E

0.1 0.1 0.5 helix 0.1 0.1 0.5 extend 0.1 0.5

20 40 40 20 40 40 40 40

µx (e Å)

µy (e Å)

µz (e Å)

R (106 e2 Å2/kBT)

-4 ( 24 -4 ( 12 -9 ( 24 -4.18 -1 ( 14 7 ( 10 8 ( 23 0.0 2(6 5 ( 22

8 ( 19 8(9 4 ( 23 -25.26 15 ( 17 20 ( 10 22 ( 36 -46.4 21 ( 13 16 ( 41

-2 ( 23 -9 ( 11 -4 ( 31 3.49 -2 ( 18 4(8 13 ( 27 -0.47 5(9 -5 ( 22

12.8 ( 0.3 9.5 ( 0.2 120 ( 4 12.7 ( 0.4 9.4 ( 0.2 123 ( 4 7.0 ( 0.2 92 ( 4

neous collective dipole moment µ(t) of the ions at time t is calculated as

µ(t) )

1

∑i eziri

N

(8)

Here the sum runs over the total number of ions N in the system, zi is the valence of the i-th ion, and ri is the coordinate vector of ion i. The ionic polarizability of the solution can be obtained from the fluctuation of the collective dipole moment as

R)

〈µ2〉 - 〈µ〉2 kBT

(9)

where µ is the collective dipole moment and 〈 〉 denotes an ensemble average. Table 3 shows the calculated collective dipole moments and the ionic polarizabilities from our simulations. Even with the long simulation carried out here, we find that it is difficult to get good statistics in calculating the collective dipole moments.

However, the ionic polarizabilities can be calculated significantly more precisely. From Table 3, one can see the large effects of truncating the Coulombic interactions on the simulation results. Truncating the Coulombic interactions at short distances overestimates the ionic polarizabilities. The R-helical polyalanine does not affect the ionic polarizability much, whereas the extended polyalanine appears to diminish the ionic polarizability of the system, at both 0.1 and 0.5 M sodium chloride concentrations. Since the ionic polarizability is sensitive to conformational changes of polypeptides, and possibly proteins, it may be possible to deduce structural information of polypeptides and proteins by measuring the static dielectric constants of these molecules in salt solutions. In DNAs,7,8 it has been shown that more detailed structural information can be obtained by measuring the frequencydependent dielectric properties of DNAs in salt solutions. To examine whether similar experiments may be used to probe the conformation of polypeptides and proteins, we have analyzed the time autocorrelation function of the collective dipole moment of the ions, which determines the frequency-dependent dielectric properties of an ionic solution. We first examine the sensitivity of the time autocorrelation function of the collective dipole moment to the choice of coulumbic cutoff. In Figure 7, the time autocorrelation functions of the components of the collective dipole moment 〈µi(0)µi(t)〉 (where i ) x, y, and z) are calculated using two different values of the Coulombic cutoff: 20 and 40 Å. The results are obtained from trajectories of 1.0 µs duration. Consistent with a previous dynamical simulation,19 the use of a short Coulombic cutoff overestimates the relaxation time of the collective dipole moment, again suggesting the importance of using large coulombic cutoff in estimating such a dynamical property. Figures 8 and 9 compare the time autocorrelation functions of the components of the collective dipole moment of ions in

15288 J. Phys. Chem., Vol. 100, No. 37, 1996

Yu et al.

Figure 8. Time autocorrelation functions 〈µi(0)µi(t)〉/〈µi(0)µi(0)〉 of the components of the collective dipole moment of ions in the x, y, and z directions for a pure 0.1 M sodium chloride solution and a similar solution with R-helical polyalanines immersed. The Coulombic potential is truncated at 40 Å.

Figure 9. Time autocorrelation functions 〈µi(0)µi(t)〉/〈µi(0)µi(0)〉 of the components of the collective dipole moment of ions in the x, y, and z directions for a pure 0.1 M sodium chloride solution and a similar solution with extended polyalanines immersed. The Coulombic potential is truncated at 40 Å.

Brownian Dynamics Simulations of Polyalanine solutions of polyalanine in the R-helical and the extended conformation. Again, one can see that the time autocorrelation function is not affected much by the R-helical polyalanine. The extended polyalanine, on the other hand, alters the behavior of the time autocorrelation function. The decay of the time autocorrelation function is slower along the direction of the longaxis of the extended polyalanine but faster along the other two orthogonal directions. Since this time autocorrelation function is sensitive to the conformation of the polypeptide, frequencydependent dielectric measurements may also yield useful insights into the conformation of polypeptides or proteins, as in the case of DNAs. IV. Discussion and Conclusions The simulations presented here suggest that the Brownian dynamics simulation method can be a useful alternative tool for studying the ion distributions around biomolecules. Besides providing an equilibrium description of the ion distributions around a biomolecule, it can also provide dynamical information. We have examined an important issue concerning the simulation of salt solutions, that is, the sensitivity of simulation results to the choice of long-ranged nonbonded cutoff. We find that the equilibrium, fluctuation, and dynamics properties of a salt solution are quite sensitive to the choice of long-ranged nonbonded cutoff. It seems that a nonbonded cutoff >30 Å is essential for obtaining reliable simulation results. This value is much larger than the nonbonded cutoffs typically used in molecular dynamics simulations of biomolecules. To improve molecular dynamics simulations, it is not sufficient to simply add counterions to the simulation models but it also seems necessary to use a large enough nonbonded cutoff for the longranged electrostatic interactions. We have also observed that polyalanine, and probably proteins in general, have differential binding affinity to different types of ions. The simulations carried out here, for example, suggest that polyalanine preferentially bind sodium over chloride ions. The preferential binding is dependent upon the conformation of the polypeptide and the ion concentrations. This finding is consistent with the experimental observations of proteins by Timasheff and co-workers.1,2 Timasheff and co-workers1,2 also suggest that there is a correlation between the ability of an ion to stabilize a protein and the distribution of the ion around the protein. An interesting question to address with Brownian dynamics simulations is whether the distributions of ions around proteins obtained from these simulations can be correlated with the ability of the ions to stabilize proteins, consistent with the

J. Phys. Chem., Vol. 100, No. 37, 1996 15289 Hofmeister’s series. Since our simulations also show that the ion distributions are dependent upon the conformation of a polypeptide, one may be able to suggest which conformations of a protein may be preferentially stabilized by a specific type of ions by simply analyzing the distribution of this type of ions around the protein, if there is indeed a correlation between ion distributions and protein stability. A suggestive study using the extended RISM theory20 has already indicated that a salt can selectively stabilize certain conformations of a polypeptide but destabilize others. This ability of salts to limit the accessible conformational space of a polypeptide has been suggested as a possible mechanism for facilitating protein folding.20 We have also illustrated that similar to DNA solutions, the frequency-dependent dielectric property of a polypeptide or protein can be sensitive to the conformation of the polypeptide. This finding suggests that one may also use dielectric measurements to study conformational changes of polypeptides and proteins. Acknowledgment. This work has been supported in part by the NIH (C.F.W.) and the NSF(J.Z.). References and Notes (1) Arakawa, T.; Bhat, R.; Timasheff, S. N. Biochemistry 1990, 29, 1914. (2) Arakawa, T.; Bhat, R.; Timasheff, S. N. Biochemistry 1990, 29, 1924. (3) Honig, B.; Sharp, K.; Yang, A.-S. T. J. Phys. Chem. 1993, 97, 1101. (4) Davis, M. E.; McCammon, J. A. Chem. ReV. 1990, 90, 509. (5) Warshel, A.; A° quist, J. Annu. ReV. Biophys. Biophys. Chem. 1991, 20, 267. (6) See, for example: Anderson, C. F.; Record, M. T., Jr. Annu. ReV. Biophys. Biophys. Chem. 1990, 19, 423, and references therein. (7) Saif, B.; Mohr, R. K.; Montrose, C. J.; Litovitz, T. A. Biopolymers 1991, 31, 1171. (8) Bone, S.; Small, C. A. Biochim. Biophys. Acta 1995, 1260, 85. (9) Manning, G. S. J. Chem. Phys. 1969, 51, 924. (10) Manning, G. S. Q. ReV. Biophys. 1978, 11, 179. (11) Van der Touw, F.; Mandel, M. Biophys. Chem. 1974, 2, 218. (12) Mandel, M. Ann. NY Acad. Sci. 1977, 303, 74. (13) Rossky, P. J.; Dudowicz, J. B.; Tembe, B. L.; Friedman, H. L. J. Chem. Phys. 1980, 73, 3372. (14) Friedman, H. L. Annu. ReV. Phys. Chem. 1981, 32, 179. (15) van Gunsteren, W. F. GROMOS 1987. (16) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. (17) Weast, R. C., Ed. Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1975; Vol. 56. (18) See, for example: Rahman, A. Phys. ReV. A 1964, A136, 405. (19) Beglov, D. B.; Lipanov, A. A. J. Biomol. Struct. Dyn. 1991, 9, 205. (20) Perkyns, J. S.; Pettitt, B. M. J. Phys. Chem. 1995, 99, 1.

JP960124R