Protonation of the Biliverdin IXα Chromophore in the Red and Far-Red Photoactive States of a Bacteriophytochrome.

The tetrapyrrole chromophore biliverdin IXα (BV) in the bacteriophytochrome from Deinococcus radiodurans (DrBphP) is usually assumed to be fully protonated, but this assumption has not been systematically validated by experiments or extensive computations. Here, we use force field molecular dynamics simulations and quantum mechanics/molecular mechanics calculations with density functional theory and XMCQDPT2 methods to investigate the effect of the five most probable protonation forms of BV on structural stability, binding pocket interactions, and absorption spectra in the two photochromic states of DrBphP. While agreement with X-ray structural data and measured UV/vis spectra suggest that in both states the protonated form of the chromophore dominates, we also find that a minor population with a deprotonated D-ring could contribute to the red-shifted tail in the absorption spectra.

Introduction

Light-sensitive signaling cascades in various plants, fungi, and bacteria are controlled by a photochemical mechanism involving photoactive proteins.[1,2] Phytochromes belong to the class of photoreceptors that exist in nature with the ability to reversibly interconvert on the timescale of μs to ms between two structurally distinct states: Pr-dark and Pfr-illuminated. These states respond to red light (Pr) and far-red light (Pfr) and regulate numerous cellular functions, such as photomorphogenesis, flowering, shade avoidance in plants, and pigment synthesis in bacteria.[3,4] The red/far-red light family of phytochromes exist as dimers. Upon red light (700 nm) absorption, the Pr state undergoes large conformational changes to photoconvert via various intermediate states into the Pfr state. The reverse photoconversion to Pr from the Pfr state is initiated either by far-red light (750 nm) absorption or by dark thermal relaxation.[5] There are also bacterial phytochromes, for which the Pr-to-Pfr photoconversion occurs in darkness.[6,7] These phytochromes with the reversed dark thermal relaxation mechanism are termed as bathy phytochromes.

Recently, X-ray structures of the photosensory unit of phytochromes from Deinococcus radiodurans (DrBphP) have become available,[5,8,9] which are of great help in beginning to understand the molecular mechanism responsible for the photoresponse of this protein. The structural unit of the core module consists of PAS (period, aryl hydrocarbon receptor nuclear translocator, and single-minded protein), GAF (cGMP-specific phosphodiesterase, adenylate cyclase, and Fhl proteins), and PHY (photochromic phytochrome) domains, while the full length phytochrome has also a histidine kinase (HK) domain attached to the PHY domain at the C-terminal end. The phytochromes can be classified based on the type and site of the covalently bound tetrapyrrole bilin chromophore that varies between organisms: phycocyanobilin (PΦB) in cyanobacterial phytochromes and phytochromobilin (PCB) in plant phytochromes have a covalent bond between the C31 atom (CAA atom in Figure a) and cysteine residue Cys259 in the GAF domain,[4,10,11] while biliverdin IXα (BV) found in bacteriophytochromes is bound to the protein via the C32 carbon atom (CBA atom in Figure a) to cysteine residue Cys24 in the PAS domain.[8,12]

Figure 1

Schematic representation of the BV IXα chromophore in the dark (left) and the illuminated state (right) phytochrome found in D. radiodurans (DrBphP). The atom naming scheme in (a) is used throughout the text. The figures were generated using MolView v2.2.

In DrBphP, the chromophore is located in the PAS-GAF domain cavity and bound to the protein via a thioether linkage with Cys24 and a hydrogen-bonding network involving residues Asp207, Arg254, His260, Ser272, Ser274, and His290 (the residue numbering is similar to the PDB structures 4O01 and 4O0P).[5] In the red light-absorbing Pr state of DrBphP, the chromophore adopts a ZZZssa geometry, while in the far-red light-absorbing Pfr state, it adopts a ZZEssa geometry with an E, anticonformation for the isomerized D-ring (Figure b). Because of these differences in the chromophore conformation, it is widely assumed that photon absorption triggers photoisomerization of the D-ring along the CHD=C1D methyl bridge.[4,13,14] In contrast, nonadiabatic QM dynamics simulations[15] and calculations with different quantum chemical methods suggest that without protein environment, the BV isomerization can occur around the CHC=C1C methyl bridge between the B and C rings.[16,17]

