Temperature Effects on Protein Motions: A Molecular Dynamics Study

The dynamics of the enzyme ribonuclease-Sa as a function of temperature has been explored through a series of molecular dynamics simulations. Long-ran...
0 downloads 0 Views 233KB Size
Download PDF
J. Phys. Chem. B 2000, 104, 10387-10397

10387

Temperature Effects on Protein Motions: A Molecular Dynamics Study of RNase-Sa Radovan Dvorsky and Josef Sevcik Institute of Molecular Biology, DubraVska cesta 21, 842 51 BratislaVa, SloVak Republic

Leo S. D. Caves, Roderick E. Hubbard, and Chandra S. Verma* Structural Biology Laboratory, Department of Chemistry, UniVersity of York, YORK YO10 5DD, UK ReceiVed: May 26, 2000; In Final Form: September 14, 2000

The dynamics of the enzyme ribonuclease-Sa as a function of temperature has been explored through a series of molecular dynamics simulations. Long-range expansion and short-range contraction of the structure lends the protein a solidlike core and a liquidlike exterior as the temperature is increased. The trend in magnitudes of fluctuations of atoms are biphasic across the 150-200K region and are increasingly non-Brownian in character as the temperature is increased. The mobility of solvent molecules is much higher than the protein atoms, even though the solvent mobility displays behavior which is dampened relative to behavior in bulk water. The region of the active site that binds the base of the nucleotide ligand shows low plasticity relative to regions that interact with the sugar-phosphate part. This suggests that the enzyme is preorganized dynamically: regions with low plasticity confer specificity while the more flexible regions have the fluidity to facilitate energetically inexpensive conformational rearrangements such as those required to achieve the transition state. Below 200K the dynamics are characterized by low amplitude harmonic motions that involve concerted motions that involve small groups of atoms. Above 200K, the dynamics are dominated by large amplitude anharmonic motions which involve long-range correlations including breathing and twisting modes such as those required for ligand binding/release/activation. The temperature-dependent transition in the character of the dynamics at ∼200K reflects the ease with which the system hops among barriers giving rise to enhanced diffusion across phase space. This enhanced plasticity is catalyzed by a significant increase in the mobility of solvent water molecules and the associated increase in frequency of different hydrogen bond arrangements and may facilitate the onset of significantly enhanced functionality above the transition temperature.

I. Introduction Many proteins (e.g., enzymes, ion pumps, transport proteins) require flexibility for their biological function. This flexibility arises from a balance between the intra and intermolecular forces and constant buffeting from solvent molecules1,2 which makes it temperature (and/or viscosity) dependent. The interplay between the large number of degrees of freedom inevitably leads to a wide variety of highly complex relaxation properties.3,4 One of these properties is the behavior of proteins as a function of temperature. It has been demonstrated in several experiments that there exists a temperature around 180-220 K below which the protein exhibits a marked reduction or complete loss of functionality.5,6 The general consensus is that below this “protein-transition” temperature, the protein conformations are kinetically trapped in distinct conformational states and hence cannot interconvert as rapidly as above this temperature. Increasingly, the evidence seems to suggest that it is the increased viscosity of the solvent which traps the protein in minima.7,8 Recent theoretical studies have shown that the presence of a mobile solvent is necessary for the onset of this phenomenon.9,10 Details of this feature have emerged from crystallographically derived B values11 or from variants of neutron scattering12 and more recently from sophisticated * Corresponding author. E-mail: [email protected]. Fax: 441904-410519. Tel: 44-1904-434532.

spectroscopic methods.5,13-17 However, the exact mechanisms underlying this phenomenon are still not clear.18 Computational studies, while limited in time scales, can often complement experimental data and provide an elegant method to investigate details of these processes at the atomic level. While this approach is still controversial owing to the approximations involved, encouraging agreements with experiments are often reported.4,9,10,19-24 To investigate the onset of this proteintransition in the crystalline state of proteins, we have performed a series of simulations using ribonuclease-Sa (RNase-Sa) as a test case. RNase-Sa is a small, 96 amino acid globular protein isolated from the bacterial strain Streptomyces aureofaciens and cleaves single-stranded RNA molecules at their 3′,5′ linkage. The protein is active in the pH range from 4.3 to 9.7 with a maximum at 7.0.25 Crystal structures have been determined of the protein in its native and complexed states26,27 revealing aspects of its functionality. In this report, we examine the dynamics of this enzyme in its pseudo-crystalline environment as a function of temperature. We performed six simulations, one each at 50, 100, 150, 200, 250 and 300 K; each simulation was 1.2 ns in length. We first examine some overall properties to judge the stability of the simulations. This is followed by a comparison of the computed equilibrium properties with those determined from X-ray crystallographic studies of the effects of temperature on RNase-Sa27 and on RNase-A.28 We then focus on some of the

