Molecular Dynamics Simulations of the Chromophore Binding Site of

J. Phys. Chem. B 2009, 113, 945–958

945

Molecular Dynamics Simulations of the Chromophore Binding Site of Deinococcus radiodurans Bacteriophytochrome Using New Force Field Parameters for the Phytochromobilin Chromophore Steve Kaminski, Grazia Daminelli, 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 29, 2008; ReVised Manuscript ReceiVed: October 13, 2008

The conformational flexibility of the tetrapyrrolic phytochromobilin (PΦB) chromophore of the bacteriophytochrome Deinococcus radiodurans (DrCBD) in the Pr state has been investigated by molecular dynamics simulations. Because these simulations require accurate force field parameters for the prosthetic group, in the present work we developed new empirical force field parameters for the PΦB molecule that are compatible with the CHARMM22 force field for proteins. For this reason, the new force field parameters for the nonbonded (partial atomic charges) and bonded (bonds, angles, dihedrals, improper) energy terms were derived by reproducing ab initio target data following the methodology used in the development of the CHARMM22 force field. This new set of parameters was employed to analyze structural and dynamical features of PΦB inside DrCBD. The 45 ns all-atom molecular dynamics (MD) simulation reveals the existence of two stable conformational states of the chromophore characterized by distinct torsional angles around the C-C bond at the methine bridge connecting rings A and B of the tetrapyrrole. This result supports experimental observations derived from NMR and resonance Raman spectroscopy. Furthermore, statistical analysis of H-bonding events allowed us to identify (a) important H-bonds between the propionic side chains of the chromophore and the apoprotein which may be relevant for the signal transduction step during the photoinduced cycle and (b) a network of eight water molecules which remain in the vicinity of the chromophore during the entire 45 ns production run. 1. Introduction Phytochromes (Phy) constitute a ubiquitous family of sensory photoreceptors, found in plants, bacteria and fungi.1,2 These photoreceptors regulate a variety of processes, ranging from seed germination, formation of leaves and control of flowering time in higher plants to phototaxis and pigmentation in bacteria. Phys are homodimeric complexes with each polypeptide containing three main domains: an N-terminal chromophore binding domain (CBD) with photosensing activity, the phytochrome domain (PHY), which is required for the formation and stability of the physiologically active state, and the C-terminal domain which is essential for signal transduction processes.1,2 Through a unique interaction between the chromophore, a linear-methine bridge tetrapyrrole (bilin), and the protein matrix (mainly the CBD and PHY domains), phytochromes can interconvert, upon light absorption, between two stable states: a physiologically inactive red-absorbing (Pr) form and an active far-red-absorbing (Pfr) form.3 The underlying molecular events controlling this process are still not fully understood, and despite the numerous studies in this field, the elucidation of the photointerconversion mechanism is still a matter of debate. For this reason, the determination of the three-dimensional structure of the Pr state of a bacteriophytochrome deletion mutant from Deinococcus radiodurans is of particular importance. The 1.45 Å resolution structure of the chromophore binding domain (DrCBD) reveals that the biliverdin chromophore binds to the protein as a 2(R),3(E)-phytochromobilin (PΦB) chromophore with a ZZZssa configuration and confor* Corresponding author.

mation at the methine bridges, with a nearly coplanar arrangement of the A, B and C pyrrole rings (see Figure 1).4,5 Furthermore, resonance Raman spectroscopic studies on DrCBD crystals6,7 suggest the coexistence of two or more PΦB conformers which differ, most likely, in the torsional angles at the A-B methine bridge. The same phenomenon has been detected for the Agrobacterium tumefaciens phytochrome (Agp1) crystals,6 plant phytochrome (Phy A),8 and cyanobacterial phytochrome in solution (Cph1).9 This structural heterogeneity of the bilin chromophore seems to be an intrinsic property of phytochromes in general. Unfortunately, it cannot be captured by crystallographic techniques. Therefore, identification and characterization of possible chromophore conformers in the Pr state and their influence on the photoactivation mechanism remains an open question. We address this question from a computational point of view, by analyzing the conformational fluctuations of the bilin chromophore inside the protein binding pocket. Full-atoms molecular dynamics (MD) simulations constitute a powerful tool for performing this task in an efficient way. The accuracy of this computational approach depends on the definition of the empirical potential energy function and on the reliability of the corresponding force field parameters. The widely used all-atom CHARMM force field has been already successfully parametrized for proteins10 and can be therefore applied for modeling the natural amino acids of the phytochrome apoprotein. The bilin chromophore, however, does not belong to these building blocks. To our knowledge, no empirical force field parameters have been developed so far for linear-methine bridge tetrapyrrole molecules. Thus, the first step toward a proper description of

10.1021/jp8047532 CCC: $40.75  2009 American Chemical Society Published on Web 01/05/2009

946 J. Phys. Chem. B, Vol. 113, No. 4, 2009

Kaminski et al.

Figure 1. 3-D structure of DrCBD (left) and structure of model compound 2(R),3(E)-PΦB and fragmentation scheme for optimization of torsion angle parameters (right). New atom types are written in bold letters.

the structural and dynamical properties of the entire phytochrome protein is the development of adequate force field parameters for the bilin chromophore which are compatible with the already existing CHARMM all-atom force field.10,11 Starting from the available X-ray structure of DrCBD in the Pr state, we employed empirical molecular dynamics to identify and characterize possible conformers of the PΦB chromophore, which may be relevant for the interpretation of the experimental data.6,7 The present paper is separated in two parts: the first part (section 2) describes the development and validation of new force field parameters for the PΦB chromophore, whereas the second part (section 3) reports the results of molecular dynamic simulation of DrCBD in solution. Special emphasis is given in elucidating the conformational flexibility of the chromophore and the stability of the hydrogen bond and water network in the protein binding pocket during the MD simulation. 2. Development and Validation of Force Field Parameters The optimization and evaluation of new force field parameters for phytochromobilin was performed following the iterative procedure used for developing the CHARMM22 force field.10 This procedure is based on the reproduction of target data at a molecular mechanics level. The target data include structural data, energies (e.g., rotational barriers) and intermolecular interaction energies of model compounds estimated from ab initio calculations.10,12,13The model compounds are simply isolated molecules containing the structural units (bond lengths, bond angles and dihedral angles) associated with those parameters that are to be optimized. The parametrization strategy is then to maximize the agreement between the target data and the corresponding molecular properties of the model compounds computed using the empirical force field. The use of this parametrization approach guaranties consistency with the available force fields CHARMM2210 and CHARMM2711 employed to model the phytochrome apoprotein. Details about the parametrization procedure are described in the paragraphs below.

Ab initio calculations required for the generation of target data were performed with the Gaussian 03 program,14 whereas molecular mechanics calculations and molecular dynamics simulations were done using the CHARMM15 version 32b2 and NAMD 2.616 codes. 2.1. Definition of the Potential Energy Function. The CHARMM potential energy function is defined as the sum of various internal (bonded) and nonbonded terms in the following way:

∑ Kb(b - b0)2 + ∑ KUB(S - S0)2 +

U(rj) )

∑

bonds

Kθ(θ - θ0)2 +

angle

∑

improper

∑

UB

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

dihedral

Kφ(φ - φ0)2 +

∑

nonbond

[( ) ( ) ]

εij

Rij rij

12

-2

Rij rij

6

+

