INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Electron transfer from cytochrome f to plastocyanin is a central reaction in the photosynthesis of all plants and bacteria (Bendall et al., 1999). The blue copper protein plastocyanin (pc) (Redinbo et al., 1994) is a mobile electron carrier protein with an eight-stranded flattened Greek-key -barrel fold. It contains a type-I copper atom coordinated by two histidines, a cysteine, and a methionine (Adman, 1992; Baker, 1994; Sykes, 1994). Its electron-donor, cytochrome f (cytf) (Carrell et al., 1999; Martinez et al., 1994), is a c-type heme protein constituent of the cytochrome b6/f complex that is anchored to the membrane by a single C-terminal transmembrane helix. The N-terminal part of cytf has two -sheet domains arranged to form an elongated structure; the heme is bound to the larger domain with the N-terminal -amino group as a ligand to the heme iron.

Although a significant number of crystal and solution structures of pc, along with a few cytf structures and an nuclear magnetic resonance (NMR)-based structure of the pc/cytf complex, have been determined, the details of the interaction mechanism of the pc/cytf redox system are not known. Long range molecular recognition processes and electrostatic effects are, however, thought to play important roles in the kinetics of the reaction, at least in vitro (Hope, 2000). The increasing interest in unraveling the details of the interaction between pc and its electron donor, cytf, is reflected in the scientific papers and reviews published recently on this topic (Hope, 2000; Gong et al. 2000a,b; Hirota et al., 2000; Illerhaus et al., 2000). The generally accepted hypothesis of electrostatic steering between pc and cytf is supported by significant local electrostatic complementarity of the molecular surfaces of the physiological partners (Bendall et al., 1999; Hope, 2000; Illerhaus et al., 2000).

In this respect, particularly worthy of note is the mutual binding complementarity between pcs and their electron donors in eukaryotic systems (Bendall et al., 1999). Higher-plant pcs are strongly acidic proteins with negatively charged residues concentrated on the so-called "eastern site" of the proteins in which Tyr-83 (spinach pc numbering), one of the putative entry/exit points for electrons, is situated. This acidic site consists of two separated clusters: the "small acidic patch" of residues 59-61 (spinach pc numbering), located near the Cu-site, and the "large acidic patch" of residues 42-45 (spinach pc numbering), located in the 5-6 loop region. Eukaryotic cytochromes f, although overall weakly negatively charged, have a ridge of positively charged residues located in the region where the heme is bound. The acidic area in eukaryotic pcs is likely to be necessary for recognition of the basic patch in cytf. The NMR-based structures of the spinach pc/turnip cytf complex obtained by Ubbink et al. (1998; Ejdeback et al., 2000) in fact show that the pc eastern site acidic residues interact with arginine and lysine residues of the electron donor partner.

It is possible that pc/cytf binding is mainly due to nondirectional Coulombic forces that provide slippery surfaces and allow a dynamic interaction between the partners with a large number of complex configurations in fast exchange (Bendall et al., 1999). Experimental studies, such as kinetic measurements on cross-linked and chemically modified systems (Gross et al., 1990; Qin and Kostic, 1993; Lee et al., 1995; Hibino et al., 1996; Hope, 2000), indicate that the initial electrostatic complex is not optimal for electron transfer to take place and that rearrangement of the complex is required. Thus, it appears that the key steps in the redox reaction between pc and cytf are: 1) the formation of an initial complex of a purely electrostatic nature; 2) the rearrangement necessary to attain a configuration of the complex optimal for electron transfer; and 3) the electron transfer itself.

The overall reaction

(1)

can be explained using the following kinetic model (Hope, 2000):

(2)

in which k_on and k_off are, respectively, the on and off rate constants for association of pc_ox and cytf_red, which might be different from those (k'_on and k'_off) for dissociation of the pc_red/cytf_ox complex; k_rearr and k_der are the forward and reverse rate constants for the rearrangement of the complex; and k_et and k_ret are the forward and reverse electron transfer rate constants, respectively.

Neglecting the rearrangement step, the overall kinetic constant k₂ can be written as (Kannt et al., 1996):

(3)

in which K_A = k_on/k_off is the equilibrium binding constant and K_F = k_et/k_ret is the equilibrium constant for the electron transfer step.

In the case of an association-controlled reaction (k_et k_off), k₂ = k_on and k₂ can be treated as a lower limit for k_on; for an activation-controlled reaction (k_et k_off), k₂ = K_Ak_et and k₂ can be considered to set a lower limit for the electron transfer rate k_et.

If the mutual rearrangement of the two proteins in the complex has to be considered and is the limiting step, interpretation of k₂ is not straightforward (Kannt et al., 1996).

Whereas the second order rate constant k₂ in Eq. 1 does not provide direct information about mechanisms, it is often used in comparative studies for estimating the effects of mutations in pc and/or cytf on the overall electron transfer rate (Hope, 2000). Experimental values of the overall kinetic constant k₂ have been determined by Kannt et al. (1996) for the in vitro electron transfer reaction between wild-type (wt) and mutant pcs from spinach and the soluble part of cytf from turnip. Although these proteins are not physiological redox partners, the reaction is still effective and can be chosen as a model system. Seven single and multiple point mutations (D42N, E43N, E43K, E43Q/D44N, E59K/E60Q, E59K/E60Q/E43N, and Q88E) were introduced into the acidic eastern site of spinach pc, and the effects on overall electron transfer from turnip cytf were analyzed. The results suggested that, in this set of mutants, the overall electron transfer rate constant is modulated by the association binding constant, whereas the electron transfer step itself appears to be almost unaffected by the mutations introduced (Kannt et al., 1996). Further support for this hypothesis is found in the recent paper by Illerhaus et al. (2000), who pointed out the role of the negatively charged residues at the eastern site of spinach pc in molecular recognition of the intact physiological partner, the spinach cytochrome bf complex.

The aim of the present work is to obtain insights into the role of the negative patch on the eukaryotic pc eastern site in the association and, as a consequence, in the overall electron transfer process with cytf, by means of computational simulation and comparison methods. The availability of the solution structure of the spinach pc/turnip cytf complex (Ubbink et al., 1998; Ejdeback et al., 2000) renders these computations possible.