10.1021/jp001933k CCC: $19.00 © 2000 American Chemical Society Published on Web 10/19/2000

10388 J. Phys. Chem. B, Vol. 104, No. 44, 2000

Dvorsky et al.

Figure 1. Stereoview of the secondary structure of ribonuclease-Sa. The protein is made up of one R-helix (residues 14-22) and five β-strands (residues 3-8, 52-57, 68-73, 78-82, 89-93). The active site residues that interact with the nucleotide are highlighted with their numbers. The figure was drawn using Molscript.95

dynamical properties of the protein and comment on the relevance of our results to enzyme function, often taking cues from studies that have been carried out on glassy substances other than proteins.29 II. Methods The crystal structure of RNase-Sa at room temperature27 (Protein Data Bank accession code 1rgg; resolution 1.2 Å; Figure 1) was used as the starting model for our calculations. The crystalline environment was generated by applying the symmetry operations of space group P212121.27 The available computational resources allowed us to include only residues within 14 Å of the central molecule, as opposed to generating the unit cell and using periodic boundaries to replicate the system.9,11 Our setup also included 657 crystallographically observed water molecules and 4 sulfate ions. Additionally, 701 water molecules were added from a preequilibrated water sphere to fill the voids between the symmetry generated molecules since waters of hydration are thought to be central to the temperature-dependent transition observed in the motions of proteins.6,9,10,15 The final system consisted of the central molecule (96 residues) surrounded by 220 residues originating from 13 symmetry-related neighboring protein molecules and 1358 water molecules, leading to ∼412 water molecules per protein. The CHARMM polar hydrogen force field30 was employed to represent the system, along with a modified TIP3 model31,32 for water. The initial structure was assumed to contain neutral His (with the proton on atom Nδ1), fully charged Asp, Glu, Arg, Lys and the two terminii and dianionic sulfate anions. The total system consisted of 7030 atoms. The nonbonded electrostatic interactions were shifted to zero at 14 Å with a dielectric constant of 1 and the van der Waals interactions were truncated with a switching potential operating between 8 and 12 Å.33 The nonbonded lists were generated in each step during the minimizations and updated every 0.05 ps during the molecular dynamics (MD) simulations. The crystal structure was initially subjected to 100 steps of steepest descent and 200 steps of adopted basis NewtonRaphson energy minimizations30 to remove unfavorable contacts. MD simulations were started from this minimized structure by assigning random velocities to the atoms from a Gaussian distribution, corresponding to mean temperatures of 50-300 K. To maintain a constant temperature the system was coupled to an external bath34 with a coupling constant of 1 ps. The equations of motion were integrated using a time step of 0.001 ps for a total time period of 1.2 ns in each system and data was collected every 0.2 ps for analysis. All the calculations were performed on a Pentium II 266 MHz CPU running under the Linux system.

Figure 2. Time series of RMS positional deviations from the initial structure of all heavy atoms of RNase-Sa for the simulations at different temperatures.

