Probing microhydration effect on the electronic structure of the GFP chromophore anion: Photoelectron spectroscopy and theoretical investigations

The discovery of green fluorescent proteins (GFPs) by Shimomura¹ has transformed biology through its use as a marker for in vivo monitoring of biological processes,^2–5 as a reporter gene, fusion tag, a cell marker,⁶ and as a monitor gene expression in different organisms.⁷ This has led to the synthesis of new GFPs with emission wavelengths spanning the visible range^8–10 that recently has produced a new class of photoactivatable fluorescent proteins¹¹ for use in ultra-resolution imaging^12,13 and data storage.^14,15

The GFP chromophore is based on p-hydroxybenzylidene-2,3-dimethylimidazolinone (HBDI, Figure 1) that resides in the β-barrel of the protein.¹⁶ There exist two forms of the chromophore, one that is neutral and the other an anionic (deprotonated) form, which is highly fluorescent (quantum yield Φ_fl = 0.79) that is generated by an S₁←S₀ excitation. Accordingly, GFP has two absorption bands peaked at 395 and 480 nm, ascribed to the neutral and deprotonated anionic forms of the chromophore, respectively.¹⁷ Excitation of either form results in green fluorescence at 508 nm, originating from the anionic chromophore as the neutral form deprotonates upon excitation.^17,18

The photophysics of the GFP chromophore is critically dependent on the environment. The fluorescence quantum yield decreases over 1000-fold from the protein (Φ_fl = 0.8) to liquid solutions (Φ_fl < 10⁻³) at room temperatures but becomes strong again in solution when cooled to 77 K.^19–23 It is nonfluorescent in the gas phase.²⁴ To obtain a fundamental understanding of the observed phenomena requires consideration of the electronic structure and excited state dynamics. To ensure reliable description of the structure and underlying electronic structure transformations, such a problem is best approached through a combined experiment and simulation approach. Photoelectron spectroscopy studies have determined the vertical detachment energy (VDE)^25–27 as well as the adiabatic detachment energy (ADE)²⁸ of the HBDI⁻ anion. The gas-phase absorption spectra for both HBDI neutral^29–31 and the HBDI⁻ anion^24,32,33 have been obtained, and compared to the absorptions in solutions and in proteins. The dynamics of the S₁ and higher excited states have been investigated using time-resolved photoelectron and action spectroscopy,^34–36 ultrafast fluorescence spectroscopy,²² and theoretical calculations,^37–40 revealing that the S₁ state evolves from the Franck-Condon (FC) region to the fluorescence state (FS) geometry on an ultrafast time scale, followed by a rotation around the bridging C–C–C bonds between phenol and imidazolinone rings. This leads to formation of a twisted intermediate geometry that subsequently undergoes internal conversion to the anionic ground state S₀. Therefore, the dramatic change in fluorescence quantum yield and excited-state lifetime is attributed to the restricted range of motion of the chromophore in the protein, where the outer β-barrel in GFP provides the chromophore a hydrophobic binding pocket, and keeps the chromophore in a preferred configuration through its hydrogen-bond network.

Despite extensive aforementioned experimental and computational studies, there still remain many open questions regarding the key fundamental variables that govern the GFP chromophore photophysics and details of radiationless relaxation in protein and solutions.^41–46 One outstanding problem is the role of autoionization as a possible relaxation pathway of the excited state under different environmental conditions. Even for the isolated HBDI⁻ anion, the energy level of the excited state S₁ relative to the electron detachment continuum D₀ of the neutral HBDI^• is still a subject of debate. A recent vibrationally resolved photoelectron spectroscopy study reported that D₀ state lies at 2.73 ± 0.01 eV above the ground state of HBDI⁻, S₀.²⁸ Several gas-phase absorption measurements have shown S₁←S₀ peaked at 2.57 eV (482.5 nm),²⁴ 2.59 eV (479 nm),³² essentially identical to the same absorption maximum of the protein in its anionic form, which would indict a bound S₁ while some other studies advocate that the S₁←S₀ peak of the isolated HBDI⁻ is blue-shifted by 0.23 eV relative to that in the protein to 2.84 eV,^33,41,47 which would suggest a unbound S₁ state. Since the fluorescence (back to S₀) and autodetachment (to D₀) both originate from the S₁ state, these two processes are competing with each other, and the relative energy level of S₁ and D₀ is expected to play a crucial role in determining the photophysics of the GFP chromophore.

