This Article Abstract Full Text (PDF) Supplemental Material Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Giachetti, A. Articles by Banci, L. Search for Related Content PubMed PubMed Citation Articles by Giachetti, A. Articles by Banci, L.

Modeling the Backbone Dynamics of Reduced and Oxidized Solvated Rat Microsomal Cytochrome b₅

Andrea Giachetti ^*, Giovanni La Penna ^* , Angelo Perico  and Lucia Banci ^*

^* Magnetic Resonance Center, University of Florence, Florence, Italy; and National Research Council, Institute for Macromolecular Studies, Genoa, Italy

Correspondence: Address reprint requests to Lucia Banci, Magnetic Resonance Center, University of Florence, via L. Sacconi 6, 50019 Sesto Fiorentino (Florence), Italy. E-mail: banci{at}cerm.unifi.it.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS SUPPLEMENTARY MATERIAL REFERENCES

 
In this article, a description of the statistics and dynamics of cytochrome b₅ in both reduced and oxidized forms is given. Results of molecular dynamics computer simulations in the explicit solvent have been combined with mode-coupling diffusion models including and neglecting the molecule-solvent correlations. R₁ and R₁ nuclear magnetic relaxation parameters of ¹⁵N in the protein backbone have been calculated and compared with experiments. Slight changes in charge density in the heme upon oxidation produces a cascade of changes in charge distributions from heme propionates up to charged residues 1.5 nm from Fe. These changes in charge distributions modify the molecular surface and the water shell surrounding the protein. The statistical changes upon oxidation can be included in diffusive models that physically explain the upper and lower limits of R₁ relaxation parameters at high off-resonance fields.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS SUPPLEMENTARY MATERIAL REFERENCES

 
Electron transfer (ET, hereafter) proteins are involved in a variety of important biological processes. These proteins exist in different oxidation states and, in most cases, they contain one or more metal ions which are donor or acceptor of electrons (Bertini et al., 2001, 1994; Kadish et al., 1999; Scott and Mauk, 1996). The different structural and dynamical properties of the two oxidation states have been the subject of a huge number of studies. In particular, the mobility of the redox forms of these proteins might be different and this difference can be related to protein function (Banci et al., 1998b; Brooks et al., 1988; McCammon and Harvey, 1987), molecular recognition and reorganization energy (Canters and Van de Kamp, 1992; Corin et al., 1991; Hake et al., 1992; Peerey et al., 1991; Walker and Tollin, 1991; Zhou and Kostic, 1992). Experimental, theoretical and computational tools are all required to monitor and understand the structural and mobility differences between redox forms at an atomic level.

Cytochrome b₅ (Cyt b₅, hereafter) is a largely widespread ET protein found in a variety of mammalian and avian species, and it has been extensively studied (Kadish et al., 1999; Lederer, 1994; Mathews, 1985). The protein is mainly found as a membrane-bound protein. However, a hydrophilic domain was identified and isolated, and it still retains its activity (Ito and Sato, 1968; Spatz and Strittmatter, 1971; von Bodman et al., 1986). Because of the moderate dimensions and availability, this domain represents an ideal system for addressing in vitro problems related to how the ET process occurs and how energy is transferred and stored in living organisms (Canters and Van de Kamp, 1992; Moore et al., 1986; Moore and Pettigrew, 1990; Scott and Mauk, 1996). Structures have been determined for some isoforms through x-ray only for the oxidized form, whereas solution structures of both oxidized (Arnesano et al., 1998a) and reduced (Banci et al., 1997a; Dangi et al., 1998b) forms of rat microsomal Cyt b₅ have become recently available, in addition to the crystal structure of the oxidized bovine liver protein (Durley and Mathews, 1996).

The comparison of the solution structures of Cyt b₅ in the two oxidation states revealed small differences, indicating rearrangements that could be relevant in determining the stability of the protein in the two redox states as well as in the ET process. Nevertheless, the structural differences were not significant enough to produce differences in physicochemical properties like the solvent accessibility, backbone mobility, and solvation entropy (Banci and Presenti, 2000; Dangi et al., 1998a,b).

Local mobility has often been thought to be important in controlling ET processes. Heteronuclear spin relaxation experiments, measured through nuclear magnetic resonance (NMR) spectroscopy, can provide direct information on local mobility (Palmer, 2001; Wagner et al., 1993; Wagner, 1993). In particular, rotating frame relaxation experiments have been used to reveal motions with characteristic times in the range of micro to milliseconds (Desvaux et al., 1995; Palmer et al., 2001). However, the interpretation of these NMR data is complicated by the high degree of cooperativity of molecular segments' motions in the relevant relaxation modes.

Two kinds of effects are expected with protein oxidation change. First, the change in electron structure, localized on the metal ion, affects the metal coordination geometry and the active site structure. Second, the change in the charge density in the active site modifies both the structure and mobility of the protein matrix and, eventually, of the protein-solvent interface. Models that are able to include in a unique frame local motions, as they are contained in the high frequency relaxation data and in the rotating frame data, are therefore required to connect molecular dynamics to NMR experiments. Computer simulations, and especially molecular dynamics (MD) simulations (Brooks et al., 1988; Kollman and Merz, 1990; van Gunsteren, 1993), allow calculations of NMR relaxation parameters taking into account the necessary structural atomic details of biological macromolecules (Case, 2002). Up to date, only a few MD and stochastic simulations are long enough to give a good statistical description of molecular rotational tumbling (Peter et al., 2001; Stocker and van Gunsteren, 2000). Usually, the measured NMR relaxation rates are analyzed within a model-free approach (Clore et al., 1990; Lipari and Szabo, 1982) which describes the involved time correlation functions by separating overall rotation and internal motions. The parameters involved in the fitting procedure, restricted or extended according to the sophistication of the fitting equations, are essentially order parameters and rotational and internal correlation times.