III. Results The root-mean-square (RMS) deviation of the atomic positions of the minimized structure relative to the crystal structure was 0.42 Å for the backbone atoms and 0.76 Å for the side chain atoms. The largest changes were in solvent exposed regions. This parallels the difference between the crystal structures determined at room temperature and at 100 K where the corresponding RMS differences are 0.5 Å and 1.4 Å respectively. During the MD simulations, the atomic positions deviate from the minimized structure, stabilizing at ∼0.7-0.94 Å (Figure 2). The time series of the radius of gyration (data not shown), an indicator of global motion, also shows a trend similar to that in Figure 2; the values stabilize at ∼12.2-12.5 Å (Table 1). In both analyses, we see that the fluctuations increase with temperature, particularly beyond 200 K, indicative of increasingly higher degrees of conformational sampling. Together, the behavior of these two properties indicate stable simulations. Coordinate Averages and Fluctuations. The RMS deviation of the MD-average structure from the starting structure (Table 1) range from 0.6 to 1.0 Å for the backbone heavy atoms and 0.7-1.2 Å for all heavy atoms. These values are similar to values reported in other simulations carried out in the crystal environment.9,11,35,36 The average values of the radius of gyration from the simulations compare well with values calculated directly from the crystal structures of Rnase-Sa (Table 1) and the changes are similar to those reported for RNase-A.28 This increase in RMS deviations and radius of gyration with increasing temperatures suggests that the structure expands with increasing temperature. This expansion can be seen in the differences between the distribution of interatomic distances at different temperatures (Figure 3a,b). The change between 100 and 300 K is very similar to that observed for RNase-Sa (Figure 3a) and for RNase-A.28 While the trend of these differences for

Temperature Effects on Protein Motions

J. Phys. Chem. B, Vol. 104, No. 44, 2000 10389

TABLE 1: Averages Derived from the Simulations of RNase-Sa at Different Temperatures temperature (K) RMS deviation (Å)a radius of gyration (Å)b

50 100 0.64 0.67 12.26 12.27 (12.4) 〈x〉 backbone atoms (Å)c 0.15 0.19 〈x〉 side chain atoms (Å)c 0.19 0.25 Lindemann parameterd 0.04 0.05 0.034 0.044 backbone 〈φ〉 (deg) 4.74 6.55 backbone 〈Ψ〉 (deg) 4.49 6.19 skewnesse 0.113 0.117 covariancef 0.22 0.26 0.32 0.45 〈x〉 water oxygens (Å)c,g D water (Å2 ps-1)h 0.003 0.007 187 185 intraprotein hbonds (no.)i average lifetimej 298 204 protein-water H bonds (no.)i 157 100 194 75 average lifetimej

150 200 250 0.63 0.65 0.75 12.35 12.35 12.4 0.23 0.30 0.06 0.051 7.37 6.93 0.119 0.29 0.65 0.009 163 144 114 44

0.28 0.36 0.07 0.062 8.63 8.24 0.161 0.58 0.90 0.019 175 99 82 19

0.31 0.42 0.089 0.068 9.43 8.99 0.168 0.68 1.29 0.043 167 59 65 15

300 0.85 12.4 (12.6) 0.34 0.46 0.094 0.069 10.00 9.50 0.182 1.0 1.76 0.084 151 54 60 11

a Average structures from simulations relative to the initial starting structure. b Values in parentheses have been calculated from the crystal structures. c 〈x〉 is the average root-mean-squared fluctuation. d The values in the top row have been calculated for all atoms; the values in the bottom row have been calculated for atoms within 4.5 Å of the centroid of the protein. e Averages of absolute values of the third moment of the distribution of atomic fluctuations. f Maximum of the absolute values of the covariances. g Computed for waters within a 6 Å shell of the central protein molecule. h Diffusion coefficient of water. i Number of hydrogen bonds with lifetime greater than 100 ps. j Average lifetime of hydrogen bond computed at 2 ps intervals.