Another important issue is related to the role of the solvent on photophysical properties of GFP. In the protein environment, crystallographic water molecules are found hydrogen-bonded to the phenolate oxygen and in the immediate vicinity of the chromophore.^16,48 Krylov and co-workers conducted theoretical investigations and found that solvent water molecules significantly increase the detachment energies (D₀←S₀), but only slightly blue-shift the excitation energies (S₁←S₀),^45,46 so that even one single water molecule was predicted to lead to a firmly bound excited state S₁. In this work, we performed negative ion photoelectron spectroscopy (NIPES) and ab initio coupled cluster based simulations on size-selected isolated (n = 0) and hydrated HBDI⁻(H₂O)_n (n = 1, 2) clusters. Microhydrated HBDI⁻ clusters also represent the first step toward modeling HBDI⁻ in aqueous solutions. Our investigations have revealed that water molecules prefer to interact with HBDI⁻ on the phenolate end, forming H-bonds with the phenolate oxygen, which shifts the highest occupied molecular orbitals (HOMOs) down in energy relative to the n = 0 case, leading to stepwise increase to larger detachment energies, while marginally influencing the excitation energies. The high level ab initio methods we benchmarked with the experiments are of high accuracy, being able to reproduce the experimental detachment energies within 1 kcal/mol. Therefore, this type of calculation combined with NIPES can not only provide a reliable description of the electronic structures of the GFP chromophore under microhydrated environments but also can serve as a structure probe to precisely identify the structure of the experimental system via comparison of the calculated detachment energies of different structures with the experimental measurements.

The NIPES experiments were conducted on an apparatus consisting of an electrospray ionization (ESI) source, a low-temperature ion trap, a time-of-flight (TOF) mass spectrometer, and a magnetic bottle TOF photoelectron spectrometer.⁴⁹ HBDI was synthesized according to published methods.⁵⁰ HBDI⁻(H₂O)_n anions were produced via spraying ∼0.1 mM HBDI in a mixture of methanol and water (3:1), titrated with a small amount of NaOH aqueous solution. The produced anions were transported by a RF-only quadrupole ion guide, then guided by a 90° ion bender into a cryogenically controlled 3D Paul trap, where the anions were accumulated and cooled via collisions with ∼0.1 mTorr cold buffer gas (20% H₂ balanced with He) for a period of 20–100 ms with the ion trap temperature set to 20 K. The cooled ions were ejected out of the Paul trap at a 10 Hz repetition rate to the extraction zone of a TOF mass spectrometer for mass analysis. The HBDI⁻(H₂O)_n anions were mass-selected and decelerated before being intercepted by a probe laser beam in the photodetachment zone. In the current study, photon energy of 266 nm (4.661 eV) from a Nd:YAG laser was used. The laser was operated at a 20 Hz repetition rate with the ion beam off at alternating laser shots, thus allowing for shot-by-shot background subtraction. The photoelectrons were collected at nearly 100% efficiency by a magnetic bottle and analyzed in a 5.2 m long flight tube. The TOF photoelectron spectra were collected and converted to kinetic energy spectra, calibrated with the known spectra of I⁻, ClO$2−$⁠, and Cu(CN)$2−$⁠. The electron binding energy spectra presented in this paper were obtained by subtracting the kinetic energy spectra from the photon energy. The energy resolution was 20 meV (full-width at half maximum, FWHM) for 1 eV kinetic energy electrons.

To calculate vertical detachment energies of the HBDI⁻(H₂O)_n complexes, we employed accurate equation-of-motion coupled-cluster formalism (EOMCC),^51,52 which utilizes the following parameterization of the wavefunction for the K-th state:

