Molecular Dynamics of Phycocyanobilin Binding Bacteriophytochromes

Dec 2, 2010 - Technische UniVersität Berlin, Institut für Chemie, Max-Volmer-Laboratorium, Sekr. PC 14, Strasse des 17. Juni 135, D-10623 Berlin, Ge...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 16677–16686

16677

Molecular Dynamics of Phycocyanobilin Binding Bacteriophytochromes: A Detailed Study of Structural and Dynamic Properties Steve Kaminski and Maria Andrea Mroginski* Technische UniVersita¨t Berlin, Institut fu¨r Chemie, Max-Volmer-Laboratorium, Sekr. PC 14, Strasse des 17. Juni 135, D-10623 Berlin, Germany ReceiVed: May 28, 2010; ReVised Manuscript ReceiVed: NoVember 8, 2010

The structural stability and conformational flexibility of two phycocyanobilin (PCB) binding bacteriaphytochromes, namely, Cph1 from Synechocystis and the GAF domain of SyB from Synechococcus, were studied using all-atoms molecular dynamics simulations techniques. In order to involve the tetrapyrrole cofactor in the simulation, new empirical force field parameters were developed for PCB which are compatible with the CHARMM22 force field for proteins. Special emphasis was made in understanding the conflicting NMRbased structures recently obtained for the two parent states, Pr and Pfr, of SyB(GAF) regarding the highly distorted cofactor conformation which is in contrast to all crystallographic measurements on other phytochrome species. The 25 ns all-atoms MD simulation of Cph1 in the Pr state and SyB(GAF) in its two parent states supports the picture of a relatively planar conformation of rings A, B, and C and ring D slightly twisted out of the ABC plane, in good agreement with crystallographic and spectroscopic experiments. SyB(GAF) converges to two very similar structures for the Pr state and also for the Pfr state with a ZZZssa conformation in disagreement with the ZZE configuration observed experimentally for other phytochromes. The failure in the prediction of the SyB(GAF) Pfr geometry is a consequence of deformed initial structure. Furthermore, in contrast to the results obtained for the SyB(GAF), the MD simulations showed a very stable Cph1 structure where an important hydrogen bond and water network in the chromophore binding pocket could be identified. We did not find evidence for structural heterogeneity either for Cph1 or SyB(GAF) at least on the nanosecond time scale. 1. Introduction Phytochromes constitute a family of red light photosensory receptors which respond to external light conditions. At first discovered in higher plants,1 they were later also found in bacteria, algae, and fungi.2,3 Phytochromes control various physiological processes like seed germination, flowering time, shade avoidance and the formation of leaves in higher plants, as well as phototaxis and pigmentation in bacteria. Through specific interactions between the polypeptide and an open-chain tetrapyrrole cofactor, phytochromes reversibly interconvert upon light absorption between two stable states. These are the physiologically inactive red-light absorbing Pr state and the activated far-red-light absorbing Pfr state.4 To understand the microscopic processes which occur during the photocycle, numerous experimental studies have been carried out over the last decades for various members of this protein family from different species. However, most of the molecular events which are involved in the photoreaction are still a matter of debate. A breakthrough in phytochrome research was the resolution of the three-dimensional structures of several phytochrome fragments from X-ray crystallography5-7 and NMR spectroscopy,8,9 which revealed a detailed insight into the structure at an atomic level. However, among these structural models, the recent NMRbased structures of the chromophore binding GAF domain8,9 of Synechococcus OSB’ bacterialphytochrome SyB(GAF) in its Pr and Pfr states do not agree with the crystallographic models obtained from other bacteriophytochromes and are also in conflict with spectroscopic data, including those of SyB(GAF) itself. The NMR-based structures suggest highly distorted minimum energy conformations of the phycocyanobilin (PCB)

chromophore and a very high structural flexibility in both parent states, in contrast to available crystal structures and results from molecular simulations of similar PΦB binding phytochromes.10 In addition, resonance Raman (RR) spectra point to far-reaching structural similarities of the chromophore structures in SyB(GAF) as compared to other PCB-binding phytochromes.11 Instead of the ZZZ/ssa conformation/configuration (see ref 12 for details of the Z/E s/a nomenclature of open chain tetrapyrroles) of the methine bridges, shown for phytochrome cofactors in the Pr state from crystal structures and supported by MD simulations,10 the NMR structure model displays an approximately perpendicular orientation of ring A with respect to the remainder of the tetrapyrrole. This distortion is removed in the Pfr model such that the authors suggest an isomerization around the A-B methine bridge as the photoinduced step in the Pr-to-Pfr transition. According to these NMR structures, the C-D entity remains largely unchanged which is in contrast to the crystallographic data and RR spectroscopic results that indicate a Z-to-E isomerization of the C-D methine bridge.7,13 In addition, neither the crystallographic nor the NMR structural models can provide a key for understanding the heterogeneity at the cofactor binding site inferred from various spectroscopic studies on biliverdin- and PCB-binding phytochromes.14-16 Specifically, it has been suggested that the origin of the heterogeneity lies in conformational differences of the A-B methine bridge.10 Thus, the aim of this work is to address this conflict between different experimental studies from a computational point of view, by analyzing the structural stability and conformational flexibility of phycocyanobilin (PCB) cofactors inside the protein

10.1021/jp104903u  2010 American Chemical Society Published on Web 12/02/2010

16678

J. Phys. Chem. B, Vol. 114, No. 50, 2010