qiqj (1) 4πεrij where K denotes force constants for bonds (b), Urey-Bradley 1,3-distances (S), bond angles (θ), dihedral angles (χ) and improper dihedral angles (φ) and b0, S0, θ0, χ0 and φ0 are the corresponding equilibrium values of the individual terms. The nonbonded interactions are described as the sum of LennardJones 6-12 and Coulomb terms. The first depends on the Lennard-Jones parameters, εij and Rij, for the atom pair i and j and the distance rij between them. The Coulomb terms are a function of the atomic partial charges qi and qj and their relative distance rij. 2.2. Model Compounds. The model compound for the protein-bound chromophore consists of an isolated protonated 2(R),3(E)-phytochromobilin (PΦB) molecule in a ZZZssa conformation/configuration of the methine bridges as found in the X-ray structure of DrCBD shown in Figure 1. The propionic side chains were protonated (neutral) and the linkage of the chromophore to the protein was saturated by a hydrogen atom.

MD Simulations of the Chromophore Binding Site of DrCBD

J. Phys. Chem. B, Vol. 113, No. 4, 2009 947

TABLE 1: Atom Types, Labels and Optimized Partial Charges of the PΦB Moleculea atom

atom type

partial charge

atom

atom type

partial charge

atom

atom type

partial charge

CAC HAC C1C H2C N_C H_C C4C C3C C2C O_C CHD HHD C1D N_D H_D C4D C3D C2D C1A N_A H_A C2A C3A C4A C4B N_B H_B C1B C2B

CPM HA C HA NR1 H CA CPY1 CT1 O CPY3 HA CPA NR1 H CPA CPB CPB CPA NR1 H CPY4 CPB CPA C NR1 H CA CPY5

-0.082 0.123 0.289 0.057 -0.555 0.332 0.312 -0.160 0.315 -0.445 -0.501 0.214 0.375 -0.543 0.306 0.364 -0.091 -0.319 0.119 -0.654 0.355 -0.038 -0.155 0.487 0.476 -0.560 0.425 0.452 -0.115

C3B O_B CHB HHB CAB HAB CBB HV1 HV2 CBC HL1 HL2 CMC HE1 HE2 HE3 CMD HD1 HD2 HD3 CMA HA1 HA2 HA3 CMB HB1 HB2 HB3 CBD

CPY6 O CPY2 HA CE1 HE1 CE2 HE1 HE2 CT2 HA HA CT3 HA HA HA CT3 HA HA HA CT3 HA HA HA CT3 HA HA HA CT2

-0.091 -0.407 -0.448 0.219 -0.332 0.307 -0.269 0.189 0.153 0.043 0.058 0.058 -0.119 0.038 0.038 0.038 -0.024 0.022 0.022 0.022 -0.075 0.050 0.050 0.050 -0.018 0.049 0.049 0.049 -0.28

HO3 HO4 CGD O2D O1D CBA HO1 HO2 CGA O2A O1A CAA HO5 HO6 CAD HO7 HO8

HA HA CC OC OC CT2 HA HA CC OC OC CT2 HA HA CT2 HA HA

0.09 0.09 0.62 -0.76 -0.76 -0.28 0.09 0.09 0.62 -0.76 -0.76 -0.18 0.09 0.09 -0.18 0.09 0.09

a

Values related to the propionate chains are written in italic letters. New atom types are in bold letters.

To account for electron delocalization effects, the entire PΦB molecule was used for the optimization of the internal bond lengths, bond angles and improper dihedral parameters and the nonbonded electrostatic parameters. Figure 1 (right) shows the chromophore structure together with the atomic labels for the heavy atoms. The list of all atoms belonging to the PΦB chromophore and their corresponding atom type is given in Table 1. Six new atom types were specifically created for the PΦB molecule to be able to define new torsion potentials and consequently improve the accordance with the internal properties of the target structure, as will be shown later. For the evaluation of torsion parameters for the chromophore, the PΦB molecule was divided into smaller fragments (see Figure 1 (right)) to reduce computational cost. The hexamethylpyrromethene (HMPM) and pentamethylpyrrolone (PMPO) molecules were used to model the torsions at the methine bridge between inner rings B and C and the methine bridges between inner rings B/C and the outer ring A/D, respectively. The torsion parameters at the chromophore-protein linkage site were evaluated using the [(aminosulfanyl)oxo]pyrrolidine acid (ASOP) and for the parameter associated with the vinyl function at ring D we used a methylenevinylpyrrolone (MVPO) molecule a as model compound. In this way all torsion angle parameters which are relevant for the chromophore flexibility were accurately and efficiently determined. 2.3. Target Data. The target data for optimizing internal force field parameters such as bond, angle and improper torsion parameters, was derived from the DFT optimized geometries of the entire PΦB chromophore. The same structure has been used for the evaluation of partial charge. The use of a complete PΦB molecule as the target structure guaranties the correct description of the electron delocalization effects. However, to reduce computational cost, the dihedral parameters were developed using the molecular fragments mentioned in the

previous section. All ab initio calculations were performed in vacuo using Becke’s three parameters functional (B3LYP)18 and the 6-31G(d) basis set. 2.4. Nonbonded Parameter. The optimization of nonbonded force field parameters involves electrostatic and van der Waals interactions. The main goal is to obtain a set of parameters that can reproduce, in a reliable way, the water-solute interaction, where water is described through a TIP3P model.17 Following MacKerell’s procedure,10 the partial atomic charges for the PΦB chromophore were obtained by computing, at a quantum mechanical level, minimum interaction energies and distances between the model compound and water molecules. We performed the calculations using fourteen water molecules located at various sites as shown in Figure 2. The PΦB geometry was previously optimized in vacuo using the hybrid density functional method B3LYP18 with the 6-31G(d) basis set whereas the water molecules were built according to the TIP3P water model.17 The minimum interaction distances were estimated through single point energy calculations of the various PΦB-water complexes at the HF/6-31G(d) level of theory. These complexes were built by systematically varying the interaction distances in steps of 0.01 Å. Interaction energies were then estimated by subtracting the sum of the monomer energies from the energy of the entire complex. Because the chromophore is a charged compound, the target HF/6-31G(d) interaction energies where not scaled.10 Initial partial atomic charges were obtained from fitting the electrostatic potential using the CHELPG19 method implemented in Gaussian03 while simultaneously reproducing the overall dipole moment of the chromophore. The partial charges of atoms belonging to the propionate side chains were not optimized. These values were taken form the already existing parameters for the heme cofactor20 without further adjustment.


Kaminski et al. TABLE 2: Optimized Parameters for Newly Defined Dihedral Angles of the PΦB Chromophore n

χ0 (deg)

Linking Fragment 1 (ASOP) CT2-SE-CT2-CPM 0.7 CT2-SE-CT2-CPM 1.1

3 1

10.0 88.0

Linking Fragment 2 (ASOP) SE-CT2-CPM-CPY1 1.4

2

45

Linking Fragment 3 (ASOP) CT2-CPM-CPY1-CA 2.4 CT2-CPM-CPY1-CA 16.4

2 1

180.0 4.0

dihedral parameter

Kχ (kcal/mol)

Double Bond Pyrrolinone Fragment Ring C-D (PMPO) CPA-CPY2-CA-CPY5 11.8 1 0.0 CPA-CPY2-CA-CPY5 3.0 2 190.0

Figure 2. Relative orientations of fourteen water molecules with respect to the PΦB molecule for the evaluation of partial charges.