where T and R_K represent cluster and state-specific excitation operators, respectively; |Φ〉 corresponds to the so-called reference function, which is usually chosen as a Hartree-Fock determinant. With the appropriate definition of the excitation operator R_K one can not only describe the vertical excitation energies but also study electron-affinities/ionization potentials (EA/IP-EOMCC formalisms). The choice of the IP-EOMCC methodology for calculating VDEs is dictated by its efficiency in encapsulating the electron correlation effects which has been demonstrated in numerous studies.^53–59 In the EA/IP-EOMCCSD methods (EA/IP-EOMCC with singles and doubles) employed in present studies the IP/EA excitation operators $RkIPCCSD$⁠, $RKEACCSD$ do not preserve the number of particles and are used to parameterize wavefunctions for N − 1 (IP) and N + 1 (EA) electron systems (assuming that the reference function describe N-electron system).⁵⁷ For the accurate determination of the vertical excitation energies (VEEs), we employed the completely renormalized EOMCCSD(T) approach (CR-EOMCCSD(T))⁶⁰ which corrects the EOMCCSD VEEs by adding non-iterative corrections due to triple excitations,

where due-to-triples corrections can be expressed in terms of trial wavefunction (〈Ψ_K|) and triply excited moments (M_K,3) of the EOMCCSD equations

The CR-EOMCCSD(T) method has already been used in the previous studies of the GFP chromophores providing highly accurate estimates of VEEs of low-lying excited states.⁶¹ 

In the EA/IP-EOMCCSD and CR-EOMCCSD(T) calculations, we used parallel implementation of these methodologies available in the 6.5 version of NWChem computational chemistry package.^57,62 In our studies IP-EOMCCSD maug-cc-pVTZ basis set was employed.⁶³ Since maug-cc-pVTZ basis set contains d and f functions centered on the C, O, N atoms, the CC calculations are expensive. For example, the accurate maug-cc-pVTZ basis set for HBDI⁻(H₂O)₂ is composed of 822 basis functions. The geometries of the HBDI⁻(H₂O)_n systems have been optimized at the DFT/B3LYP level using maug-cc-pVDZ basis set.

Figure 2 shows the 20 K NIPE spectra of HBDI⁻(H₂O)_n (n = 0, 1, 2) at 266 nm. As reported before,²⁸ the first band with electron binding energy (EBE) = 2.7–3.5 eV of HBDI⁻ corresponds to the transition from the ground state of the anion S₀ to the ground state of the HBDI^• neutral D₀ with associated vibrational excitation in the D₀ state. The 2nd band with EBE = 3.5–4.5 eV corresponds to the transition of S₀ to D₁, the first excited state of the neutral. In general, the spectra of the hydrated cluster anions show similar spectral patterns with respect to the bare HBDI⁻, but the spectra shift to higher electron binding energy with increasing number of water molecules. The experimental VDE for each anion is obtained from the maximum of the first spectral band from each spectrum and it is equal to 2.73 ± 0.01, 3.15 ± 0.05, and 3.50 ± 0.05 eV for n = 0, 1, and 2, respectively (Table I). No vibrational structures are resolved in the hydrated clusters because the addition of water molecules introduces additional intermolecular coordinates that might be excited in photodetachment; and these spectra were taken under less optimal instrumental resolution due to the weak mass intensities. The experimental VDE of each anion species measures the energy required for the vertical transition from the anion to the neutral without geometry relaxation and corresponds to the energy difference between the neutral radical and anion, both calculated at optimized anion geometry.

In our investigation we have analyzed several possible structures of microhydrated HBDI⁻ clusters, involving both cis and trans solute conformers. The two configurations differ in their intramolecular hydrogen bonding –C–H…–N for cis (Figure 1) and C–H…–O for trans-conformer (see the supplementary material⁶⁴ for the structure of trans-versus cis-conformer). In accordance with previous studies⁴⁵ cis configuration is the most stable one (2.78 kcal/mol below trans) (Table I). Calculations indicate that in the absence of waters the two conformers have essentially the same electron binding energies. The stability of cis conformer persists for the microhydrated structures. Among the different hydration configurations, the solvation of phenol oxygen (O_p) leads to the lowest energy structures for both cis and trans-conformers. For both 1 and 2 water structures, the cis conformer hydrated at the phenolic oxygen is therefore identified as the lowest energy structure. The addition of the 1st water raises EBE by 0.37 eV to 3.156 eV. The addition of the 2nd water results in further stabilization by 0.32 eV with final EBE at 3.479 eV. Overall the EBE of hydrated cis-confomer is in excellent agreement with the experiment, which identifies it as species generated in the experimental measurements. The next lowest energy structure is a trans-conformer solvated also at the phenol oxygen. While its electron binding energy also matches well with the experiment, it has a relatively high energy (+2.79 kcal/mol for 1 water, +2.99 kcal/mol for 2 water). Other higher energy stable solvated structures can be generated by shifting solvation from phenol to imidazolone ring, either at oxygen or nitrogen position. Given the presence of intramolecular hydrogen bond at the nitrogen position in cis conformer, nitrogen solvated structures for single water clusters are slightly higher in energy than their oxygen counterparts. The opposite situation is observed for trans, where the presence of intramolecular hydrogen bond at the oxygen position disfavors its solvation.