each temperature increment is much more complex (Figure 3b), the general pattern is one of contraction at small distances (originating in increased local packing density) and expansion at larger distances. An interesting feature of the 150-200 K region is that it is relatively flat indicative of a state where some sort of phase transition is occurring. Differences between these distance difference matrices reveal an expansion in ∼60% of the protein while the contraction is distributed over ∼30% of the protein as the temperature is increased. The mean fluctuations of the backbone atoms and side chain heavy atoms increase with temperature (Figure 4a) with the latter, on average, 27-35% larger than the former. This temperature dependence is biphasic and is similar to that reported in experimental28 and theoretical studies.9,37-40 The increase is not as sharp as in solution studies9 due to the pseudocrystalline environment used in our work. The biphasic behavior is more pronounced for the side chain atoms relative to the backbone atoms with the change in slopes of the two phases increasing by 33% for the backbone atoms and 60% for the side chain atoms. The distributions of fluctuations becomes broader and asymmetric with the maxima shifting to higher values with temperature (Figure 4b). A similar behavior was reported from studies on RNase-A28 and lysozyme.41 The fluctuations of the solvent (Table 1 entry for water oxygen atoms) are 2-3-fold larger and the increase is much more nonlinear than for the protein atoms. Again, the increase is steeper beyond 150 K than below it. The diffusion coefficients computed for water (Table 1) show a steep increase beyond 200 K. The positional fluctuations of the atoms are higher at the surface of the protein relative to the interior, as one would expect. In fact an inverse relationship was found (data not shown) between the packing as judged from the solvent accessible surface area of an atom and its mobility or its positional fluctuation.42 However, this distribution becomes more diffuse as the temperature increases, suggestive of greater fluidity of atomic motion at elevated temperatures and at

Figure 3. Distributions of interatomic distances: (a) difference between 300 and 100 K simulations and crystals; (b) difference between sequential temperature simulations.

increasing distances from the center of the system. This was further confirmed by computing the Lindemann disorder parameter, which has been used to demonstrate the existence of a solid core and liquidlike exterior in proteins43 (the Lindemann parameter is calculated as x∑i〈∆ri2〉/4.52N where N is 907, the number of atoms in our system and ri is the position of atom i such that ri2 ) (ri - 〈ri〉)2 with 〈 〉 denoting configurational averages). This parameter undergoes a transition between 200 and 250 K (Table 1). The smaller values in our study, compared to the solution studies43 may reflect the constraints of the crystal environment. The small difference between the overall and the internal Lindemann parameter for temperatures below 200 K suggests that below the transition temperature, the protein’s overall character, particularly in the constrained environment of crystals, is largely solidlike. The larger difference at 250 and 300 K suggests that the increased

10390 J. Phys. Chem. B, Vol. 104, No. 44, 2000

Figure 4. (a) Mean squared fluctuations of the backbone and side chain heavy atom positions during the simulations as a function of temperature. (b) Distribution of mean squared fluctuations of the heavy atoms as a function of temperature during the simulations.

liquid behavior at higher temperature arises from increased mobility among the residues that are far from the centroid. This is in accord with the view of the protein transition being linked to the mobility of the solvent9,10,44 since it is the residues distant from the centroid that are, in general, most exposed to solvent. The distribution of the atomic fluctuations across the protein derived from the simulations are in good agreement with those derived from the crystallographic B values (the relationship between the crystallographic B value and the mean square fluctuation of an atom is B value ) (8Π2〈∆r2〉)/3 where 〈∆r2〉 is the mean-square dynamical fluctuation of the atomic position vector r). There is a remarkable agreement in the change in fluctuations between the 100 and 300 K crystal structures as compared to the values derived from the simulation (Figure 5a), particularly in the regions of secondary structure. The largest differences occur toward the C-terminal end of the helix and loop regions and in regions that are involved in nucleotide binding; this inhomogeneous increase in plasticity might lead to the temperature dependent onset of functionality.5,6 An examination of the fluctuations as a function of temperature (Figure 5b) shows that large increases in fluctuations with temperature occur in the loop regions and are also largest in regions which interact with the sugar-phosphate region of the nucleotide relative to the regions that interact with the base: these are particularly clear in the N-terminal half of the loop between R and β2 and in the loops between β2 and β3 and between β4 and β5. There is also a pattern of separation between 50 and 150 K and 200 and 300 K, with the latter showing much higher mobility.