Biliverdin photoisomerization is followed by a series of structural changes in the chromophore binding pocket that ultimately lead to larger structural changes in the protein complex. These major changes include: (i) a secondary structure transformation from an anti-parallel β-sheet to an α-helix in the tongue region (residue range 461–490) between the GAF and PHY domains; and (ii) a large opening of more than ∼3 nm between the two PHY domains of the dimer in the Pfr state (Figure b).[5] The light-induced changes in the PHY domains are speculated to reposition the HK domain as well as the phosphate-accepting histidine residues in the ATP-binding site within that HK domain.[18] These structural rearrangements would eventually promote a phosphotransfer process within the phytochrome dimer and thereby regulate the downstream signal transmission.[19]

Figure 2

Cartoon representation for the crystal structure of PAS-GAF-PHY domains of DrBphP in the (a) Pr-dark state (PDB: 4O0P(5)) that absorbs red light (700 nm) to trigger the secondary structure rearrangement from β-sheets to α-helix (orange) and a PHY domain opening to form the (b) Pfr-illuminated state (PDB: 4O01(5)) that absorbs far-red light (750 nm) to revert back to the Pr state; (c,d) experimental UV/vis absorption spectrum reported by Takala et al., which exhibits the Q-band and Soret band for Pr and Pfr states, respectively.[5]

The Pr and Pfr X-ray structures (Figure a,b) provide glimpses of structural changes that follow photoexcitation but do not address the crucial question of how the protein achieves high signaling efficiency. In addition, there is as yet no consensus about the precise sequence of events in what is essentially a dynamic process. To better understand the molecular mechanisms that drive this reversible photoactivation pathway, experimental data should be complemented with excited state molecular dynamics (MD) simulation, as in previous work on photoreceptors.[20−22] However, performing such simulations requires an accurate structural model with detailed information about the electrostatic interactions in the chromophore binding pocket, and hence the protonation form of BV.

Various authors have proposed a fully protonated BV in Pr and Pfr states based on H-bond analysis from X-ray crystallography data,[8,9] NMR experiments,[23,24] and RR spectra.[25,26] In contrast, others have shown that the Pr → Pfr photoactivation includes proton release/uptake[9,27−33] by the chromophore in the intermediate states with a higher amplitude for the proton release than uptake.[27,31−33] Although resonance Raman studies that are sensitive to N–H in-plane vibrations in BV have suggested that the C pyrrole ring nitrogen remains protonated in both Pr and Pfr states of plant and cyanobacterial phytochromes,[25,34,35] ambiguity remains as to which of the other rings could be involved in these protonation changes. Nevertheless, these observations suggest that in addition to the photoisomerization, protonation changes may be relevant for photoactivation. Furthermore, because phytochromes from different species vary in (i) their chromophore linkage sites, (ii) orientation of pyrrole rings A/D in both Pr and Pfr state, and (iii) isomerization direction of the CHD=C1D double bond,[36] it remains challenging to reconcile the aforementioned experimental findings into a consistent picture that is valid for phytochromes in general and DrBphP in particular.

To resolve the chromophore protonation state for DrBphP and prepare an accurate atomistic model for computationally costly excited state MD simulations,[37] we have used classical MD simulations in combination with quantum chemistry and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations to systematically investigate the effect of the chromophore protonation state on the structural and spectroscopic characteristics of the DrBphP photosensory domain in both Pr and Pfr conformations. Only the five most probable protonation patterns were investigated:

all nitrogen atoms protonated (all-prot);

nitrogen in ring A deprotonated and rings B, C, and D protonated (deprot-A);

nitrogen in ring B deprotonated and rings A, C, and D protonated (deprot-B);

nitrogen in ring C deprotonated and rings A, B, and D protonated (deprot-C);

nitrogen in ring D deprotonated and rings A, B, and C protonated (deprot-D).

To the best of our knowledge, this is the first study to incorporate the PHY domain along with CBD in MD simulations and QM/MM calculations of DrBphP, aimed at systematically exploring the effects of the chromophore protonation state on the stability and absorptivity of this protein in solution.