The partial atomic charges obtained for the PΦB molecule are listed in Table 1. The quality of the optimized partial charges was estimated by reproducing minimum interaction energies and waterchromophore distances derived at the ab initio level of theory (see Table 4). Distances were deliberately underestimated in about 0.2 Å compared to the target data following the established parametrization procedure.10 Because an accurate reproduction of both interaction energies and distances was impossible we always tried to find in all cases a good compromise between them. As an additional criterion for evaluating the quality of the partial charges we simultaneously tried to reproduce the overall dipole moment. To account for the missing explicit inclusion of polarizability in classical force fields, the dipole moment should be overestimated10 to a certain extent compared to the values obtained with QM calculations. Therefore, to be consistent with CHARMM force field parameters for proteins, the partial charges of the chromophore were adjusted so as to obtain an empirical dipole moment 10% larger than the ab initio target values. Empirical- and QM-calculated dipole moments as well as interaction energies and distances for the PΦB molecule are given in Table 4. Lennard-Jones parameters ij and Rij were calculated by means of the Lorentz-Berthelot mixing rules using the CHARMM27Lennard-Jones parameters of atoms types. Because these atomicLennard-Jones parameters are transferable between atom types which are chemically similar, those parameters required for the PΦB chromophore were simply collected from other compounds such as heme. Therefore, in this work we did not optimize any Lennard-Jones atomic parameters for PΦB. 2.5. Internal Parameters. The parametrization of the internal force field parameters of PΦB was done iteratively starting from a set of reasonable initial values from already parametrized analogue structures. Most of the initial internal parameters (bond, angles, dihedral, improper) and atom types correspond to the CHARMM parameters previously optimized for the heme molecule20 (toppar_all22_prot_heme.str). Among all internal force field parameters of the PΦB chromophore defined in eq 1, in the present work we optimized the corresponding equilibrium bond lengths (b0), bond angles (θ0), improper (φ0) and dihedral angles (χ0) as well as harmonic force constant Kφ, Kχ and periodicity and phase parameter (n, δ) for the dihedral torsions. The harmonic force constants Kb, KS, Kθ together with the all parameters of the propionate side chain were kept fixed during the optimization procedure.

Single Bond Pyrrolinone Fragment Ring C-D (PMPO) CPB-CPA-CPY2-CA 2.0 3 0.0 CPB-CPA-CPY2-CA 4.4 2 138.0 CPB-CPA-CPY2--CA 2.1 1 243.0 Pyrrol Fragment CPY4-CPA-CPM-CPA CPY4-CPA-CPM-CPA CPY4-CPA-CPM-CPA

1 (HMPM) 9.2 4.5 2.2

2 1 3

183.0 2.0 190.0

Pyrrol Fragment 2 (HMPM) CPB-CPA-CPM-CPA 8.9 CPB-CPA-CPM-CPA 1.5 CPB-CPA-CPM-CPA 1.9

1 2 3

184.0 0.0 196.0

Double Bond Pyrrolinone Fragment Ring A-B (PMPO) CPY1-CA-CPY3-CPA 1.9 1 0.0 CPY1-CA-CPY3-CPA 6.4 2 178.0 CPY1-CA-CPY3-CPA 1.8 3 160.0 Single Bond Pyrrolinone Fragment Ring A-B (PMPO) CA-CPY3-CPA-CPB 1.2 3 62.0 CA-CPY3-CPA-CPB 4.4 2 162.0 CA-CPY3-CPA-CPB 2.5 1 252.0 Vinyl Fragment CPY5-CPY6-CE1-CE2 CPY5-CPY6-CE1-CE2 CPY5-CPY6-CE1-CE2

(MVPO) 0.9 2.35 0.5

1 2 3

0.0 180.0 0.0

The target data consists of the DFT optimized geometries of the model compound (for evaluation of bond, angle and improper torsion parameters) and of rotational energy profiles of its fragments (for the evaluation of dihedral torsion parameters). The quantum mechanical geometry optimizations were carried out under tight convergence criteria to a final gradient of 10-5 au, using the X-ray crystallographic structure as input geometry. All geometry optimizations in CHARMM involving the entire chromophore and fragments were performed using the adopted basis Newton-Raphson method (ABNR) followed by several steps using the Newton-Raphson (NRAPH) minimizer. To ensure a reliable conformational behavior of the chromophore, ten new harmonic torsion potentials were defined (see Table 2) to describe rotations at the methine bridges, the linkage to the protein (CYS 24) and at the vinyl functionality on ring D. The corresponding torsional force field parameters were determined by fitting the torsional profiles of the various molecules (Figure 3) obtained by quantum mechanical calculations with the corresponding molecular mechanics potential energy function. This function was constructed according to eq 1 and employing, in a first stage, the initial internal parameters chosen for the entire PΦB molecule and the corresponding point charges calculated as described in section 2.4. The QM and MM potential energy curves were scanned in steps of 10° by



Figure 3. Potential energy surfaces for partially constrained model compounds. The surfaces were scanned in steps of 10°. The black curves denote the target data derived using ab initio methods (DFT/B3LYP/6-31G(d)), and the gray dashed-curves result from optimized torsion parameters (see Table 2). Potential curves were shifted by setting the lowest energy conformations for each curve equal to zero.

performing constraint energy minimizations with a rms force criterion of 10-4 for the ab initio calculations and 10-6 for

analogue molecular mechanics calculations. In each step only the torsional angle in question was kept fixed to a certain value


Kaminski et al.

TABLE 3: Optimized Parameters for All Newly Defined Internal Coordinates (Except Torsions) of PΦB bond parameter

Kb [(kcal/mol)/Å2]

b0 (Å)

bond parameter

Kb [(kcal/mol)/Å2]

b0 (Å)

NR1-CPA NR1-C CT2-CPM NR1-CA CPY1-CA CT1-CPY1 CPM-CPY1 CPY2-CA HA-CPY2 CPY2-CPA

377.2 260.0 230.0 377.2 305.0 230.0 360.0 360.0 367.6 360.0

1.3817 1.4190 1.4830 1.4117 1.4772 1.4907 1.3266 1.3800 1.0902 1.4455

CPY3-CA CPY3-CPA CPY4-CPA CPY4-CPB CT2-CPY4 CPY5-CA CT3-CPY5 CPY6-C CPY6-CE1 CPY6-CPY5

360.0 360.0 299.8 340.7 230.0 305.0 230.0 250.0 450.0 305.0

1.3676 1.4225 1.4052 1.4018 1.4946 1.4620 1.4850 1.4711 1.4290 1.3470

angle parameter

Kθ [(kcal/mol)/rad2]

θ0 (deg)

angle parameter

Kθ [(kcal/mol)/rad2]

θ0 (deg)

NR1-CPA-CPB CPA-NR1-CPA CPM-CPA-NR1 NR1-C-O HA-CT2-CPM C-NR1-CA HA-CPM-CPY1 CT2-CPM-CPY1 NR1-CA-CPY1 CT1-CPY1-CPM CT1-CPY1-CA CPY1-CT1-C CPM-CPY1-CA CPY2-CPA-NR1 CPY2-CA-NR1 CA-CPY2-CPA

122.0 139.3 88.0 80.0 49.3 50.0 12.7 45.8 122.0 45.8 45.8 52.0 61.6 88.0 88.0 94.2

110.0 116.3 131.8 118.0 107.5 133.5 117.44 117.49 115.0 116.6 128.0 113.7 124.1 112.39 129.0 127.0

CPY1-CA-CPY3 CA-CPY3-CPA CPY3-CPA-CPB CPM-CPA-CPY4 CPB-CPY4-CPA CPY4-CPB-CPA CT2-CT2-CPY4 NR1-CA-CPY5 CPY5-CA-CPY2 CT3-CPY5-CA NR1-C-CPY6 C-CPY6-CPY5 C-CPY6-CE1 CE2-CE1-CPY6 CA-CPY5-CPY6

61.6 94.2 61.6 61.6 30.8 30.8 70.0 122.0 61.60 45.8 20.0 52.0 70.0 40.0 40.0

127.57 122.80 125.07 132.50 136.01 145.01 114.7 112.4 124.97 115.9 109.5 116.0 114.49 117.6 116.5

improper parameter