Prompers and Brüschweiler (2002) introduced a different approach to the interpretation of the NMR dynamics of proteins, named, in its more advanced form, isotropic reorientational eigenmode dynamics. This method gives first a statistical estimate of the isotropically averaged covariance matrix of the second-rank orientations of the network of relaxing vectors. Then, it derives the dynamics from an MD simulation by describing the time evolution of the projections, onto the static eigenvectors of the covariance matrix, of the second-rank components related to the relaxing vectors. The resulting approximate dynamics normally produces an overestimate of the relaxation rates that are more relevant in determining the NMR relaxation parameters. This problem is overcome by a fitting to the NMR relaxation data to decrease the lower relaxation rates.

A more basic approach couples the diffusion theory, to describe the time evolution of protein configuration, with computer simulations, used to estimate the time-independent statistics. This method has been applied both on MD (Fausti et al., 2000; La Penna et al., 2003; Perico and Pratolongo, 1997) and Monte Carlo (La Penna et al., 2003) simulations. In this article, we report a model of nuclear spin relaxation parameters of Cyt b₅ obtained by combining detailed MD computer simulations of the protein in the explicit solvent with diffusion theory. This tool allows the calculation of time correlation functions relevant to NMR relaxation without the need of the information on time evolution of configurational variables contained in the MD computer simulation. The structural and statistical properties obtained by MD are summarized and the application of diffusion theory to the calculation of molecular dynamics are presented. The theory provides a model for the upper and lower limits of the relaxation parameters, based on the statistical results of the computer simulation of the atomic model and without the need of parameter fitting or a priori assumptions on timescale separation.

The biological relevance of these findings as well as perspectives of this approach are discussed.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS SUPPLEMENTARY MATERIAL REFERENCES

 
Molecular dynamics simulation setup
The calculations to obtain the point charges and the MD statistics of Cyt b₅ in explicit water have been performed with the package NWCHEM (HPCC group, 2002). In protoporphyrin IX, or heme type b, only the Fe ion is covalently bonded to the protein. In Cyt b₅, the heme can be bonded in two different orientations related by a 180° rotation around the CHC-CHA meso direction. This possibility gives rise to two protein conformations called A and B, the ratio of which depends on the origin of Cyt b₅, being that conformation A is always the most populated.

The point charges of the heme, the Fe-bonded histidines and the iron ion (Fe(II) and Fe(III) in the reduced and oxidized states, respectively), have been computed ab initio by using the RESP procedure (Bayly et al., 2002) as it has been implemented in the NWCHEM package. For the oxidized form, the Hartree-Fock calculations have been restricted to a single unpaired electron, as it is experimentally found for the Fe⁺³ ion in this complex.

First, a 6-31G basis-set has been used to optimize the geometry of the heme only, where all methyl, vinyl, and propionate substituents in protoporphyrin IX have been replaced with hydrogen atoms. The initial D_4h symmetry of the heme group has not been imposed (Rovira et al., 2001). Second, two imidazole molecules have been added to the resulting heme configurations with the imidazole planes translated and oriented as in the Cyt b₅ crystallographic structure of recombinant trypsin-solubilized fragment of Cyt b₅ (PDB entry 1EHB) (Wu et al., 2000). In the A conformation, the imidazole planes are almost contained in the plane formed by the two segments CHB-CHD of heme and N2(ligand 1)–N2(ligand 2).

A detailed study of the electron energy obtained with Hartree-Fock calculations using the 6-31G^* basis-set has been performed. By changing the orientation of one or two imidazole planes with respect to the CHB-CHD segment, several local minima have been identified: the two imidazole planes can be parallel to each other and parallel either to the C_meso-Fe-C_meso line or the N(heme)-Fe-N(heme) line (two conformations, named I and II, respectively); the two imidazole planes can be orthogonal to each other, parallel either to the two orthogonal C_meso-Fe-C_meso lines, or parallel to the two orthogonal N-Fe-N lines (two conformations, named III and IV, respectively); finally, the imidazole can have planes oriented at 45°, one parallel to the C_meso-Fe-C_meso line and the other parallel to the N-Fe-N line (one conformation, named V). In both the oxidation states, conformation II has been found to be the most stable and the other conformations have therefore been discarded for the successive calculations.

Finally, the vinyl, methyl and propionate side chains have been added in a standard geometry to the heme conformation II to construct a model of protoporphyrin IX-(imidazole)₂ complex. The geometry of this model has been optimized with a 6-31G^* basis-set in both the oxidation states. The resulting molecular orbitals have been used for the subsequent RESP calculations.

The choice of equivalent atoms on which constraining identical point charges in the RESP calculation has been done is according to the AMBER force field (Cornell et al., 1995), as included in the NWCHEM package. All the atoms having the same point charge in the above database have been kept equivalent in the RESP calculations. The excess of charge obtained by linking the Cß of His in the -protonated form to the imidazole has been equally spread over the Cß and Hß atoms of His-39 and His-63.

By comparing the results of RESP calculations for conformation II of reduced and oxidized forms, the major delocalization of the +1 charge excess is on Fe, on the four N heme atoms, and on the H1 and H2 atoms of the imidazole ligands. Therefore, only those point charges localized on these atoms have been changed in the oxidized form with respect to the reduced form, adding +0.4 charge to Fe, +0.1 to each N heme atom, and +0.05 to each of the H1 and H2 atoms of His-39 and His-63. The point charges of the other atoms of heme and His ligands have been kept identical to the reduced form. Point charges in reduced and oxidized forms are provided as Supplementary Material.