In addition, we also assessed the accuracy with which various density functional theory (DFT) functionals predict the absorption maxima by comparing the DFT spectra to spectra computed at the correlated multiconfigurational level of ab initio theory. In line with the assumptions in previous computational studies,[37−44] our results support a fully protonated biliverdin chromophore in both Pr and Pfr states of DrBphP, but also suggest that a very minor population, in which the D-ring is deprotonated, might be responsible for a long red tail in the absorption spectrum of this protein.

Methods

Protein Structures

The atomic coordinates of the core photosensory module (PAS-GAF-PHY: 523 amino acids) in dark and illuminated state of the photoreceptor protein reported by Takala et al., were taken from PDB models 4O0P and 4O01, with resolutions of 3.80 and 3.24 Å, respectively.[5] Residues missing from these refined X-ray structures (see Table S1 in the Supporting Information for an overview of missing residues) were added based on the ab initio loop-modeling approach in the MODELLER program[45,46] that finds the optimal structural model based on a potential energy function. We used the Dowser program to place internal water molecules inside hydrophobic cavities of the protein using protein–water and water–water interaction energy algorithms.[47] With a probe radius of 0.2 Å and new DOWSER parameters for the BV chromophore molecule (Table S2), water molecules were inserted inside protein cavities, in which the interaction energy of the water was below a threshold value of −12 kcal/mol (Figure S1).

The protonation states and tautomeric forms of amino acids with ionizable side chains were determined by pKa estimation with H++ Server[48] and PROPKA 3.1[49,50] at pH = 7. The rotameric states of these residues were inspected and modified to a new rotamer if that would increase the number of hydrogen bonds. For amino acids with nonionizable side chains, the H-atoms were added using force field parameters corresponding to their protonation form at pH = 7. As in previous works,[38−41] we assume the same protonation state of the protein residues for all protonation states of the chromophore. While a change in protonation of the chromophore might affect the pKa’s of the surrounding protein residues, the number of possible states that would need to be sampled is too large to handle with standard MD simulations and was therefore not attempted here. The structural protein models of the Pr and Pfr states after equilibration are available in the Supporting Information (Pr.pdb and Pfr.pdb).

Force Field Model

The AMBER03 force field was used to model interactions between the atoms in our simulations.[51] To model the intra- and intermolecular interactions of the chromophore, we derived force field parameters by means of ab initio calculations for a QM model that consisted of BV and the covalently attached side chain of Cys24. For all protonation states, the geometry of the QM model was optimized in the gas phase at the CAMB3LYP/6-31G(dp) level of DFT theory. The partial charges for the different protonation forms of BV were estimated by performing a two-stage restricted electrostatic potential (RESP) fit to the electron density obtained at the HF/6-31G(d) level of theory, using constraints on symmetry-related atoms in the second step (Tables S4 and S5).[52] The atomic charges on the main chain of Cys24 were described following the approach proposed by Altun and co-workers (Tables S6 and S7).[53] Additional atom types were created specifically for the chromophore to define new dihedrals and improve structural description. The complete force field for use with the Gromacs MD program[54−56] is available in the Supporting Information. The DFT and ab initio calculations for the RESP fit were carried out with the Gaussian09 program,[57] while the actual RESP fittings were performed with the AmberTools15 package.[58]

Classical MD Simulation

For MD simulations, the protein was solvated in a cubic periodic box (Pr: 12 × 12 × 12 nm and Pfr: 13 × 13 × 13 nm), which was filled with TIP3P water molecules.[59] The large ∼3 nm opening of the PHY domain in the Pfr state necessitated the larger box size for this system. The total charge on the protein systems was neutralized by replacing solvent water molecules with Na+ and Cl– ions until a physiological salt concentration of 0.15 M was reached. The neutral solvated model contained 81 water molecules added by DOWSER (Figure S1) and another 50 454 and 68 294 water molecules in the Pr and Pfr models, respectively.

The short-range attractive and repulsive dispersive interactions were described by a Lennard-Jones potential with a cut off of 1.0 nm. Electrostatic interactions were calculated at each time step using the particle mesh Ewald method[60] with a grid spacing of 0.12 nm. The LINCS algorithm was used to constrain bond lengths in the protein,[61] while SETTLE[62] was used to constrain the internal degrees of freedom of the TIP3P water molecules. These constraints allowed us to perform simulations with a time step of 2 fs.