Kφ [(kcal/mol)/rad2]

φ0 (deg)

improper parameter

Kφ [(kcal/mol)/rad2]

φ0 (deg)

CPA-CPB-NR1-CPM CPA-CPB-NR1CPY2 CPA-CPB-NR1CPY3 CPY1-NR1-CACPY3

200.0 140.0 145.0 140.0

0.0 0.0 0.0 180.0

CPA-CPM-NR1CPY4 CPY4-CPA-CPB-CT2 CPY5-NR1-CA-CPY2 CPY6-CE1-C-CPY5

61.0 64.0 7.0 180.0

0.0 0.0 180.0 0.0

while the rest of the structure was allowed to relax. Because of the chemical similarity between the outer rings A-B and C-D (see Figure 3C) one set of target data (here the structure of ring C-D was chosen) was sufficient to model the corresponding methine bridge torsions. Due to electron conjugation at the central methine bridge, the ab initio potentials for the two bond torsions are practically identical. We now can summarize the iterative procedure for the parametrization of the nonbonded parameters. At first, partial charges and internal parameters (no torsions) have been derived for the entire chromophore. These parameters have been transferred to the smaller fragments to evaluate torsion parameters that were transferred back to the chromophore to refine the residual internal parameters and partial charges until selfconsistency within all force field parameters for the chromophore was achieved. The optimized bond length, bond angle and improper torsion parameters for the protein bound chromophore are listed in Table 3 whereas the corresponding optimized dihedral parameters are given in Table 2. Furthermore, it is important to mention that, for the MM in vacuo calculations, the force field parameters of the PΦBpropionic side chain in its neutral form were taken from the acetic acid included in the CHARMM force field (see toppar_all22_prot_model.str/RESI ACEH). In the condensed phase, molecular dynamics simulations of PΦB with charged propionate side chain were performed using the force field parameters derived for the heme cofactor. During the parametrization of the dihedral torsions emphasis was placed on the accurate reproduction of the local minima,

heights of rotational barrier as well as the overall shape of the potential surfaces (Figure 3). Most problematic in this respect was the potential describing single bond rotations of the linking fragment to the protein (Figure 3f,g), because of the high flexibility of this structure. Because we are not interested in predicting Z/E isomerization energy barriers at double bonds, only a certain region of the potential energy curve around a local minimum was scanned. Within these limited dihedral-angle intervals the torsion potentials were in general very well reproduced. However, for the rotation around the CdC bond at the AB methine bridge (Figure 3a) although the overall curvature and local minimum at around 195° could be well predicted, the energy barrier and the second local minimum slightly differ from the QM target. 2.6. Validation of Force Field Parameters. Validation of the final set of all internal and nonbonded parameters optimized for the PΦB cofactor was done by comparing the energy minimized structures of the isolated PΦB molecule and the model fragments (see section 2.2) obtained at the QM and MM levels. Energy minimizations of the PΦB molecule were performed starting from the crystal structure (the protein-bound PΦB in DrCBD). To avoid the formation of hydrogen bonds between the propionate side chains and the PΦB backbone during minimizations in vacuo, which may distort the validation of the new force field parameters, the PΦB-acidic groups were removed. Figure 4 (bottom) shows the resulting MM- and QMminimized structure of PΦB and the corresponding root-meansquare deviations (rmsd) for all heavy atoms after optimal



TABLE 4: Interaction Energies and Distances for Chromophore-Water Complexes Evaluated at the ab Initio (HF/6-31G(d)) and Empirical Levels interaction

HF/6-31G(d) Emin (kcal/mol)

empirical/CHARMM Rmin (Å)

Emin (kcal/mol)

Rmin (Å)

1 -2.549 3.39 -2.762 3.38 2 -3.789 3.27 -3.883 3.43 3 -5.418 2.63 -5.386 2.58 4 -6.154 2.60 -6.036 2.55 5 -3.059 3.39 -3.0427 3.39 6 -3.095 4.09 -4.479 3.40 7 -2.329 3.33 -2.497 3.33 8 -2.672 3.80 -2.583 3.51 9 -2.909 2.06 -3.887 1.82 10 -3.692 2.61 -4.298 2.54 11 -2.450 3.41 -2.597 3.34 12 -0.739 3.25 -0.977 3.26 13 -5.489 2.51 -5.648 2.36 14 -4.966 3.20 -5.019 3.14 average differential energy: 0.307 kcal/mol average difference distance: 0.184 Å dipole moment (Debye)/HF-6-31G(d): total ) 13.370, x ) -13.369, y ) 0.104, z ) -0.162 dipole moment (Debye)/empirical: total ) 14.646, x ) -14.556, y ) -1.426, z ) 0.772

deviations ∆E (kcal/mol)

∆R (Å)

-0.213 -0.094 0.032 0.118 0.016 -1.384 -0.168 0.082 -0.98 -0.61 -0.147 -0.238 -0.159 -0.053

-0.01 0.16 -0.05 -0.05 0.0 -0.69 0.0 -0.29 -0.24 -0.07 -0.07 0.01 -0.15 -0.06

the torsional angles toward all other internal and nonbonded parameters. The average deviation for all 15 defined improper angles is 0.665°, which is again satisfactory. In addition to the PΦB chromophore, we evaluated the overall structural agreement between the QM and MM optimized structures of the small model fragments shown in Figure 4. As starting structures for the unconstrained geometry optimizations we chose geometries of the fragments lying in the vicinity of important conformational local minima. These unconstrained geometry minimizations we performed using the same conditions as described above for minimizations of the entire chromophore. Taking as a reference the equilibrium structure of PΦB in the Pr state, we chose a trans-like conformation of the fragment mimicking rings C and D and a cis-like conformation of the fragments involving rings A and B as well as rings B and C. For the fragment describing the attachment of the cofactor to the protein we used a starting geometry where all three considered torsion angles were near their local minima. Concerning the fragment with the vinyl functionality, both cisand trans-conformations were taken into account. The resulting rms deviations between the optimized QM and MM structures listed in Figure 4 are satisfactory.

Figure 4. Optimized ab initio (DFT/B3LYP/6-31g(d)) (dark gray) and empirical (light gray) structures of the PΦB chromophore without propionate chains and the small model fragments used in the parametrization proccedure. rmsd are given for non-hydrogen atoms after optimal superposition.

superposition of the two structures. In the case of the PΦB molecule there is an overall very good agreement between the MM-optimized and the ab initio structures as indicated by the rmsd of only 0.282 Å. The average deviation for 39 bond lengths between QM and MM optimized structure including all heavy atoms is about 0.009 Å, which is satisfactory for our needs as well as the deviation of 0.941° for 40 bond angles. Significantly larger is the average difference of 3.393° involving all methine bridge, vinyl and link-fragment torsion angles. This large deviation can be explained in terms of a higher sensitivity of

3. Molecular Dynamics Simulations 3.1. Computational Details. To check the performance of the new force field parameters for PΦB in a protein environment, we performed molecular dynamics simulations of the entire CBD from Drph1 phytochrome photoreceptor. Initial coordinates for the protein were taken from the Protein Data Bank (reference: 2O9C/conformation A/resolution: 1.45 Å). Missing heavy atoms of several residue side chains (ASP 4, ARG 70, GLU 193, HIS 196, ARG 310) were added using the CHARMM software, subsequently hydrogen atoms were incorporated to the structure by means of the HBUILD routine.21 The assignment of a proton to titratable groups was done on the basis of visual inspection of the environment of the charged amino acids and histidine residues. In this way, all lysine (2) and arginine (19) residues found in the crystal structure were protonated and all aspartate (11) and glutamate residues (19) were considered to be deprotonated. His 260 and His 290 were modeled with a hydrogen at the δ-nitrogen whereas the remaining histidine residues were protonated at the -nitrogen.