Dvorsky et al. To probe the origin of the larger increase in fluctuations above 200 K than below, we examined the anharmonicity associated with the motions. This was carried out by computing the skewness parameter. This is the third moment of the atomic displacement distribution function and is a measure of deviation from a Gaussian distribution.45 At temperatures below 200 K (Table 1) we find that the motions are predominantly harmonic, with values similar to those reported by others.9,38-40 The change is characterized by a small linear increase in anharmonicity of about 5%. There is a sudden increase of about 35% between 150 and 200 K after which the increase to 300 K is highly nonlinear. The smaller values of skewness at higher temperatures in our study relative to studies of proteins in solution are most probably due to the pseudo-crystalline phase, particularly as we know that the protein transition is modulated by solvent molecules.6,9,10,44 The spread of skewness across the protein backbone (Figure 5c) reveals two features: (a) The anharmonicity is higher in the loop regions and termini of secondary structural elements. (b) Low temperature motion is characterized by little anharmonicity with a sparse, inhomogeneous distribution of the skewness across the protein; in contrast, at temperatures above 200 K, the skewness is not only larger but is also exhibited by many more residues and is spread across the protein, particularly in regions that are in contact with the ligand. The effect of this increasing anharmonicity on the overall dynamics was examined by computing the RMS difference between subsequent structures along each trajectory as a function of the time lag between structures (Figure 5d). This analysis captures the essential diffusive motion of a system and is used to extract the diffusion coefficients.46 Two main features are evident: (a) beyond ∼100 ps, the system behaves in a Brownian manner (fluctuation is proportional to xt where t is time, in the long-term limit (with diffusion coefficients increasing with temperature); (b) in the sub-100 ps time scales, anomalous (nonBrownian) diffusion, indicated by the transition from linear to a power law behavior,3,47,48 sets in from 200 K onward. Behavior of Dihedral Angles. The overall fluctuations of the backbone dihedral angles gradually increase with temperature (Table 1), indicative that increasingly large volumes of phase space are sampled. Although the magnitudes of the difference between the average value of backbone dihedrals during the simulations are larger than the corresponding dihedrals in the crystal structure (Figure 6a), large differences are observed in the same regions of the protein. The mobility of the dihedrals as measured through their transitions between minima,19 increases with temperature, particularly from 200 K onward. McCammon and Karplus49 found that the motions of backbone atoms were primarily responsible for the creation of voids that gave rise to side chain plasticity. In our study, a cross over point occurs around the putative transition temperature where the number of residues undergoing backbone dihedral transitions exceeds the number of residues undergoing side chain dihedral transitions (Figure 6b). While we have not come across any studies that have reported this behavior, this suggests that the emergent functionality at temperatures greater than the transition temperature may be primarily linked to motions of the backbone atoms. Davis and Agard50 have found from an NMR study that the specificity of alpha-lytic protease might be modulated by the motions of the backbone atoms. Examination of the distribution of the number of dihedral transitions across the protein in our study (Figure 6c) shows that among the backbone dihedrals, it is only from 200 K onward that a sizable number of transitions is seen. This is particularly evident in the R-β2 loop, the β2-β3 loop and the β4-β5 loop; all

Temperature Effects on Protein Motions J. Phys. Chem. B, Vol. 104, No. 44, 2000 10391

Figure 5. Distribution across the enzyme of the mean square fluctuations of the backbone atoms: (a) change between 100 and 300 K separately during the simulations and from the crystal; (b) absolute values during the simulations at each temperature during the simulation; (c) distribution of skewness of the atomic motion averaged over the backbone atoms for each residue as a function of temperature; (d) RMS positional differences between subsequent structures as a function of the time lag between subsequent structures during the simulations.

10392 J. Phys. Chem. B, Vol. 104, No. 44, 2000

Dvorsky et al.

Figure 6. (a) Change in backbone dihedral angles between 100 and 300 K from the simulation and from the crystal structure. (b) Number of residues undergoing dihedral angle transitions in the backbone or side chain during the simulations. (c) Distribution of the total number of dihedral angle transitions (backbone and side chain) during the simulations.

these regions are in intimate contact with the nucleotide.51 A similar pattern is seen among the side chain transitions (data not shown). Some dihedral transitions occur below 200 K (in the loops between β2 and β3 and between β4 and β5), which is of interest in the light of the finding that even below the transition temperature, functional motions may be enabled.14 We have seen that both, the amplitudes of motions as well as their spread across the protein, increases with increasing temperature. So we now turn our attention to collective motions of the protein as a whole. For this, in the following sections, we compute the correlations between the characteristic motions of individual residues. Collective Motions. Equal Time Atomic Displacement Correlations. The importance of large scale collective motions in protein function is increasingly being recognized and characterized both experimentally and theoretically.10,52-60 This has involved concerted motion at both, the small and the large scale.61-66 These were examined through equal-time covariances between the positional fluctuations of CR atoms67 (the covariance between two atoms i and j is calculated as 〈(ri - jri)(rj jrj)〉 where ri is the position vector of atom i, jri is the mean position of atom i and 〈 〉 denotes ensemble averages). The covariances were normalized with respect to the 300 K data as it contained the largest amplitude. The radial distribution of the covariances is shown in Figure 7a for the 6 temperatures. We note: (a) the magnitude of covariances are barely above 0.1 for 50-150 K; for temperatures

above that, the covariances increase significantly in magnitude, particularly for 250 and 300 K (also see Table 1); (b) the initial decays of covariances are near-linear/power-law for 50-100 K, increasing to exponential by 300 K. The distributions decay to zero between 2 and 4 Å for 50-150 K and ∼4.5-6 Å for 200-300 K. The correlation distance of ∼6 Å at 300 K is close to both experimental68,69 and theoretical53,67 observations. An analysis of the probability distribution of covariances (normalized for each temperature separately) as a function of distance (data not shown) shows that around 6-8 Å (which would be about the average separation between the centroids of residues that are not covalently linked), two types of behavior are seen: for temperatures below 200 K, the number of correlations stronger than an arbitrarily chosen value of 0.2 increases as the temperature is lowered and contrastingly, for temperatures g200 K, the number increases as the temperature is increased. The behavior in the 100 ps) motion is Brownian in character and is increasingly diffusive as the temperature is increased, with large variance beyond 200 K; motion at sub-100 ps time scales is increasingly non-Brownian (fluctuation is not proportional to xtime) as the temperature is increased, with significant anharmonicity increasing in amplitude and spreading across the protein. 6. The mobility of the backbone dihedrals substantially increases from 200 K upward and is localized to regions that are in close contact with the ligand. Transitions in the side chain dihedral angles of residues in contact with the sugar-phosphate moiety are seen only from 200 K onward and may be coupled to motions that bring the ligand to its reactive conformation.