binding pocket of phytochromes via molecular dynamics (MD) simulations. Classical MD simulations are in particular a well established tool which can, with a certain degree of accuracy and at a relatively low computational cost, reproduce conformational states of macromolecules. In a former work,10 we performed MD simulations using a suitable set of empirical force field parameters for the bilin cofactor to identify a possible origin for structural heterogeneity of the phytochromobilin (PΦB) chromophore in DrCBD (chromophore binding domain of Deinococcus radioduras) observed experimentally. Up to now, only a few computational studies16-20 were done on phytochrome systems by other research groups. To our knowledge, none of them dealt in detail with structural and dynamic properties of the cofactor molecule in the protein binding pocket on the nanosecond time scale. In the present work, we extended the computational studies of phytochrome photoreceptors to two PCB binding phytochromes, by performing MD simulations of the entire photosensory core (PAS-GAF-PHY) of Cph1 from Synechocystis in its Pr state6 and of the chromophore binding domain (GAF) in the Pr and Pfr states of SyB bacterialphytochrome from Synechococcus.8,9 In order to address the open questions concerning structural and dynamic properties at the PCB binding site and to include the PCB cofactor in the MD simulations, a set of molecular mechanics force field parameters were developed for PCB. 2. Computational Methods The Potential Energy Function. MD simulations were performed using the all-atoms CHARMM21 potential energy function. The overall energy U is given as a sum of pair interactions via various bonded and nonbonded terms:

U(rj) )

∑

Kb(b - b0)2 +

bond

UB

∑

Kθ(θ - θ0) +

∑

Kφ(φ - φ0)2 +

angle

improper

∑ KUB(S - S0)2 +

2

∑

Kχ(1 + cos(nχ - δ)) +

dihedral

∑

nonbond

[ [( ) ( ) ] ] εij

Rij rij

12

Rij 6 + rij qiqj (1) 4πεrij

-2

Here, K denotes the force constants for bonds (b), Urey-Bradley 1,3-distances (S), bond angles (θ), dihedral angles (χ), and improper dihedral angles (φ), whereas b0, S0, θ0, χ0, and φ0 are the respective values of the individual terms at equilibrium. Electrostatic and dispersion interactions (nonbonded) are described as a sum of Lennard-Jones and Coulomb terms. The first depends on the parameters εij and Rij, for the atom pair i and j and the distance rij between them. The Coulomb interactions are calculated between pairs of atoms at a relative distance rij, containing partial charges qi and qj. Molecular Dynamics Simulations. Protein structural models derived from X-ray and NMR measurements were used as starting geometries for simulations of the macromolecules. For Cph1, the initial set of coordinates for the protein was obtained from the Protein Data Bank (code, 2VEA; resolution, 2.45 Å). Model building was carried out as described previously.22 In all MD simulations, the cofactor was covalently attached to the protein via a thioether bridge. In that study on Cph1, also the investigation of different protonation states of titratable histidine residues in the chromophore binding pocket

Kaminski and Mroginski were included. For the comparison with the structures obtained from NMR solution measurements, where HIS139 and HIS169 are protonated at the δ-nitrogen, we decided to perform simulations of Cph1 using the same protonation state for the respective histidines, namely, HSD260 and HSD290. With respect to SyB(GAF) phytochrome, three-dimensional low energy structures of both parent states, Pr and Pfr, derived from NMR spectroscopy were kindly provided by Andrew Ulijasz.9 The protein structures were solvated in a 90 × 90 × 90 Å3 cubic box of TIP3P water.23 For both protein systems, SyB(GAF) and Cph1, several equivalent residues especially in the binding pocket exist which will be subject to further discussions. For the sake of clarity, a sequence alignment, shown in the Supporting Information (Figure S5), has been performed using clustalX.24 All tables and figures included in the Supporting Information are denoted with a prefix S. In order to neutralize the net charge of the system and to ensure protein stability throughout the MD simulations,25 sodium and chloride ions were added to both systems by initially placing a counterion in the vicinity (3 Å distance) of each charged residue not involved in salt bridges. The remaining net charge was removed by randomly placing counterions in the water box. Computational details on molecular modeling of Cph1 and SyB(GAF) including minimization, heating, and equilibration were identical to those of former simulations on other phytochromes10,22 and are, therefore, described in more detail in the Supporting Information. For each of the three phytochromes investigated in this work, Cph1\Pr, SyB(GAF)\Pr, and SyB(GAF)\Pfr, production runs of 25 ns were carried out. MD simulations were performed with the NAMD26 program (version 2.6) in combination with the CHARMM2221 force field. All molecular structures were drawn using the VMD 1.8.6 visualization software.27 Development and Validation of Force Field Parameters for PCB. Classical MD simulations require a set of force field parameters, suitable for the molecular system of interest. While these parameters are available for amino acids from a variety of different force fields, this is usually not the case for protein bound cofactor molecules. To make the PCB cofactor participate in the dynamics of Cph1 and SyB(GAF), a set of new force field parameters for PCB were developed. We followed Mackerell’s21 procedure in order to ensure consistency with the standard CHARMM22 force field. A detailed discussion on the parametrization of tetrapyrroles has already been made in a related work10 for phytochromobilin, and will be only briefly reviewed here. Details of the performance of the derived new parameters for PCB are given as Supporting Information. In general, the empirical force field parameters for bonded and nonbonded interactions were optimized to reproduce molecular properties, mainly obtained from calculations at higher levels (ab initio) of theory. This so-called target data includes, e.g., structural information from optimized geometries, rotational energy profiles, and intermolecular interaction energies of model compounds with water molecules in close contact. For the parametrization procedure, model compounds of tetrapyrroles were taken as isolated molecules for simplicity and computational feasibility. Optimization of internal force field parameters for bonds, angles, and improper torsions were performed on the complete PCB structure in order to account for electron delocalization effects. The corresponding target data were obtained from density functional (DFT) calculations at the B3LYP/6-31G* level of theory. Average deviations for various internal coor-

Molecular Dynamics of PCB Binding Bacteriophytochromes

J. Phys. Chem. B, Vol. 114, No. 50, 2010 16679