Kaminski et al.

Figure 5. Root-mean-square deviations (rmsd) for the protein (light gray), its backbone atoms (dark gray) and the chromophore heavy atoms (black) and as a function of time.

This gives a strongly negatively charged protein with a total charge of -10 e. Afterward, the protein was solvated in a cubic box of water with an initial volume of 884736 Å3 (26307 solvent water molecules). The shortest initial distance between protein and water box walls was approximately 15 Å. This distance is deliberately larger than the chosen van der Waals long-range cutoff (12 Å) to avoid long-range interactions between the protein and its images. The TIP3P potential was used for the water molecules.17 The protein charge was neutralized by inserting chloride and sodium counter-ions in the following manner: each charged residue not involved in salt bridges was individually neutralized by placing a counterion within a 3 Å radius from it, whereas the remaining ions were placed randomly in the water box. All subsequent calculations concerning heating, equilibration and production were performed with the NAMD program (version 2.6)16 using the CHARMM27 force field. At first the solvent waters were heated to 300 K and equilibrated for 40 ps while keeping all other parts of the system fixed. Subsequently, several energy minimization runs (2000 steps of conjugate gradient) were performed in which the initial harmonic restraints (10 kcal/(mol Å2)) applied to the chromophore and all heavy backbone atoms of the protein were gradually released until the entire system was free. These harmonic positional restrains were described with a quadratic potential energy function. Afterward, the system was heated up to 300 K during 60 ps using Langevin dynamics with restraint chromophore and protein backbone atoms. These steps were followed by 100 ps molecular dynamics simulation under constant pressure, constant temperature (NPT) conditions using a combination of the Langevin Piston Nose-Hoover method,22,23 as implemented in NAMD. Here again, the harmonic quadratic restraints were stepwise released until the system was totally unrestrained. After the system was equilibrated, we proceeded with the production run consisting in 50 ns molecular dynamics simulation using a reduced Langevin damping factor (varying from 5.0 to 1.0) to more closely approximate free dynamics (NPT conditions at 300 K). All simulations were done under periodic boundary conditions with the long-range electrostatic interaction described via the particle-mesh-Ewald method.24 For the van der Waals interactions a cutoff (12 Å) was used in combination with a switching function. To use a 2 fs time step, all bond lengths between heavy atoms and hydrogen have been con-

strained to their minimum energy values by applying the SHAKE algorithm.25 All molecular structures were drawn using the visual Molecular Dynamics Software VMD 1.8.6.26 3.2. Dynamic Properties of the Protein. In Figure 5 the root-mean-square deviations (rmsd) for all heavy atoms of the protein and chromophore were plotted over the simulation time using the first frame of the simulation as the reference structure (conformations were stored each 1 ps). Because the rmsd increased progressively during the first 5 ns of the MD simulation, they were considered as part of the equilibration phase and were not used for further analysis. The rmsd during the remaining 45 ns of the simulation, on the other hand, remained stable and were therefore used for the statistical evaluation of the trajectory. The average of the rmsd evaluated over the 45 ns of the production run yield 2.4 Å for the entire protein, 1.9 Å for the protein backbone and 1.3 Å for the chromophore (excluding hydrogen atoms). These values are similar to those previously calculated for cytochrome c oxidase.13 The slightly higher average rmsd for the protein results from the larger mobility of the protein side chains. During the simulation the overall shape of the protein is conserved as indicated by the small fluctuations (0.12 Å) of its calculated average radius of gyration of 20.33 Å. The rmsd involving all heavy atoms of the average protein structure resulting from the MD simulation and the crystal structure of DrCDB is 1.82 Å. The mayor structural differences between these two structures are found in the unstructured coil region near the N-terminus (residues 4-11) as well as in the three R helices defined by the residue sequences 58-62, 66-69 and 81-88. Structurally very conserved on the other hand, are the six β strands surrounding the chromophore. To analyze the mobility of the molecular system, the atomic root-mean-square fluctuations (rmsf) were calculated from the MD simulation. Under the assumption of isotropic atomic fluctuations, the rmsf are related to the experimental B-factors (Debye-Waller or temperature factor) from the X-ray structure via the following expression27

〈∆ri2〉 )

3 Bi 8π2

(2)

where Bi and the square root of 〈∆ri2〉 represent the temperature factor and the root-mean-square fluctuations associated with



Figure 6. Per residue averaged root-mean-square fluctuations over all atoms derived from of experimental B-factor (black) and calculated from MD simulation (red).

atom i, respectively. The 〈∆ri2〉 of each atom was calculated with respect to its average position during the MD simulation. Figure 6 shows the experimental and the calculated rmsf averaged over the individual amino acids. The theoretical values are systematically larger than the experimental ones by almost a factor of 2. Therefore, one should compare these values in a qualitative manner28 to determine whether the relative mobilities predicted by the MD simulation for different regions of the system are in accordance with the experiment. In fact, the results plotted in Figure 6 indicate that for DrCBD that is, in general, the case. For example, we see that for one of the most flexible parts in the protein, a 3-10 helix (residues 130-134), simulation and experiment are in a qualitative good agreement. However, large deviation in flexibility can be observed for the loop region between two β-strands comprising residues 279-282 as well as for the loop region between residues 104-108. In particular, the predicted mobility of residues lying in the vicinity of the chromophore, such as CYS24, ASP207, MET 259, HIS 260, SER 272, SER 274 and HIS290, agrees very well with the corresponding B-factors. 3.3. Dynamic Properties of the Chromophore. When analyzing the 45 ns production run of the MD trajectory, we noticed a particularly interesting behavior of the chromophore and its interaction with residue HIS 260. During the simulation, the system seems to jump between two states characterized by the presence (state 1) or absence (state 2) of a water molecule localized between rings B and C, the so-called pyrrole water illustrated in Figure 11. In state 1, the pyrrole water is hydrogen bonded to HIS260 and to the NH group on ring A as observed in the crystal structure. In state 2, on the other hand, this pyrrole water moves away from the binding pocket and HIS260 comes closer to ring A strengthening its interaction with the chromophore. Each state is therefore characterized by distinct N_C (ring A in PΦB)-ND1 (H260) distances as shown in Figure 7. To evaluate in more detail the structural properties of the chromophore in each of these states, average structures were

calculated over the time intervals (each 4 ns long) indicated in Figure 7 with dashed boxes. The results are presented in Figure 8, where the average structures for state 1 (left) and for state 2 (right) are superimposed with the crystallographic structure. For state 1, we obtain an overall good agreement concerning the chromophore and linkage to the protein (rmsd for PΦB of 0.58 Å). The distance between ring A of the chromophore and HIS 260 (5.23 Å) is comparable to the X-ray structure (4.99 Å). The average structure of state 2 reveals, on the other hand, significant distortions of the chromophore structure (rmsd for PΦB: 0.92 Å) that mainly involve ring A and ring C. Unlike the crystal structure, ring A moves out of the plane formed by the two coplanar inner rings. The dihedral angles at the A-B methine bridge of the average structures predicted for these two chromophore states are 13.41° (single bond)/18.44° (double bond) in state 1 and 18.48° (single bond)/17.21° (double bond) in state 2. These torsional angles are significantly larger compared to those measured in the X-ray structure: 7.55° (single bond)/6.773° (double bond). In addition, from an energetic point of view, we cannot distinguish one state from the other. The energy plots along the simulation show that the protein and the chromophore remain stable. This would be in agreement with the fact that at room temperature and during 45 ns simulation both states are almost equally occupied (see Figure 7). Other regions of the chromophore are structurally similar for both states, although slight deviations also exist in the vicinity of ring D. The conformational flexibility of the chromophore was further investigated by analyzing the statistical distribution of all dihedral angles at the methine bridges. Except for the single bond torsion between rings A and B (A-B methine bridge) all distributions could be perfectly fitted with a single Gaussian function with widths oscillating between 8° to 10°. As an example we plotted in Figure 9 (bottom, left) the statistical distribution of the dihedral angle at the double bond of the A-B methine bridge. In this case a single Gaussian function with


