J. Phys. Chem. 1996, 100, 19241-19250
19241
ARTICLES Molecular Dynamics of Microperoxidases in Aqueous and Nonaqueous Solutions S. Melchionna,*,†,‡ M. Barteri,† and G. Ciccotti‡ Dipartimento di Chimica and INFM-Dipartimento di Fisica, UniVersita` “La Sapienza”, 00185 Roma, Italy ReceiVed: March 8, 1996; In Final Form: July 1, 1996X
We report an extended molecular dynamics study of the heme-containing peptide “microperoxidase” (MP), a model system of several heme-proteins. The most commonly studied forms of MP, the octa- and the undecapeptide, are studied both in pure water conditions and in mixed solvent (80% methanol and 20% water) conditions that stabilize the monomeric form. The equilibrium structures and their stability are visualized by graphical tools and analyzed in terms of the molecular shape, the vibrational amplitudes, and torsional angles along the backbone. The obtained structures are used to confirm the experimental data on aggregation of Urry and Pettegrew7,8 by using simple geometric arguments based on the possible assemblies of several monomeric units in space. We have investigated in detail the H-bond network and the influence on the heme conformational properties. We have found that heme interacts with the peptide through H-bonds formed by the propionates, the axial histidine, and different partners of the peptidic chain. These interactions are regulated by the degree of exposure of HIS to the solvent. These results permit us to predict that the undecapeptide and the octapeptide can be good model systems for horseradish peroxidase and cytochrome c′, respectively.
1. Introduction Microperoxidases (MP) are small heme peptides obtained by hydrolysis of cytochrome c with pepsin and trypsin.1,2 The two most commonly studied species of microperoxidases are the undecapeptide (MP11), obtained by peptic digestion of cytochrome (cyt) c at pH ) 2.5, and the octapeptide (MP8), the subsequent product obtained by hydrolysis of MP11 with trypsin. MP8 and MP11 exhibit the heme group bound respectively to 11 and 8 aminoacidic residues of the original primary sequence of cytochrome c, from Val11 to Glu21 (Figures 1 and 2).3 The prosthetic group is attached to the peptidic chain by two thioether linkages to the R-carbons of the saturated vinyl groups of Cys14 and Cys17 on the two pyrrole moieties. Microperoxidases have been widely studied in the past as model systems for the active sites of hemoproteins.4,5 In particular they have been used as models for several hemebased enzymes, as catalysts of oxygen transfer (cytochrome P-450) and redox reactions (peroxidases), offering the possibility of studying pentacoordinated complexes as well as sixcoordinated complexes.6 Unfortunately we still ignore the detailed structure of these molecules since it is not yet possible to obtain suitable crystal samples to perform X-ray diffraction measurements on these systems due to the high flexibility of the peptide. In contrast, the hydrophobic character of the porphyrinic ring produces aggregates with large molecular weight in aqueous solution at alkaline and neutral pH.7,8 The different coordinations of Fe in polymeric MP and the fact that in the aggregates it is found that the iron atoms form a mixture of different spin states are some of the limitations in obtaining precise experimental data using different spectroscopic techniques. * To whom correspondence should be addressed. † Dipartimento di Chimica. ‡ INFM-Dipartimento di Fisica. X Abstract published in AdVance ACS Abstracts, October 1, 1996.
S0022-3654(96)00727-7 CCC: $12.00
Figure 1. Schematic representation of microperoxidases. The amino acid sequence from 11 to 21 of cytochrome c is shown, indicating the covalent linkages of Cys14 and Cys17 to the heme group, via the two vinylic groups (CH-CH3), and the His18 to iron attachment. In brackets are shown the residues present in the undecapeptide only. In the position opposed to the peptide (sixth coordination position of Fe) an imidazole molecule is coordinated to the iron atom. The porphyrin core is composed of four pyrrolic rings (hexagons) and four β-carbons rings (pentagons). In the inset is represented the generic aminoacid block unit. Bold line: backbone linkages. R: lateral group specific of the amino acid. φ,ψ: torsion angles around the N-CR and CR-C bonds, respectively. φ is the dihedral angle formed by the two planes passing through the triplets of atoms C(i-1) - N(i) - CR(i) and N(i) - CR(i) C(i), while ψ is the dihedral angle formed by the planes passing through N(i) - CR(i) - C(i) and CR(i) - C(i) - NR(i+1), respectively, where italic indexes refer to the aminoacidic sequence.
This paper describes the structure and stability of microperoxidases by use of molecular dynamics (MD), one of the most important theoretical tools to investigate the structural and dynamical properties of microscopic biological systems.9,10 We © 1996 American Chemical Society
19242 J. Phys. Chem., Vol. 100, No. 50, 1996
Figure 2. 3-D representation of the heme plane and peptidic chain. In the heme plane are in evidence the pyrrolic nitrogens, in the text often referred to as Npyrr, the four propionic oxygens, the imidazole (upper), and the histidylimidazole (lower) with its (NH) group (referred to in the text as (NH)His).
Figure 3. Sketch of possible MP monomer alignment. Left: irregular head-to-tail arrangment. Right: regular stacked arrangment. The thioether and histidine attachments to iron are represented by a solid segment. The dashed segments represent the linkages between different monomic units through coordination in the sixth position to the iron atom.
have focused our study on the monomeric form of microperoxidases since a clear knowledge of the aggregation patterns is still lacking. We have studied the conformational behavior of both the undecapeptide and octapeptide, and we have performed a comparative study of microperoxidases in pure water (MP11W, MP8W) and in the presence of mixed solvent (20% water and 80% methanol, obtaining the systems MP11M and MP8M). It is well-known that in the presence of organic solvent the effective attraction between different heme planes is reduced due to the lower dielectric constant so that, if the methanol concentration is high enough, the monomeric form is recovered. In our study we were interested in describing the behavior of microperoxidases in organic solvent since this is the typical experimental condition. One of the early questions about MP was the interpretation of the experimental data of aggregation experiments performed by Urry and Pettegrew in the 1960s7,8 by using circular dichroism. These authors concluded that the undecapeptide exhibits preferentially an irregular alignment (so-called headto-tail) of the heme planes, while the octapeptide aggregates in a more regular stacked or oblique conformation (see Figure 3). Their conclusions were drawn by guessing the structures of the monomers in pure water with the help of geometrical models (based on excluded volume or the hydrophobic content of the residues) and by observing that the Fe atom could be coordinated only to the available aminic groups of the backbone side groups. With this simple method those authors have predicted the possible assemblies in space of several monomers. In this paper we have used the detailed structures obtained with MD to check the early hypothesis of Urry and Pettegrew, and we have confirmed that the undecapeptide seems to aggregate in an
Melchionna et al. irregular fashion, while the octapeptide is well suited to pack in a stacked fashion. A detailed description of the physiological behavior of heme in very different heme-proteins is still lacking. Microperoxidases are often used as model systems of heme-proteins, but the electrochemical and magnetic properties, as measured by cyclic voltammetry11,12 and resonant Raman measurements,15 of these relatively simple molecules are still unexplained. The catalytic properties of the porphyrin are strictly dependent on the distribution of the electronic charge among the heme plane, the axial ligands, and possibly other peptidic partners. The charge arrangement around the iron atom depends on the degree of engagement of the proton of the (N-H)His group in H-bonds with solvent or peptidic partners. In fact a proton reduction of the axial histidine, i.e. a negative charge surplus on the iron atom, would stabilize the lower oxidation state of Fe and lower the reactivity of heme. As an example, microperoxidases have shown to be charge carriers by undergoing direct electron transfer at glassy carbon electrodes of cyclic voltammetry cells, and measurements of the redox potential of MP11 as a function of pH exhibits reversible redox behavior at neutral pH range, while at alkaline and acid pH range the reactivity of MP11 is lowered.11,12 In this paper we have faced the problem of obtaining qualitative indications on the redox behavior of these molecules by using the statistical data on the classical system, i.e. without considering electronic degrees of freedom or polarization effects, as obtained by standard MD. From our data we have observed that MP8 and MP11 show a very different degree of H-bonding of (N-H)His with peptidic or solvent H-bond acceptors. For MP11 the axial histidine is weakly engaged in H-bonds, while for MP8 stable H-bonds stabilize the lower oxidation state of Fe, suggesting that MP11 is the most reactive species. The formation of H-bonds involving the (N-H)His group was previously observed by use of resonance Raman spectroscopy13 and in our previous MD simulation of MP11 in vacuo and in pure water solvent.14 Moreover, resonant Raman studies of monomeric N-acetyl-MP8 showed the presence of a mixture of intermediate and high-spin states,15 i.e. with the Fe atom lying slightly out of the heme plane or completely out-of-plane. The mixture of different spin states is again usually referred to the formation of the same H-bond involving the histidylimidazole. This H-bond formation weakens the Fe-His bond and therefore favors the high-spin state. In this paper we find that MP8 is favored in a high-spin state, due to charge rearrangement after H-bond formation. In contrast, the spin state of MP11 depends both on the heme distortion, due to the cysteines’ covalent linkages and attraction of the propionates toward the peptidic chain, and on the stronger Fe-His bond. The result of these two competing forces on the position of the iron in MP11 is, on the basis of our classical model, unpredictable. In section 2 we present the microscopic model and summarize the technical details used to perform our simulations. In section 3 we report and discuss the results of simulations separated in four different subsections, characterization of the monomer conformation and stability, analysis of the aggregation patterns resulting from the obtained monomer structures, description of the H-bond network observed during the trajectories, and heme conformational properties related to H-bond formations. In section 4 we collect a few concluding remarks. 2. Computational Model and Details The molecular dynamics simulations were performed using the software package GROMOS developed by van Gunsteren and Berendsen (a detailed description of the GROMOS treat-
Molecular Dynamics of Microperoxidases ment of biological molecules may be found in ref 16). Our microscopic model of microperoxidases is purely classical. Electronic degrees of freedom are not taken into account, nor is charge redistribution, due to polarization effects, considered. The complex electronic charge arrangement in the heme plane and neighboring sites is not well-known, in particular for coordination of heme with different axial ligands, as is the case of imidazole. In our model, restricted to studying a system with frozen electronic distribution, as well as frozen bond lengths, the redox behavior of the molecule can be inferred indirectly through simple charge redistribution mechanisms. Since in our study we were interested in detecting the intramolecular H-bond network, we have adopted the GROMOS force field, which is well suited for modeling such an interaction. In fact, although no explicit hydrogen-bonding potential is introduced, in previous investigations it has been shown that the Coulombic potential and properly rescaled van der Waals potential parameters, chosen by the GROMOS force field parameter set, effectively represent correctly the H-bond interaction network.17 We refer to ref 14 for a detailed description of the most important approximations introduced in the computational model, while in the following we recall the main technical features for the sake of completeness. The long-range electrostatic interactions were treated using the “neutral charge group” method. Contrary to standard methods that compute the charge-charge interaction terms directly, as in the Ewald sum method, in the neutral charge group method the usual Coulombic force exerted between two pointcharged atoms is replaced by a dipole-dipole interaction between groups of atoms, neutral as a whole.10 This choice is consistent with the periodic boundary conditions and minimum image convention applied to the simulated system since the longrange dependence of the Coulombic potential (∝ 1/r) is replaced by the faster dipole-dipole interaction term (∝ 1/r3). The neutral charge group approach is widely used in the simulation of biopolymers since for peptides we deal mostly with localized neutral charge groups.10 The hydrogen atoms attached to carbons were not considered explicitly, but we used the modified set of force field parameters consistent with this choice.10,16 The water molecules were represented by the single point charge (SPC) model,18 where the molecules are modeled as rigid bodies with three force field centers. The methanol molecules were modeled in the same manner, by replacing a proton with the heavier CH3 united atom. We used the charged version of the GROMOS force field parameters to model the nonbonded interactions in both aqueous and nonaqueous solvent. This empirical set was used for the atoms of the backbone and of the porphyrin without any modification. The cutoff of the forces for the nonbonded interactions was taken at 8 Å for all the simulations without using a neighbor list. To avoid discontinuities in the Hamiltonian, and consequent spurious heating of the systems during the simulations, we have used a polynomial switching function for the nonbonded interactions. The nonbonded interactions have been computed up to a cutoff of 10 Å, where the switching function goes smoothly from 1 at 8 Å to 0 at 10 Å. The chosen cutoff and switching function are compatible with the shortrange form of the van der Waals interaction and with the dipole-dipole Coulombic interactions. With this choice we have not observed any drift in the total energy during all the simulations. Since the physical system allows us to use the neutral charge group approach instead of the complete form of the Coulombic potential, in our computational model we implicitly neglect the
J. Phys. Chem., Vol. 100, No. 50, 1996 19243 effects of long-range electrostatic interactions on the equilibrium properties of microperoxidases. Therefore we have not performed a systematic study of the dependence of the results with the cutoff for the equilibrium conformations, but we have chosen the cutoff that furnished an optimal computational efficiency together with reliable results. The simulation cell was built with a parallelepiped shape and with a size such that one has enough solvent molecules so that, when applying the periodic boundary conditions, the solute molecule weakly interacts with its images (in all the simulations the box edges were always superior to 20 Å). We kept all bond lengths fixed in time by use of the SHAKE algorithm.19 This tool allows saving of CPU time while it is known that the equilibrium properties are not modified.20 All the simulations were performed at constant energy and volume. The so-called constant pressure-constant temperature algorithms present in GROMOS21 were used only to set up the initial temperature (300 K) and internal pressure (1 atm) of the system before beginning the equilibrium runs. We took as a starting configuration of MP11 and MP8 in water the crystallographic coordinates of the corresponding residues and prosthetic group of reduced bonito cytochrome c (entry pdb1cyc of the Brookhaven Data Bank31). Energy minimization and subsequent temperature rescaling were carried out to obtain a relaxed system. For the starting structure of the solute in mixed solvent conditions, we took the equilibrium structures obtained from the last step of simulation in pure water and changed the force field parameters of 80% of randomly selected water molecules to obtain the mixture of methanol and water. Then we performed energy minimization and equilibrium runs. For MP11 and MP8 we used 539 and 459 solvent molecules, respectively. The equilibrium box volumes were found to be 18.3 and 14.1 nm3 for MP11 and MP8 in water and 30.6 and 44.5 nm3 for the MP11 and MP8 in organic solvent, respectively. The equations of motion were integrated in time with a time step of 2 fs, for a total of 100 ps. This time span is usually too short to sample properly the configurational space of peptides, but for microperoxidases the three covalent linkages between the peptide and the heme group (Figure 2) infer a substantial rigidity to the solute molecule so that with this trajectory length a proper sampling can be achieved. With this length of the trajectories we have always observed during the equilibrium simulations the stationary behavior of relevant mechanical properties of the system such as the potential and kinetic energy and the virial and the configurational distance (compare eq 1 and Figure 5 in the following discussion). The conformational evolution of the peptides has been monitored in time in order to determine whether the stability of the mechanical observables corresponded to an effective conformational rigidity. In Figure 4 a sequence of conformations of MP11W and MP8W taken at regular intervals during the equilibrium runs for the peptidic backbone and heme group are reported. From this picture we observe that MP11W and MP8W oscillate around the peculiar equilibrium conformations that will be discussed in the subsequent section, in particular relative to the central residues that are the main argument of the following discussion. For instance, during the 100 ps run MP11W shows regularly the characteristic helical arrangement of the peptidic backbone and orientation with respect to the heme plane. These plots show that the conformations have a stationary behavior during the runs, and the same pattern has been observed for MP11M and MP8M. The quantities reported in this paper were computed over 500 sets of data saved at a time interval of 0.2 ps of the total trajectory. A first 10 ps equilibration MD run was performed
19244 J. Phys. Chem., Vol. 100, No. 50, 1996
Melchionna et al. TABLE 1: Diagonalized Gyration Tensor (T), Gyration Radius (Rg), and End-to-End Distance (∆) Computed for the Atoms Belonging to the Peptidic Backbone of the Simulated Systems (MP11W, MP11 in Pure Water; MP11M, MP11 in Mixed Methanol+Water Solvent; MP8W, MP8 in Pure Water; MP8M, MP8 in Mixed Methanol+Water Solvent) and for the Crystallographic Data of the Corresponding Residues of Bonito Cytochrome ca cyt c MP11W MP11M MP8W MP8M
T11
T22
T33
Rg
∆
39.72 21.25 36.13 17.38 16.62
9.89 13.15 9.72 7.12 6.10
4.90 5.75 6.94 4.26 4.55
7.38 6.33 (0.24) 7.27 (0.48) 5.35 (0.14) 5.22 (0.09)
15.46 11.5 (1.3) 14.5 (1.4) 8.5 (0.4) 8.1 (0.5)
a The gyration tensor is defined: T˜ Rβ ) 〈Σi mi/Σjmj(riR - R0R)(riβ R0β)〉, where R0 is the position vector of the center of mass of the system and the greek letters refer to the three components x, y, z. The gyration radius is Rg ) [∑RTRR]1/2, and the end-to-end distance is ∆ ) 〈|r1 rN|〉. In parentheses for Rg and ∆ are reported the relative standard deviations.
TABLE 2: Root Mean Square Displacement of the Cr Atoms rmsd
Figure 4. Sequence of five different conformations of MP11W (left column) and MP8W (right column) taken at regular intervals during the equilibrium runs. The plots refer to the peptidic backbone and heme group.
on every simulated system before beginning to collect the data to be analyzed later. 3. Results and Discussion Microperoxidases are generally found in aggregates with high molecular mass (>10 kDa), while the monomers are experimentally obtained by working at high dilution and in the presence of strong ligands, such as imidazole, to coordinate the sixth position of the iron atom.6 Our results concern the analysis of the monomer structure, coordinated with imidazole in the sixth position of Fe (see Figures 1 and 2). Different aggregation patterns have been proposed for the aggregates, with either an irregular or regular alignment. In Figure 3 two possible aggregation patterns, the irregular headto-tail alignment (a) and the more regular stacked alignment (b), are schematically drawn. By altering the solvent conditions other possible patterns have been proposed.28 In the following we will only infer the possible regularity of the aggregates that can be built using monomer conformation. Monomer Characterization. The global structures of the peptidic chain of MP11W, MP11M, MP8W, and MP8M are monitored by calculating the gyration tensor (in its diagonalized
CR
MP11W
MP11M
MP8W
MP8M
11 12 13 14 15 16 17 18 19 20 21
1.38 1.31 1.25 1.01 0.89 1.18 1.37 1.09 1.17 1.62 1.83
1.90 1.90 1.94 1.64 1.77 1.83 1.92 1.82 2.00 2.03 2.45
1.71 1.77 1.80 1.89 1.79 1.60 1.70 1.83
1.98 1.88 1.74 1.67 1.89 2.00 1.88 2.29
form), the gyration radius, and the end-to-end distance, as defined and reported in Table 1. The data on the gyration radius show that a global modification is expected in the MP11 molecule switching from pure water to mixed solvent conditions rather than in the MP8 molecule, since in MP8W and MP8M the gyration tensor shows the same kind of anisotropy and size. In MP11M the size has increased considerably (the linear dimensions of MP11M and MP11W have a ratio of 1.15). For the undecapeptide the size fluctuations increase of a factor of 2 when going from pure to mixed solvent, despite the expectation that the less polar solvent stabilizes the molecule. This feature is instead observed in the MP8 molecule, which exhibits very small fluctuations in nonaqueous solvent. The end-to-end data exhibit a broader structure for MP11M than the MP11W. The fluctuations of the end-to-end distance of MP11 are globally 3 times larger than MP8, as indicated in parentheses in Table 1. The flexibility of the backbone chain, i.e. of the atoms belonging to the peptide backbone, is studied by computing the root mean square displacements of the single CR atoms from the average values, as reported in Table 2. The MP11M and MP8M systems show very similar behavior with the largest rms among the different simulation conditions. The smallest rms are found for the MP11W system, while MP8 has the same trend but smaller differences between pure solvent and mixed solvent conditions. The values of MP11W are close in amplitude to the experimental and simulation results obtained for proteins,9 while MP11M, MP8W, and MP8M have systematically larger values (about a factor of 2.2 with respect to typical X-ray data values). The peptide conformation is better understood by looking at the set of dihedral angles {φ,ψ} of torsion around the N-CR and CR-C backbone bonds, respectively (inset of Figure 1).
Molecular Dynamics of Microperoxidases
J. Phys. Chem., Vol. 100, No. 50, 1996 19245
TABLE 3: {O,ψ} Dihedral Angles and Relative Fluctuations for the Central Residues CR
cyt c
MP11W
MP11M
MP8W
MP8M
14 15 16 17 18
-104 -8 -58 -11 -35
-61 (11) -58 (11) -71 (13) -65 (12) -115 (14)
φ (deg) -65 (14) -156 (16) -5 (32) -60 (13) -75 (17)
13 (42) -100 (21) 64 (14) -121 (22)
-58 (11) -67 (13) 73 (9) -139 (15)
14 15 16 17 18
-63 -80 -3 -120 157
-34 (11) -34 (13) -52 (14) -42 (13) 119 (27)
ψ (deg) -110 (14) -60 (16) -74 (22) -44 (17) -46 (18)
-32 (39) -71 (22) 173 (14) -85 (21) 158 (11)
163 (12) -43 (15) 142 (12) -42 (15) -105 (20)
The φ angle is the dihedral angle formed by the two planes formed by the triplets of atoms of the backbone sequence (i) (i) (i) C(i-1) - N(i) - C(i) R and N - CR - C , respectively (the italic index refers to the aminoacidic sequence). The ψ angle is the dihedral angle formed by the two planes containing the triplets (i) (i) of atoms of the backbone sequence N(i) - C(i) R - C and CR (i) (i+1) C - N , respectively. The analysis of the angles {φ,ψ} is here reported only for the central residues, from Cys14 to His18, which give the dominant contributions to the whole structure. A direct calculation of the time evolution of the dihedral angles furnishes a good picture of the stability of the peptide by direct evidence of the dihedral torsions taking place along the backbone of the peptide. These calculations monitor better the stability of the peptide than for instance the rms displacements of the CR atoms. In fact, the relevant structural transitions of the peptide are due to torsional movements around the backbone bonds, so that the knowledge of the backbone torsional angles and relative fluctuations permits us to monitor directly the structural stability of the peptidic chain. In contrast, the rms displacements are superpositions of both the local thermal motion of the single atoms and the rigid body motions, and they do not represent effectively the peptidic chain stability. In Table 3 are reported the observed values of {φ,ψ} and the relative fluctuations. The values show that the four systems adopt in general different conformations, and the relative structures need to be analyzed separately. The structures appear to be quite rigid and characterized by a quite well-defined set of dihedral angles in all the solvent conditions, which is not observed in our previous calculations of MP1114 and MP822 in vacuum conditions (the fluctuations increase systematically up to a factor of 2). The original conformation of the cyt c corresponding residues seems to be definitively lost in all the systems. As we will see later on by looking at the configurational distance, in reality only MP11W maintains some memory of the original cyt c structure. The MP11W system shows the most stable behavior, as observed from the analysis of the rms displacements. This system keeps a helical arrangement, as first observed by Ehrenberg and Theorell in 1955,1 confirmed in a subsequent NMR study by Kimura et al.23 and predicted by Jehanly et al. using the method of Chou and Fasman,24 since the values correspond, within the fluctuations, to the ideal values of a righthanded R-helix (φ ) -57°, ψ ) -47°). This arrangement is lost in the MP11M, due to a dihedral transition at the level of GLN16. In comparison to MP11W, MP8W is much more mobile and adopts a random coil conformation. Moreover, we observe that the octapeptide in mixed solvent has a conformation similar to that of MP8W at the level of residues Ala15 to Cys17. The fluctuations reveals a slightly more stable structure for the octapeptide in mixed solvent than in pure water.
The structure and stability of MP seem to be effectively analyzed by computing the configurational difference between the simulated structures and the reference structure of cyt c. To this aim we use the definition of the square configurational distance D(ν) Γ ,
DΓ(ν)(t) )
1
N(ν)
B(R,β,γ)ri(ν) - r°i(ν)|2] min [ ∑ |R
N(ν) (R,β,γ)
(1)
i)1
where {r(ν) i } is the set of coordinates of the subunit ν of the molecule under study relative to its center-of-mass position, and {r°i(ν)} are the corresponding reference crystallographic coordinates of cyt c.25 The minimum is calculated over the three Euler angles (R,β,γ) that by a global rotation around the centerof-mass position, obtained by the application of the orthogonal matrix B R(R,β,γ), optimizes the superposition between the two structures and then eliminates spurious contributions. The calculation of D(ν) Γ is very efficient by use of the formalism of quaternions, as described by Kneller,26 so that the minimization over rotations may be solved as an eigenvalue problem of a (ν) four-dimensional matrix, as a function of {r(ν) i } and {r° i }. The minimum eigenvalue of this matrix corresponds exactly to the wanted distance. The configurational distance is a direct measure of how the new structures have gone far away from the original structure of cyt c. Since it is a sum of the configurational distance atom by atom, it is not a local measure of the peptidic arrangement, as for the set of {φ,ψ} angles or for the rms displacements of single atoms, but it is a good parameter for comparing effectively structures of strongly inhomogeneous systems, as in the case of proteins. The time evolutions of D(ν) Γ referred to the CR atoms of the backbone of our systems are plotted in Figure 5. Figure 5 shows that in general the average value of D(ν) Γ of MP11W is closer to the original structure of cyt c than the others. Its value in time is very close to the evolution of D(ν) Γ for MP11M. In contrast, the structures of both MP8W and MP8M are altered, resulting in a net difference in the conformations from the X-ray data. From the plot it is also evident that the fluctuations of D(ν) Γ , which may be thought again as a measure of the stability of the molecule, are comparable. Moreover, in Table 4 are reported the average components of the configurational distance of the single CR atoms. These values represent the average distances of the atoms of the dynamical structure from the corresponding atoms of cyt c, once the sum of eq 1 has been minimized with respect to rototranslations. The single atoms’ components have a local character, permitting us to distinguish which part of the peptide has moved farther from the original structure. In Figure 6a,b,c,d are shown the most recurrent configurations of the studied systems. The most recurrent structures are those structures chosen because they have the closest configuration to the average configuration computed over the whole trajectory. This means that given the set of trajectories as the output of the MD run, we first compute the average positions of all the atoms and then we compute which frame has the minimum configurational distance from the average structure itself; we call this set of coordinates the most recurrent structure. The ensemble-average structure is not significant because, when rigid bonds are used to simulate the molecule, the averaged coordinates over the trajectory result in a “strange” conformation that violates the constraints on the bond lengths. In fact, the average distance between a pair of atoms is not in general equal to the distance between the average coordinates of the pair. Figures 6a-d represents schematically the backbones of the studied systems together with the His18 residues and the heme group.
19246 J. Phys. Chem., Vol. 100, No. 50, 1996
Melchionna et al.
Figure 5. Square configurational distance (Å2) as a function of time (plotted every 0.2 ps of the total MD trajectory).
TABLE 4: Single Average Components of the Configurational Distance Referred Only to the Cr Atoms (Distance between the Position of the Selected Atom and the Position of the Corresponding Atom of cyt c Once the Structures of the Whole Cr Sequence Are Optimally Superposed; That Is the Sum of Eq 1 Has Been Minimized with Respect to Rototranslations; See Text for Details) CR
MP11W
MP11M
MP8W
MP8M
11 12 13 14 15 16 17 18 19 20 21
3.42 2.85 2.65 1.02 1.38 1.23 1.32 0.65 0.99 0.95 4.85
5.50 2.30 4.88 3.39 1.56 3.04 3.08 1.91 2.89 3.76 2.46
3.81 2.05 0.91 1.49 2.09 1.75 3.66 3.84
4.94 3.50 1.90 1.30 1.88 1.21 1.87 5.76
As discussed above, we observe that only the MP11W maintains a helical-like conformation. The smaller MP8W shows a similar conformation that is particularly distorted in the terminal residues. The systems in the mixed solvent exhibit a quite dramatical change, in particular for the undecapeptide. In general the octapeptide keeps a closed and packed conformation, with a tendency to form a loop (its conformation resembles a doughnut), while the undecapeptide shows a more linear elongated conformation. From these snapshots the reader may have an idea also of the orientation of the heme group with respect to the backbone. In accordance with Kimura et al.23 we find that the peptide chain folds up in a way to cover only one side of the heme plane. Aggregation Patterns. The most recurrent structures that we have obtained in pure water conditions are the starting point to interpret the data on the aggregates, which are generally thought to be arranged as linear polymers. Experimental studies have shown that in the aggregates the Fe atom is preferentially coordinated symmetrically with six nitrogen atoms in octahedral geometry.17 Therefore in the aggregate the possible coordination
of the Fe atom in the sixth available position is with one of the available aminic groups of the peptidic chain. Using simple geometrical arguments, we have assembled together two or more monomeric units of MP (assembling together only peptides of the same size), coordinating the Fe atom with the aminic groups of the peptidic chain that are available in the two peptides separately. In the undecapeptide, where Fe may coordinate both with R-NH2 of Val11 and -NH2 of Lys13 (here we use the standard nomenclature of increasing greek indexes when moving from the backbone atoms to the side chain atoms), the rodlike structure of the backbone and the angle formed by the backbone average axis and the heme plane normal do not favor a regular piled polymeric structure. Recently by use of energy minimization and MD simulation of differently constructed dimers of MP11 in water conditions, the same pattern has been observed for the dimer.27 The results indicated that the heme planes of two different microperoxidases are oriented with an angle between the two heme normals of 130° with a Fe-Fe distance of 14 Å. By using these values, our polymeric assembly results in an irregular head-to-tail conformation, as first postulated by Urry.7 For instance, a possible geometry obtained in this way is the octamer structure that has been postulated in ref 28. In MP8 the only donor available for coordination with Fe is the R-amino group of the terminal cysteine. In this case the symmetry of the backbone is no longer rodlike: the backbone resembles more a torus that lies on a plane almost parallel to the heme plane (Figure 6c,d). This structure favors a regular stacked conformation. In this case two subsequent heme planes have a smaller distance than for MP11 (