10396 J. Phys. Chem. B, Vol. 104, No. 44, 2000 The number of residues exhibiting transitions in the backbone exceeds those with transitions in the side chains at temperatures higher than 200 K. This suggests that the large scale functional motions are driven by motions in secondary structural elements as a whole.49,50,76 These motions may modulate ligand binding and overall specificity while side chain plasticity is involved in the fine-tuning. 7. Correlations of atomic motion increase with temperature particularly from 200 K upward. Below 200 K, the dynamics are characterized by small amplitude harmonic motions, coupling small groups of atoms. Above 200 K, the dynamics are characterized by long range collective motions. The 300 K simulation exhibits strong correlations between regions of the protein involved in ligand interactions and display breathing modes such as those that regulate ligand binding/release and perhaps activation. Of course the functional modes involved in binding/release of ligands may not necessarily be the same as those involved in catalysis;83 however there may exist some coupling between the two. Indeed they may well be altered due to the presence of intermolecular contacts in the crystalline state.11 This issue remains unresolved as there is contrasting evidence of inhibition of activity in a crystal of ribonuclease5 on one hand, and large scale conformational changes (the T to R transition of insulin) in its crystalline84 state on the other. We note that proper sampling of the phase space at different temperatures is yet another unresolved issue,85 particularly as the starting structure used for our studies was resolved at 300 K. Studies exploring this are currently in progress. 8. Principal coordinate analyses show increased amplitudes of motions and sampling of larger amounts of phase space as the temperature is increased. The latter, particularly evident from 200 K onward, is characterized by anharmonic motion which involves liquidlike diffusion spanning several basins of attraction,76-79,86 a feature that characterizes nonlinear dynamical systems.56,77 Motions that may well be functionally important3 involving regions of the protein in direct contact with the ligand as described by PC1 and PC2 arise primarily beyond 200 K; however some motion is seen at 150 K. Such motions might be linked to residual enzyme activity at sub-transition temperatures such as that seen in the enzyme glutamate dehydrogenase (DGH) by Daniel et al.14 Interestingly, in their neutron scattering study they find no evidence for (