Kaminski et al.

Figure 7. Time dependent distance between ring A (N_C) of the chromophore and HIS 260 (ND1). Dashed boxes indicate the time interval from which average structures of state 1 and state 2 were determined.

Figure 8. Chromophore binding pocket: Crystallographic structure (dark gray) superimposed with average structures (light gray) resulting from a 4 ns MD simulation of state 1 (left) and state 2 (right). rmsd values for the heavy atoms of PΦB are 0.58 Å (state 1) and 0.92 Å (state 2).

maximum at 18.7° and a bandwidth of 14° is sufficient for a satisfactory description of the dihedral angle distribution. The distribution of the A-B torsion around the single bond, however, displays an asymmetric shape (see Figure 9 bottom right), which can only be fitted if at least two Gaussian components are considered. The maxima of these two components are located at 9.5° and 16.2° with bandwidths of 12° and 9°, respectively, which are slightly larger than those required for fitting the angular distribution of the other torsions. Because the band widths of the distribution curves reflect the degree of fluctuation associated with each torsional angle, we conclude that the A-B methine bridge is more flexible than the B-C and C-D methine bridges. This conclusion is in agreement with the

experimental X-ray B-factors for the chromophore which are larger for ring A than for the rest of the PΦB chromophore. In addition, these results would also support experimental observations derived from NMR9 and resonance Raman spectroscopy.6 By plotting the N_C-ND1(HIS260) distance as a function of the dihedral angles at the A-B methine bridge (Figure 9, top left and right) we could extract some information about how these quantities correlate. The scatterplots for the two dihedral angles show basically two cluster of points located at the average N_C-ND1 distances characterizing the protein’s state 1 and state 2 described above. In this respect, the upper cluster corresponds to state 1 whereas the lower cluster stands for state 2. For the scatter plot associated with the double bond torsion



Figure 9. Dihedral angle distributions (bottom) for the double (left) and single (right) bond torsion between rings A and B estimated over the 45 ns MD simulation. Scatter plots (top) show the correlation between these dihedral angles and the N_C (ring A)-ND1 (HIS 260) distance distribution.

(Figure 9 top left), both clusters of points are arranged along a line corresponding to the maximum of the Gaussian curve describing the distribution of double-bond dihedral angle at the A-B methine bridge. This indicates that the fluctuations at the A-B double bond torsions are practically not affected by the position of the HIS260 residue. The opposite is found for the A-B single-bond torsion. In this case, in the corresponding scatterplot (Figure 9 top right) the upper points cluster is not only more dispersed than the lower one but also shifted to lower angles, suggesting that the position of the HIS260 influences the average value and fluctuations of the A-B single bond torsion. More precisely, we see that in state 2 (short distance) the torsional angle at the A-B methine bridge is larger than that for state 1 (large distance), indicating that the displacement of HIS60 toward the chromophore brings ring A out of the plane formed by rings B and C. In fact, a N-H · · · N hydrogen bond is built between HIS260 and ring A in PΦB. The presence of the pyrrole water in the protein-binding pocket (state 1) on the other hand, allows for a larger flexibility of the A-B methine bridge. Summarizing, the simulation predict two stable states for the PΦB chromophore in the binding pocket characterized by (a) presence/absence of the pyrrole water, (b) two distinct interaction distances between PΦB and HIS260 and (c) the out of plane movement of the PΦB-ring A. 3.4. Hydrogen Network between the Chromophore and Protein Pocket. To find out which residues in the vicinity the chromophore play an important role in terms of hydrogen bonding, average distances between hydrogen bond donor and acceptor between the chromophore and its surrounding amino acids were estimated. The most important hydrogen bond donor and acceptor pairs involving the chromophore and the protein (H-bond with donor-acceptor distance in the X-ray structure

Figure 10. Hydrogen-bonding network between the PΦB chromophore and surrounding residues in the DrCBD protein pocket. Related statistics are given in Tables 5 and 6.

or during the MD simulation shorter that 3.4 Å31) are illustrated in Figure 10. The corresponding average donor-acceptor distances, evaluated over the 45 ns of the trajectory, are listed in Table 5 together with the experimental values extracted from the crystallographic structure for comparison. The agreement between measured distances from the X-ray and the average values obtained from the MD simulation is in general good, in particular for interactions 1, 2 and 3. For these


Kaminski et al. TABLE 6: Statistics of H-Bond Events between the Chromophore and Surrounding Residues with Occupancies Larger Than 20%a atom 1 (protein)

atom 2 (PΦB)

average lifetime (ps)

events

occupancy (%)

ASP 207 O ASP 207 O ARG 254 HH11 ARG 254 HH11 ARG 254_HH21 ARG 254_HH21 SER 257 HN HSE 260 HE2 SER 272 HG1 SER 274 HG1

H_D H_A O2D O1D O2D O1D O1D O1A O2A O1A

11.4 8.8 60.8 90.9 32.9 26.4 146.3 405.6 24.2 512.1

1436 1016 541 169 930 1011 201 90 1028 60

36 21 73 34 68 59 65 81 55 68

a

Statistical parameters were evaluated over 45 ns MD simulation.

TABLE 7: Statistics of H-Bond Events (with Occupancies Larger Than 20%) between Important Water Molecules and Residues in the Vicinity of the PΦB Chromophorea water

Figure 11. Stable water network in the chromophore binding pocket of Drph1. The structure refers to an average over the last 15 ns of the simulation. Dotted lines indicate acceptor-hydrogen distances less than or equal than 2.4 Å. Related statistics are given in Tables 7 and 8.

TABLE 5: Distances (Å) between Hydrogen Bond Donors and Acceptors Participating in the Interactions Illustrated in Figure 10a H-bond

X-ray

MD

∆R

1 2 3 4 5 6 7 8 9 10

2.85 2.80 2.99 2.67 2.64 2.68 2.80 4.99 2.96 2.77

2.96 (0.42) 3.00 (0.33) 2.89 (0.21) 4.21 (0.71) 3.31 (0.89) 3.29 (0.51) 3.20 (0.80) 4.49 (0.82) 3.44 (0.50) 3.28 (0.33)

0.11 0.20 -0.10 1.54 0.67 0.48 0.40 -0.50 0.48 0.51

a Average distances and fluctuations (in parentheses) are calculated over the 45 ns of the simulation and compared to the crystallographic structure.

three interactions fluctuations of less than 0.40 Å are predicted, indicating that the corresponding H-bonds remain stable during the simulation. These H-bonds formed between the chromophore and ARG254 and SER257 are very important because they are responsible for the rigidity of the propionate side chain at ring B. Furthermore, the two H-bond interactions with the Arginine residue are conserved in other phytochrome species, such as in Cph129 and in RpBpphP3,30 and its rigidity may be of relevance for the signal transduction mechanism between chromophore and protein taking place during the photoinduced reaction cycle. The average H-bond donor-acceptor distances for the interactions 5-10 are slightly larger than the experimental ones (∆ ∼ 0.5 Å). In addition, most of these interactions (5-9) are associated with relatively high standard deviations (0.5-0.8 Å), reflecting a significant mobility of the residue side chains involved in the corresponding interaction. This is in particular the case for the three H-bond interactions involving the propionic