The 10 phytochrome systems were initially subjected to an energy minimization for 10 000 steps using the steepest descent algorithm, followed by multiple equilibration steps: 200 ps simulation with position restrains on all protein atoms and two 500 ps simulations on the position-restrained system using NVT and NPT ensembles for temperature (300 K) and pressure (1 bar) equilibration, respectively. After minimization and equilibration, all phytochrome systems were simulated for 100 ns. All classical MD simulations were performed with the GROMACS 4.6.5 MD package.[54−56]

Spectra Calculations

For each chromophore protonation state, a snapshot was taken at every nanosecond from the 100 ns trajectories. Single-point vertical excitation energies were calculated for these 100 structures at the QM/MM level and superposed into a spectrum.[86,87] For these QMMM calculations, a smaller BV chromophore QM subsystem (QM-small, Figure S3) was selected, in which the methyl groups of rings A, B, and D and the propionate side chains were replaced by H atoms, and the Cys24 attachment was excluded by truncating the chromophore system at the CAA-CBA bond (Figure a). The rest of the system was described at the AMBER03 force field level.[51] In these calculations, electrostatic interactions with the rest of the chromophore, protein, and solvent were added to account for the polarization of the QM region by the environment. The choice for a QM subsystem that is smaller than the actual chromophore was validated by comparing the spectra computed for this small QM subsystem to those computed with the complete chromophore plus pyrrole water molecule inside the QM subsystem (Supporting Information, Figure S6).

The excitation energies (ΔE) of the snapshots i were computed both at the level of time-dependent DFT (TD-DFT)[63−69] and at the correlated multiconfigurational level of ab initio theory with the extended multiconfiguration quasidegenerate perturbation theory (XMCQDPT2) method.[70] For the TD-DFT method, we used various functionals (BLYP,[71,72] B3LYP,[73−75] PBE0,[76] and CAMB3LYP[77]) and the 6-31G(d,p) basis set. In the XMCQDPT2 computations, a complete active space self-consistent field (CASSCF) wave function[78] based on an active space of 12 electrons in 12 orbitals, averaged over first 5 singlet states and expanded in the cc-pVDZ basis set,[79] was used as a reference for the second-order perturbation theory calculation (i.e., XMCQDPT2/SA(5)-CASSCF(12,12)/cc-pVDZ). The active space orbitals used in these calculations are shown in Figure S10 of the Supporting Information. The selection of this active space was a trade-off between accuracy and computational efficiency and validated by computing energy profiles for the most relevant chromophore distortions (Figure S11).

The excitation energies were converted into a spectrum by convoluting the energy gaps with Gaussianswhere I is the intensity as a function of excitation energy (E), me is the electron mass, M is the number of snapshots included in the analysis, ΔE is the excitation energy in snapshot i, and μ is the transition dipole moment of that excitation. A width of σ = 0.02 eV was chosen for the convolution. The QM calculations were performed with the Gaussian09 package[57] (DFT and HF) and Firefly[80] (XMCQDPT2).

We note that a similar approach has recently been used by Polyakov et al.[42] However, the major difference between their work and ours is that we have computed absorption spectra of ensembles of the protein in solution at room temperature which can account for possible heterogeneities that may affect the interpretation of the spectra. In addition, while they focused on the fully protonated chromophore, we also investigated alternative protonation states.

Results and Discussion

MD Simulation

The stability of the chromophore pocket with different protonation forms of BV was assessed by analyzing the network of non-covalent interactions between the chromophore and the rest of the system in the MD trajectories and comparing these interactions to those in the X-ray structures (Figure ).

Figure 3

Stick representation of the conserved residues and the H-bonding network surrounding BV. Carbon atoms are colored in cyan and orange for BV and protein, respectively, with nitrogen in blue, oxygen in red, and water molecules as red spheres. The hydrogen bonding network for both models of DrBphP is shown as black dashed lines.

Pr—Dark State Simulations

The His260 residue stacked below the plane of the BV B-ring and C-ring is a stable H-bond acceptor for the pyrrole ring N–H atoms in the all-prot and deprot-D models. In contrast, in the deprot-A, deprot-B, and deprot-C models, this interaction breaks off within 25 ns and does not reform again (Figure S4). Furthermore, the absence of a proton on the A, B, or C ring of the chromophore disrupts the electrostatic interaction network involving the BV pyrrole N-atoms, Asp207 backbone carbonyl, and pyrrole water.