Molecular dynamics (MD) trajectories have been simulated for the reduced and oxidized forms of Cyt b₅. Both simulations started from the minimized averaged NMR structures (Arnesano et al., 1998a; Banci et al., 1997a) of each oxidation state (1AQA and 1AW3 PDB entries, respectively) with the missing starting and final residues added according to standard geometries. The heme moiety is taken in the A conformation with respect to the protein matrix in both the oxidation states. The net charges of residues have been chosen as they are at pH = 7, i.e., Glu, Asp, and the heme propionates carry a net –1 charge (in electron units of charge) and Lys and Arg +1. These starting structures have been energy-minimized in the vacuum. The force field has been taken from the AMBER database included in the NWCHEM package, except for the point charges modified according to the procedure summarized above. The stretching and bending parameters for the heme-His linkage have been taken identical to the C-N-Fe and N-Fe-N environments in the heme (Arnesano et al., 1998a; Banci et al., 1997b), whereas no energy restraint is applied on the orientation of the His imidazole planes with respect to the heme. The latter choice has been motivated by the need for introducing as little bias as possible in defining the heme orientation with respect to the protein matrix.

A number of 10 and 9 sodium ions have been added to the starting structures of reduced and oxidized forms, respectively, to balance the net negative charges on the two oxidation states. The ions have been located close to the most negative regions for the electrostatic potential computed on the molecular surface. The resulting structures have been merged into cubic boxes of TIP3P water molecules (Jorgensen et al., 1983). The side of the cubic boxes have been set to have a distance of 1 nm between the solute atom farthest from the solute center and the closest box side. All water molecules with the oxygen atom closer than 0.25 nm from any solute atom have been removed. An amount of 4371 and 3690 water molecules have been added to the reduced and oxidized forms, respectively, and the box sides were 5.25 and 4.97 nm.

The simulation boxes with the protein, sodium ions and water molecules have been energy-minimized. In the following MD simulations a thermostat and a barostat (Berendsen et al., 1984) have been used to keep the statistics close to the NTP ensemble corresponding to the given pressure and temperature. The pressure was always set to 0.1 MPa, and the compressibility kept isotropic to keep the simulation box cubic. The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain the lengths of bonds involving hydrogen atoms to their equilibrium values, with a time step of 2 fs. A cutoff of 1 nm was used for the Lennard-Jones and direct-space part of electrostatic interactions, and the smooth particle-mesh Ewald algorithm (Essmann et al., 1995) was used for evaluating electrostatic interactions. The solvent molecules were initially simulated at the temperature of 300 K, then all systems were gradually heated from 0 to 300 K, with 5.8 ns of MD simulations at 300 K performed for both reduced and oxidized systems.

The analysis of the time evolution of the radius of gyration (data not shown) suggested that a rather long period of the trajectory is needed for the system to reach equilibration in the given force field. The first 1.3 ns of the trajectories have been excluded for the following analysis and the last 4.5 ns (45,000 configurations) have been kept. Moreover, we have computed all the following averages including 9000 configurations in the statistics, since the distribution of square radius of gyration and the first 10 first-rank relaxation modes (see below) did not change significantly reducing the number of configurations by a factor 5.

To analyze in detail the charge density around the iron ion, we make use of a suitable function emphasizing the changes of the radial part of the charge density with the oxidation state. Once the distributions P_i(r) of distances between charged atoms and Fe have been computed for each kind i of charged atom or group, each distribution can be multiplied by the charge carried by the related atom or group and the resulting charge density can be properly normalized, as

(1)

where M is the number of charge types summarizing the charge density, z_i is the point charge carried by each type of charge, and P_id(r) is the ideal distance distribution as it is obtained for a uniform density of charge carriers in the same shell at distance r from Fe.

We shall compute the above g_c function using both the full set of atomic point charges and a reduced set of point charges, assuming that the major contribution to the charge density is due to side chains carrying net charges and sodium ions. Four kinds of charges will be assumed to represent the entire charge density:

1. The negative charge of carboxylate groups, located on the two carboxylate oxygens of Glu and Asp residues () and on the heme propionate oxygen atoms (46 atoms).
   
2. The positive charge of guanidinic groups of Arg, located on the two N atoms () (6 atoms).
   
3. The positive charge of ammonium groups of Lys, located on the N atom (z₃ = 1) (10 atoms).
   
4. The charge of the sodium ions, which carry the charge z₄ = 1 (10 and 9 atoms in the reduced and oxidized forms, respectively).
   

The distance distributions calculated for distance pairs involving Fe and these groups have been averaged over all the equivalent groups in the protein and divided by the same distribution of the corresponding ions as uniformly distributed in the average volume of the simulation box.