W1 W1 W2 W2 W2 W3 W3 W5 W5 W6 W6 W7 W7 W8 W8 a

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

atom 2


events

occupancy (%)

VAL 252 O LEU 223 HN ARG 254 HH12 ARG 254 O THR 256 OG1 TYR 216 HH PΦB O2D THR 224 HG1 PΦB O2A PΦB O1A PΦB CGA PΦB O1A SER 272 OG PΦB O2A PΦB H_B

138.8 43.0 13.6 28.9 24.6 386.3 261.4 1302.2 49.2 22.7 9.5 69.6 22.2 18.6 225.9

292 823 1301 127 1538 93 117 23 569 1499 2350 496 1273 1228 150

90 79 39 79 84 80 67 67 62 76 50 77 63 51 75

Statistical parameters were evaluated over 45 ns MD simulation.

side chain at ring C and the neighboring SER272, SER274 and HIS260 residues (interactions 5, 6 and 7, respectively), suggesting that this propionic side chain is looser than the one attached on ring B, and probably less efficient for signal transduction. Furthermore, these three interactions become weaker during the simulation as indicated by the slight increase of the corresponding donor acceptor distances compared to the X-ray values. Large fluctuations are predicted for the interactions 7 and 8 involving ring A in PΦB with HIS260 which may be related with the structural flexibility in this region and the movement of the pyrrole water as previously discussed in section 3.3. The most pronounced discrepancy between the MD result and the X-ray H-bond donor-acceptor distance (∆R ) 1.54 Å) is found for interaction 4 defined between the TYR216 and the propionic side chain on ring B. The large average distance of 4.21 Å indicates that the H-bond is broken during the MD simulation. In fact, a closer look at the trajectory shows that for certain time spans a water molecule (W3 in Figure 11) places itself between the TYR216 and the carboxylic group of the propionic side chain at ring B. Additional statistical information concerning the formation of H-bonds in the chromophore binding pocket and their lifetimes is summarized in Table 6. The statistical analysis was performed over the 45 ns MD simulation with a time resolution of 5 ps. The criteria for defining a hydrogen bond was that the distance between hydrogen and the acceptor should be less or

MD Simulations of the Chromophore Binding Site of DrCBD TABLE 8: Statistics of H-Bond Events between Important Water Molecules Surrounding the Chromophore water

water


events

occupancy (%)

W1 W2 W3 W4 W5 W7

W2 W4 W4 W5 W6 W8

39.1 6.9 93.1 7.6 33.0 16.8

905 2077 291 1816 738 1626

79 32 60 31 54 61

equal than a 2.4 Å.31 The stability of H-bonds formed between the protein and the chromophore can be efficiently described in terms of H-bond occupancy as the product between average lifetime and number of events divided by the temporal length of the MD simulation. The results presented in Table 6 show that stable H-bonds are those formed between the ARG254 and SER274 residues and the oxygen atoms at the propionic side chain bound to ring B, with occupancies higher than 65%. This salt bridge between the ARG254 and PΦB seem to be responsible for rigidity of chromophore’s B ring. In a similar way, the propionic side chain attached to ring C forms a stable H-bond with the HIS260 and with the SER274, characterized by long average lifetime of 405 and 512 ps, respectively, and low number of events. The H-bonds formed between the ASP207 and the NH groups of rings B and C, on the other hand, do not seem to be very stable. They are constantly broken and formed during the simulation, as indicated by the large number of events, but each time they are formed they only exist for a very short time (lifetime ∼ 10 ps). Thus the predicted occupancies are fairly low (36% and 21%). 3.5. Water Network in the Chromophore Binding Pocket. To investigate the formation of a water network in the chromophore binding pocket, we selected eight water molecules which remained during the production run within a 5 Å distance from the chromophore. An illustration of the most stable water network we identified is given in Figure 11, showing the geometric average over the last 15 ns of the MD simulation. Among these water molecules, only W3 is already observed in the X-ray structure; W5, W7 and W8 are solvent waters that at an early stage of the simulation occupy the position of a crystallographic water, and the remaining water molecules are solvent water that find a stable position in the vicinity of the chromophore. As described in the previous section, the dynamical behavior of this group of water molecules was analyzed by performing a statistical analysis of the H-bonding events involving the selected waters and calculating the so-called H-bond occupancy. The results of the statistical analysis are listed in Tables 7 and 8. Very stable interactions are the H-bonds formed between W2 and the ARG 254 and THR256. These H-bonds are characterized by a large number of events (1127 and 1538) with average lifetimes of ca. 25 ps and, consequently, large occupation factors (79% and 84%). Similar results are found for the H-bond interactions formed between W6 and the PΦB-propionic side chain on ring C and between W7 and the SER272 residue. An also very important interaction is that between THR224 and W5. This water forms only few H-Bonds with the threonine, but each of them is stable for 1.3 ns. Furthermore, water W8 seems to be an important water bridge connecting the chromophore’s ring D with the propionic side chain attached to ring C. The H-bonds involving W8 and the PΦB chromophore are very stable as indicated by the high occupancy values (75% for the W8-PΦB(H_B) hydrogen bond and 51% for the W8-PΦB(O2A) hydrogen bond). This water

J. Phys. Chem. B, Vol. 113, No. 4, 2009 957 bridge reduces the mobility of ring D in the Pr state. In addition, W3, W5, W6 and W7 establish relatively stable H-bonds with the two PΦB propionate side chains. The corresponding hydrogen bond occupancies are larger than 50%. Although water W4 is not H-bonded neither to the chromophore nor to the apoprotein, its presence seems to be of fundamental importance for the stability of the entire water network. This water is located between the propionate side chains and interacts with W3, W2 and W5 with hydrogen bond occupancies larger than 30%. These H-bonds are constantly formed and broken during the entire production run; however, the relative positions of these water molecules do not significantly change. The most stable water-water interaction is that predicted between W1 and W2 with 79% H-bond occupancy. As shown in Figure 11, these two water molecules do not directly interact with the chromophore. Noticeable is the dynamical behavior of the crystallographic pyrrole water. This water is located in the center of the chromophore cavity forming H-bonds with the NH groups of rings A, B and C. Our molecular dynamics simulations show that these H-bonds are not strong enough to trap the pyrrole water in the cavity for a long time. In particular, we observe that the pyrrole water moves away from its initial position, and after a while the empty space is filled up with a solvent water. This process is repeated several times during the 45 ns production run. As mentioned in the previous section, the movement of the pyrrole water out of the chromophore cavity significantly alters the chromophore structure. 4. Conclusions In the present work, we could successfully extend the CHARMM force field by developing parameters for phytochomobilin (PΦB), a member of the biologically important family of linear tetrapyrroles. The set of parameters reproduces in an accurate manner the structural properties, rotational barriers and intermolecular interactions of the chromophore as compared to ab initio target data. In addition, the quality of the new force field parameters was tested by performing molecular dynamics simulations of the entire bacteriophytochrome DrCDB. As a result, both structural (average structure) and dynamic properties (root-mean-squarefluctuations) of the chromophore are in satisfactory agreement with the available experimental data. In contrast to the X-ray structure where the chromophore’s rings A, B and C lie in a coplanar arrangement, the MD simulations favors a twisted conformation of the A-B methine bridge. Furthermore, the MD simulations reveal the existence of structural heterogeneity of PΦB chromophore in DrCBD. We could identify two conformational states of the chromophore characterized by distinct torsional angles around the C-C single bond at the A-B methine bridge: 13.8° for state 1 and 18.5° for state 2. The formation of these two states is associated with the formation of hydrogen bonds between the residue HIS260 and the PΦB cofactor, as well as with the presence/ absence or the pyrrole water in the chromophore cavity. The existence of conformational heterogeneity involving ring A is in qualitative agreement with NMR and RR spectroscopic data. A statistical analysis of H-bonding events was performed to study the dynamical behavior of H-bonds in the chromophore binding pocket. Among the ten H-bond donor and acceptor pairs considered in our analysis, those formed between the PΦB chromophore and HIS260, ARG254 and SER274 residues are in particular very stable. These H-bonds are responsible for the rigidity of the propionic side chains at rings B and C and may therefore of relevance for understanding the signal transduction