Figure 1. Calculated rmsf values by averaging over all non-hydrogen atoms per amino acid residue of Cph1, derived from experimental B-factors (black) and calculated from the MD simulation (blue). Colored regions denote secondary structure motives: R helix (in dark gray) and β sheet (in light gray).

dinates of optimized geometries at the QM and MM level are given in Table S3 of the Supporting Information. To reduce the computational effort for the development of dihedral parameters, however, rotational energy profiles were calculated using smaller molecular fragments of PCB. Such profiles, as obtained from optimized empirical torsion parameters (see Table S4, Supporting Information), for the fragments are illustrated in Figure S3 of the Supporting Information. The optimization of partial atomic charges of PCB was performed considering 12 chromophore-water interactions, for which single point calculations at fixed geometries were sequentially (one water molecule + cofactor) carried out, as shown in Figure S1 of the Supporting Information. A set of empirical charges (Table S2, Supporting Information) were obtained by reproducing minimum interaction distances and energies, as well as the overall dipole moment of PCB (Table S1, Supporting Information), as calculated at the HF/6-31G* level of theory. Because the dipole moment of a charged molecule (+1e for PCB) depends on its position with respect to the origin of the coordinate system, in order to compare the ab initio and molecular mechanics calculations, the relative orientation of the water-cofactor complexes to the origin were identical for both of them. Furthermore, since the protein-cofactor carries a net charge (-1), the interaction energies were not scaled, in consistency with the procedure derived by MacKerell.21 The final set of force field parameters for PCB was validated by comparing energy minimized structures of the isolated cofactor molecule and the small fragments obtained at a quantum mechanical (QM) and molecular mechanical (MM) level (see Figure S2, Supporting Information). For this purpose, the cofactor conformation/configuration of the Cph1 X-ray structure (ZZZssa) was used as the initial geometry with the acidic groups at the propionate chains removed, to avoid the formation of hydrogen bonds between the chains and the backbone of PCB during in Vacuo minimizations. Root-mean-square deviations (rmsd’s) between the MM- and QM-minimized structures of PCB (0.27 Å) and the corresponding small fragments (less than 0.2 Å) are shown in Figure S2 of the Supporting Information. An overall good performance could be achieved with the set of