The calculated EOMCCSD ionization potentials are shown in Table I. In all IP-EOMCCSD calculations core electrons were not correlated. To estimate basis set effects we performed IP-EOMCCSD calculations using maug-cc-pVDZ and maug-cc-pVTZ basis sets. The IP-EOMCCSD maug-cc-pVDZ results are consistently smaller by 0.2 eV than that of maug-cc-pVTZ for all systems considered in Table I (vide infra). Using more accurate maug-cc-pVTZ basis set we were able to reproduce the experimental values within corresponding margins of error for the HBDI⁻(H₂O)_n (n = 1, 2) complexes. For the HBDI⁻ anion the IP-EOMCCSD values of 2.788 (cis) and 2.786 (trans) eV are in a very good agreement with the experiment (2.73 eV) with errors not exceeding 0.06 eV for both configurations. For the HBDI⁻(H₂O)_n system, the IP-EOMCCSD result for the cis-O_p configuration (3.156 eV) is in a good agreement with experimentally inferred value of 3.15 eV, while the trans-configuration gives slightly lower value of 3.145 eV. The IP-EOMCCSD calculations for remaining conformers yield the results, which are significantly different from the experimental value. Overall our results suggest that it is cis-O_p configuration that gives rise to experimentally observed photoelectron spectrum. Analogously, the IP-EOMCCSD results for the HBDI⁻(H₂O)₂ complex suggest that only cis-O_p–O_p and trans-O_p–O_p are the experimentally observed structures, but the latter is ∼3 kcal/mol higher in energy. In both cases the IP-EOMCCSD results (3.479 and 3.467 eV, respectively) are within experimental error bar result (3.50 ± 0.05 eV). It should be noted that our values differ from the previous calculations,⁴⁵ where ω B97X-D/6-311++G(2df, 2pd) electron detachment energies for HBDI⁻(H₂O)_n (n = 0, 1, 2), even corrected by the DFT/EOM-IP-CCSD using a smaller basis set, appear to be underestimated by ∼0.24–0.30 eV compared to our experimental numbers.

The analysis of leading amplitudes defining IP-EOMCCSD cluster operators (Eq. (1)) shows that for all experimentally observed HBDI⁻(H₂O)_n (n = 0, 1, 2) complexes the leading excitations correspond to the removal of the electron from the HOMO orbital. Although these amplitudes are dominant for the N − 1 particle systems their role slightly decreases with the water molecules being added. Therefore, we expect increasing role of the correlation effects in describing larger HBDI⁻(H₂O)_n complexes with n > 2. It is also instructive to explore the basis set effects on the quality of the calculated IP-EOMCCSD VDEs. For this purpose, we have performed IP-EOMCCSD calculations using smaller maug-cc-pVDZ basis set for all experimentally observed conformers (HBDI⁻ (cis and trans), HBDI⁻(H₂O)₁ (cis-O_p, trans-O_p), and HBDI⁻(H₂O)₂ (cis-O_p–O_p, trans-O_p–O_p)) which showed consistent 0.2 eV differences between maug-cc-pVDZ and maug-cc-pVTZ results. For example, the maug-cc-pVDZ and maug-cc-pVTZ IP-EOMCCSD results for the cis-O_p–O_p system amount to 3.277 and 3.479 eV, respectively, which demonstrates the importance of basis set quality in computational studies of VDEs.