958 J. Phys. Chem. B, Vol. 113, No. 4, 2009 mechanism between chromophore and protein at an intermediate stage of the photocycle. Finally, the MD simulations allowed us to identify a stable water network comprising at least eight water molecules in the close vicinity of the chromophore. The presence of this water network has a strong influence on the chromophore’s structural properties. In particular, the propionic chains and ring D gain rigidity from this network. Additional Information The parameter and topology files for PΦB are available from the authors upon request. Acknowledgment. We acknowledge the help from the Konrad-Zuse-Zentrum in Berlin for supplying us with computational resources and we would like to thank Dr. Bernd Kallies for technical support, as well as Prof. Peter Hildebrandt for helpful discussions. The financial support by the Deutsche Forschungsgemeinschaft, SfB498 to M.A.M and G.D. and the financial support by the Investitionsbank Berlin and the Excellence Cluster “Unifying concepts in catalysis” to S.K. and M.A.M. is gratefully acknowledged. Abbreviations RR: resonance Raman DrCB: Deinococcus radiodurans chromophore binding domain PΦB: Phytochromobilin QM: quantum mechanical MM: molecular mechanical MD: molecular dynamics References and Notes (1) Rockwell, N. C.; Su, Y. S.; Lagarias, J. C. Phytochrome structure and signaling mechanisms. Annu. ReV. Plant Biol. 2006, 57, 837–858. (2) Quail, P. H. Phytochrome photosensory signalling networks. Nature ReV. Mol. Cell Biol. 2002, 3 (2), 85–93. (3) Briggs, W. R.; Spudich, J. L. Handbook of Photosensory Receptors; Wiley Verlag: Weinheim, Germany, 2005. (4) Wagner, J. R.; Brunzelle, J. S.; Forest, K. T.; Vierstra, R. D. A light-sensing knot revealed by the structure of the chromophore-binding domain of phytochrome. Nature 2005, 438 (7066), 325–331. (5) Wagner, J. R.; Zhang, J. R.; Brunzelle, J. S.; Vierstra, R. D.; Forest, K. T. High resolution structure of Deinococcus bacteriophytochrome yields new insights into phytochrome architecture and evolution. J. Biol. Chem. 2007, 282 (16), 12298–12309. (6) von Stetten, D.; Gu¨nter, M.; Scheerer, P.; Murgida, D. H.; Mroginski, M. A.; Krauβ, N.; Lamparter, T.; Zhang, J.; Anstrom, D. M.; Vierstra, R. D.; Forest, K. T.; Hildebrandt, P. Chromophore heterogeneity and the photoconversion in phytochrome crystals and solution studied by resonance Raman spectroscopy. Angew. Chem. Int. Ed. 2008, 47, 1–4. (7) Wagner, J. R.; Zhang, J.; von Stetten, D.; Gu¨nter, M.; Murgida, D. H.; Mroginski, M. A.; Walker, J. M.; Forest, K. T.; Hildebrandt, P.; Vierstra, R. D. Mutational analysis of Deinococcus Radiodurans Bacteriophytochrome reveals key amino acids necessary for structural integrity and the proton exchange cycle during phoroconversion. J. Biol. Chem. 2008, 283, 12212–12226. (8) Schmidt, P.; Gensch, T.; Remberg, A.; Gartner, W.; Braslavsky, S. E.; Schaffner, K. The complexity of the P-r to P-fr phototransformation kinetics is an intrinsic property of native phytochrome. Photochem. Photobiol. 1998, 68 (5), 754–761. (9) van Thor, J. J.; Mackeen, M.; Kuprov, I.; Dwek, R. A.; Wormald, M. R. Chromophore structure in the photocycle of the cyanobacterial phytochrome Cph1. Biophys. J. 2006, 91 (5), 1811–1822.

Kaminski et al. (10) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102 (18), 3586–3616. (11) Foloppe, N.; MacKerell, A. D. All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 2000, 21 (2), 86–104. (12) Reuter, N.; Lin, R.; Thiel, W. Green fluorescent proteins: Empirical force field for the neutral and deprotonated forms of the chromophore. Molecular dynamics simulation’s of the wild type and S65T mutant. J. Phys. Chem. B 2002, 106 (24), 6310–6321. (13) Olkhova, E.; Hutter, M. C.; Lill, M. A.; Helms, V.; Michel, H. Dynamic water networks in cytochrome c oxidase from Paracoccus denitrificans investigated by molecular dynamics simulations. Biophys. J. 2004, 86 (4), 1873–1889. (14) Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (15) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. Charmm - A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4 (2), 187–217. (16) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781–1802. (17) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926–935. (18) Becke, A. D. A New Mixing of Hartree-Fock and Local DensityFunctional Theories. J. Chem. Phys. 1993, 98 (2), 1372–1377. (19) Breneman, C. M.; Wiberg, K. B. Determining Atom-Centered Monopoles from Molecular Electrostatic Potentials - the Need for High Sampling Density in Formamide Conformational-Analysis. J. Comput. Chem. 1990, 11 (3), 361–373. (20) Autenrieth, F.; Tajkhorshid, E.; Baudry, J.; Luthey-Schulten, Z. Classical force field parameters for the heme prosthetic group of cytochrome c. J. Comput. Chem. 2004, 25 (13), 1613–1622. (21) Brunger, A. T.; Karplus, M. Polar Hydrogen Positions in Proteins - Empirical Energy Placement and Neutron-Diffraction Comparison. ProteinsStruct. Funct. Genet. 1988, 4 (2), 148–156. (22) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. ConstantPressure Molecular-Dynamics Simulation - the Langevin Piston Method. J. Chem. Phys. 1995, 103 (11), 4613–4621. (23) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101 (5), 4177–4189. (24) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - An N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089–10092. (25) van Gunsteren, W. F.; Berendsen, H. J. C. Algorithms for Macromolecular Dynamics and Constraint Dynamics, 34th ed.; 1977; pp 1311-1327. (26) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14 (1), 33. (27) Willis, B. T. M.; Prior, A. W. Thermal Vibrations in Crystallography; Cambridge University Press: New York, 1975. (28) Brooks, C. L.; Karplus, M.; Pettitt, B. M. A. Theoretical Perspective Of Dynamics, Structure, and Thermodynamics. AdV. Chem. Phys. 1988, 71, 1–200. (29) Essen, L. O.; Hughes, J.; Mailliet, J.; The structure of a complete phytochrome sensory module in the Pr ground stateProc. Natl. Acad. Sci. U.S.A. 2008, 105, 14709–14714. (30) Yang, X.; Stojkovic, E. A.; Kuk, J.; Moffatt, K. Crystal structure of the chromophore binding domain of an unusual bacteriophytochrome, RpBphP3, reveals residues that modulate photoconversion. Proc. Natl. Acad. Sci. U.S.A. 2007, 104 (30), 12571–12576. (31) Deloof, H.; Nilsson, L.; Rigler, R. Molecular-Dynamics Simulation of Galanin in Aqueous and Nonaqueous Solution. J. Am. Chem. Soc. 1992, 114 (11), 4028–4035.

JP8047532

Molecular Dynamics Simulations of the Chromophore Binding Site of

Recommend Documents