The distance distribution of the salt bridge between the Asp207–Arg466 side-chains (Figure S5) shows least disruptions over 100 ns in the all-prot and deprot-D models. In the deprot-A, deprot-B, and deprot-C Pr models, this salt bridge ruptures within the initial 20 ns and does not reform in the remainder of the 100 ns simulation time.

Hydrophobic residues Phe203, Tyr176, and Tyr263 form a pocket around the D-ring of BV, which allows the D-ring to isomerize without steric clashes between the residue side chains and methyl group on the D-ring.[5,8,29] Altering the interactions involving Tyr176 and Tyr263 causes photobleaching in proteins where these amino acids are mutated.[30] In our Pr simulations with A-ring, B-ring, or C-ring deprotonated, the rapid loss of the aforementioned salt bridge between Asp207 and Arg466 also disrupts the hydrogen bond between Asp207 and Tyr263 (Figure S5). This rupture causes Tyr263 to move away from Tyr176 and Phe203 and thus destabilizes the whole hydrophobic cavity. In contrast, the cavity remains intact in simulations where the chromophore is fully protonated (all-prot) and deprotonated at the D-ring (deprot-D).

The chromophore B-ring propionate side-chain is stabilized through a salt bridge with Arg254 and accepts a hydrogen bond from the Tyr216 hydroxyl group. The C-ring propionate side-chain makes polar contacts with S274 and a water-mediated contact with His260.[4,5,8] These propionate side-chain interactions remain intact for all Pr models which suggests a negligible impact of the chromophore protonation on the side-chain interactions (Figure S5).

Pfr—Illuminated State Simulations

In the illuminated state, the residue configuration inside the chromophore binding pocket changes because of (i) photoisomerization, (ii) secondary structure changes in the conserved tongue region in contact with the chromophore, and (iii) a global opening of the PHY domain.[5]

The His260, D207, pyrrole water, and pyrrole N-atoms of BV form a polar contact network similar to the Pr state in all simulations. H260 maintains a stable contact with BV pyrrole-ring N-atoms by either acting as a hydrogen bond acceptor (all-prot and deprot-D) or donor (deprot-A, deprot-B, and deprot-C). However, the conserved interactions of Asp207-BV pyrrole N-atoms, Asp207–Ser468 and Ser468–Tyr263, are lost (Figures S6 and S7) in the Pfr state with BV deprotonated at the A, B, or C ring. Instead of re-establishing these hydrogen bonds later in the simulation, Asp207 moves out of the chromophore-binding pocket to form new interactions with the solvent molecules. In contrast, when the chromophore is fully protonated or deprotonated at the D-ring, the hydrogen bonds involving Asp207, Ser468, and Tyr263 remain intact. These observations support both a fully protonated chromophore and a chromophore with a deprotonated D ring, rather than the other chromophore protonation states, in which the hydrogen bonds are less stable.

Similar to the Pr state simulations, the B-ring and C-ring propionate side-chain interactions with residues Arg254, Tyr216, and Ser274 are stable in Pfr models regardless of the BV protonation state (Figure S7).

On the basis of the analysis of the conserved interactions between the chromophore and the protein in our MD simulations, we inferred that the A-ring, B-ring, and C-ring are protonated in Pr and Pfr but could not determine whether the D-ring is also protonated, as the chromophore binding pocket was stable in both situations (all-prot and D-deprot models). To differentiate further between these states, we also computed absorption spectra that can be directly compared to experiment.

Effect of DFT Functional

In Figure , we compare the UV/vis spectra for Pr with an all-protonated chromophore, calculated at various levels of theory. Whereas the TD-DFT spectra contain the Soret band around 400 nm and the Q-band around 700 nm, the XMCQDPT2 spectrum lacks a clear Soret band. The latter is due to the limited number of states (nstates = 5) in the CASSCF calculations. With 10 states, the Soret band is predicted (Figure S8, Supporting Information), but this is computationally too expensive and also not needed as the Q-band absorption near 700 nm is much more sensitive to the structural changes associated with photoactivation from Pr to Pfr.[4,5,9]

Figure 4