Diffusion theory and diffusive models
Nuclear relaxation rates of Cyt b₅ in both the oxidation states were experimentally determined to understand the extent of configurational fluctuations in the two redox states (Banci et al., 1998a; Dangi et al., 1998a). Particularly relevant NMR relaxation parameters are the relaxation rates of the z-component of the magnetization due to ¹⁵N nuclei in the static magnetic field (R(N_z), or R₁, hereafter) and in the rotating frame (R(N_z'), or R₁ hereafter). The R₁ parameter is more sensitive to slow motions affecting the protein backbone: its frequency dependence can be experimentally determined, when affected by characteristic times in the range of milli- and microseconds, due to the setup of the NMR experiments.

In this work we shall ignore the effects of paramagnetic atoms on the above relaxation rates. These effects have been estimated to be at maximum 0.14 and 0.17 Hz for R(N_z) and R(N_xy) (or R₂), respectively, over experimental average values of 2 and 9 Hz for R₁ and R₁, respectively, in the oxidized form (Banci et al., 1998a; Dangi et al., 1998a). Therefore, the paramagnetic effect accounts for a maximum variation of <10% in the measured data of the oxidized form.

The relationships between the above NMR relaxation rates and spectral densities are given by the following equations (Peng and Wagner, 1992):

(2)

(3)

The meaning of the symbols contained in the above equations is

(4)

where _i is the gyromagnetic ratio of nucleus i and r is the modulus of the given N-H vector that, in our model, has been kept constant (see previous subsection),

(5)

where is the chemical shift anisotropy of each nucleus N,

(6)

where ₁ is half of the radio-frequency amplitude in the plane perpendicular to the static magnetic field and is the chemical shift in rad/s (Larmor frequency) of each N nucleus, and finally,

(7)

In the above equations, we have assumed that ¹⁵N nuclear spin relaxation is governed by the modulation of the dipolar coupling with its bonded ¹H and of the chemical shift anisotropy. The time dependence of both these interactions is determined by the same process: the spectral densities J in the equations above contain the stochastic motions of the N-H bond, the movement of this vector also being assumed to describe the motions of the principal axis of the chemical shift tensor. We have assumed = 160 ppm for all the N nuclei in the protein. The contribution of the chemical shift exchange or modulation (Desvaux et al., 1995) can be, in principle, included. It is here neglected because of the need of computationally expensive models for including the functional dependence of ¹⁵N chemical shift from the atomic configuration (Xu and Case, 2002).

The relationship between J and the time correlation function (TCF, hereafter) statistically describing the contribution of the molecular dynamics on spin transitions is given by Cavanagh et al. (1996), as

(8)

Finally, the TCF is given by

(9)

where are irreducible spherical tensors (Rose, 1957) and is the direction of the given N-H vector. P₂ is the Legendre polynomial of order 2, and ß is here the angle that the N-H vector spans in time t.

A quantity that summarizes the mobility of the N-H vector is the integral of the orientational part of the TCF,

(10)

The quantity is called correlation time of the given internuclear unit vector.

The mode-coupling diffusion (MCD) theory of the dynamics of a biological macromolecule in solution is adopted here for the computation of TCFs. The diffusion theory treats the solvent hydrodynamically and uses a detailed molecular model for the macromolecule in terms of beads (atoms or groups of atoms) connected by real or effective bonds diffusing in an interatomic potential (the same used in the MD simulation). The beads are represented as points of coordinates r_i, with i running over the N_a beads in the frictional model, and friction coefficients _i = 6 a_i, with as the solvent viscosity. The Stokes' radii a_i of the beads are calculated here by using the solvent-accessible surface area (SASA) method (Pastor and Karplus, 1988) with a zero probe radius, by summing the surfaces of each constituent group (La Penna et al., 2000b).

In this work, we shall use two kinds of bead representations for the macromolecule. The macromolecular model with beads located on solute atoms (the "bare" model, hereafter) uses one bead for each amino acid (positioned on the C atom) and one bead for the heme (positioned on Fe; La Penna et al., 2000a). In the second model, the beads are positioned on the water molecules surrounding the macromolecule (the "surface" model, hereafter) and all the beads have the same Stokes' radius of 0.14 nm (i.e., the usual Stokes' radius of the water molecule).

These bead frictional models, both with and without the inclusion of water layers, have been widely used for rigid macromolecules (de la Torre et al., 2000; Fernandes et al., 2002). The MCD approach in the form summarized below constitutes the generalization of diffusion theory to flexible macromolecules that gives the exact rotational diffusion in the limit of rigid macromolecules (La Penna et al., 1999).

The MCD approach (La Penna et al., 1999; Perico and Pratolongo, 1997), can be briefly summarized as follows. The N_a beads of the given frictional model can be connected by N_b bonds (l_i, i = 1, ..., N_b). The variables l_i can be organized in the 3 x N_b dimensional array l containing all the bond vectors l_i. This array entirely describes the model configuration, and its dynamics is regulated by the operator L, adjoint to the diffusion Smoluchowski operator D,

(11)

where U is the potential energy of the beads as a function of the bead coordinates, k_B is the Boltzmann constant, and T is the absolute temperature.

The diffusion tensor D is given by

(12)

(13)

with H and T representing the hydrodynamic interaction matrix and the Rotne-Prager (Rotne and Prager, 1969) tensor, respectively, and D_i = k_BT/_i is the diffusion coefficient of each bead.

By expanding the conditional probability (solution to the Smoluchowski equation) in a complete set of eigenfunctions of L, the TCF of any coordinate-dependent dynamic variable with zero average f(t) may be expressed in the standard form

(14)

where –_i and _i are, respectively, the eigenvalues and the normalized eigenfunctions of the operator L,

(15)

This eigenvalue equation is solved by expansion of the eigenfunctions in a suitable basis-set. Once a basis-set has been chosen, the eigenvalue equation for the diffusion operator becomes a generalized eigenvalue equation in matrix form that can be solved with standard methods. By using the RM2-II basis-set (i.e., a second-order approximation) in MCD theory (La Penna et al., 1999), the second-rank TCFs are sum of exponential functions, with the constant _i the rates of the relaxation modes. Note that the isotropic reorientational eigenmode dynamics method, mentioned in the previous section, gives an alternative derivation of these relaxation rates: the time correlation functions for the projection on the static eigenvectors of the covariance matrix of the second-rank tensor components related to relaxing vectors (f in Eq. 14) are derived directly from the MD simulation.

Once the second-rank TCFs have been computed with the RM2-II basis-set, they can be analytically Fourier-transformed to the spectral densities (Eq. 8) and analytically integrated to give correlation times (Eq. 10). In the MCD approach, the incomplete conformational sampling of MD affects the rates and the coefficients in the multiexponential expansion of the TCFs of interest (Eq. 14), introducing errors that increase with the rate.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS SUPPLEMENTARY MATERIAL REFERENCES

 
Configurational statistics
The Fe atom in the heme group has been covalently bonded to the N2 atoms of His-39 and His-63 to mimic the Fe axial coordination in Cyt b₅, but no torsional potential contributions related to all dihedral angles centered on the Fe-N2 bonds have been included. Therefore, it is meaningful to analyze the deviations of these dihedral angles from the NMR structures both in the reduced and oxidized states. In the reduced form, the dihedral angles NC-Fe-N2-C1 in His-39 and His-63 are –43 ± 9 and 32 ± 13°, respectively, averaged over the MD trajectory, which compare with –45 and 25° in the most refined available solution structure (1B5A PDB entry) (Dangi et al., 1998a).

For the oxidized form, the averages over the MD trajectory of the same dihedral angles are –33 ± 10 and 33 ± 11°, respectively, compared to –39 and 21° of the available crystallographic structure (1CYO PDB entry) (Durley and Mathews, 1996) and –45 and 24° of the NMR-minimized average structure (1AW3 PDB entry) (Arnesano et al., 1998a). For both the redox forms, only moderately larger average dihedral angles are observed compared to NMR data. On average, the MD trajectories of both the redox forms display the two imidazole planes parallel to each other and only slightly rotated from the CB-Fe-CD meso direction (conformation II in Methods, corresponding to –45 and +45°) toward the NC-Fe-NA direction. Conformation II has been found the most stable by ab initio calculations (see Methods) and both MD-averaged dihedral angles are consistent with this conformation within statistical errors. Moreover, in the oxidized form the MD-averaged dihedral angles are consistent with the values derived by paramagnetic effects on NMR parameters (see Fig. 5 in Banci et al., 2002).

View larger version (27K): [in this window] [in a new window]   FIGURE 5  Configurations of reduced (left) and oxidized (right) forms at 3 ns, corresponding to about the maximum increase of molecular SASA in the oxidized form. The heme is displayed with all the bonds connecting heavy atoms and Fe is displayed as a sphere. Arrows point to the C-terminus. Figure prepared with the MolMol program (Koradi et al., 1996).

The effect of the different oxidation state on the secondary structure can be monitored through the time evolution of the state of each residue along the MD trajectories (Kabsch and Sander, 1983). In Fig. 1, every 45 ps for each residue the helical state is identified by a solid bar, with the ß-sheet state by a shaded bar, and compared to -helices in the NMR structures, shown as vertical bars on the left. It can be observed that h1, h3, h4, and h5 are all maintained during the simulation, showing a larger stability of h3 and h4 in the oxidized state with respect to the reduced state. The region corresponding to h2, a short helical portion in both NMR structures, is disordered after MD equilibration, whereas h6 is partially lost in the oxidized state only. In both redox forms, the ß-sheet structure is well-maintained during the whole MD simulation. Except for h2 and h6, the -helical and ß-sheet populations obtained by the MD trajectories represent the wobbling of the secondary structure elements in both reduced and oxidized states solution structures.

View larger version (78K): [in this window] [in a new window]   FIGURE 1  Time evolution of secondary structure of reduced (A) and oxidized (B) forms along the MD simulations: helices are identified in solid bar, ß-sheets in shaded bar; all the other motives are not displayed. Time resolution is 45 ps. Residues in -helix in the NMR structures are displayed as vertical bars on the left.

Because of the slow decay through space of electrostatic interactions, the modification of charge density on the heme group upon oxidation is expected to produce long-range effects both on the protein matrix and on the solvent and ions surrounding the protein. This may be particularly important in Cyt b₅ because the heme is close to flexible portions of the protein and to the bulk polar environment (water and counterions), where charged atoms are allowed to move upon change in the electric field. The residues more sensitive to a change of charge density in the heme are those containing a net charge in the side chains, i.e., Glu–, Asp–, Lys+, and Arg+. Any structural change affecting these charged residues is also expected to affect, more or less directly, the solvent and ionic protein environment. This long-range effect can significantly contribute to the variation of entropy upon heme oxidation.

In Cyt b₅ at pH 7 there are 21 negatively charged (9 Asp and 12 Glu) and 13 positively charged (10 Lys and 3 Arg) residues. The further addition of the –2 and –1 charges of the heme in the reduced and oxidized forms, respectively, produces a net negative charge of –10 and –9 in the two oxidation states, respectively. In our model the negative charge in the molecule has been balanced by 10 and 9 sodium ions in the reduced and oxidized forms, respectively. To summarize the different behavior of the negative and positive charge locations in the two redox forms of the protein, we have computed a set of charge density distribution functions for the full atomic point charges used in the MD simulation and for simplified representations of net charges.

Assuming that the heme group did not significantly change its position within its pocket, as it is revealed by the very minor changes in the geometrical parameters with respect to the His ligands, we analyzed the radial part of the charge density distribution (i.e., g_c(r) of Eq. 1) in terms of distances between Fe and any other charged point in the protein and sodium counterions.

The g_c(r) function is plotted in Fig. 2 for the two redox forms of the protein and using the atomic point charges. It can be noticed that the function displays a rugged behavior in the range 0–1 nm, which is a direct effect of the charge density modification within the heme pocket. The peaks within 0.5 and 1 nm contain the contributions of the heme propionates (see below) and the smoother behavior of the charge density beyond 1 nm is due to residues in the protein matrix not directly bonded to the heme. Therefore, in this latter region the changes in charge density with oxidation is due to the structural change of the protein matrix and of the counterion distribution around the heme pocket. It can be observed that a significant change in charge density occurs at 1.5 nm from Fe, where a negative contribution in the reduced form is compensated by some positive contribution in the oxidized form.

View larger version (13K): [in this window] [in a new window]   FIGURE 2  Charge radial distribution function of Eq. 1 in the text, calculated with the atomic point charges. Reduced form (solid line) and oxidized form (dotted line).

To separate the various contributions to charge density, we have calculated the charge density contributions of those classes of groups expected to be more relevant for understanding the charge density modifications, i.e., the groups of residues carrying a net negative or positive charge and the sodium counterions (see Methods). In Fig. 3, the g_c(r) function for this reduced charge-set is plotted for the two redox forms of Cyt b₅. Different effects of the redox state of Fe can be observed in the three r ranges 0–1 nm, 1–2 nm, and beyond 2 nm. In the first range the closest negative density in the oxidized form is due to the bending of one of the two heme propionates. In the reduced form one of the 10 sodium ions makes a strong salt bridge with the A propionate (data not shown), but no sodium ion can come too close to Fe in the oxidized form because of the strong hydrogen bonds between the propionate and surrounding atoms, i.e., the backbone amidic H atoms of Ser-64 and His-63. In the second range, 1 < r < 2 nm, the positive charge density significantly changes upon heme oxidation beyond 1.2 nm. In this region, a decomposition of g_c in terms of positive and negative contributions shows that negative charges move farther from Fe and positive charges become closer to Fe in the oxidized state, i.e., the opposite of what is expected. The charge movement inverts the density oscillation in the 1–2 nm region and this behavior accounts for the change in charge density in the model observed at distances 1.5 nm (Fig. 2). Beyond 2 nm, the effect of the change in oxidation state becomes almost negligible.

View larger version (14K): [in this window] [in a new window]   FIGURE 3  Total charge radial distribution functions calculated with simplified point charges (see text). Reduced form (solid line) and oxidized form (dotted line).

The molecular surface can be indirectly changed by the electrostatic interactions because of possible change in the distribution of hydrophobic and hydrophilic groups on the surface accompanying the charge density modification. This effect can be observed by looking at the time evolution of the SASA computed for the whole molecule in the two different oxidation states (Fig. 4) (Eisenhaber et al., 1995). It can be observed that the molecular surface in the oxidized state tends to increase during the MD trajectory. Starting from a value very similar to that of the reduced form, it becomes 10% larger after 1 ns and back to the reduced form values at the end of the trajectory. This behavior shows that the molecule in the oxidized form is changing shape during the trajectory whereas in the reduced state the starting shape is more stable. A comparison with the time evolution of the secondary structure (see Fig. 1) shows that the molecular expansion in the oxidized form occurs when h2 is not populated.

View larger version (42K): [in this window] [in a new window]   FIGURE 4  Time evolution of molecular SASA for the reduced (solid line) and oxidized (dotted line) forms.

All the structural modifications described above for the oxidized form are within 0.3 nm (estimated over C of residues 1–88) from the x-ray bovine structure (Durley and Mathews, 1996) and from rat microsomal NMR structure (Arnesano et al., 1998a). This range of variability is commonly observed as a consequence of the actual solution environment at room conditions.

To conclude the analysis of conformational changes occurring in Cyt b₅ upon oxidation, we computed the solute-solvent radial distribution function (RDF, hereafter), g(r), following a reported procedure (La Penna et al., 2003). This function contains relevant information on the structure of the water shells surrounding the protein. In Fig. 6, the RDF function of distance pairs involving any of the protein atoms and any of the water oxygen atoms have been plotted for the reduced and oxidized forms. As already observed (La Penna et al., 2003), the first peak of the function contains pairs involving heavy protein atoms and the water oxygen in the first solvent shell. In both cases, for large r-values the RDF approaches values >1 because of the higher density of water as obtained in the simulation box with respect to bulk water at the same temperature and pressure conditions. The most significant difference in the RDF between the two redox forms is, however, concentrated in the first shell: in the reduced form, a well-structured water shell can be observed with water molecules ordered within 0.3 nm; in the oxidized form, this water shell is looser, spreading up to 0.35 nm far from the solute, and the distinction is no longer possible between the first and second solvation shells. The less-defined first water shell in the oxidized form is related to the more fluctuating molecular surface. As observed for other proteins and nucleic acids, the radial organization of the first water shell is mainly due to hydrophilic and charged groups exposed to the solvent by the molecular surface. In the oxidized form of Cyt b₅, these groups are more mobile, as is shown by the wider distribution of both the Fe-negative and Fe-positive radial charge distributions displayed in Fig. 3. Therefore, a water shell has fewer chances to be created in the oxidized form. All these observations suggest that the decrease in water organization, as it occurs in the MD simulation in explicit solvent, can be relevant for explaining the positive entropy change upon oxidation of Cyt b₅ (Dangi et al., 1998a). The implications of these observations can be consistently included in the following diffusive model.

View larger version (18K): [in this window] [in a new window]   FIGURE 6  Solute-solvent radial distribution function (RDF) for the reduced (solid line) and oxidized (dotted line) forms.

In the reduced form, the first peak in RDF is within 0.30 nm, and 360 ± 13 water molecules have been counted on average in this range. In the oxidized form, the first peak is approximately in the range 0–0.35 nm, where 530 ± 14 molecules are located. Therefore, the surface model for the reduced form is made of 360 friction points located on the surface, whereas in the oxidized form it is made of 530 friction points. Once the water density has been explicitly calculated from the MD trajectory, friction locations are chosen randomly on the molecular SASA, as already described in La Penna et al. (2003).

By analyzing the first-rank rates obtained by the diffusive bare and surface models, it is possible to derive estimates of the diffusion tensor principal components of a molecule with an average shape derived by the MD trajectory. The protein in both redox states has ₁ < _{2 3}, thus representing a diffusion tensor with eigenvalues in the order D_zz > D_xx D_yy and a rod-like diffusion tensor. By using the relationships in Table 2 of La Penna et al. (1999), we obtain the results in Table 1 (this article). The ratio D_||/D for both models is comparable with the estimates in the literature (Dangi et al., 1998a), even if a slightly less anisotropic diffusion tensor is found here for the oxidized form (1.30 compared to 1.35). As expected, the diffusion tensor eigenvalues decrease by using the solvated diffusive models compared to the bare models, because of the larger size of the surface bead aggregates. However, the ratio does not change with the model and the rotational diffusion anisotropies of the two redox forms are not significantly different.

View larger version (23K): [in this window] [in a new window]   FIGURE 7  Correlation times (Eq. 10 in text) of H-N bonds in the reduced (solid line) and oxidized (dotted line) forms calculated with the bare model. Residues in -helix in the NMR structures are displayed as horizontal bars for reduced (top) and oxidized (bottom) forms.

A direct comparison of R₁ ¹⁵N relaxivities can be made between the values calculated through the two diffusive models and experiments. The R₁ values for reduced and oxidized forms are shown in Fig. 8, A and B, respectively. It can be observed that the choice of the diffusive model produces mainly a shift of the data because of the global change of relaxation rates. In both redox states the experimental data show a smaller variation along the protein chain compared to the calculated values, particularly in the regions that have been found less ordered in the simulations. The diffusive method allows the calculation of the P₂(t) TCF as the sum of many exponential functions, each decaying with one of the eigenvalues (rates) of Langevin equation in its matrix representation. The simplification of such TCF in terms of two exponential functions is possible and the amplitude of the function with the lowest rate is the so-called Lipari-Szabo order parameter (Fausti et al., 1999). For this protein, this simplification gives almost the same NMR relaxivities of using the entire set of rates (data not shown), thus showing that for these molecular statistics the information on the mobility pattern is all contained in the order parameters and in the five lowest second-rank relaxation modes. For both the models, the differences between experiments and calculations arise from the low convergence of the MD trajectories.

View larger version (29K): [in this window] [in a new window]   FIGURE 8  R1 NMR relaxation rates of 15N at (1H) = 600 MHz for the reduced (A) and oxidized (B) forms: experiments (squares with error bars), results of the bare model (solid line), and results of the surface model (dotted line).

In Fig. 9, the order parameters obtained with and without block-averaging have been plotted for the two redox forms and compared with the order parameters obtained by fitting R₁ relaxation rates (Dangi et al., 1998a). It can be observed that rarely sampled conformational transitions in the less ordered regions dramatically affect the order parameters and, therefore, the mobility and the relaxation rates, but the systematic removal of these statistical errors can partially recover the experimental information and flatten the R₁ variation within the molecular chain.

View larger version (23K): [in this window] [in a new window]   FIGURE 9  Lipari-Szabo order parameters for the reduced (A) and oxidized (B) forms: results obtained by fitting R1 NMR experiments (squares), results from averaging the MD trajectory (solid line), and from averages with blocks of 300 ps each (dotted line).

1

1

1

View larger version (29K): [in this window] [in a new window]   FIGURE 10  Same as Fig. 8 for the R1 NMR relaxation rates of 15N at (1H) = 600 MHz and e = 1870 Hz for the reduced (A) and oxidized (B) forms.

1

1

1

The extent of this effect can be physically related to the molecular solvation and surface. We have previously observed that these two latter parameters display by far the most evident differences between the statistics of the reduced and oxidized forms of Cyt b₅: the different charge density on the heme produces structural changes that modify the molecular surface: this produces a molecular expansion (Fig. 4) and a looser solvation shell in the oxidized form (Fig. 6).

The bare and surface diffusive models can be considered as lower and upper limits, respectively, for the R₁ relaxation rates or, in general, for the low frequency contributions to spectral densities when the effects of the chemical shift modulation are ignored. Any frequency dependence of the NMR relaxation rates can then be related to a slow exchange between the leading terms in the different diffusion models. This kind of behavior is expected, for instance, when a molecular statistics allows the equilibrium between a basin of configurations A, with a first organized water shell, and another basin of configurations B, in which the first shell of water molecules is less defined. This is exactly the behavior suggested by Fig. 4, where, in the oxidized form, the molecule seems to return in a more compact form after the expansion.

Typical rotating frame experiments are performed in the fast exchange limit. Assuming k_a as the unique conversion rate from the basin of configurations A to the basin of configurations B, and k_b as the rate of the inverse process, the NMR spectrum displays a single peak for each ¹⁵N nucleus when the amplitude of chemical shift modulation within the different conformations is smaller than k = k_a + k_b. Within this simplification, when ₁ is larger than k and, therefore, the contribution from these conformational transitions to the R₁ parameters of each ¹⁵N nucleus is zero, the relaxation parameters tend to the population average for each nucleus. The MD simulations for the two forms of Cyt b₅ suggest that when the population of basin A is 100%, the asymptotic behavior of R₁ would be the result of the bare model, whereas in the opposite case, when the population of basin A is zero, R₁ would be the result of the surface model. These different ₁ asymptotic behaviors of R₁ can then be included in more complete models of R₁ when chemical exchange events are significant, both in the fast and slow exchange limits (Trott and Palmer, 2002).

The MCD theory, in the present form, is not yet able to derive the protein dynamics related to infrequent events coupled with the molecular wobbling. Therefore, the RM2-II MCD calculations cannot obtain the timescale 1/k for the exchange between conformations revealed by the frequency dependence of R₁ in the experiments. Nevertheless, the appearance of collective conformational exchange events in the oxidized form accounts for the larger number of ¹⁵N nuclei whose R₁ relaxation rates are frequency-dependent in the kHz region in the oxidized form with respect to the reduced form. Within this frame, the populations and the rate constant k within the two basins are the same for all the nuclei. Therefore, the contribution to R₁ parameters due to chemical shift modulation depends only from the amplitude of this modulation, i.e., the difference in chemical shift of the nucleus in the two basins of conformations. This difference strongly depends on the average position and hydration of the nuclei within each basin.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS SUPPLEMENTARY MATERIAL REFERENCES

 
Molecular dynamics simulations of 4.5 ns has been performed for the Fe(II) and Fe(III) redox forms of Cyt b₅, using as starting configurations the structures independently determined for each oxidation state. The two systems also differ for the point charges of iron and a few atoms of its ligands.

The configurational statistics obtained by MD simulations have been analyzed in detail. The most populated conformation of the histidine ligands has been found with the two imidazole planes parallel to each other and almost parallel to the B-D meso direction of the heme for both the oxidation states. The average dihedral angles describing this conformation are consistent with the solution NMR structures and with the pattern of chemical shifts due to the unpaired electron in the oxidized form.

One of the heme propionates is bent toward the iron ion in the oxidized form, whereas in the reduced form it is screened by a sodium ion. The bent conformation of the propionate in the oxidized form is strongly stabilized by hydrogen bond interactions with the backbone amidic H of Ser-64 and His-63. These interactions prevent the approach of sodium ions to the heme site and the screening of the heme propionate negative charge.

The distribution of secondary structure motives reveals that in both oxidation states the ß-sheet structure is maintained, as well as the core of -helices h3–h5. Two NMR-derived -helices are found unstable. Helix 2 in both redox forms is rarely populated, thus showing its large propensity to unfolding as already suggested by NMR data in GdmCl 2 M solution. The mobility of this region is important because it contains His-39 bonded to Fe. Helix 6 is less stable in the oxidized form because of the breaking of electrostatic interactions between the positively charged groups of the C-terminal and h6 region (residues 83–90) and negatively charged residues in h4 region (residues 53–60). In the oxidized form, the breaking of these interactions increases the stability of helix 4 and the solvent exposure of the region 53–60. This region has been identified as one of the negative patches in the recognition of Cyt c and modulation upon oxidation of the mobility and solvent accessibility has been already suggested by analyzing MD simulations of the bovine liver oxidized form.

The breaking of these electrostatic interactions and the related increase in molecular surface are indirect consequences of the larger negative charge close to Fe in the oxidized state: the larger positive charge in the oxidized state of iron is efficiently delocalized within the coordination site, whereas the negative charge of the bent propionate and the extrusion of a positive sodium ion from the active site has a stronger effect on the charge distribution in the protein. Negatively charged side chains at 1.5 nm from Fe change conformation and positively charged side chains, which are in the reduced form involved in electrostatic interactions, are therefore affected.

The change in size and nature of the molecular surface occurring in the oxidized state has significant effects on the structure of the first solvation shell: in the reduced form the first water shell is relatively tightly bound to the protein, whereas in the oxidized state it moves with the molecular surface and it is, therefore, more mobile. The consequence is that, on average, there are 50% more water molecules in the oxidized form than in the reduced one in the first solvation shell, but most of these molecules are less correlated to the solute.

The different structure of the solvation shell strongly affects the hydrodynamic properties of the molecule as suggested by reduction potential and NMR experiments. Mode-coupling diffusion theory has then been applied to calculate the molecular dynamics of H-N bonds in the protein backbone. The Smoluchowski diffusion equation, represented as the eigenvalue equation for the adjoint of the diffusion operator, has been solved by computing averages obtained along with the MD statistics. Two different diffusive models of the solute molecule have been used: a first bare model ignores the contribution of the solvation shell to the molecular friction, whereas a second surface model considers only the contribution of the first shell of water molecules correlated to the solute. The number of water molecules in this shell has been computed through the analysis of the solute-solvent radial distribution function and is not an adjustable parameter.

The results of the bare and surface diffusive models nicely match the lower and upper limits, respectively, of ¹⁵N NMR R₁ relaxation measurements on both the redox forms of the protein when chemical shift modulation contribution to NMR relaxation is not effective. This result suggests that a complete description of these relaxation parameters can be given in terms of an exchange process between conformational basins related to different molecular shapes. The MD statistics shows that the basin containing expanded molecules is more populated in the oxidized state because of the sampling of conformations where arrays of electrostatic interactions are broken. Solvation analysis shows that the sampling of these conformations produces a larger and looser first solvation shell in th