The problem of correctly reproducing the association of pc with cytf or cytochrome c has been addressed previously by various computational techniques. Roberts et al. (1991) used a computer-graphics alignment of protein electrostatic potentials followed by a systematic orientational search of intermolecular electrostatic energies at different protein-protein separation distances to build a complex of pc with cytochrome c. Pearson and colleagues used both a manual docking procedure (Pearson et al., 1996) and a Brownian dynamics (BD) simulation (Pearson and Gross, 1998) to dock poplar pc and turnip cytf. Soriano et al. (1997) also used a manual docking procedure to generate three docked structures of the complex of pc from poplar and cytf from turnip, which proved to be consistent with the structures obtained previously by Pearson et al. (1996). Ullmann et al. (1997) used a combination of Monte Carlo and molecular dynamics (MD) simulations to generate putative docked complexes of French bean pc and turnip cytf for each of which they analyzed electronic coupling between the copper and the heme sites.

Two different approaches are used in our study. First, a comparative analysis of the computed molecular electrostatic potentials (MEPs) of the wt pc and the seven mutants from Kannt et al. (1996) was carried out using our PIPSA (protein interaction property similarity analysis) procedure (Blomberg et al., 1999). The Hodgkin similarity index (SI) (Hodgkin and Richards, 1987) is used as a quantitative descriptor of the similarities and differences between the electrostatic potentials of the proteins. Blomberg at al. (1999) used Hodgkin SIs for a principle component analysis of the electrostatic potentials of modeled PH domains to classify them according to their binding properties. We recently used this approach to compare the molecular interaction properties of different blue copper protein subfamilies in a semiquantitative way (De Rienzo et al., 2000). In the present work, SIs are used as quantitative descriptors to estimate relative electrostatic interactions between pc mutants and cytf relevant for their association. As far as we are aware, this is the first time that SIs calculated for a macromolecular system have been used as quantitative descriptors of the kinetic properties of the system itself in a quantitative structure-activity relationship analysis.

Second, the effects on the pc/cytf association rate of mutations of residues in the pc eastern site were studied by Brownian dynamics simulation of the association of the wild-type and mutant pcs with cytf. The "Simulation of Diffusional Association of proteins" (SDA) program (Gabdoulline and Wade, 1998) used for this was previously applied to study the association of several protein-protein systems (Gabdoulline and Wade, 1997, 2001; Elcock et al., 1999).

These two different approaches to address the problem of long range electrostatic interactions in protein-protein association were applied in a concerted manner to provide complementary insights into the association processes of redox systems.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Construction of structures of mutant pcs

The crystal structure of spinach pc (Xue et al., 1998), used for MEP SI calculations, and the NMR-based structure of the spinach pc/turnip cytf complex (Ubbink et al., 1998; Ejdeback et al., 2000), used in the Brownian dynamics simulations, were obtained from the protein databank (pdb codes 1ag6 and 2pcf, respectively).

SI calculation

The crystal structure 1ag6 was mutated (D8G) to the spinach pc wt sequence (as given in the SWISSPROT protein sequence database). Mutants were then obtained by substituting the appropriate residues (D42N, E43K, E43N, Q88E, E43Q/D44N, E59K/E60Q, E59K/E60Q/E43N). All mutations were performed with the biopolymer module of the InsightII software package (MSI, San Diego, CA, 1995). The WHATIF program (v. 19990307-1124) (Vriend, 1990; Hooft et al., 1996) was used to add polar hydrogen atoms to the wt and mutant structures.

BD simulation

The 10 structures of the pc/cytf complex deposited in the pdb file 2pcf were obtained by using constraints derived from NMR experiments with labeled pc to perform classical MD simulations (Ubbink et al., 1998; Ejdeback et al., 2000). As a consequence of the experimental technique used, the pc/cytf structures are best defined in the region close to the heme site of cytf, whereas the interactions in the other regions at the complex interface are less accurately defined.

After a detailed analysis of these 10 structures, structure number 9 was selected and used for the Brownian dynamics simulations. In this structure, pc shows the greatest conformational deviation from the ensemble mean (root mean square deviation (rms) = 2.94 Å as reported in the 2pcf file). Two different sets, A and B, of wt and mutant pc/cytf complexes were built. Set A consists of the wt and seven pc mutants (D42N, E43K, E43N, Q88E, E43Q/D44N, E59K/E60Q, E59K/E60Q/E43N) obtained by mutating the selected structure with the biopolymer module of InsightII (MSI, 1995). The WHATIF program (v. 19990307-1124) (Vriend, 1990; Hooft et al., 1996) was then used to add polar hydrogen atoms to the wt and mutant structures. Set B was built by mutating the appropriate residues in the wt complex after having performed energy minimization and MD on it, using the program CHARMM22 (Brooks et. al., 1983). This was done to optimize the interactions at the pc/cytf interface complex. The minimization procedure consisted of 100 steps of steepest descent, followed by a conjugate gradient minimization until the rms gradient of the potential energy was less than 0.0001 kcal/mol·Å. The united-atom CHARMM22 force-field parameters, a 12-Å nonbonded cutoff, and a dielectric constant of 80 were used. The energy-minimized coordinates were used as a starting point for MD simulation. During the MD simulation, the lengths of the bonds involving hydrogen atoms were constrained according to the SHAKE algorithm (van Gunsteren and Berendsen, 1977), allowing an integration time step of 1 fs. The secondary structure was maintained by applying harmonic distance constraints with a force constant of 1.0 kcal/mol·Å² and a scaling factor of 10, between the backbone oxygen and nitrogen atoms of residues in opposite strands of -sheets and between the hydrogen-bonding backbone oxygen and nitrogen atoms in -helices (between residues i and i + 4). The coordination geometry of the Cu-site in pc was maintained by applying weak harmonic distance constraints between the Cu atom and its ligand atoms His-37-ND1, His-87-ND1, Cys-84-SG, and Met-92-SD. The coordination geometry of the heme-site in cytf was maintained by applying weak harmonic distance constraints between the Fe atom and its axial atom ligands His-25-NE2 and Tyr-1-N, and between the Fe atom and the pyrrole ring nitrogen atoms, which are the equatorial ligands. The formal charges of +2e on the Cu and the Fe ions were redistributed between the metal ions and their ligands as described in the following section (MEP calculations). The structure of the complex was thermalized to 300 K with a 5-K increment per 6000 steps implemented by randomly assigning individual atomic velocities from a Gaussian (Brooks et. al., 1983) distribution. After heating, the system was allowed to equilibrate until the potential energy was approximately stable (20 ps). Velocities were scaled by a single factor every 0.050 ps to maintain constant temperature. An additional 10 ps of equilibration were run without velocity scaling. Data were collected every 0.5 ps over the last 40 ps of each simulation. The time-averaged structure was subjected to minimization by the same procedure as used before the MD simulation. The mutant series (set B) was produced by making mutations in the final minimized structure with the biopolymer module of the InsightII package (MSI, 1995) and then adding polar hydrogen atoms with the WHATIF program (v. 19990307-1124) (Vriend, 1990).