QM/MM UV/vis-absorption spectra of DrBphP evaluated at the TDDFT level with various DFT functionals and at the XMCQDPT2 level of theory in the Pr state. The experimental spectrum is shown in dashed orange lines, with a Q-band peak at 700 nm. The best agreement is obtained with the BLYP functional and XMCQDPT2 method with an absorption maxima at 692 nm (black line) and 708 nm (blue line), respectively.

Compared to the experimental spectrum with a maximum of the Q-band at 700 nm in the Pr state,[5] the TDDFT spectra are blue-shifted, while the XMCQDPT2 spectrum is red-shifted. The hybrid functionals with 20–25% exact exchange, B3LYP and PBE0, deviated by 35 nm (0.09 eV) and 43 nm (0.116 eV), respectively, while the long-range-corrected metahybrid CAMB3LYP functional, which was developed to yield more accurate results for charge-transfer excitations, shows the largest blue shift of 56 nm (0.154 eV). While hybrid functionals and metahybrid functionals normally improve the accuracy of electronic excitations involving Rydberg or charge-transfer states, we here observe that the absorption maximum (692 nm/1.79 eV) obtained with the Generalized Gradient Approximation (GGA) BLYP functional is closest to both the experiment (700 nm/1.77 eV) and the XMCQDPT2 result (708 nm/1.75 eV). The performance of the different DFT functionals for the CBD-PHY dimer is in line with the results of a recent benchmark study on the CBD monomer.[40] Because of the better agreement of the TD-BLYP/6-31G(d) and the XMCQDPT2/cc-pVDZ results for Pr, we only consider these methods in what follows.

TDDFT/MM Spectra of Various Protonation States

On the basis of the analysis of the conserved interactions between the chromophore and the protein in our MD simulations, we inferred that the A-ring, B-ring, and C-ring are protonated, but could not determine whether also the D-ring is protonated, as the chromophore binding pocket was stable in both situations (all-prot and D-deprot models). Therefore, we compare the computed absorption spectra of the five chromophore protonation states to experiments. Figure shows these absorption spectra for both Pr and Pfr protein conformations. While in experiments the photoactivation of DrBphP is accompanied by a red shift in the Q-band absorption of 0.12 eV (from 700 to 750 nm), our TD-DFT/MM calculations show much smaller red shifts (Table ). The largest red shift (0.04 eV) was obtained when the chromophore D-ring is deprotonated, in which case the absorption maximum changes from 714 nm in Pr to 731 nm in Pfr. With the fully protonated chromophore, the red shift is much smaller: only 8 nm or 0.02 eV. The maximum shift, however, is obtained if we assume that somewhere along the Pr to Pfr transition, the D ring becomes deprotonated. In this scenario, the Q-band absorption maximum moves from 688 nm in the Pr state to 731 nm in the Pfr state, which corresponds to a red shift of 0.11 eV (43 nm), in reasonable agreement with the experiment. Previously, such loss of a proton has been suggested to occur during the Pr-to-Pfr conversion in cyanobacterial phytochrome Cph1 and bacterial phytochrome Agp1 on the basis of flash-induced transient absorption measurements and pH dependence of the absorption spectra.[27−30,33,34] On the basis of the MD snapshot in Figure S13, which highlights a proton wire from the chromophore to solvent, we speculate that the proton could leave the pyrrole D-ring of BV via that wire. Alternatively, also His260 might be involved, as it could potentially shuttle the proton onto the pyrrole water molecule, which exchanges frequently with bulk solvent in both Pr and Pfr simulations.

Figure 5

Computed UV/Vis absorption spectra of the photoactive state models Pr and Pfr with different protonation forms of the BV chromophore computed using the TD-DFT method, BLYP/6-31G(d,p).

Table 1

Q-Band and Soret Band Peak Computed from Single-Point Excited State Energies Using BLYP/6-31G(d,p) for Five Different Protonation Forms of BV in Both Pr and Pfr States

	Pr—dark state (nm)		Pfr—illuminated state (nm)
BV protonation form	Q-band	Soret band	Q-band	Soret band
all-prot	696	409	704	476
deprot-A	708	421	720	433
deprot-B	696	439	708	413
deprot-C	712	416	716	424
deprot-D	716	430	742	423
Exp.	700	394	750	409