optimized empirical parameters employed. This was also confirmed by a statistical analysis of internal coordinates for the entire minimized PCB molecule. In addition, the results presented in Table S3 of the Supporting Information indicate that MM and QM optimized structures are in satisfactory agreement. 3. Results and Discussion Structural and Dynamic Properties of the Polypeptides. For both protein systems, Cph1 and SyB(GAF), the statistical analysis of the MD simulations was performed using the last 20 ns of the calculated trajectories. The average rmsd of the polypeptide of Cph1 yielded ∼3.2 ( 0.3 Å, which is very similar to the one for SyB(GAF)/Pr with ∼3.2 ( 0.2 Å. For the Pfr state of SyB(GAF), however, a significantly larger average value of ∼4.2 ( 0.4 Å was obtained, indicating that during the simulation SyB(GAF)/Pfr deviates more strongly from its initial structure compared to the others. The radius of gyration of the polypeptide of Cph1 was calculated to be ∼30.0 ( 0.3 Å. The respective values for SyB(GAF) are ∼16.4 ( 0.1 Å for Pr and ∼16.4 ( 0.2 Å for Pfr. The small fluctuations in each case indicate that no significant changes in the overall shape of the protein occurred during the production run for any of the three systems. Furthermore, the deviations between the average gyration radii from the simulations and the ones obtained from the initial experimental structures are also relatively small with ∼0.5, ∼0.4, and ∼0.9 Å for (SyB(GAF)/Pr, SyB(GAF)/Pfr, and Cph1, respectively. To monitor the flexibility of the polypeptide, root-meansquare fluctuations (rmsf’s) per residue were computed from the MD trajectories of the three systems and plotted over the entire range of residues in Figures 1 and 2 for Cph1 and SyB(GAF), respectively. For Cph1, temperature B-factors, as obtained from X-ray analysis, were used to calculate the experimental rmsf (black curve in Figure 1) via

〈∆ri2〉 )

3 Bi 8π2

16680

J. Phys. Chem. B, Vol. 114, No. 50, 2010

Kaminski and Mroginski

Figure 2. Calculated rmsf values by averaging over all non-hydrogen atoms per amino acid residue of the Pr (black) and Pfr (blue) states of SyB(GAF). Colored regions denote secondary structure motives: R helix (in dark gray) and β sheet (in light gray).

Figure 3. Superimposition of NMR (dark gray) and MD average structures (white) for the Pr (left) and Pfr (right) states of the PCB cofactor.

for a qualitative comparison. The comparison between experimental and calculated rmsf values, however, should only be done qualitatively by analyzing the relative flexibility of different regions in the polypeptide chain. For the first 250 residues of Cph1 (top graph of Figure 1), which are located in the PAS and in the GAF domains, the protein flexibility computed out of the MD simulations is in overall good agreement with that estimated from the experiments. In both cases, the largest fluctuations occur in the region near residue 80. For the PHY domain of the protein (bottom graph of Figure 1), more evenly distributed rmsf values are obtained by the MD simulations, whereas those derived from the experiments reveal a significantly larger flexibility. In the large R helix connecting the GAF and the PHY domain of the protein (residues ∼300-345), the experimental rmsf values are small in the middle of the helix and increase significantly to the edges. The fluctuations obtained from the simulations are in this region rather small, indicating an overall more rigid helical structure at room temperature. In the case of SyB(GAF), a comparison between the two curves in Figure 2 representing the rmsf for the Pr (black curve) and Pfr (blue curve) states reveals very similar flexibility of the residues over almost the entire sequence. An exception, however, is the region between residues 70-83 and 108-112, where deviations between Pr and Pfr of more than 1 Å are noted. The strongest protein movements occur in the region between residues 110 and 130. Here, the fluctuations are more pronounced for the Pfr state. The respective residues belong to a large unstructured loop region that in SyB(GAF-PAS) interacts with the PAS domain. This result is in qualitatively good

agreement with the experimental NMR results, for which a superimposition of the 20 lowest energy structures yielded the largest variability in the same region, for both parent states. Overall, the magnitude of the protein flexibility for both parent states of SyB(GAF) is comparable to that obtained from the MD simulation of Cph1. Structural and Dynamic Properties of the PCB Cofactor and Its Binding Pocket. In the case of Cph1, the calculated rmsd of the average PCB-cofactor structure obtained from the MD simulations with respect to the corresponding experimental structure is quite low (rmsd ∼0.8 Å), indicating that both structures differ little from each other. This is especially true for rings A and B, while ring C is almost coplanar to rings A and B in the average structure from the MD simulation but displaced more out of plane in the crystal structure. Large deviations are observed at the methine bridge connecting rings C and D. This methine bridge exhibits a rather planar geometry in the X-ray structure (single bond torsion angle, 19.0°; double bond torsion angle, 2.3°), whereas the MD simulation average reveals a much more distorted conformation (single bond torsion angle, 45.4°; double bond torsion angle, 22.3°). Large conformational differences are also observed for both propionate chains of PCB. An observation related to this phenomenon is shown in the right picture of Figure 4. Here, residue ARG254 which forms a salt-bridge to the propionate chain of ring B in the crystal structure is replaced by a neighboring arginine residue (ARG222) during the MD simulation. In contrast to Cph1 phytochrome, the PCB structures in both parent states of SyB(GAF) deviate significantly from the initial NMR structures, as indicated by the rmsd of ca. 2.1 Å for the

Molecular Dynamics of PCB Binding Bacteriophytochromes

J. Phys. Chem. B, Vol. 114, No. 50, 2010 16681

Figure 4. Left: Average structure of the chromophore binding pocket calculated from the MD simulations (white) and superimposed with the crystallographic structure (dark gray) of Cph1. Right: Comparison of PCB-protein interactions in the binding pocket. ARG254 is exchanged by ARG222 for interactions with a propionate chain of PCB during the simulation (white/average structure from the MD simulation) compared to the initial X-ray structure (dark gray).

Figure 5. Visualization of calculated rmsf values for non-hydrogen atoms (in Å) of the PCB cofactor in Cph1 (left) and in both parent states Pr (middle) and Pfr (right) of SyB(GAF).

Pr state and of ca.1.8 Å for the Pfr state, as illustrated in Figure 3. Here, the lowest energy NMR structures of PCB (dark gray) are superimposed with the corresponding average structures obtained from the MD simulations (white) in the Pr (left) and Pfr (right) states. The MD simulations indicate an almost coplanar structure of ring A with respect to the inner rings B-C in the Pr state, while the lowest energy NMR structure shows ring A in an almost 90° out-of-plane conformation. Also, the highly distorted methine bridge between the inner rings B-C as derived from the experiments is contrasted by a more coplanar structure in the MD simulations. Ring D evolves rapidly from an almost perpendicular orientation with respect to the inner rings in the starting (NMR) geometry to an anti configuration at the CD methine bridge in the simulations of both Pr and Pfr states. Interestingly, as a result of the MD simulations, the PCB cofactor in the two parent states converged to very similar structures with a more planar ZZZssa-like conformation/ configuration (see structures in white color in Figure 3). This structural behavior of PCB in SyB(GAF) Pr/Pfr states has been verified by several shorter (three, each 2.5 ns) independent MD trajectories. While the observations for Pr are in agreement with related crystallographic structures of Cph1 and DrCBD, the results for the Pfr state (ZZZssa) are, however, in conflict with structural models for the Pfr state of biliverdin-binding phytochromes.7 The conformational fluctuations of the cofactor in the three systems throughout the simulations are plotted in Figure 5. Here, the calculated rmsf per atom, as illustrated by a color code, allow for a spatially resolved analysis of the flexibility of the PCB cofactor. The absolute values, in Å, are given within color bars for each system. For Cph1, rings B, C, and D correspond to the more rigid part of the cofactor, whereas the propionate chains and ring A are subject to larger fluctuations.

The most striking observation for the SyB(GAF) species is the difference in absolute rmsf values between the two parent states. Figure 5 shows that the PCB cofactor in the Pr state is subject to substantially larger fluctuations compared to Pfr. With rmsf values up to 2.5 Å and more, the cofactor in the Pr state is significantly more flexible compared to most of the surrounding residues (with exception of ASP86) in the binding pocket. In addition, the relative flexibility of PCB in the two states is different. The highest mobility is noted for rings C and D in Pr, while rings A and B are most flexible in the Pfr state. Interestingly, comparison with the rmsf of PCB in the Pr state of Cph1 reveals similarities to the SyB(GAF)/Pfr rather than to its Pr state. This is true not only for the absolute rmsf values but also for the relative molecular flexibilities. For the PΦB chromophore in DrCBD, the alternating variation of the distance between ring A and His260 during the MD simulation may be related to the existence of structural heterogeneity of the cofactor.10 This alternating distance behavior was not observed for Cph1 or for SyB(GAF). In fact, His260 stays in both proteins outside the hydrogen bond distance from ring A of PCB during the entire simulation time. Also, the single bond dihedral angle in the methine bridge between rings A-B only fluctuates around a single value, in contrast to PΦB in DrCDB, where two states could be identified. To summarize, MD simulations on Cph1 and SyB(GAF) did not give, at least on a nanosecond time scale, evidence for structural heterogeneity of the tetrapyrrole cofactor. Hydrogen Network between the Chromophore and the Surrounding Residues. The most important hydrogen bond interactions involving the PCB chromophore and the protein matrix are defined by a donor-acceptor distance shorter than 3.4 Å in the X-ray structure or during the MD simulation.28 In this way, for Cph1, a total of eight hydrogen bonds could be clearly identified as illustrated by dotted lines in Figure 6, where

16682

J. Phys. Chem. B, Vol. 114, No. 50, 2010

Kaminski and Mroginski TABLE 2: Statistics of Hydrogen Bond Events between the Chromophore and Surrounding Amino Acid Residues with Occupancies Larger than 20%a average occupancy atom 1 (protein) atom 2 (PCB) lifetime (ps) events (%)

Figure 6. Hydrogen bonding network between the PCB chromophore and surrounding residues in the protein binding pocket. The structure was obtained as an average from the MD simulation of Cph1. Related statistics are given in Tables 1 and 2.

TABLE 1: Distances (in Å) between Hydrogen Bond Donors and Acceptors in the Chromophore Binding Pocket of Cph1, Participating in the Interactions Illustrated in Figure 6 H-bond atom 1 (protein) atom 2 (PCB) crystal 1

N (Y257)

2

NH2 (R222)

3 4 5 6 7 8 X

NE (R222) OG1 (T274) NE2 (H260) O (D207) O (D207) NE2 (H290) ND1 (H260)

O1D (prop. chain ring O1D (prop. chain ring O2D (prop. chain ring O1A (prop. chain ring O1A (prop. chain ring N (ring B) N (ring B) O (ring D) N (ring A)

MDb

Diff.

3.06

3.65 (0.64) 0.59

2.67a

2.87 (0.30) 0.20

3.33a

2.98 (0.42) 0.35

2.74

3.96 (0.46) 1.22

3.31

3.02 (0.36) 0.29

2.72 3.05 2.80 5.36

3.16 3.27 3.95 5.17

B) B) B) C) C) (0.34) (0.25) (0.44) (0.33)

0.44 0.22 1.15 0.19

a ARG254 instead of ARG222. b Average distances and fluctuations (in parentheses) are calculated over the last 20 ns of the MD simulation.

an average structure of the binding pocket from the MD simulation is shown. Tables 1 and 2 show the corresponding statistics of hydrogen bond events evaluated over the last 20 ns of the MD trajectory. In Table 1, average donor-acceptor distances and fluctuations (in parentheses) from the MD simulation are compared to the respective values from the crystal structure. The interactions of PCB with ARG222 (number 2 and 3) during the simulation are compared with those between PCB and ARG254 from the X-ray structure, because of the positional exchange shown in Figure 4 (right). For most of the interactions (1, 2, 3, 5, and 7) that were monitored, the distances measured in the crystal structure are within the fluctuations obtained from the MD simulations and are thus in good agreement. The largest deviations (>1 Å) are noted for interaction numbers 4 and 8, where both THR274 and HIS290 depart from the chromophore site during the simulation. SER272 forms a very weak hydrogen bond with the propionic side chain of ring C, which seems to be further weakened and even disrupted during the simulation. The corresponding donor-acceptor distance increases from 3.71 to 4.47 Å. A similar effect has been previously predicted for DrCBD.10 In Table 2, the stability of hydrogen bonds between residues in the binding pocket and PCB is presented in terms of their

O/ASP 207 H/ring A O/ASP 207 H/ring B HH21/ARG 222 O1D/prop. chain ring HE/ARG 222 O2D/prop. chain ring HH21/ARG 222 O1D/prop. chain ring HE/ARG 222 O2D/prop. chain ring HN/TYR 257 O1D/prop. chain ring HE2/HIS 260 O1A/prop. chain ring HE2/HIS 260 O1A/prop. chain ring

28.5 58.6 395.3

183 186 30

26 54 59

141.6

38

27

85.3

133

57

140.4

87

68

59.3

84

25

38.7

124

24

164.9

94

77

B B B B B C C

a Statistical parameters were evaluated over a 20 ns MD simulation of Cph1. The positions of heavy atoms related to PCB are explained in Figure S3 in the Supporting Information.

frequency (events) and duration (lifetime) throughout the simulation. The product between these values divided by the simulation time yields the occupancy, representing an estimate of the overall stability of a hydrogen bond. The largest occupancies for Cph1 could be observed for the interactions between HIS260 and the propionate chain at ring C of PCB on the one hand and those between ARG222 and the propionate chain at ring B of PCB on the other hand. The observation, that the propionate chains are involved in most of the interactions listed in Table 2, is in agreement with the results obtained for DrCBD, where the propionates were also tightly bound to surrounding residues via hydrogen bonds. While an overall structurally very stable protein binding pocket is estimated from the MD simulation of Cph1, both parent states of SyB(GAF) display a completely different behavior. In Figure 7A, the position of important residues in the chromophore binding pocket is presented before and after the MD simulation for the Pr (left) and Pfr (right) states. In each case, the structures were aligned with respect to the cofactor, to obtain the position of the residues relative to those of PCB. The movement of residues is illustrated by black arrows. For clarity, only the cofactor of the initial structure in both parent states is shown. Pronounced positional changes relative to the initial structure can be observed for TYR142 in Pfr and TYR54, HIS169 in Pr, as well as PHE82 in both states. Figure 7B shows a segment of the chromophore binding pocket in SyB(GAF), namely, the cofactor together with a short R helix (residues 137-146) which are connected to each other via a cysteine residue. The superimposition of the initial NMR structures (dark gray) and the average structures from the MD simulation (white) reveals, for both parent states, a change of the relative position of the cofactors inside the binding pocket, which is more pronounced for the Pr state than for the Pfr state. For the Pfr state, significant conformational differences, compared to the NMR structure, are in particular noted at the cysteine residue. As a result of the cofactor movements in the binding pocket, the only stable protein-cofactor interaction observed in the Pr state is the one between ARG101 and the propionate chain of PCB on ring B. No other residues stay in close contact with the chromophore during the simulation, whereas, for the Pfr

Molecular Dynamics of PCB Binding Bacteriophytochromes

J. Phys. Chem. B, Vol. 114, No. 50, 2010 16683

Figure 7. (A) Superimposition of NMR (dark gray) and average MD structures (white) for the Pr (left) and Pfr (right) states. The segment of the protein that is shown illustrates the movement of the cofactor in the binding pocket throughout the simulation, relative to the initial structure. (B) Positions of important residues in the binding pocket before (shiny color) and after (normal color) a 25 ns MD simulation for the Pr (left) and Pfr (right) states. Black lines denote the direction of motion during the simulation for each residue. (C) Superimposition of NMR (dark gray) and average MD structures (white) for the Pr (left) and Pfr (right) states. The segment of the protein that is shown illustrates the different interactions of ARG101 in the two parent states.

state, the only residue in close contact with the cofactor is HIS139, which forms a stable hydrogen bond with ring A of the chromophore. This observation is clearly in conflict with the NMR-based structural models, where HIS139 behaves in the opposite way by turning away from the chromophore in the Pfr state.9 The reduced number of chromophore-protein interactions which remain throughout the MD simulations of SyB(GAF) is in contrast to the results from MD simulations of Cph1 and DrCBD10 where a larger number (9-10) of amino acids interact with the chromophore through a stable hydrogen bonding network. Although not clearly visible in Figure 7A, the most striking positional change in SyB(GAF) refers to the ASP86 in both parent states, which moves far away from its initial position in between the inner pyrrole rings B and C. Furthermore, ASP86 is subject to a large rmsf (Pr, 2.3 Å; Pfr, 1.7 Å), indicating that

its position remains rather unstable during the simulation, in contrast to other residues in the binding pocket. This observation is in accordance neither with NMR solution structures of SyB(GAF) nor with results obtained from MD simulations of Cph1 (and DrCBD). To confirm that this behavior has statistical significance, three independent MD simulations of 2.5 ns each were performed for both parent states. The time-dependent distance plots between ASP86 (backbone oxygen) and ring B (pyrrole hydrogen) of the chromophore for Pr (left) and Pfr (right) resulting from these simulations are shown in Figure 8. In each of the six MD simulations, the ASP86 leaves its initial position, which happens especially in the Pfr state, rather suddenly within the first 400 ps of the simulation. In the Pr state, the distance between the inner ring B and the ASP86 increases from 3.3 Å up to 13 Å. For comparison, the red curve in the left picture of Figure 8 shows in addition the time-

16684

J. Phys. Chem. B, Vol. 114, No. 50, 2010

Kaminski and Mroginski

Figure 8. Time-dependent distance plot between ASP86 (backbone oxygen) and ring B (pyrrole hydrogen) of the chromophore in SyB(GAF) from three independent MD simulations in the Pr (left) and Pfr (right) states. The additional red curve in the left picture shows the distance behavior as observed from MD simulations of Cph1.

dependent ASP86-ring B distance behavior obtained for Cph1 in the Pr state. In this respect, it is notable that in the Pr state already the initial (lowest energy NMR structure) distances between ASP86 (backbone oxygen) and the three pyrrole rings A/B/C (pyrrole hydrogen) are beyond the hydrogen bond distance (ASP86-ring A, 3.30 Å; ASP86-ring B, 4.49 Å; ASP86-ring C, 3.77 Å) and therefore larger compared to the respective distances in the crystallographic structures of Cph1 and DrCBD. In the latter structures, such large movements of ASP86 during the MD simulations have not been observed. For SyB(GAF), the NMR structures show an interesting behavior for ARG101, which is a highly conserved residue in the chromophore binding pocket. As already known from several crystal structures (Cph1, DrCBD), this residue forms a salt bridge with the propionate chain of the cofactors ring B in the Pr state. The structure of the Pfr state from NMR measurements reveals a high mobility of this residue and the formation of a salt bridge to residue GLN185. Figure 7C illustrates the large movement of ARG101 and shows the distances it overcomes upon photoconversion (dotted lines) for interactions with the cofactor (2) and GLN185 (1) in the Pr (left picture) and Pfr (right picture) states, respectively. A superimposition of the NMR structures (dark gray) with the average structures (white) confirms the stability of ARG101 in both states throughout the MD simulation. In particular for the Pfr state, we observed that the salt bridge between ARG101 and GLN185 becomes stronger, as reflected by the corresponding donor-acceptor distance (1) which decreases from 3.56 Å (as revealed by the NMR structure) to 2.78 Å in the MD average structure. The distances of ARG101 to the propionate chains (2) in the Pr state differ by only 0.07 Å when comparing theory and experiment. Water Network in the PCB Binding Pocket. The structural stability of the chromophore binding pocket is also influenced by water-cofactor and water-protein interactions. From the MD simulation of Cph1, seven water molecules could be identified which form stable hydrogen bonds with the cofactor or with amino acid residues in the binding pocket. Figure 9 shows an average structure of the chromophore binding pocket built over the last 15 ns of the MD simulation containing the spatially conserved water molecules. Dotted lines indicate donor-acceptor distances within 3.4 Å (hydrogen bond distance). The water network predicted in the work for Cph1 is consistent with the water network observed in DrCBD.10 Not only are a similar number of water molecules involved (Cph1, 7; DrCBD, 8) but also the spatial arrangement of water molecules in the binding pocket is comparable. Like in DrCBD, the propionate chains of the cofactor play a dominant role in the network, because of the large number of contacts to water residues.

Figure 9. Stable water network in the chromophore binding pocket of Cph1. The structure refers to an average over the last 15 ns of the simulation. Dotted lines indicate donor-acceptor distances less or equal than 3.4 Å.

In both Cph1 and DrCBD, the position of one water molecule (Cph1, W6; DrCBD, W8) remains stable, forming a bridge between the propionate chain at ring C of the cofactor and ring D. The stability of this water molecule, as it is also shown in Tables 3 and 4, plays most probably an important role for the structural rigidity of the cofactor at ring D in the Pr state. The hydrogen bond occupancy factors for the interaction between W6 and rings C and D are 83 and 74%, respectively. The overall similarity of the water network is rather remarkable, since several residues in the binding pocket are not conserved between Cph1 and DrCBD. The largest deviation concerning the two proteins is related to the pyrrole water (W7 in Figure 9). The position of W7 is more stable within the inner rings B and C of the PCB cofactor in Cph1 than in DrCBD. This water remains in its position throughout the entire MD simulation, as indicated by the relatively large occupation factors estimated for the hydrogen bond interaction between W7 and rings A, B and C (65, 53, and 76%). In DrCBD, unlike Cph1, the chromophore cavity is only partially occupied by water.10 The movement of the pyrrole water out of the chromophore cavity has been considered as the origin for the structural heterogeneity of the chromophore in DrCBD. This effect has not been reproduced for Cph1, at least not on the nanosecond time scale.

Molecular Dynamics of PCB Binding Bacteriophytochromes TABLE 3: Statistics of Hydrogen Bond Events with Occupancies Larger than 20% between Important Water Molecules and Residues in the Binding Pocketa water W1 W1 W1 W2 W2 W3 W3 W4 W4 W5 W5 W5 W6 W6 W7 W7 W7 W7

O H H O H O H H H O H H O H H O O O

atom 2

average lifetime (ps)

events

occupancy (%)

HG1/SER 272 O2D/PCB O1A/O2A/PCB HG1/SER 274 O2A/PCB HG1/SER 272 O1A/PCB O2A/PCB O2D/PCB O1G/SER 274 HE2/SER 290 O1A/PCB H_B/PCB O1A/PCB ND1/HIS 260 H_C/PCB H_D/PCB H_A/PCB

157.7 155.0 260.4 323.3 129.6 152.8 192.5 1214.9 1735.9 51.4 115.5 54.7 487.1 39.6 148.1 86.5 51.9 201.3

43 36 54 31 76 94 88 17 7 336 125 220 34 374 111 150 203 76

34 28 71 51 49 72 85 100 61 86 72 60 83 74 82 65 53 76

a

Statistical parameters were evaluated from a 20 ns MD simulation of Cph1.

TABLE 4: Statistics of Hydrogen Bond Events with Occupancies Larger than 20% between Important Water Molecules Surrounding the Chromophorea water

water

average lifetime (ps)

events

occupancy (%)

W5 W2 W1

W6 W3 W3

31.5 32.9 188.9

339 121 53

53 20 50

a

Statistical parameters were evaluated from a 20 ns MD simulation of Cph1.

For SyB(GAF), a water network comparable to the one described for Cph1 has not been found, neither for the Pr nor for the Pfr state. The observed exchange of water molecules in the chromophore binding pocket was significantly higher compared to Cph1 (and DrCBD). The absence of amino acid residues in a stable close contact to the cofactor throughout the simulation is most probably the reason for this observation. 4. Conclusion In the present work, MD simulations of two PCB-binding members of the phytochrome family, Cph1 and SyB(GAF), were carried out using a new set of force field parameters derived for the PCB cofactor. Emphasis was laid on analyzing the structural stability and dynamical properties of the bilin chromophore and its immediate protein environment in order to understand the conflicting results derived from X-ray crystallography and NMR spectroscopy regarding the conformation and the flexibility of the protein bound PCB chromophore. The force field parameters developed for PCB were shown to reproduce various kinds of target data, calculated at higher levels of theory in Vacuo, in an overall satisfactory manner. Furthermore, the quality of these parameters for describing structural and dynamic properties of tetrapyrrole cofactors in a complex protein environment was tested by performing MD simulations of Cph1 phytochrome. We have shown that for this protein the structural properties (average structure) as well as its flexibility (rmsfs) are in good agreement with the experimental data. The rmsd between the average structure for PCB derived from the MD simulation and that from the crystal

J. Phys. Chem. B, Vol. 114, No. 50, 2010 16685 structure yield only 0.8 Å. The largest deviations are predicted at the CD methine bridge which is slightly more twisted in the theoretically predicted structure compared to the experimental one. In terms of structural heterogeneity, as observed from resonance Raman and NMR spectroscopic studies of various phytochromes, a possible origin has been identified in a previous work on the PΦB binding phytochrome in DrCBD. The interaction of the cofactor with a histidine residue in close contact leads to two well-defined PΦB conformations differing in the torsional angle around the C-C single bond at the AB methine bridge, formed on the nanosecond time scale during the simulations. A similar phenomenon has not been observed in the present work for PCB binding systems, Cph1 and SyB(GAF) Pr/Pfr. Hence, no evidence for heterogeneity in the PCB cofactor was found. However, for the Pfr state of SyB(GAF), HIS260 formed a very stable hydrogen bond to ring A of the chromophore similar to the one predicted for state 2 in DrCBD.10 From the MD simulations of Cph1 and DrCBD (former work), stable hydrogen bonds and water networks could be identified. Although the two systems differ in several aspects (different protein size, cofactor, and residue types in the binding pocket), the predicted water networks are very similar with respect to the number of water residues and their spatial arrangement in the binding pocket. In these networks, the propionate chains as well as residues in close contact to rings C and D of the cofactors are strongly involved, stabilizing a largely rigid structure in this region of the bilin. In addition, this network may also play an important role in the transient deprotonation/reprotonation during the photocycle. MD simulations of SyB(GAF) revealed structural and dynamic properties, especially in the binding pocket, which differ from those predicted for Cph1 and DrCBD. The highly distorted initial low energy PCB geometries derived from NMR data were found to be unstable within the simulations both in the Pr and in the Pfr state. Also, several adjacent residues were not structurally stable, especially the ASP86, as it was shown from several independent trajectories in both parent states of SyB(GAF). Thus, the present study does not support the structural models and the molecular mechanism of the photoinduced Prto-Pfr conversion recently suggested for SyB(GAF).9 Acknowledgment. Computational resources from the Norddeutscher Verbund fu¨r Hoch und Ho¨chstleistungsrechnen (HLRN) are acknowledged. We also thank the DFG Cluster of Excellence “UniCat” for the financial support as well as very helpful discussions with Prof. Hildebrandt and David von Stetten. Abbreviations PCB: phycocyanobilin GAF: cGMP phosphodiesterase/adenyl cyclase/FhlA SyB(GAF): GAF domain of cyanobacterial phytochrome from Synochococcus OSB’ Cph1: cyanobacterial phytochrome MM: molecular mechanics MD: molecular dynamics DRCBD: chromophore binding domain of Deinococcus radiodurans Supporting Information Available: Further details concerning the protocols for the MD simulations as well as the validation of force field parameters for phycocyanobilin (the parameters and topology files for phycocyanobilin are available

16686

J. Phys. Chem. B, Vol. 114, No. 50, 2010

from the authors upon request). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Butler, W. L.; Norris, K. H.; Siegelman, H. W.; Hendricks, S. B. Proc. Natl. Acad. Sci. U.S.A. 1959, 45, 1703–1708. (2) Quail, P. H. Nat. ReV. Mol. Cell Biol. 2002, 3, 85–93. (3) Rockwell, N. C.; Su, Y. S.; Lagarias, J. C. Annu. ReV. Plant Biol. 2006, 57, 837–858. (4) Briggs, W. R. Spudich, J. L. Handbook of Photosensory Receptors; Wiley Verlag: Weinheim, Germany, 2005. (5) Wagner, J. R.; Brunzelle, J. S.; Forest, K. T.; Vierstra, R. D. Nature 2005, 438, 325–331. (6) Essen, L. O.; Hughes, J.; Mailliet, J. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 14709–14714. (7) Yang, X.; Kuk, J.; Moffat, K. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 14715–14720. (8) Cornilescu, G.; Ulijasz, A. T.; Cornilescu, C. C.; Markley, J. L.; Vierstra, R. D. J. Mol. Biol. 2008, 383, 403–413. (9) Ulijasz, A. T.; Cornilescu, G.; Cornilescu, C. C.; Zhang, J. R.; Rivera, M.; Markley, J. L.; Vierstra, R. D. Nature 2010, 463, 250–256. (10) Kaminski, S.; Daminelli, G.; Mroginski, M. A. J. Phys. Chem. B 2009, 113, 945–958. (11) Ulijasz, A. T.; Cornilescu, G.; von Stetten, D.; Kaminski, S.; Mroginski, M. A.; Zhang, J. R.; Bhaya, D.; Hildebrandt, P.; Vierstra, R. D. J. Biol. Chem. 2008, 283, 21251–21266. (12) McNaught, A. D. Wilkinson, A. Compendium of Chemical Terminology; Blackwell Science: 1997. (13) Mroginski, M. A.; Murgida, D. H.; Hildebrandt, P. Acc. Chem. Res. 2007, 40, 258–266. (14) von Stetten, D.; Gunther, M.; Scheerer, P.; Murgida, D. H.; Mroginski, M. A.; Krauss, N.; Lamparter, T.; Zhang, J.; Anstrom, D. M.; Vierstra, R. D.; Forest, K. T.; Hildebrandt, P. Angew. Chem., Int. Ed. 2008, 47, 4753–4755.

Kaminski and Mroginski (15) Schmidt, P.; Gensch, T.; Remberg, A.; Ga¨rtner, W.; Braslavsky, S. E.; Schaffner, K. Photochem. Photobiol. 1998, 68, 754–761. (16) van Thor, J. J.; Mackeen, M.; Kuprov, I.; Dwek, R. A.; Wormald, M. R. Biophys. J. 2006, 91, 1811–1822. (17) Borg, O. A.; Durbeej, B. J. Phys. Chem. B 2007, 111, 11554– 11565. (18) Matute, R. A.; Contreras, R.; Gonzalez, L. J. Phys. Chem. Lett. 2010, 1, 796–801. (19) Bornschlogl, T.; Anstrom, D. M.; Mey, E.; Dzubiella, J.; Rief, M.; Forest, K. T. Biophys. J. 2009, 96, 1508–1514. (20) Gabriel, J. L.; Hoober, J. K. J. Theor. Biol. 1991, 151, 541–556. (21) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E. I.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586–3616. (22) Mroginski, M. A.; von Stetten, D.; Escobar, F. V.; Strauss, H. M.; Kaminski, S.; Scheerer, P.; Gunther, M.; Murgida, D. H.; Schmieder, P.; Bongards, C.; Gartner, W.; Mailliet, J.; Hughes, J.; Essen, L. O.; Hildebrandt, P. Biophys. J. 2009, 96, 4153–4163. (23) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (24) Chenna, R.; Sugawara, H.; Koike, T.; Lopez, R.; Gibson, T. J.; Higgins, D. G.; Thompson, J. D. Nucleic Acids Res. 2003, 31, 3497. (25) Ibragimova, G. T.; Wade, R. C. Biophys. J. 1998, 74, 2906–2911. (26) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781–1802. (27) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. (28) Deloof, H.; Nilsson, L.; Rigler, R. J. Am. Chem. Soc. 1992, 114, 4028–4035.

JP104903U