MEP calculations

The University of Houston Brownian Dynamics program, version 6.1 (Madura et al., 1995), was used to compute the MEPs of the pcs and cytf.

The N- and C-terminal residues were treated as charged. The ionizable residues were considered in their usual protonation states at pH 7. The histidine residues were considered to be either in their singly or doubly protonated forms, depending on their hydrogen-bonding environment. Atomic charges and radii were assigned to the protein residues from the OPLS (Jorgensen and Tirado-Rives, 1988) nonbonded parameter sets. The iron radius was set to 1.3 Å and the copper radius to 1.2 Å, as in the Amber (v. 4.1) parameter set (Weiner et al. 1986).

The formal charge of +2e for the Cu⁺² state in pc was partially redistributed over the copper ligands so that a charge of +1.5e was assigned to the copper in its oxidized state and the remaining charge of +0.5e was evenly distributed among the Cu ligands. Thus, in pc, a charge of +0.125e was added to the OPLS default values for each of the two His-ND1 Cu ligands and for the Met-SD and the Cys-SG Cu ligands (De Rienzo et al., 2000).

Similarly, the formal charge of +2e for the Fe⁺² state in cytf was partially redistributed over the iron ligands, so that a charge of +0.4e was assigned to the reduced iron and the remaining charge of +1.6e was distributed among its ligands. Thus, in cytf a charge of +0.267e was added to the OPLS default values for each of the four pyrrole nitrogens, which are the equatorial ligands to the iron, and a charge of +0.266e was added to each of the two iron axial ligands, His-25-NE2 and Tyr-1-N. Cys-14-SG and Cys-17-SG, which are covalently bound to the heme, were given the same charge values as the methionine sulfur atom.

A grid dimension of 110 × 110 × 110 Å³ was assigned, together with a 1.0-Å grid spacing, for computing the electrostatic potential by finite difference solution of the linearized Poisson-Boltzmann equation. The grid was centered on the global center of mass of the superimposed structures.

The dielectric constants of the solvent and the protein were set to 78 and 2, respectively. The dielectric boundary was determined from the van der Waals surface of the protein, and dielectric boundary smoothing (Davis and McCammon, 1991) was implemented.

The MEP grids were computed at 100 mM ionic strength, assuming a Boltzmann distribution corresponding to a temperature of 300 K.

Calculation of the SIs

For every pc mutant a, the MEP, _a(i, j, k), computed on a three-dimensional grid, was compared with the electrostatic potential of wt pc, _WT(i, j, k), by calculating the Hodgkin SI (Hodgkin and Richards, 1987). This index provides a measure of the similarity between the magnitudes and distributions of the MEPs of the two proteins compared.

(4)

in which the scalar products (M_a,M_WT), M_a², M_WT² represent the sum of products of values calculated for a given MEP M on grid points (i, j, k):

(5)

in which the summation is over the grid points within a defined region of interest. This region, called the "skin," is chosen to be at a distance  from the van der Waals surface of each protein and to have a thickness . Thus, the SIs for comparison of two proteins are calculated for grid points within the intersection of their skins. In our previous work (Wade et al., 2001; De Rienzo et al., 2000), we tested different values for  and , because these parameters influence the SI values obtained and the computational effort. The principal observation was that the choice of  and  mainly depends on the type of potential (i.e., electrostatic or hydrophobic) considered and analyzed. Therefore, we decided to set the two parameters  and  so that the region in which the potentials are compared is not so close to the protein surfaces as to be highly sensitive to small changes in protein structure and not so far that only the major components of the potentials (e.g., electrostatic monopoles) can be detected. In the present study, values of  = 3 Å and  = 4 Å were used to define the skins for the electrostatic potential analysis. The SI values lie in the range 1 to 1 with values near to 1 implying that the proteins have highly similar potentials on the overlapping skins, and values near to 1 implying that the potentials have inverted distributions.

Similarly, the MEPs of the wt and mutant pcs were compared with the MEP of the pc Q88E mutant. This comparison was carried out because, within the estimated experimental errors, this mutation does not alter the rate from that of wt pc: at pH 7.5, its overall rate constant k₂ is a little smaller than that of the wt, whereas at pH 6.0 it is a little larger. Thus, in this second set of calculations, the structure of the Q88E mutant, and not the structure of wt pc, was chosen as the reference.

The analysis of the molecular potential similarity was performed for two different regions: 1) the potentials were compared in the intersecting regions of the complete skins of the proteins and 2) the analysis was restricted to the so-called eastern site. For the latter, a vector was chosen having its origin at the center of mass of the superimposed structures and pointing toward the hydroxyl oxygen of Tyr-83, which lies in the middle of the eastern site. Only the intersection of those parts of the skins inside a conical region centered on the vector and with a 30° angular extent was considered. The details of the method have been described by Blomberg et al. (1999).

The index

(6)

in which M_b is either M_WT or M_Q88E (the MEP of the Q88E mutant), was computed from the SI values and used as a molecular descriptor of the difference between the MEPs of two proteins when evaluating correlations with k₂ rate constants. As the difference between MEPs of any two proteins (a and b) in the set is small due to the fact that only a few mutations were introduced in the wt primary structure, it is possible to assume that