At this point, the combined MD simulations and TDDFT/MM calculations suggest a protonated chromophore in the Pr state of DrBphP and a chromophore that is deprotonated at the D-ring in the Pfr state. However, because the systematic error in TDDFT of about 0.25 eV[81] is beyond the measurable red-shift, further confirmation is essential. We therefore also computed the absorption spectra at the correlated XMCQDPT2 level of theory, which has a higher accuracy with an error of about 0.1 eV.[70]

XMCQDPT2 Spectra

In previous works, the effect of deprotonation at the A, B, and C rings on the optical response of the chromophore has been computed,[14,37,41,43] but not at the D-ring. Because, we rule out deprotonation of the A, B, and C rings in both the Pr and Pfr states of the photosensory domain based on the MD simulations, we have not computed the excitation spectra for these species at XMCQDPT2 level of theory. The XMCQDPT2/SA (5)-CASSCF(12,12)/cc-pVDZ//Amber03 QM/MM spectra are shown in Figure . While for the fully protonated chromophore, the spectra agree reasonably with the TD-BLYP results (Figures and 6), the spectrum of the protein with a chromophore that is deprotonated at the D-ring, shows additional peaks. In particular, there is a broad optical transition beyond 800 nm that was not seen in the TD-DFT computations. Inspection of the molecular orbitals involved in these excitations reveals that the lowest energy S0–S1 transition in the deprotonated chromophore at 850 nm corresponds to a single-electron highest occupied molecular orbital–lowest unoccupied molecular orbital (HOMO-LUMO) transition, while the excitation around 700 nm is a two-electron or double excitation. Because of the adiabatic approximation in the exchange–correlation response kernel (ALDA), these double excitations are notoriously difficult to describe with TD-DFT, which is why they are missing here. The apparent agreement between the TD-DFT results and experiment must thus have been fortuitous and therefore cannot be used to support that the D-ring of the chromophore is deprotonated in either Pfr or Pr. Consistent with the TD-DFT results, the red shift between the Pr and Pfr states with the fully protonated chromophore is underestimated also at the XMCQDPT2/cc-pVDZ level (0.032 eV), but the deviation with respect to the experimental absorption maximum (0.12 eV) is closer to the anticipated error of the method (∼0.1 eV).[70]

Figure 6

Calculated UV/vis absorption spectra for the fully protonated BV model [(a)-Pr, (d)-Pfr] and model with the D-ring deprotonated [(b)-Pr, (e)-Pfr] of DrBphP at the correlated XMCQDPT2 level of theory. Plots (c,f) represent the unnormalized spectra for the Pr and Pfr states, respectively, where the experimental spectra were scaled to the peak of the fully protonated BV model for better comparison.

However, before concluding that the chromophore is fully protonated in both Pr and Pfr states of the protein, we notice that because of the much lower oscillator strength of the chromophore with the deprotonated D-ring (Figure ), the spectrum would be dominated by absorption of the protonated chromophore, masking the deprotonated species. Furthermore, the spectra available in the literature for this protein were only recorded up until 850 nm and sometimes even baseline corrected.[4−7] Therefore, we speculate that, in particular, the red tail of the absorption might contain a minor contribution from proteins, in which the D-ring is not protonated. An indication in support of our speculation is the spectra of phytochromes with covalently modified biliverdin chromophores, which show a weak but a visible rise in absorption beyond 850 nm.[82−84] Nevertheless, to experimentally verify whether the deprotonated chromophore indeed absorbs beyond 850 nm, we suggest recording the spectra over a larger range of wavelengths and pH values.

Conclusions

We have performed MD simulations of the photosensory domain of the D. radiodurans phytochrome protein with five different biliverdin protonation states in an attempt to determine the most likely protonation state. An analysis of the trajectories in terms of chromophore–protein interactions as well as absorption spectra suggests that the biliverdin chromophore in DrBphP is fully protonated, in line with previous assumptions. However, on the basis of the QM/MM spectra calculations at the correlated level of theory, we cannot rule out a minor population of proteins containing a chromophore with the deprotonated D-ring. Therefore, in our follow-up work, aimed at unraveling the photoisomerization mechanism, we will not only perform excited state dynamics simulations of the protein with a fully protonated chromophore, but also include simulations in which the D-ring of the chromophore is deprotonated. Because we found that the optical S0–S1 transition of the latter has a double excitation character, we will need a suitable multiconfigurational method rather than a TD-DFT model.