In addition to IP calculations, excited state energy calculation was performed for all cis and trans lowest energy conformers. For this purpose, we utilized the CR-EOMCCSD(T) approach, which is known to be very efficient in dealing with excited-state correlation effects. The CR-EOMCCSD(T), in analogy to the ground-state CCSD(T) approach,⁶⁵ is characterized by the N⁷ numerical scaling (where N symbolically designates system size). To alleviate the computational complexity calculations were performed using maug-cc-pVDZ basis. To estimate error introduced by the basis set size, additional maug-cc-pVTZ calculations were performed for cis and trans isolated HBDI⁻ structures. Our results indicate that both cis- and trans-conformers are characterized by nearly the same excited state energies. The addition of the first water molecule has also very little impact on excited state energies. The addition of the second water molecule results in slightly higher excited state energies (+0.07 eV). Calculations with higher basis set (maug-cc-pVTZ) decrease the excited state energies by ∼0.05 eV. This basis set error is significantly smaller compared to that observed in electron binding energy calculations. Our calculations show a very weak dependence of VEE as number of water, VEE = 2.73, 2.72, and 2.79 eV for cis-(n = 0), cis-O_p (n = 1), and cis-O_p–O_p (n = 2) structures, respectively, consistent with a previous SOS-CIS(D)/cc-pVTZ calculation, which reported VEEs of 2.61, 2.64, and 2.71 eV for the corresponding conformers.⁴⁵ Overall our excited energy calculations clearly indicate that for small solvated clusters the excited state energies lie below the corresponding electron detachment energies, and this difference increases with cluster size, i.e., V DE − V EE = 0.10, 0.44, and 0.69 eV for n = 0, 1, and 2, respectively. Thus, the first few water molecules are shown to progressively stabilize the excited state of the GFP chromophore anion and therefore reduce or completely shut off the electron detachment process as one of possible deactivation processes upon photoexcitation of GFP.

In summary, we have demonstrated that the EOMCC methodologies can be used to analyze the photoelectron spectra for HBDI⁻(H₂O)_n (n = 0, 1, 2). We demonstrate that high-quality basis set was pivotal for obtaining agreement between theoretical and experimental predictions. This justified the use of scalable codes capable of performing large-scale calculations involving over 800 basis functions. The IP-EOMCCSD / maug-cc-pVTZ calculations reproduce experimental data with 1 kcal/mol accuracy and thus allowed us to identify experimentally observed conformers.

NIPE spectra of HBDI⁻(H₂O)_n (n = 0, 1, 2) indicate solvent water molecules appreciably and steadily stabilize the ground state of the chromophore S₀ relative to the electron-detached neutral continuum D₀, resulting in incremental increase of VDE from 2.73 eV for the isolated chromophore to 3.15 and 3.50 eV for the mono- and di-hydrated clusters. IP-EOMCCSD/maug-cc-pVTZ calculated VDEs based on the lowest conformers, in which solvent waters bind to the phenol O end in cis-configuration, reproduced the experimental data within 1 kcal/mol; while the VDEs calculated at the same level based on other structures differ significantly with the experimental numbers. Thus, NIPES combined with IP-EOMCCSD / maug-cc-pVTZ calculations can be used as a structural probe to identify the exact conformers of clusters existed in the experiments. In addition, our CR-EOMCCSD(T) excited state calculations reveal that the solvent water molecules have much smaller influences in VEE, compared to that of VDE. Consequently, the first few water molecules, such as those existed in GFP as crystallographic waters have significant impacts in stabilizing the photoexcited state relative to the photodetached neutral state, and play crucial roles in influencing UV excited-state photoresponse of GFP and related biochromophore ions.

This work was supported by U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, and performed using EMSL, a national scientific user facility sponsored by DOE’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory, which is operated by Battelle Memorial Institute for the DOE. The calculations have been performed using the EMSL and PNNL Institutional Computing both resources located at PNNL, which is sponsored by the Department of Energy Office of Biological and Environmental Research. This material is based upon work supported by the U.S. Department of Energy under EPSCoR Grant No. DE-SC0012432 with additional support from the Louisiana Board of Regents (W.A.S.). We thank Professor De-Qing Zhang of Institute of Chemistry, Chinese Academy of Sciences for providing us the HBDI sample.