(6`)

and therefore the index sqrt(2  2SI) is a measure of the differences in potential between the two proteins:

(6'')

The sqrt(2  2SI) values range from 0 to 2 with values near to 0 implying that the protein MEPs are highly similar and values near to 2 implying that MEPs are anticorrelated. When the index equals 1, the proteins have orthogonal MEPs.

Variations in the electrostatic potentials can provide an approximation to variations in the binding free energy among the mutants, because this is due to the product of each charge on a protein (a and b) with the electrostatic potential generated by the other protein (b or a) at a given position in the space. Thus, descriptors computed from SIs are directly related to the association equilibrium constants (K_A), which, in this set of mutants, are correlated to k₂.

Simulation of diffusional association by Brownian dynamics

The SDA program (Gabdoulline and Wade, 1998) was used to simulate the diffusion of pc toward cytf. The Brownian dynamics method implemented in the program allows one to simulate protein-protein encounter, compute the bimolecular diffusional association rates by solving the diffusion equation (Ermak and McCammon, 1978), and examine their dependence on protein mutations and on the nature of the physical environment. The theoretical background and the details of the methods have been described previously (Gabdoulline and Wade, 1996, 1998, 1999; Elcock et al., 1999). The methodology and parameters used in this work are identical to those in (Gabdoulline and Wade, 2001), in which they are described in full. It is important to note that these model and parameters can be applied, without modification, to compute bimolecular diffusional association rates for any given pair of proteins, under the primary assumption that the proteins may be treated as rigid bodies. We report a brief outline of the simulations and the most important parameters used.

Cytf and pc were treated as translating and rotating rigid bodies. The translational motions were simulated for pc relative to cytf, and therefore, the center of mass of cytf was positioned at the center of the simulation space. The translational and rotational diffusion constants of pc and cytf were calculated for the proteins (including their polar hydrogens), using a standard set of atomic masses implemented in the Macrodox program (version 3.2.1) (Northrup et al., 1987a,b). The translational and rotational diffusion coefficients of pc were computed to be 0.137 × 10¹ Ų/ps and 0.407 × 10⁴ ps¹, respectively, whereas the corresponding values for cytf were 0.127 × 10¹ Ų/ps and 0.327 × 10⁴ ps¹, respectively.

A continuum model of the solvent was adopted, and intermolecular forces and torques were given as the sum of repulsion and electrostatic forces. Short-range attractive interactions were not modeled because of their lesser importance for encounter complex formation. Short-range repulsions were treated by an exclusion volume that prohibits van der Waals overlap of the proteins. Electrostatic forces were computed by an approximation to the Poisson-Boltzmann forces (Gabdoulline and Wade, 1996). This approximation is as computationally efficient as a test charge model but permits the low dielectric regions of both proteins to be accounted for. Effective charges (Gabdoulline and Wade, 1996) were derived for both proteins by fitting to reproduce their Poisson-Boltzmann electrostatic potential in a homogeneous solvent dielectric. Then interactions were computed for the effective charges of each protein moving on the electrostatic potential grid of the other protein.

The electrostatic forces were approximated as the sum of a charge interaction (q) term and a charge desolvation term (Elcock et al., 1999) due to the approach of the low dielectric cavity of one protein to the effective charges of the other protein. A numerical desolvation factor of 1.7 was used for the desolvation term. An ionic strength of 100 mM was used for the simulations of all the mutants and the wt complexes, to reproduce the experimental parameters used by Kannt et al. (1996) to determine the overall rate constant k₂. The association of wt pc with cytf was also simulated at 200, 300, 400, 500, and 700 mM ionic strength to reproduce the experimental ionic strength conditions (Kannt et al., 1996).

At the beginning of each trajectory, pc was positioned with a randomly chosen orientation at a randomly chosen point on the surface of a sphere of radius b = 100 Å centered on the center of mass of cytf. Brownian dynamics simulations were then performed until pc diffused outside a sphere of radius q = 500 Å. Radius b was chosen so that at this distance the forces between the pc and cytf were centrosymmetric and isotropic; radius q was chosen so that the diffusive flux of pc was centrosymmetric.

During a trajectory, it is necessary to check for the formation of the diffusional encounter complex, which is the endpoint of the diffusional association phase (from this point, the proteins would rearrange into a bound complex). This is achieved by defining reaction criteria based on contact formation. The contacts monitored were those established between hydrogen-bond donors and acceptors, which are within a distance of 5 Å at the interface of the pc/cytf complex structure. During the simulations, the formation of any two independent contacts out of all the interprotein atomic contacts previously defined was monitored.

Four thousand trajectories were run for the wt and the mutant complex forms to have ~200 to 500 reaction events for the formation of the pc/cytf electrostatic complex, which means association constants computed with a 5 to 10% accuracy. In the case of the mutants E59K/E60Q and E59K/E60Q/E43N, 20,000 runs were necessary to obtain a similar number of reaction events during the simulations, because the electrostatic steering was much weaker. A variable time step was used: a time step of 1.0 ps was set when the proteins were positioned up to 50 Å apart, whereas at larger distances, the time step increased with increasing protein separation distance.

The SDA calculations were performed for structures in sets A and B. The SDA package distribution can be obtained at the internet address http://www-z.embl-heidelberg.de:8080/ExternalInfo//wade/pub/soft/SDA where a detailed explanation of the theoretical background of the average Boltzmann factors (computed to draw Fig. 1) can also be found.

View larger version (48K): [in this window] [in a new window]   FIGURE 1   (a) MEPs computed at an ionic strength of 100 mM for the wt and mutant spinach pcs. The isopotential contours are at 0.9 (red) and +0.9 kcal/mol/e (blue). The picture was generated with the InsightII program. (b) Electrostatic steering for pc and cytf. The contours enclose energetically favorable regions for each protein to reside in around the other, as derived from computation of the Boltzmann factors (BFs) of the interaction energies, defining the distribution of the center of one protein (2) around another protein (1) and vice versa. Plastocyanin is represented by a yellow ribbon, and the copper ion is represented by the sphere. Cytf is represented by a blue ribbon, and the heme is shown in a red stick representation. The pink contour shows the distribution of pc around cytf, and the green contour shows the distribution of cytf around pc. Contour levels for both pc and cytf contours are at 10.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

MEP SI analysis and comparison

Important differences in MEPs were found between the wt and mutant pcs (Fig. 1 a). The double and triple mutations (E59K/E60Q and E59K/E60Q/E43N) introduced in the small acidic patch significantly reduce the negative potential of that region. The mutants in the large acidic patch do not produce such a devastating effect. In fact, single and double mutations in this area appear to have less impact, probably because of the greater extent of the large patch with respect to the small patch. The Q88E mutation enlarges the negative potential without heavily modifying its shape.

The SI values calculated between each mutant and the wt pc are reported in Table 1 as are the experimental k₂ data from Kannt et al. (1996). In Fig. 2, the logarithmic values of the experimental rate constants (k₂) obtained at pH 7.5 and pH 6.0 for the mutants with respect to the wt (ln(k₂/k₂__WT)) versus the molecular descriptor sqrt(2  2SI), derived from SIs computed over the complete skins of the proteins (Fig. 2, a and b) and focusing on the acidic eastern site (Fig. 2, c and d), are reported. In all cases, good linear regressions are obtained (r² ~ 0.9). The trends indicate that the more similar the MEPs of wt and mutant pcs, the more similar the overall experimental rate constants (k₂). Moreover, the high correlation between the SI and SIes descriptors (r² = 0.99) indicates that almost all of the variance in the overall MEPs of the pc mutants is explained by changes in surface charge distribution occurring in a very restricted area of the eastern site centered on Tyr-83, although mutations were introduced in the more extended total acidic area of the eastern site (both the small and large patches).

View this table: [in this window] [in a new window]   TABLE 1   Experimental k2 values (Kannt et al., 1996) and In(k2/k2_WT) values at pH 7.5 and pH 6.0; similarity index values computed comparing the MEPs of the mutants relative to the wt pc on the complete skin of the proteins (SI and sqrt(2  2SI)) and on the "eastern site" (Sles and sqrt(2  2Sles))

View larger version (27K): [in this window] [in a new window]   FIGURE 2   The experimental k2 values (Kannt et al. 1996 for the mutants relative to the wt pc (ln(k2/k2_WT)) are plotted against the electrostatic potential SI, expressed as sqrt(2  2SI) and computed between the wt and mutant pcs over the complete skins of the proteins at (a) pH 7.5 with linear regression: ln(k2/k2_WT_pH7.5) = 4.29sqrt(2  2SI), r2 = 0.897, s = 0.424, n = 8, and omitting Q88E: ln(k2/k2_WT_pH7.5) = 3.95sqrt(2  2SI)  0.26, r2 = 0.972, s = 0.210, n = 7; and (b) pH 6.0 with linear regression: ln(k2/k2_WT_pH6.0) = 4.41sqrt(2  2SI) + 0.20, r2 = 0.902, s = 0.424, n = 8, and omitting Q88E: ln(k2/k2_WT_pH6.0) = 4.05sqrt(2  2SI)  0.07, r2 = 0.987, s = 0.142. The same relative rates are plotted against values of sqrt(2  2SIes) computed between the wt and mutant pcs on the eastern site at c, pH 7.5, with linear regression: ln(k2/k2_WT_pH7.5) = 3.88sqrt(2  2SIes) + 0.01, r2 = 0.903, s = 0.412, n = 8, and omitting Q88E: ln(k2/k2_WT_pH7.5) = 3.57sqrt(2  2SIes)  0.25, r2 = 0.974, s = 0.201, n = 7; and (d) pH 6.0, with linear regression: ln(k2/k2_WT_pH6.0) = 3.99sqrt(2  2SIes) + 0.22, r2 = 0.908, s = 0.412, n = 8, and omitting Q88E: ln(k2/k2_WT_pH6.0) = 3.66sqrt(2  2SIes)  0.06, r2 = 0.988, s = 0.136, n = 7. The experimental errors of k2 values are reported in Table 1.

The linear correlation coefficients in Fig. 2 are notably increased by omitting the Q88E mutant, which is underestimated by the regression models obtained. This is due to the fact that the wt and the Q88E mutant have very similar experimental rate constants (differences being within the estimated experimental errors) for their reaction with cytf at both pH 7.5 and pH 6.0 (Table 1). This implies that the addition of a negative charge in the proximity of the eastern site neither speeds up nor slows down the reaction in vitro. On the other side, the additional negative charge introduced in the large acidic patch affects, although slightly, the computed MEP of the protein. This perturbation is quantified by an SI value, which is not 1 (more precisely SI < 1). This means that we predict the Q88E mutant to be less reactive than the wt pc, whereas actually (experimentally) it is not.

A second set of descriptors was computed by selecting the Q88E mutant as the reference protein and comparing the MEPs of the wt and mutant pcs with respect to the MEP of the Q88E mutant. The computed similarity indices are reported in Table 2. In Fig. 3, the logarithmic values of the experimental rate constants (k₂) obtained at 1) pH 7.5 and 2) pH 6.0 for the wt and mutant proteins relative to the Q88E mutant (ln(k₂/k₂__Q88E)) versus sqrt(2  2SI) computed over the complete skins of the proteins are reported. The correlation at pH 6.0 (r² = 0.96, s = 0.26) is slightly better than that at pH 7.5 (r² = 0.95, s = 0.31) and both correlations show better statistics than the corresponding correlations in Fig. 2, a and b, obtained using the wt pc as the reference. The SI values computed with respect to the Q88E mutant by focusing on the eastern site (SIes and sqrt(SIes)) are also reported in Table 2. However, the correlation plots are not shown because poorer linear square correlation coefficients for ln(k₂/k_{2_Q88E}) versus sqrt(SIes) are obtained (r² = 0.77 at both pH 7.5 and pH 6.0).

View this table: [in this window] [in a new window]   TABLE 2   Experimental k2 values (Kannt et al., 1996) and In(k2/k2_Q88E) values at pH 7.5 and pH 6.0; similarity index values computed comparing the MEPs of the mutants and the wt pc relative to the Q88E mutant on the complete skin of the proteins (SI and sqrt(2  2SI)) and on the "eastern site" (Sles and sqrt(2  2Sles))

View larger version (13K): [in this window] [in a new window]   FIGURE 3   The experimental k2 values (Kannt et al., 1996) for the proteins relative to the Q88E mutant (ln(k2/k2_Q88E)) are plotted against values of sqrt(2  2SI) computed between the Q88E mutant and the other pc variants over the complete skins of the proteins at (a) pH 7.5 with linear regression ln(k2/k2_Q88E_pH7.5) = 4.12sqrt(2  2SI) + 0.38, r2 = 0.947, s = 0.305, n = 8; (b) pH 6.0 with linear regression ln(k2/k2_Q88E_pH6.0) = 4.26sqrt(2  2SI) + 0.38, r2 = 0.964, s = 0.258, n = 8. The experimental errors of k2 values are reported in Table 2.

Analysis of the structures of the pc/cytf complex

The 10 structures constituting the NMR-based ensemble deposited as 2pcf in the pdb were analyzed to select the most suitable structure for the Brownian dynamics simulations. With the Brownian dynamics method implemented in the SDA program, the association rate constant k_on is calculated by simulating the formation of the diffusional encounter complex, in which the electrostatic interactions are expected to be optimized. Therefore, by "suitable structure," we mean the structure that can best represent the electrostatic complex that is formed in the first step of the overall electron transfer reaction between pc and cytf. Experimental data (Hope, 2000; Illerhaus et al., 2000; Bendall et al., 1999; Gong et al., 2000b; Redinbo et al., 1994) and computational studies (Roberts et al., 1991; Pearson et al., 1996; Pearson and Gross, 1998; Ullmann et al., 1997) show that both the small and the large acidic patches on the pc eastern site are involved in intermolecular interactions in the electrostatic complex. From this perspective, the most suitable NMR-based structure appears to be that deposited as number 9 in the 2pcf pdb coordinate file, in which the highest number of electrostatic and hydrogen-bonding interactions, involving both the large and the small acidic patches, stabilizes the pc/cytf complex. Structure number 9 was thus selected and used to build a set of mutant structures of pc/cytf complexes (sets A), as described in the Materials and Methods section. To optimize further the interactions at the pc/cytf interface, structure 9 underwent an MD run. The resulting structure was then used to build a second set (set B) of mutant pc/cytf structures (see Materials and Methods for details).

From analysis of the MD simulation trajectories, it is evident that the most flexible regions in the structure of the complex are the regions of spinach pc around the large (42-45) and the small (59-61) acidic patches and around Asp-9. The overall RMSD computed by superposing the complete C traces of the wt structures of cytf and pc in the complex in set B onto the C traces of the original structure of the complex (wt in set A) is 4.01 Å. The internal RMSD of the C atoms of pc in the two structures is 1.79 Å, and that of cytf is 3.95 Å. The larger value for cytf is mainly due to changes in the relative orientation of its two domains. Furthermore, comparison of the original structure of the complex and that from MD showed that the conformational differences cause differences in side chain packing at the pc/cytf interface of the two proteins. In fact, a higher number of interactions stabilizes the MD complex, thus causing tighter packing of pc and cytf than in the original complex structure.

Brownian dynamics simulations

Two sets of BD simulations were run with the SDA program using the structures in sets A and B for comparison purposes. The calculated k_on (k_C) values for the mutant and wt pcs are reported in Table 3. The k_C values computed for the encounter complexes, defined by two contacts at a distance of 6.0 Å between the donor and acceptor groups at the interface in the complex, were compared with the experimental k₂ values.

View this table: [in this window] [in a new window]   TABLE 3   Experimental k2 values (Kannt et al., 1996) and In(k2/k2_WT) values at pH 7.5 and pH 6.0; computed association rate constants kC and In(kC/kC_WT) for complexes in set A (kCA) and in set B (kCB) and their values relative to the experimental k2 at pH 7.5 (kC_rel_7.5_A,B = kC_A,B/k2_7.5) and pH 6.0 (kC_rel_6.0_A,B = kC_A,B/k2_6.0). kC values reported here were computed for the encounter complexes defined by two contacts at a distance between donor and acceptor groups of 6.0 Å

The logarithmic values of k_C obtained for the structures in sets A (Fig. 4, a and c) and B (Fig. 4, b and d) relative to the k_C value obtained for the wt (lnk_C/k_{C_WT}) were compared with the logarithmic values of the experimental k₂ relative to the k₂ value determined for the wt (lnk₂/k_{2_WT}) at pH 7.5 (Fig. 4, a and b) and pH 6.0 (Fig. 4, c and d). The linear regressions reported in Fig. 4, a and b, for pH 7.5 show good correlation coefficients (r² = 0.82 and r² = 0.89, respectively). However, the square correlation coefficients improve notably for the correlations reported in Fig. 4, c and d, (r² = 0.86 and r² = 0.92, respectively), in which the experimental k₂ values at pH 6.0 are used for comparison with the calculated association rate constants. This is due to the fact that at pH 6.0 both the k_C and k₂ values of the Q88E mutant are higher than that of the wt, whereas at pH 7.5, the k_C value of the Q88E mutant is higher but the k₂ value is lower than the corresponding values of the wt protein.

View larger version (26K): [in this window] [in a new window]   FIGURE 4   Comparison of relative overall rate constants determined experimentally with the association rate constants computed by Brownian dynamics simulation. The experimental k2 values (Kannt et al., 1996) of mutant pcs relative to the wt protein (ln(k2/k2_WT)) are plotted against corresponding relative computed association rates (ln(kC/kC_WT)). The calculated kC values were obtained (with encounter complexes defined by two contacts at a 6.0 Å or smaller separation) for structures in (a) set A at pH 7.5: ln(k2/k2_WT_pH7.5) = 0.72ln(kC/kC_WT_A)  0.74, r2 = 0.823, s = 0.538, n = 8; (b) set B at pH 7.5: ln(k2/k2_WT_pH7.5) = 1.12ln(kC/kC_WT_B)  0.62, r2 = 0.894, s = 0.431, n = 8; (c) set A at pH 6.0: ln(k2/k2_WT_pH6.0) = 0.76ln(kC/kC_WT_A)  0.53, r2 = 0.864, s = 0.479, n = 8; (d) set B at pH 6.0: ln(k2/k2_WT_pH6.0) = 1.17ln(kC/kC_WT_B)  0.42, r2 = 0.925, s = 0.372, n = 8. The experimental errors of k2 values and the estimated errors of kC are reported in Table 3.

Although good estimates of relative association rates were obtained, the calculated association rates are generally overestimates of the experimentally determined overall rate constants for the wt pc/cytf complex (see Table 3). The overestimation factor for the wt pc is ~15 in both sets A and B. The estimated error on this factor, computed from a large series of Brownian dynamics simulations, is approximately ±5. Taking this into account and looking at the series of data, it is evident that not all the mutants are similarly overestimated. In fact, the overestimation factor varies from 8 to 49 for the different datasets and mutants. Furthermore, in structure set A, the mutants in the large acidic patch are clearly more overestimated than those in the small acidic patch. On the other hand, there is not such a clear difference between overestimations for the small and large acidic patch mutants in set B (Table 3).

This difference between sets A and B is probably due to the different relative importance of hydrogen-bonding interactions in the small and in the large acidic patches in stabilizing the complexes of the two sets (Table 4). Considering the complete set of SDA simulations, hydrogen-bonding interactions in the small acidic patch appear many more times than those in the large acidic patch in the encounter complexes formed by structures in set A. On the other hand, the interactions in the large acidic patch appear to be much more relevant for encounter complex formation in set B (see Table 4). This means that mutating a residue in the small acidic patch in structure B of the wt complex has consequences for the computed k_on values that are less drastic than the corresponding mutation in structure A of the wt complex. In contrast, mutating a residue in the large acidic patch in structure B of the wt complex has consequences of greater impact than the corresponding mutation in structure A of the wt complex.

For comparison, Brownian dynamics simulations were also run using a third set of structures (set C). These structures were built by selecting the first structure deposited in 2pcf and mutating the appropriate residues to obtain its mutant structures in the same way as for set A. Structure number 1 was selected because it seemed a good representative of all other structures in the 2pcf pdb file with the exception of structure number 9. The linear regressions obtained with set C (ln(k₂/k_{2_WT_pH7.5}) = 0.91ln(k_C/k_{C_WT})  0.78 and ln(k₂/k_{2_WT_pH6.0}) = 0.96ln(k_C/k_{C_WT})  0.56) have square correlation coefficients (r² = 0.79 and r² = 0.84, respectively) that are worse than the corresponding ones obtained for structure set A (Fig. 4, a and c). Moreover, the factor by which k_C is overestimated with respect to k₂ is larger (16-69) than that computed for set A structures, and the large acidic patch mutants show k_C values that are much more overestimated than for set A (data not reported).

What appears to be common to these three sets of structures (A, B, and C) is the difficulty of reproducing, with calculations, the experimental rates of the mutants in the large acidic patch. A possible explanation for this is that the pc large acidic site is quite far from the redox site (Cu site) and therefore, this region is modeled in the pc/cytf complexes without the aid of well-defined NMR-based constraints. This implies that the structure at the pc large patch interface is not as reliable as that at the small patch interface with cytf. A second explanation might be found in the MEP distribution of the large patch mutants. As we mentioned previously, introducing mutations in the large acidic patch does not disturb the wt pc MEP as much as introducing mutations in the small acidic patch, perhaps because of the different extent of the two regions (Fig. 1 a).

Fig. 1 b shows the electrostatically and sterically most probable positions occupied by the center of mass of pc around cytf (represented by a pink solid contour), when cytf is kept fixed, and by the center of mass of cytf around pc (green solid contour), when pc is kept fixed. The contoured values are average Boltzmann factors (BFs) computed with the SDA program. BFs can be used to represent the electrostatic steering because it has been shown that the diffusional association rate enhancement due to electrostatic interactions is proportional to the average BF in the region of the transition state for binding (Zhou et al., 1998). BFs were computed only for the wt complex in set A. It is evident that it is most favorable for wt pc to be around the cytf region where the heme site and the basic area are found. Similarly, cytf mainly interacts with the highly negative region found in pc and commonly known as the "eastern site." It is possible to deduce that there is a directional and reciprocal local steering of the cytf basic area and the pc acidic area toward each other, which supports the hypothesis of formation of an initial electrostatic complex involving both the large and small acidic regions.

Finally, we computed the k_on values for the wt complex in set A at six different values of the ionic strength I (with I  100 mM), and we compared the variation of the computed k_on and the experimental k₂ rate constants (Kannt et al., 1996) as a function of I. Good agreement in the ionic strength dependence of the computed and experimental values is found as shown in Fig. 5.

View larger version (10K): [in this window] [in a new window]   FIGURE 5   Ionic strength dependence of experimental overall rates and computed association rates. Experimental lnk2 values at pH 6.0 (Kannt et al., 1996) determined at six different values of ionic strength () and lnkC values computed for the structure of the wt complex in set A at the same six values of ionic strength () are plotted against the square root of the ionic strength (sqrtI). The linear correlations are, respectively: lnk2 = 5.88sqrtI + 20.43, r2 = 0.910, s = 0.384, n = 6; lnkC = 5.77sqrtI + 23.50, r2 = 0.998, s = 0.005, n = 6.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

MEP comparative analysis

The MEPs of spinach pc and its mutant forms, Q88E, D42N, E43N, E43K, E43Q/D44N, E59K/E60Q, and E59K/E60Q/E43N, were compared by computing MEP SIs. The sqrt(2  2SI) value is a quantitative measure of the difference in the magnitude of the electrostatic interaction potential between the two proteins compared. Assuming that both the proteins compared interact with the same redox partner with overall similar features, this descriptor will relate to binding energy and, via the Boltzmann factor in the region of the transition state (Zhou et al., 1998) to the electrostatic enhancement of association rates. This relation is born out in our results (see Fig. 2), which show that the more the wt pc MEP distribution is perturbed by mutations, the more the reaction with cytf is hampered. This means that the overall rate constants for the studied reactions of wt and mutant pcs with cytf are mainly affected by electrostatically steered association, at least when mutations are introduced in the eastern site region. These considerations are in complete agreement with the experimental results by Kannt et al. (1996). In fact, the experimental k₂ and K_A values that they found for the wt pc and the mutants D42N, E43N, E43K, E43Q/D44N, and E59K/E60Q interacting with cytf (the value of K_A for E59K/E60Q/E43N is not given in the paper) are highly correlated (r² = 0.974). This suggests that the modulation in the trend of the overall electron transfer rate constant (k₂) between wt pc and cytf is mainly due to the mutation-affected binding constant (K_A), whereas the actual electron transfer rate constant (k_et) is almost unaffected by mutation or only affected by a constant factor. The unique exception is the mutant Q88E, whose k₂ value is highly similar to that of the wt pc, whereas its K_A is double that of the wt. This might be ascribed to a decreased electron transfer rate, probably as a consequence of the fact that the increased binding of pc and cytf hinders the rearrangement of the complex to the configuration most suitable for electron transfer. Pure electrostatic considerations are not sufficient to predict correctly the k₂ value of Q88E, and this will be discussed further on with the aid of computed association constants (k_on) for the complex cytf/pc. SIs, in fact, can be considered as approximate measures of the K_As (binding free energies), which are defined by k_on/k_off. Therefore, variation of the k_off value seems to play an important role for the Q88E mutant.

We can conclude that, when one or more negatively charged residue(s) in the pc eastern site is mutated into a positively charged or polar residue, the electrostatic properties of the eastern face change, reducing its ability to attract the cytf basic patch and consequently influencing the overall reaction rate. These calculations and experiments relate to the interaction of spinach pc with the soluble part of turnip cytf. Nevertheless, recent studies of Illerhaus et al. (2000) provide experimental support (consistent with our calculations) for the importance of the acidic patch pc residues for the overall electron transfer between spinach pc and its intact biological partner, the cytbf complex from spinach.

The Q88E mutant is an outlier in all the models based on SIs shown in Fig. 2, always being underestimated by the linear equations obtained. In fact, although the experimental k₂ values determined for the Q88E mutant are, by considering the experimental error (Kannt et al., 1996), only slightly higher than those obtained for wt pc, at both pH 7.5 and pH 6.0, its MEP is significantly affected by the addition of a negative charge in the proximity of the acidic patch and therefore SI(Q88E, wt) is <1. For these reasons, predicting the features of the Q88E mutant is not straightforward. When this mutant is omitted from the dataset, the linear correlation coefficient increases significantly in all the models (Fig. 2). Taking a different approach, we considered the Q88E mutant, and not the wt, as the reference protein. This can be justified by considering that, at pH 6.0, the k₂ value of the Q88E mutant, taking into account the value of the experimental errors, is slightly higher than that for wt. Therefore, it is the most efficient variant in the interaction with the redox partner and can be taken as a template to design new mutants. The SIs were thus computed by comparing the MEPs of the wt pc and its variants with that of the Q88E mutant, and similarly, the experimental k₂ values of the wt and the mutant pcs were compared with the Q88E mutant reaction rate constants (Table 2; Fig. 3). The graph (Fig. 3 b) obtained using the experimental values at pH 6.0 shows particularly good agreement between the molecular descriptor sqrt(2  2SI) and the experimental values. Here, the Q88E variant is well predicted by the model.

In contrast, the predictive ability of the descriptor sqrt(2  2SIes), limited to the electrostatic properties at the eastern site, is notably lower than that obtained using the wt pc as the reference protein. It is worth remembering that in the other mutants and the wt pc, the residue at position 88 is a Gln and not a Glu. This implies that comparing, for example, Q88E with the single point mutant D42N, we are actually comparing the reference protein with a double mutant of the reference protein. The consequent differences in protein MEPs are clearly more effective on SIes values, computed comparing the MEPs on the restricted eastern site where the mutations are located, than on the SI values, computed comparing the MEPs on the complete protein surface.

In summary, from a methodological point of view, it is reasonable to choose the reference protein as the protein for which the highest or lowest activity has been measured, independently of whether it is the wt or a mutant form. A drawback of this approach is that further experiments, leading to an amplification of the protein set initially used to build the model, might require the definition of a new reference mutant whose activity is now the highest (or lowest). In this case, the model should be rebuilt. This limits the predictive capabilities of the approach. In the case presented here, our main interest is to design mutants with altered function compared with the wt pc, and therefore, the use of wt as the reference protein still appears to be the best choice. As a whole, computation of MEP SIs as described here provides an efficient way to estimate approximately relative association rates when these are dominated by changes in the electrostatic potential of the mutated protein as well as a useful tool for interpretation of binding data for a set of protein mutants.

Brownian dynamics simulations

The association of wt and mutant pcs with cytf was simulated, taking into account the intermolecular long-range electrostatic forces, to compute association rate constants (k_C) of the partner proteins. Two different sets of protein structures, set A and set B, were studied. These were derived from the NMR-based structure of the pc/cytf complex (pdb file 2pcf) and are different from the structures used for the MEP analysis, which used the crystal structure of pc (pdb file 1ag6). Mutants in the two sets have different conformations due to differences in the wt structures used, as explained in the Materials and Methods section. In the wt complex of set A (wt A), the most effective interactions at the interface between the two partners are those in the pc small acidic patch. On the other hand, in the wt complex that underwent MD and energy minimization (wt B), a larger number of contacts was established at the interface of the pc/cytf complex, especially in the large acidic patch, with respect to the original wt structure (wt A). Conformational changes are due to modifications of the hydrogen-bond network, which stabilizes the complex, during the MD simulation. This can be clearly seen by comparing the lists (reported in Table 4) of hydrogen-bonding donor and acceptor atoms at the pc/cytf interface in wt A and wt B. The differences between the two structures of the complex are reflected in the total pc/cytf interaction energies (computed using CHARMM22): 108.3 kcal/mol for wt A and 285.3 kcal/mol for wt B. Furthermore, the Cu atom and the heme Fe atom are 11.4 Å apart in wt A, whereas they are further apart (12.4 Å) in wt B, implying that the conformation of this second complex structure is less suitable for electron transfer. The differences noted between the structures in the two sets affect the Brownian dynamics simulations, as shown by the estimated association rate constants k_C (Table 3) and the contacts determined for encounter complex formation (Table 4). The k_C values computed within the set B series reproduce the trend within the experimental k₂ values better than the k_C values computed using structures in set A. This is again probably due to the different conformations and better energetic stabilization of the set B complexes.

As a whole, all the models proposed in Fig. 4 and resulting from studies on both structure sets A and B allow the prediction of the overall rate constant values from the computed association rate constants for the different pc variants in interaction with wt cytf. The good predictive ability of the models, in agreement with the results obtained above with the SI descriptors, highlights the fundamental role of electrostatics not only in promoting the formation of the pc/cytf encounter complex, but also in characterizing the overall electron transfer process.

It is worth noting that the two sets of descriptors used (SIs and k_C) are correlated. The index sqrt(2  2SI) explains 91% of the variance of the index ln(k_C/k_{C_WT}) computed for set B structures, and 85% of the variance of ln(k_C/k_{C_WT}) computed for set A structures. These results suggest that, whereas in both cases electrostatic interactions are important for steering the protein partners into the correct electron-transfer orientation, they are more effective for the formation of complexes in set B than in set A.

The importance of electrostatics as steering interactions is further confirmed by the dependence of the k_C and k₂ values on the ionic strength: above 100 mM ionic strength, the overall rate constant values are reduced and so are the computed association rate constants by increasing the ionic strength of the environment (Kannt et al., 1996; Illerhaus et al., 2000). At ionic strengths below 100 mM, it is difficult to determine the overall rate constants for the wt pc/wt cytf complex experimentally, probably because of rate saturation effects (Kannt et al., 1996). This suggests that, at low ionic strength, pc and cytf may get stuck in the electrostatically favorable complex, and the transition to the electron-transfer optimal conformation becomes harder. This idea is supported by experimental evidence (Kannt et al., 1996) that overall rate saturation phenomena are not evident for the interaction with cytf of those pc variants (such as E59K/E60Q and E59K/E60Q/E43N) in which the magnitude and the distribution of the eastern site electrostatic potential are most perturbed with respect to the wt protein. It is also consistent with results obtained from kinetic analysis of the overall electron transfer process in the spinach pc/cytbf complex (Illerhaus et al., 2000), that show that removal of negative charges in the eastern site leads to a marked decrease of the overall ele