Letter pubs.acs.org/JPCL
Alternative Electron-Transfer Channels Ensure Ultrafast Deactivation of Light-Induced Excited States in Riboflavin Binding Protein Laura Zanetti-Polzi,† Massimiliano Aschi,† Andrea Amadei,*,‡ and Isabella Daidone*,† †
Department of Physical and Chemical Sciences, University of L’Aquila, via Vetoio (Coppito 1), 67010 L’Aquila, Italy Department of Chemical and Technological Sciences, University of Rome “Tor Vergata”, Via della Ricerca Scientifica, 00185 Rome, Italy
‡
S Supporting Information *
ABSTRACT: Flavoproteins, containing flavin chromophores, are enzymes capable of transferring electrons at very high speeds. The ultrafast photoinduced electron-transfer (ET) kinetics of riboflavin binding protein to the excited riboflavin was studied by femtosecond spectroscopy and found to occur within a few hundred femtoseconds [Zhong and Zewail, Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 11867−11872]. This ultrafast kinetics was attributed to the presence of two aromatic rings that could transfer the electron to riboflavin: the side chains of tryptophan 156 and tyrosine 75. However, the underlying ET mechanism remained unclear. Here, using a hybrid quantum mechanical−molecular dynamics approach, we perform ET dynamics simulations taking into account the motion of the protein and the solvent upon ET. This approach reveals that ET occurs via a major reaction channel involving tyrosine 75 (83%) and a minor one involving tryptophan 156 (17%). We also show that the protein environment is designed to ensure the fast quenching of the riboflavin excited state.
U
ltrafast electron-transfer (ET) dynamics in proteins play a crucial role in various biological systems, in particular in enzymes with redox reactions.1−3 Flavoproteins, which contain flavin-based chromophores, are examples of such enzymes and are involved in various catalytic processes.4−6 Although most of the reactions involving flavins in biological systems are not light-driven, photoinduced ET reactions of flavins in flavoproteins have been the subject of both experimental and theoretical investigations7−9 also because they provide useful models for the study of ultrafast primary photoreactions in biological systems, such as in photosynthetic reaction centers or photosensory proteins.10−13 Moreover, fast quenching of lightinduced excited states of flavins stored in flavoproteins has been proposed as a strategy to suppress uncontrolled flavin photochemical reactivity.14 The ultrafast photoinduced ET of riboflavin binding protein (RBP, see Figure 1A) to the bound cofactor riboflavin (Rbf, vitamin B2) has been previously studied by femtosecond spectroscopy.15−17 The fluorescence quenching of RBP was observed to occur on a femtosecond time scale by measurements of the excited-state flavin (Rbf*) fluorescence decay.15,16 Femtosecond−picosecond transient absorption measurements16 on RBP confirmed that the ET state formation takes place in the regime of hundreds of femtoseconds, corresponding to the time scale of the ultrafast fluorescence quenching dynamics. Ultrafast ET dynamics in RBP, which is the fastest ET process occurring in flavoproteins,7,15−17 were attributed to the tightly bound structure of the ligand into the binding site of RBF. In fact, the isoalloxazine ring of Rbf is stacked between © 2017 American Chemical Society
Figure 1. (A) Rbf binding site in RBP. The two aromatic residues in the binding site are highlighted: Trp156 in gray, and Tyr75 in green. (B) Representative structure of riboflavin (Rbf). The two moieties, the isoalloxazine ring and the ribityl tail, are highlighted.
the parallel planes of the side chains of a tryptophan (Trp156) and a tyrosine (Tyr75) with crystallographic distances between the side chains centers of mass and the isoalloxazine moiety center of mass of 0.53 and 0.35 nm, respectively18 (see Figure 1B). In particular, it was proposed on the basis of the oxidation potentials of indole (the side chain of Trp) and phenol (the side chain of Tyr)19 that the ET to the excited flavin occurs from the side chain of Trp156, rather than from the side chain of Tyr75, by a one-electron half-reduction reaction of the flavin. Received: June 20, 2017 Accepted: June 30, 2017 Published: June 30, 2017 3321
DOI: 10.1021/acs.jpclett.7b01575 J. Phys. Chem. Lett. 2017, 8, 3321−3327
Letter
The Journal of Physical Chemistry Letters
MD simulation, and their effect on the QC electronic properties (e.g., energies) are calculated. At each frame of the MD simulation the perturbed electronic Hamiltonian matrix (H̃ ) is constructed:
However, no direct evidence was provided ruling out reduction by tyrosine.7,20 The interpretation of the experimental results, and the elucidation of the underlying mechanisms, must rely on theoretical−computational approaches. Nevertheless, explicit modeling of ET reactions in a complex protein−solvent environment still represents a challenge from a theoretical− computational point of view because the involved time scales (picosecond−nanosecond) and length scales (nanometer) are well beyond the scales usually attainable by the most accurate ab initio dynamical theoretical approaches.21−23 In particular, most of the difficulties rising when dealing with very complex systems mainly depend on the following aspects: (i) the intrinsic complexity of the problem at the electronic level; (ii) the demand of using multiple trajectories to reconstruct the reactive ensemble in phase space (including both the chromophores and the environment), providing the appropriate statistical information to model the kinetics; (iii) the need for a correct and physically consistent determination of the initial conditions to be used for determining the timedependent observables. For these reasons, a complete and exhaustive model for this kind of problem needs to involve both high-level electronic calculations and a statistical− mechanical perspective typical of large molecular systems, based on atomistic molecular dynamics (MD) simulations.3,24 The aim of the present work is to provide direct evidence of the mechanism underlying the experimental results on the lightinduced ET in RBP. By means of a hybrid quantum/classical theoretical−computational approach, similar in spirit to other hybrid methods,3,21,24−30 based on MD simulations and the perturbed matrix method (PMM), hereafter termed MD-PMM, we are able to explicitly simulate the kinetics of the ET to the excited Rbf and to provide a detailed picture of the underlying mechanisms. The MD-PMM has been successfully applied for the modeling and calculation of, for example, amide I infrared IR spectra,31 ultraviolet−visible (UV−vis) absorption spectra,32 and reduction potentials in proteins.33,34 It has also been applied to model ET reactions in solution35−37 and is here for the first time applied to calculate the ET kinetics inside a protein. In the MD-PMM approach,35−37 the part of the system where the quantum processes of interest occur, the quantum center (QC), is treated at the electronic level, and the rest of the system is modeled as an atomic−molecular semiclassical subsystem exerting an electrostatic effect on the QC electronic states. The main difference with the other hybrid methods is that in the MD-PMM the whole system (including the QC) phase space is sampled by classical MD simulations based on typical atomistic empirical/semiempirical force−fields. The electrostatic perturbation of the environment on the electronic properties of the QC by the environment is included a posteriori: the electronic properties of the isolated QC (unperturbed properties) are calculated quantum-chemically in vacuum (i.e., in the gas phase) and then, for each configuration generated by all-atoms classical MD simulations of the whole system, the electrostatic effect of the instantaneous atomistic configurations of the environment is included as a perturbing term to the electronic properties of the isolated QC providing the electronic properties of the QC embedded in the perturbing environment. To include the electrostatic perturbation of the environment on the QC electronic states, the electrostatic potential and electric field that each atom of the environment exerts on the QC is evaluated at each frame of the
H̃ ≃ H̃ 0 + Iq̃ T = + Z1̃ + ΔVI ̃
(1)
[Z1̃ ]j , j ′ = −E·⟨ϕj0|μ̂ |ϕj0′ ⟩
(2)
with qT and μ̂ the QC total charge and the dipole operator, respectively; = and E are the electrostatic potential and field, respectively, exerted by the environment on the QC center of mass at each frame of the simulation. ΔV includes all the other terms treated as a simple short-range potential, and I ̃ is the identity matrix. The diagonalization of H̃ at each MD frame provides the QC quantum properties in the presence of the perturbing environment. For the specific task of modeling the ET reaction in RBP, the isoalloxazine moiety of Rbf, the side chain of Trp156, and the side chain of Tyr75 are selected as QCs. The diagonalization of the corresponding perturbed electronic Hamiltonian matrices (eq 1) at each MD frame allows the calculation at each MD configuration of the energy variation upon ET, i.e., the time evolution of the transition energy. More details on the theoretical background of the MD-PMM can be found in the Supporting Information. In the ultrafast spectroscopic experiment of Zhong and Zewail,16 after the riboflavin was photoexcited with a pump pulse fixed at ≈480 nm, the femtosecond-resolved fluorescence transient was followed to monitor the decay of the Rbf excited state, that is the result of the ET to the isoalloxazine ring from neighbor aromatic residues. To model with the MD-PMM approach the charge separation (CS) process and to investigate its kinetics, the experimental excited-state ensemble has to be properly reproduced by the calculations, i.e., the whole system configurations to be used as starting points for the modeling of the ET process have to be properly selected. To this aim, the absorption spectrum of Rbf in RBP is here calculated with the MD-PMM methodology (see the Supporting Information). The calculated spectrum, shown in Figure S3, is in good agreement with the experimental spectrum and is used to select the initial configurations for the CS process. In the experiment of Zhong and Zewail, the pump pulse was fixed at 480 nm, i.e., at the red shoulder of the steady-state absorption spectrum of Rbf after the complexation with apoRBP, in order to minimize the effect of vibrational relaxation.16 As no vibrational effects are included in the present computational approach, the best reproduction of the experimental excited-state ensemble can be achieved by selecting the whole system configurations with energy corresponding to the wavelength of the main peak. Thus, from the continuous trajectory of perturbed energies obtained by applying the MD-PMM approach (see the Supporting Information) to a 20 ns MD simulation (EqMD, see Methods and the Supporting Information), 60 configurations of the whole system were randomly selected with a perturbed frequency corresponding to the main peak wavelength (i.e., with a perturbed energy between 267 and 271 kJ/mol). These configurations were used as starting points for 60 independent MD simulations (Neq-MDs, see Methods) to be analyzed in order to reconstruct the CS kinetics. As already mentioned, two aromatic residues are present in the Rbf binding site: Trp156 and Tyr75. As the oxidation 3322
DOI: 10.1021/acs.jpclett.7b01575 J. Phys. Chem. Lett. 2017, 8, 3321−3327
Letter
The Journal of Physical Chemistry Letters
Figure 2. (A) Black curve: MD-PMM kinetic profile for the charge separation process of eq 3 as obtained by means of the distribution of the time intervals required for the first ET to occur (Δ< = 0.0). Red curve: biexponential fitting of the black curve providing τf = 176 fs and τs = 5.9 ps. (B) Black curve: MD-PMM kinetic profile for the charge separation process of eq 4 as obtained by means of the distribution of the time intervals required for the first ET to occur (Δ< = 0.0). Red curve: biexponential fitting of the black curve providing τf = 213 fs and τs = 1.28 ps. (C) Black curve: MD-PMM kinetic profile for the overall charge separation process (considering eqs 3 and 4 as alternative pathways) as obtained by means of the distribution of the time intervals required for the first ET to occur (Δ< = 0.0). Red curve: biexponential fitting of the black curve providing τf = 120 fs and τs = 530 fs. (D) Comparison between the calculated and experimental decay. Black line: MD-PMM kinetic profile for the overall charge separation process reported in panel C. Red squares: experimental points taken from Zhong and Zewail.16 The intensity of the highest experimental point has been set to 1.
mean lifetimes τf = 1/kf = 176 fs (17%) and τs = 1/kf = 5.9 ps (83%) for the fast and slow components of the CS process, respectively. The experimental transient decay of Zhong and Zewail16 provided a fast and slow decay time constant of ≈100 fs (85%) and ≈700 fs (15%), respectively. Thus, while the calculated fast time constant (176 fs) is in good agreement with the experimental one, the slow time constant is rather overestimated by the calculations (5.9 ps). In addition, the calculated probabilities for the fast and slow components (17% and 83%, respectively) are in disagreement with the experimental values. The disagreement between the experimental and the calculated results motivated us to consider the alternative ET pathway involving Tyr75:
potential for Tyr is higher than that for Trp, the experimentally observed ultrafast quenching dynamics has been hypothesized to dominantly result from the ET between the isoalloxazine ring and the aromatic side chain of Trp156 in the protein to form anionic Rbf and cationic Trp:16 Rbf* + Trp → Rbf − + Trp+
(3)
Along the Neq-MDs, the transition energy of this CS process is calculated with the MD-PMM approach (see the Supporting Information for more details). To estimate the kinetics of the CS process according to the adiabatic approximation (i.e., assuming instantaneous relaxation of the QC quantum state to the Hamiltonian eigenstate), we extracted the distribution of the time lengths needed to reach the first crossing of the diabatic energies (i.e., the energies of the QC quantum states approximating the Hamiltonian eigenstates in the reactants and products condition, respectively) from the Neq-MDs, i.e., the time necessary for each of the 60 transition-energy trajectories to change its sign (see Figure S1, in which an example of the time course of the transition energy in which an ET event occurs is shown). Note that, to properly compare the calculated and the experimental results, we excluded from the calculation the transitions occurring in less than 50 fs, which we have estimated to be the experimental temporal resolution. The obtained kinetic profile, reported in Figure 2A, is fitted by a biexponential function, as done in the experiments,16 providing the rate constants of the process, kf and ks, corresponding to
Rbf* + Tyr → Rbf − + Tyr +
(4)
In agreement with the above-mentioned remark of Zhong and Zewail on the oxidation potential for Tyr and Trp, our quantum calculations in vacuum provide a higher oxidation potential for Tyr than for Trp by ≈80 kJ/mol. However, the calculation of the free energy difference ΔAET of the ET reaction involving Tyr and Trp in the protein−flavin environment provides a different result. In fact, by calculating ΔAET from EqMD as35 ΔAET = −kBT ln⟨e−β Δ< ⟩ 3323
(5) DOI: 10.1021/acs.jpclett.7b01575 J. Phys. Chem. Lett. 2017, 8, 3321−3327
Letter
The Journal of Physical Chemistry Letters
show a relevant number of ET events, confirming that the CS reaction involving Tyr75 is a possible ET channel, and the kinetic profile shows that it is also a very fast channel. The biexponential fit of such a profile provides indeed mean lifetimes τf = 213 fs (78%) and τs = 1.28 ps (22%) for the fast and slow components of the CS process. The probabilty of the fast component is considerably higher for this ET reaction than for the one involving Trp156 and results are in very good agreement with the experimental estimate (85%), suggesting that the CS reaction involving Tyr75 might be the dominant ET process in determining the experimental decay. However, both the fast and slow components are slightly overestimated with respect to the corresponding experimental values. According to the transition energies calculated via the MDPMM approach, both the CS process involving Tyr75 and the one involving Trp156 are possible channels for the ET reaction. Thus, they should be treated as alternative ET pathways for each of the Neq-MDs. By monitoring for each trajectory the time necessary to reach the first crossing of the diabatic energies for both channels, it is possible to determine the faster ET channel for each trajectory and thus to reconstruct the kinetic profile of the overall reaction involving both Tyr75 and Trp156 as partners for the CS reaction (see Figure 2C). Such an analysis confirms that the channel involving Tyr75 is the dominant pathway (83%) and provides mean lifetimes τf = 120 fs (55%) and τs = 530 fs (45%) for the fast and slow components of the CS process. Both the fast and slow component mean lifetimes are in very good agreement with the corresponding experimental estimates (≈100 fs and ≈700 fs), and the relative probabilities are well reproduced. Even more remarkable is the excellent agreement between the calculated and experimental kinetic decay (see Figure 2D). Thus, according to our calculations, the experimentally observed decay is the result of two alternative CS channels: one involving Tyr75 and the other involving Trp156. To deepen the molecular mechanism at the basis of the CS process we also analyzed the perturbation exerted by each protein residue and by the solvent on the QCs. Such a perturbation determines the time evolution of the transition energy and thus the crossings of the two diabatic surfaces. Defining Δ=Tyr/Trp = =Tyr/Trp − =Rbf , i.e., the difference between the electrostatic potential exerted by the perturbing environment on Tyr/Trp (=Tyr/Trp) and on Rbf (=Rbf ), we analyzed the electrostatic potential difference Δ=Tyr/Trp exerted by each protein residue and by the solvent at the beginning of the 60 Neq-MDs (first 10 fs) and for transition energies around 0.0 along the Neq-MDs, i.e., for configurations near the crossing of the two diabatic surfaces. The difference ΔΔ=Tyr/Trp = Δ=Tyr/Trp(t = 0) − Δ=Tyr/Trp(Δ< = 0) between the contribution of each residue at the starting configurations and at the crossing of the two diabatic surfaces unveils which are the most determinant residues involved in the ET reaction: for both Tyr 75 (see Figure 4A) and Trp 156 (see Figure 4B), only a limited number of residues shows a relevant difference between the initial and crossing configurations. This result, together with the previous observation that neither the distance nor the relative orientation between the ET partners relevantly vary in the crossing configurations, suggests that the protein scaffold is engineered to drive a very fast ET reaction. This is also supported by the observation that for the faster and dominant reaction channel involving Tyr (panel A) the number
with Δ< the transition energy associated with the ET reactions of eqs 3 and 4 calculated with the MD-PMM approach, we obtain ΔAET ≈ −45 kJ/mol for reaction 3 and ΔAET ≈ −48 kJ/ mol for reaction 4, i.e., the free energy difference is essentially the same for the two ET reactions, which are both thermodynamically favored. In order to understand why the protein environment has a stronger stabilizing effect on the transition energy in the ET reaction involving Tyr than in the one involving Trp, thus compensating the gas-phase oxidation potential difference, we analyzed the contribution of the electrostatic potential exerted by the environment (protein and solvent) to the transition energies of both ET reactions (the electrostatic potential is indeed the main determinant to the transition energy shift between the gas phase and the condensed phase, see eq 1). In particular, the single-residue contribution to the electrostatic potential was calculated to elucidate which residues of the environment determine such a shift. The single-residue analysis is reported in Figure 3 in which the electrostatic potential
Figure 3. Electrostatic potential contribution to the transition energy, Δ
