MM-Based Calculations of Absorption and Emission Spectra of

Article pubs.acs.org/JPCB

QM/MM-Based Calculations of Absorption and Emission Spectra of LSSmOrange Variants Maike Bergeler,*,† Hideaki Mizuno,‡ Eduard Fron,§ and Jeremy N. Harvey*,† †

Department of Chemistry, Quantum Chemistry and Physical Chemistry Section ‡Department of Chemistry, Biochemistry, Molecular and Structural Biology Section, §Department of Chemistry, Molecular Imaging and Photonics, KU Leuven, Celestijnenlaan 200F, 3001 Leuven, Belgium S Supporting Information *

ABSTRACT: The goal of this computational work is to gain new insight into the photochemistry of the fluorescent protein (FP) LSSmOrange. This FP is of interest because besides exhibiting the eponymous large spectral shift (LSS) between the absorption and emission energies, it has been experimentally observed that it can also undergo a photoconversion process, which leads to a change in the absorption wavelength of the chromophore (from 437 to 553 nm). There is strong experimental evidence that this photoconversion is caused by decarboxylation of a glutamate located in the close vicinity of the chromophore. Still, the exact chemical mechanism of the decarboxylation process as well as the precise understanding of structure−property relations in the measured absorption and emission spectra is not yet fully understood. Therefore, hybrid quantum mechanics/molecular mechanics (QM/MM) calculations are performed to model the absorption and emission spectra of the original and photoconverted forms of LSSmOrange. The necessary force-field parameters of the chromophore are optimized with CGenFF and the FFToolkit. A thorough analysis of QM methods to study the excitation energies of this specific FP chromophore has been carried out. Furthermore, the influence of the size of the QM region has been investigated. We found that QM/MM calculations performed with time-dependent density functional theory (CAM-B3LYP/D3/6-31G*) and QM calculations performed with the semiempirical ZIndo/S method including a polarizable continuum model can describe the excitation energies reasonably well. Moreover, already a small QM region size seems to be sufficient for the study of the photochemistry in LSSmOrange. Especially, the calculated ZIndo spectra are in very good agreement with the experimental ones. On the basis of the spectra obtained, we could verify the experimentally assigned structures.

■

INTRODUCTION Fluorescent proteins (FPs) are widely used in noninvasive (protein) labeling experiments.1−3 One of the most famous FPs is the jellyfish green fluorescent protein (GFP),4 but a broad range of other FPs exist. All FPs have in common that a certain part of the protein serves as a chromophore, which emits light in the visible range of light. Different FPs can fluoresce at different wavelengths, leading to many different colors. The absorption and emission wavelengths can be tuned by, for example, adjusting the length of the mesomeric system, the chemical composition, and the environment of the chromophore. One FP with emission in the orange range of light is LSSmOrange. It is special because of the large shift between the absorption and emission bands, which led to its abbreviation “LSS”, meaning “large Stokes shift” (typically for gaps between absorption and emission above 100 nm). Although this synonym is commonly used within the community, it is actually better to refer to LSS in LSSmOrange as “light-induced spectral shift” (suggested by Fron et al.5) because the difference between the absorption and emission wavelengths is not only due to vibrational and environmental relaxation and is hence © 2016 American Chemical Society

not a Stokes shift. Because of the large spectral shift, LSSmOrange has many applications in multicolor labeling.6−9 The chromophore of LSSmOrange is a 2-[(5-)-2-hydroxydihydro-oxazole]-4-(p-hydroxybenzylidene)-5-imidazolinone moiety, which is embedded in a β-barrel protein structure (see Figure 1). It is covalently bound to two amino acids of the protein structure. The phenol group of the chromophore can be protonated or deprotonated depending on the protein environment and the electronic state of the system. It has been found that LSSmOrange undergoes (like many other FPs) an excited-state proton transfer (ESPT).5,10 Interestingly, Fron et al. have observed that besides the ESPT also a Kolbe-like decarboxylation of a glutamate close to the chromophore (Glu215)5,11,12 can occur (in the first excited state (ES)). This decarboxylation process leads to a change in the absorption wavelength of the chromophore (from 437 to 553 nm) and is thus also referred to as a photoconversion process (Phot. conv.). It is noteworthy that the photoconversion leads to a Received: September 30, 2016 Revised: November 18, 2016 Published: November 22, 2016 12454

DOI: 10.1021/acs.jpcb.6b09815 J. Phys. Chem. B 2016, 120, 12454−12465

Article

The Journal of Physical Chemistry B

Figure 1. Schematic drawing of the chromophore embedded in the LSSmOrange protein structure (in cartoon representation). “Phot. conv.” means photoconverted.

snapshots and then calculating the vertical emission energy from these optimized structures22 or Car-Parrinello MD has been employed to simulate the dynamics in the ES.23 The disadvantage of the first method is that only local minima on the potential energy surface are considered in the spectra, although at a given temperature, sampling from MD may mean that higher-energy structures are also accessible and thus contribute to the spectra. In the second case, the length of MD simulation is significantly limited in comparison to the one in classical MD simulations. Because of these limitations, here we investigated the possibility of obtaining reasonable emission spectra from TDDFT and ZIndo vertical emission energies calculated on classical MD simulation snapshot structures.

change in the charge of the immediate chromophore environment (by 1 elementary unit), which might influence the properties of the chromophore. The chemical mechanisms leading to the observed absorption and emission properties are not yet completely understood. By measuring, for example, the time-resolved absorption spectrum, one can get some information about the mechanisms.5 Still, there are limitations in assigning structural changes to the observed properties, and at this point, computations can provide more detailed insight into the structure−property relations. Thus, in this work, force-field (FF)-based molecular dynamics (MD) and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations have been performed to model the absorption and emission spectra of the wild-type (WT) and photoconverted forms of LSSmOrange. For MD studies, the necessary FF parameters of the chromophore in the ground state (GS) and ES have been developed. A detailed description of the procedure can be found in the Computational Methodology section. MD snapshots have been taken as input structures for QM/MM and QM + polarizable continuum model (PCM) calculations. For the QM part of the QM/MM calculation, time-dependent density functional theory (TDDFT) with the CAM-B3LYP density functional has been chosen to obtain vertical excitation and emission energies. Furthermore, QM calculations with the semiempirical ZIndo method in a polarizable continuum solvation model have been performed. With this method, reliable excitation energies have already been obtained for other systems.13 Caricato et al. observed that the ZIndo + PCM excitation energies of the systems they studied were even closer to the experimental than the TDDFT ones. Similar QM/MM-based calculations of absorption spectra have already been performed for several different fluorescent molecules.14−21 In contrast, the QM/MM calculation of emission spectra directly from MD snapshots has not yet been commonly done. This might be due to the fact that an optimization of the FF resembling the desired ES has to be carried out, which can be very different from a harmonic potential and thus might not be well reproduced by the parametrized FF. Instead, several emission spectra have been calculated by performing excited-state QM/MM geometry optimizations starting from structures drawn from GS MD

■

COMPUTATIONAL METHODOLOGY

As a starting structure for the calculations, we have chosen the crystal structure with protein data bank entry “4Q7R”. Its chromophore is called “OFM”, which is a 2-[(5-)-2-hydroxydihydro-oxazole]-4-(p-hydroxybenzylidene)-5-imidazolinone moiety. It should be noted that the phenol group of the chromophore can be protonated or deprotonated, depending on the protein environment. We studied two different forms of LSSmOrange, the WT and the photoconverted form, in which glutamic acid 215 exists in its natural form (at pH = 7 nondecarboxylated, as glutamate) and is decarboxylated, respectively. In the GS WT form of LSSmOrange, the protonated chromophore seems to be most stable (denoted “WT_0”), whereas in the photoconverted form, the deprotonated chromophore seems to be favored.5 In the ES, ESPT of the phenol proton occurs, which leads to the deprotonated WT_1 structure. FF parameters have been generated for the protonated and deprotonated chromophores in the GS and the first excited singlet state. Thus, in total, four different chromophore FFs have been parametrized. Additionally, the decarboxylated glutamic acid has been parametrized. For all other atoms, the CHARMM22 FF has been chosen. The FF for the decarboxylated glutamic acid has been obtained by a comparison of the FF parameters of glutamic acid and RCH3-containing amino acids such as leucine, isoleucine, and threonine. For the chromophore parameters, CGenFF24 has been employed. The CGenFF interface provides the first parameter set and a judgment about which of those parameters 12455


Article


have set them to 4.50 kcal/(mol rad2), multiplicity 2, and phase 180.00°, as done by Armengol et al.20 For the dihedral angle parametrization, the maximum force constant has been set to 20 kcal/mol and the energy cutoff to 20 kcal/mol. The default maximum force constant of 3 kcal/ mol is not appropriate for our system as double bonds and rings are involved, which require a significantly larger force constant. We chose the energy cutoff of 20 kcal/mol because for our study we are only interested in structures that may be formed at 310 K and we are not studying any rotational barriers. We started the optimization in the simulated annealing mode and switched to the downhill mode in the second iteration step. With this routine, we obtained a root-meansquare error (RMSE) between the DFT and FF energies of 1.225 kcal/mol for the WT_0, GS chromophore, which appears reasonably low to us. We obtained similar GS dihedral force constants as those obtained by Reuter et al. for GS GFP ϕ and τ dihedrals, namely, 2.438 kcal/mol (compared to 1.400 kcal/ mol) and 5.816 kcal/mol (compared to 6.840 kcal/mol), respectively. As Reuter et al. observed that dihedrals ϕ and τ (and only these) differ quite significantly for the neutral and the anionic chromophore, we chose the two parametrized dihedrals they obtained for the anionic GFP for our anionic system as well, that is, 2.700 kcal/mol for ϕ and 3.900 kcal/mol for τ (periodicity = 2 and phase = 180.00°). In the Supporting Information (SI) we plotted the occurrence of ϕ and τ dihedral angles between −60 and +60° for all studied systems during the MD simulation (as presented in ref 20). One can observe that the chromophore is mostly only slightly distorted from planarity, by up to 20°. A similar change in these two dihedral angles has been observed by Armengol et al. for red FPs mNeptune1 and mCardinal. For the first excited singlet state chromophore, we started the FF parametrization from the GS values. The charges were replaced by CHelpG charges.30 Adaptions have been made such that the charge of aliphatic hydrogens equals 0.09 au and the overall charge resembles that of the whole molecule. An additional important fitting criterion was that the dipole resembles the one created by the originally calculated CHelpG charges. The FFToolkit has not been employed to optimize these charges because the system size of the chromophore is too big. The bond length, angles, and dihedrals have been optimized with CGenFF similar to the optimization of the GS FFs. Therefore, the FFTK source code had to be slightly modified to read excited-state energies. After dihedral angle parametrization, we obtained a RMSE between the DFT and FF energies of 0.946 kcal/mol for the WT_1, ES chromophore. For the molecular dynamics setup, we first evaluated the protonation state of the amino acids at pH = 7 in LSSmOrange with the program H++.31−34 As a result, we obtained that no aspartic acid is protonated, no glutamic acid is protonated, no HSD exists, amino acids 17 and 75 are of HSE type, and amino acids 25, 172, 204, and 221 are of HSP type. We also checked whether cysteine residues might be able to form disulfide bridges by comparing the distances between them. In LSSmOrange, there are no such cysteine residues. The preparation of the MD simulation has been done with CHARMM (CharmmGui).35,36 The protein has been solvated in a cubic water box of 74 Å edge length. Periodic boundary conditions (PBCs) have been employed. An equilibration of 50 ps followed by a 20 ns MD simulation (NpT at T = 303.15 K and p = 1.01325 bar, with Langevin dynamics) has then been performed with NAMD.37 Shake38 has been employed for all

have to be optimized. For several parameters (mainly the two five-membered rings), CGenFF assigned estimated parameter error terms (referred to as “penalty values”) larger than 50, which means that these FF parameters should be optimized. It should be noted that in line with the concept of the CGenFF program, we use improper dihedral terms only for carbonyl groups as most other out-of-plane motions seem to be reproduced correctly by valence angle and (proper) dihedral terms alone.25 The (proper) dihedral terms have been carefully parametrized with respect to DFT energies of the corresponding scans, and according to the output of our fitting procedure, they do indeed suffice to describe our system properly. We have done the optimization with the FFToolkit26 of visual molecular dynamics (VMD) 1.9.2 and VMD 1.9.127 in combination with Gaussian 09.28 For the optimization of the charges in the neutral chromophore, the chromophore structure has been split into smaller parts. Only in the two five-membered rings, the charges were assigned with high penalties (larger than 50); thus, we optimized these charges separately (with the FFToolkit that takes into account the QM dipole of the molecule) and then combined the newly obtained charges with those having a low enough penalty. Furthermore, all hydrogen atoms connected to sp3 carbon atoms (aliphatic) were set to a charge of 0.09 au, as it is commonly done in the CHARMM FF. As the overall sum of the charges has to be equal to the molecular charge, we had to slightly adapt the charges to yield exactly the desired overall charge. The charges of the anionic chromophore have been adapted from the neutral chromophore by taking into account the existing CHARMM22 parametrized charges of phenoxide (from toppar_all22_prot_model.str). All (constrained) optimizations and Hessian calculations of the GS capped chromophore, and parts of it, have been performed with B3LYP/6-31G* including dispersion correction D3. It should be noted that the FF parameters of the GS GFP chromophore fitted against B3LYP/6-31G* QM energies already exist.29 The LSSmOrange chromophore is similar to the GFP chromophore, but the second five-membered ring including oxygen and nitrogen does not exist in GFP. This elongation of the mesomeric system in comparison to the GFP chromophore might change the electronic structure of the chromophore significantly. Thus, we perform a full FF parametrization for LSSmOrange in an analogous way as described in Reuter et al.29 We will compare our parameters to those obtained in this earlier work as they provide a very good representation of their system. We observed small differences in the assigned charges of the five-membered ring that exists in both chromophores, which probably occur due to the different groups being attached to the ring in the two chromophores. As our CGenFF optimized charges reproduce the CHelpG charges well, we kept our charges. It should be noted that Reuter et al. observed significant differences in the two dihedral force constants describing the planarity of the chromophore (τ and ϕ, as defined in ref 20 for red FPs mNeptune1 and mCardinal) for the neutral and the anionic chromophore. Thus, we have decided to choose their dihedral force constants for the two angles, τ and ϕ, for the anionic system. All other anionic dihedral angles are the same as those in the neutral system. It is noteworthy that we have adapted a few dihedral angles after the parametrization with FFTK. The dihedral force constants obtained with the FFTK parametrization for the four dihedral angles X-CG2R57-CG2R57-X seem to be too low. Thus, we 12456


Article


Figure 2. Representative scheme of the MD (with PBC) and QM/MM model systems.

Table 1. Overview of the Different Systems Studieda

C: chromophore (OFM), CH: protonated chromophore at tyrosine, B: deprotonated Asp161, BH: protonated Asp161, R: Glu215 without −CO−2 . GS: ground state, ES: excited state. WT_0 is wild-type LSSmOrange before ESPT and WT_1 after ESPT.

a

With CAM-B3LYP, we mainly employed the 6-31G* basis set for spectral calculations but also performed tests with larger basis sets (as denoted). For ZIndo calculations, the PCM with a dielectric constant similar to that of the surrounding protein (diethyl ether, ϵ = 4.24) has been applied. A nonequilibrium solvation is assumed because the structure of the environment does not adapt to the electronic state of the solute. For benchmarking purposes, RICC2 calculations have also been performed on a minimal model system with TURBOMOLE 7.052 by making use of the freeze option (with default number of orbitals being frozen). Furthermore, CASPT2 calculations have been run with MOLPRO53 after performing a CASSCF(6,6) or CASSCF(8,8) calculation. The active spaces in our CASSCF(6,6) and CASSCF(8,8) were simply chosen by including the appropriate number of frontier orbitals, and the orbitals can roughly be described as having π or π* character. The TDDFT calculations show that the most important transition has a predominant HOMO to LUMO character, with the HOMO and the LUMO corresponding to the π and π* orbitals. These calculations have been performed on a minimal model system of the chromophore. The absorption spectra have been calculated from vertical excitation energies based on GS MD structures, and the emission spectra, from vertical emission energies starting from

bonds involving hydrogen atoms; thus, a time step of 2 fs could be applied. To check whether the simulation time is long enough, we calculated the all-protein-atom root-mean-square deviations (RMSDs), excitation energies, and root-mean-square fluctuation (RMSF) windows (for the chromophore only, averaged over 1 ns each) over time. The results are shown in the SI. The excitation energies seem to be converged as no shift over time can be observed. The RMSD does not seem to be completely converged for all simulated states, but the RMSF of the chromophore itself is relatively small; thus, we assume that the simulation time is long enough for the purpose of this study. For the QM/MM and QM in a PCM calculations, snapshots have been taken from the MD run every 0.1 ns starting from 2 ns in MD simulation. In total, this leads to 180 snapshot structures per simulation. We deleted all water molecules further away than 30 Å radius from the chromophore center. Figure 2 shows the MD and QM/MM system sizes. QM/MM calculations were performed with QOMMMA.PY39 using GAUSSIAN28 for the QM part and TINKER40 for the MD part. Semiempirical ZIndo/S41−49 (ZIndo will be used synonymously here) and CAM-B3LYP50 TDDFT single-point calculations have been performed on MD structures to obtain vertical excitation energies to the first excited singlet state. ORCA51 has been chosen for the ZIndo/MM test calculations. 12457


Article


Table 2. Excitation Wavelengths (λcalc, in nm), the Corresponding Energy Differences (ΔE, in kcal/mol), and Oscillator Strengths ( f) Obtained for Different Computational Methodsa

a

As an example, the chromophore model system (“chromophore”) and the B3LYP/Def2TZVP optimized structure of a minimal model system of an arbitrarily chosen snapshot structure from the GS/WT_0 MD simulation have been studied. The two model systems are depicted next to the table entries.

“excited-state MD structures”. The systems, which have been studied here, are described in detail in Table 1. All excitation energies with an oscillator strength above 0.1 (from the first three excitation energies) have been taken into account for modeling the absorption spectrum. For the emission spectra, only the first excitation energy, which resembles the energy difference between the GS and the first excited singlet state, has been taken. Lorentzians of 10 nm widths and height according to the oscillator strength have been added up and plotted with GNUPLOT. We normalized each spectrum by dividing it by its integral. Different QM regions have been considered. The smallest model system consists of the chromophore only. All larger model systems have been built by either choosing special amino acids to be taken into account or employing an automatized script that evaluates whether an amino acid is taken into account completely or partially (e.g., by its side chain) on the basis of a spherical cutoff around the whole chromophore or by employing both together. We have chosen the criterion such that if one of the backbone atoms is within the threshold distance, the whole amino acid is taken into account; otherwise, only the side chain is taken into account. This has been done for cutoff values of 1.5−5.0 Å. It is noteworthy that as the script is run for every snapshot, the QM regions can differ in the snapshots. Therefore, a direct comparison between the QM total energies is not possible.

For all QM in PCM calculations, the amino acid backbones have been cut such that the whole peptide bond has always been taken into account and additionally the carbon Cα atom, which is saturated with hydrogen atoms. A principal component analysis of the chromophores based on the MD snapshot structures has been carried out with BIO3D.54 Molecular structures were visualized with PyMOL.55

■

RESULTS AND DISCUSSION

Comparison of QM Methods. Several benchmark studies have been published that compare the performance of quantum chemical methods to obtain excited-state energies for different organic molecules.56−61 One should note that the results were compared with either the experiment or the results obtained from high-level QM methods such as RICC2 or CASPT2. In summary, it has been observed that TDDFT can yield relatively good results with hybrid functionals that contain between 22 and 25% of exact exchange and with range-separated hybrid functionals, such as, for example, CAM-B3LYP. With this method an accuracy of 0.2−0.3 eV could be achieved (compared to the experiment). Also LC-ωPBE62−64 and ωB97XD65 provided good excitation energies. The most precise results could be achieved with coupled-cluster (CC2 and CC3) and second-order perturbation (CASPT2) theories employing large basis sets (e.g., aug-cc-pVTZ). 12458


Article


Table 3. Excitation Wavelengths (“λcalc”, with the Largest Oscillator Strength) in nm (and Excitation Energies in kcal/mol in Brackets) and Their Corresponding Oscillator Strengths (f) Obtained for Different QM Modelsa CAM-B3LYP + PCM

CAM-B3LYP/MM

ZIndo + PCM

ZIndo/MM

QM model(charge)

λcalc

f

λcalc

f

λcalc

f

λcalc

f

chromophore(0) 1.5 Å(0) 1.8 Å(−2) 2.0 Å(−1) 2.5 Å(−1) 3.0 Å(0) 4.0 Å(0) 5.0 Å(−2) C_144_bb(0) C_144(−1) C_215_sc(−1)

365.6(78.2) 364.5(78.4) 355.3(80.5) 366.0(78.1) 363.9(78.6) 363.0(78.8) 363.5(78.7) 364.9(78.4) 374.5(76.4) 378.4(75.6) 350.7(81.5)

0.93 0.82 0.78 0.94 0.80 0.91 0.91 0.49* 1.01 1.04 0.84

344.5(83.0) 346.2(82.6) 350.2(81.7) 353.1(81.0) 356.2(80.3) 358.8(79.7) 360.3(79.4) 360.6(79.3) 349.2(81.9) 349.6(81.8) 343.8(83.2)

0.81 0.73 0.83 0.84 0.78 0.76 0.69 0.68 0.93 0.95 0.78

399.7(71.5) 398.0(71.8) 391.3(73.1) 395.6(72.3) 391.4(73.0) 393.4(72.7) 392.8(72.8) 391.4(73.0) 404.0(70.8) 407.2(70.2) 388.2(73.6)

0.99 0.97 1.00 1.00 0.98 0.99 1.00 0.97 1.00 1.00 1.01

381.8(74.9) 375.7(76.1) 377.3(75.8) 377.6(75.7) 378.9(75.5) 380.8(75.1) 380.8(75.1) 381.0(75.0) 380.4(75.2) 380.5(75.1) 372.5(76.8)

0.83 0.81 0.83 0.84 0.81 0.81 0.82 0.79 0.87 0.88 0.84

a

As an example, an arbitrarily chosen snapshot structure from GS/WT_0 (H-bond with oxygen of Glu144) has been studied. For PCM a dielectric constant of 4.24 has been chosen, which approximately resembles the protein environment. Experimental maximum absorption is at 437 nm (65.4 kcal/mol). C_XXX means chromophore plus amino acid XXX; C_XXX_bb means chromophore plus backbone of amino acid XXX (only); and C_XXX_sc means chromophore plus side chain of amino acid XXX (only). CAM-B3LYP stands for CAM-B3LYP/D3/6-31G*. The asterisk denotes that in this case two very similar excitation energies with similar oscillator strengths occur.

excitation energy. It is noteworthy that it is not completely comparable to the other values as it is calculated with the PCM and not with the inclusion of MM charges. The ZIndocalculated oscillator strengths are by about 0.1−0.2 larger than the oscillator strengths calculated with TDDFT and RICC2 for the minimal and the chromophore model systems. This has a small impact on the exact shape of the spectra because the oscillator strength is chosen as the intensity of each excitation or emission. Nevertheless, the impact should not be significant if we assume that for one chosen QM method all excitations are assigned with a slightly too large or too low oscillator strength. The overall differences between the RICC2, TDDFT(CAMB3LYP), and ZIndo results are small. Therefore, it appears to be a valid approach to us to employ TDDFT(CAM-B3LYP) and ZIndo for the excited-state calculations of the LSSmOrange chromophore. As for the absorption and emission spectra, a large amount of single-point calculations is required and it would be advantageous to reduce the computational cost of each calculation. Therefore, we compared the basis set dependence of the excitation energies within the CAM-B3LYP calculations. An increase of the basis set size from 6-31G to Def2TZVP changes the excitation energy by only 2.7 kcal/mol for the chromophore model system and by 3.5 kcal/mol for the minimal model system. The oscillator strengths differ only by less than 0.05, respectively, which is negligible. These changes (especially in the complete chromophore model system) are relatively small; thus, we decided to calculate the TDDFT(CAM-B3LYP) spectra with a rather small basis set, namely, 6-31G*. Influence of QM System Size. In Table 3, ZIndo and CAM-B3LYP excitation energies and oscillator strengths for several QM regions are given for one arbitrarily chosen snapshot structure out of the 20 ns GS/WT_0 MD simulation. In this snapshot structure, a hydrogen-bonding network between the protonated OFM and the oxygen backbone atom of Glu144 is formed. As the hydrogen-bonding interaction could influence the excitation energy of the chromophore, we studied two model systems, called, “C_144” and “C_144_bb”. The first one takes into account the whole amino acid Glu144, and the second, only the backbone of amino acid Glu144 (in addition to the

In addition, several benchmark studies exist that evaluate different QM methods for the calculation of excitation energies for FPs, mainly the GFP.61,66−70 In general and also here, TDDFT with hybrid functionals, such as PBE0, B3LYP, and CAM-B3LYP, performed reasonably well in comparison to RICC2.70 However, it should be noted that differences between RICC2 results and experimental values also exist. RICC2 seems to overestimate the excitation energy by up to 0.3 eV. Moreover, these differences between RICC2 and experiment as well as between RICC2 and other QM methods also depend on the overall charge of the chromophore (for the neutral chromophore, larger deviations have been observed than for the singly negatively charged chromophore68). In this article, we studied the performance of selected quantum chemical methods (based on the literature benchmark studies) for excited-state calculations on the LSSmOrange chromophore. In Table 2, the excitation wavelengths and oscillator strengths for the whole LSSmOrange chromophore (capped after the peptide bond to the next amino acid) and for its minimal model system are reported. We have performed TDDFT calculations with CAM-B3LYP and ωB97XD and semiempirical ZIndo, RICC2 coupled-cluster, and CASPT2 calculations here. The latter two methods have been employed only on the minimal model system because they are computationally very demanding. We observed that in the gas-phase minimal model system the TDDFT results with the chosen functionals are all in good agreement with the RICC2/Def2TZVP result. With LC-ωPBE, the deviation is slightly larger (overestimation by 4.8 kcal/mol), whereas with CAM-B3LYP and ωB97XD, the results are very similar (deviations smaller than 0.5 kcal/mol). The ZIndo result is by 5.3 kcal/mol lower than the RICC2/Def2TZVP result. Given that the RICC2 energies slightly overestimate the excitation energies in comparison to the experiment, the ZIndo results might even be closer to the experimental ones than the RICC2 and TDDFT(CAM-B3LYP) results. In the QM/MM (and QM in PCM) calculations of the larger chromophore model system, the obtained TDDFT excitation energies are also in relatively good agreement with each other. Only the ZIndo excitation energy deviates more significantly, by 9.7 kcal/mol compared to the CAM-B3LYP/Def2TZVP 12459


Article


and 215 are on the opposite sides of OFM and thus might cause a dipole moment. In the chosen snapshot structure, Arg95 is present (explicitly, and not as MM charge) in all model systems larger than or equal to 2 Å, Glu144 in all model systems larger than 1.8 Å, Glu215 in all model systems larger than or equal to 1.8 Å, and Glu148 in all model systems larger than or equal to 5 Å. Thus, Glu215 and Glu144 might cause the small deviations in the excitation energies observed when the model system size was increased from 1.5 to 1.8 Å. Lys70 was present in the QM region when the radius was 3 Å and above. Glu215 is decarboxylated in the photoconverted form of the chromophore, leading to a change in the surrounding charge by one elementary unit. Additionally, we have analyzed the hydrogen bonds formed between the phenol/phenolate oxygen of the neutral/anionic chromophore and the surrounding. If certain amino acids form hydrogen bonds with the chromophore over a large simulation time, these might be important to be included in the QM region. From our analysis, we observed that the main H-bond networks of the phenol group are formed between the phenol hydrogen and Ser146, Glu144, and water (see SI). The phenolate mainly forms hydrogen bonds with water and amino acid Ser146 and to a small extent with Leu199. As, overall, the excitation energies for the chromophore model system and the larger model systems do not differ in a major way for CAM-B3LYP and ZIndo, the calculations of absorption and emission spectra are based on the computationally affordable chromophore model system as QM region. For the WT_1 (ES) and photoconverted model systems (ES), we have additionally taken into account Ser146 in the QM region (C_146 model) because this amino acid seems to be an important H-bonding partner based on our H-bonding analysis. We employed ZIndo with PCM and CAM-B3LYP with MM charges. However, as our analysis in Table 3 is based on only a single snapshot structure, we also performed excitation energy calculations with a 2.5 Å model region to find out how much the two spectra differ (see SI for details). For the ZIndo + PCM spectra, we observed only very small changes, mainly in the positions of the spectra corresponding to a shift of less than 1.5 kcal/mol. For the CAM-B3LYP/MM absorption spectrum of the photoconverted LSSmOrange, we observed a slightly larger deviation in the positions of the spectra, by +2.4 kcal/mol. We performed these larger QM system calculations also to find out whether the explicit inclusion of Glu215 in the QM region affects the spectra. This amino acid is of special interest because it is negatively charged in the WT and (after decarboxylation) neutral in the photoconverted LSSmOrange. However, we could not identify any significant influence of this single amino acid on the spectra. It should be mentioned that besides the QM region, also a larger basis set in the CAM-B3LYP calculations might slightly change the CAM-B3LYP spectra. On the basis of our studies described in section Comparison of QM Methods, the error of choosing the 6-31G* basis set in the CAM-B3LYP calculations (compared to Def2TZVP) amounts to +2.7 kcal/mol. The plus sign indicates that the approximation leads to an overestimation of the excitation energy. In total, the overestimation of the CAM-B3LYP excitation energies based on our approximations is by 4−5 kcal/mol. With ZIndo, we observed an overestimation of the excitation energies by approx. 1.5 kcal/mol due to the restricted QM region. However, we assume that the qualitative picture stays the same.

chromophore). Furthermore, we have built a model system that includes amino acid Glu215 because this amino acid is charged in the WT and neutral after the decarboxylation in the photoconverted form of LSSmOrange and thus might play an important role in understanding the differences between the spectra. Besides these QM model systems based on specific amino acids, we have chosen distance cutoff criteria from 1.5 to 5.0 Å (see Computational Methodology for details). The approximate number of atoms in the different QM models are approx. 63−64 (chromophore only), 87−88 (C_144 model system), 160 (2.0 Å model system), 245 (2.5 Å model system), 370 (3.0 Å model system), 470 (4.0 Å model system), and 650 (5.0 Å model system). We studied the difference between employing a continuum solvation model (with ϵ = 4.24) around the QM system (QM + PCM) and taking the exact positions and sizes of the MM charges into account (QM/MM) for both QM methods, TDDFT(CAM-B3LYP) and semiempirical ZIndo. First of all, one can observe that within each of the four employed methods the excitation energies do not change significantly with the QM region. Similar findings have been reported for the LSSmKate1 FP.21 From the 2.0 Å QM model on, the differences in the excitation energies compared to those of the 5.0 Å model system are lower than 1.5 kcal/mol. Interestingly, the convergence of the QM/MM excitation energies compared to the convergence of the QM + PCM excitation energies occurs from different sides. In other words, one can observe a slight decrease in the CAM-B3LYP(QM/ MM) and ZIndo(QM/MM) (with the exception of the “chromophore”) excitation energies and a slight increase in the ZIndo (QM + PCM) and CAM-B3LYP(QM + PCM) excitation energies for larger QM regions. In the case of TDDFT within QM/MM, the observed small red shift of the absorption wavelengths is in line with the observations made by others.17,71 We also investigated the influence of the QM region on the excitation energies for three additional snapshot structures taken from the WT_0 MD simulation as well as for three snapshot structures taken from the photoconverted GS and excited-state MD simulations (see Section 1 in the Supporting Information). We observed that for all three investigated systems the ZIndo + PCM and ZIndo/MM excitation energies do not significantly change with the QM model size. The fact that for both QM methods the excitation energies converge very fast with the QM region and that the QM/MM and QM + PCM excitation energies converge to almost the same excitation energy suggests that the excitation energy does not significantly depend on the environment and that in this specific FPfor a sufficient QM region of approx. 2.5 Åthe impact of the more distant MM charges is small. On comparing the ZIndo and CAM-B3LYP results with each other, we observe that the excitation energies obtained from the ZIndo(QM + PCM) calculations are approximately 6−7 kcal/ mol lower than the CAM-B3LYP(QM + PCM) excitation energies (1st vs 3rd row). The ZIndo and CAM-B3LYP QM/ MM excitation energies (2nd vs 4th row) vary by 4−8 kcal/ mol. Furthermore, we investigated the occurrence of charged amino acids in the proximity of the chromophore because these amino acids might influence the excitation energies more significantly than the neutral ones. We found five amino acids that have a charged group close to the chromophore, which are Lys70, Arg95, Glu144, Glu148, and Glu215. Amino acids 95 12460


Article


imidazoline oxygen atom are explicitly taken into account in the QM region (the Lys70 side chain is included from 3 Å on and the Arg95 side chain from 2 Å on). For the photoconverted form with an anionic phenol group, we found a significantly larger charge shift for the phenol unit, which suggests that detailed treatment of interactions between the environment and this phenol group might be important. Possible interaction partners of the phenol oxygen atoms have been identified by the hydrogen-bonding analysis described above. Again the imidazoline oxygen atom shows a significant charge shift, and thus it might be important to take into account possible interactions of this group with the protein environment. Absorption and Emission Spectra. In Figures 3 and 4, the calculated ZIndo + PCM and CAM-B3LYP/MM absorption and emission spectra are shown. The experimental wavelengths with the maximal absorption intensity of the WT_0 LSSmOrange and its photoconverted form are 437 and 553 nm, respectively. The experimental emission spectra of the WT_1 and the photoconverted form are almost identical, with maximum emission wavelengths at 573 nm (WT_1) and 569 nm (photoconverted form).72 One can observe that ZIndo gives very similar absorption spectra compared to the experimental ones (see Figure 3). The CAMB3LYP spectra are shifted by 8−13 kcal/mol to larger excitation energies (or smaller wavelengths) in comparison to the ZIndo spectra (and the measured spectra) (see Figure 4). Despite this disagreement concerning the absolute positions of the peaks, the relative positions of the two absorption spectra and the emission spectra with respect to the absorption spectra are in agreement with those of the experimental spectra. It should be noted that emission from the photoconverted form is very slightly blue shifted compared to that from the WT in the experiments. This small shift (4 nm) is not reproduced by the calculations. The difference between the CAM-B3LYP and the experimental maximum absorption and emission wavelengths is most probably due to a combination of the small basis set and QM region error (as described above) and due to the fact that CAM-B3LYP is known to overestimate excitation energies.70 To directly compare the shape of the CAM-B3LYP (QM/MM) and ZIndo (QM + PCM) spectra, we overlayed the first one with the second in Figure 5. The overlay has been done by shifting each CAM-B3LYP spectrum to the maximum absorption or emission wavelength of the corresponding ZIndo spectrum. One can see that the CAM-B3LYP spectra are all slightly narrower than the ZIndo spectra. Our calculated ZIndo and CAM-B3LYP WT_0 absorption spectra further indicate that the LSSmOrange chromophore is (in contrast to many other FPs) protonated in the GS (as we assumed in the beginning on the basis of experimental observations) because the relative positions of the spectra match the experimental ones reasonably well. Furthermore, the photoconverted absorption spectra (in which the chromophore is deprotonated) are well-separated from the WT_0 absorption spectra for both ZIndo and CAM-B3LYP. However, these spectra are not directly comparable to a (hypothetical) WT_1 absorption spectrum because Glu215 is decarboxylated, which might slightly influence the spectrum. The experimentally observed full width at half-maximum (FWHM) for the WT_0 absorption spectrum is 84 nm, and for the photoconverted form, 63 nm. The FWHMs in the calculated ZIndo/CAM-B3LYP WT_0 spectra are 49/36 nm, and in the photoconverted form, 63/37 nm. Thus, the line

We performed an analysis of the shift of the natural bond orbital (NBO) charge from the GS to the ES as it has been done by Armengol et al.20 for a randomly chosen snapshot structure of the WT_0, GS, and Phot. conv. GS systems. The results are shown in Table 4. We observed that for WT_0 the Table 4. NBO Charges (in au) in the Chromophore in the GS (Q(S0)) and First Excited (Q(S1)) States as Well as the Difference between the Two Charges (ΔQ)a moiety phenol methine bridge imidazoline oxazoline tail O4(phenol) O2(imidazoline) O1(oxazoline) O(oxazoline) N1(oxazoline) C9(oxazoline) phenol methine bridge imidazoline oxazoline tail O4(phenol) O2(imidazoline) O1(oxazoline) O(oxazoline) N1(oxazoline) C9(oxazoline) a

Q(S0)

Q(S1)

WT_0 GS +0.11053 +0.24316 +0.18176 +0.05917 −0.52896 −0.42617 −0.04788 −0.17698 +0.18813 +0.17901 −0.74550 −0.72400 −0.69014 −0.64034 −0.62754 −0.63013 −0.75230 −0.75894 −0.46797 −0.51980 +0.24564 +0.19323 Phot. Conv. GS −0.81343 −0.56384 +0.16397 −0.00293 −0.54053 −0.49336 −0.04421 −0.16223 +0.18395 +0.17304 −0.85245 −0.78230 −0.69346 −0.67280 −0.64387 −0.64778 −0.75654 −0.76318 −0.47641 −0.52330 +0.24457 +0.20133

ΔQ +0.13263 −0.12259 +0.10279 −0.1291 −0.00912 +0.02150 +0.04980 −0.00259 −0.00664 −0.05183 −0.05241 +0.24959 −0.16690 +0.04717 −0.11802 −0.01091 +0.07015 +0.02066 −0.00391 −0.00664 −0.04689 −0.04324

The atom names are the ones defined in Section 5 of the SI.

charge shifts in the defined moieties are of similar sizes compared to the ones calculated for mCardinal and mNeptune1. However, if we look at certain atoms instead of whole groups, the charge shifts on the oxygen atoms of the oxazoline and the phenol subunit are much less pronounced than for the other two FPs (the oxazoline oxygens are compared to the acylimine oxygens of mCardinal and mNeptune1). Only on the imidazoline oxygen atom, we observe a comparably large charge shift. We conclude that the protein environment around the imidazoline oxygen atom might influence the excitation energies but not the protein environment around the oxygen atoms of the oxazoline group. This part of the chromophore is also quite different in the LSSmOrange protein compared to the other two chromophores (oxazoline group instead of acylimine group). As, overall, the charge shift in the oxazoline group is not negligible, we investigated which atoms have the biggest charge shifts and observed that the sp2-hybridized carbon and nitrogen atoms gain a significant amount of negative charge. As this moiety (CN) is rather buried in the center of the chromophore, we do not expect significant changes in the excitation energies based on interactions of this moiety with the protein environment. When we compare these results with those in Tables 3 and S1, we see that the influence of the protein environment on the excitation energies is small, even if the protonated Lys70 or Arg95 amino acids close to the 12461


Article


Figure 3. Calculated (Calc.) and experimentally measured (Exp.) ZIndo + PCM absorption (Abs.) and emission (Emiss.) spectra of the WT_0, WT_1, and the photoconverted (Phot. conv.) form of LSSmOrange (in wavelengths). As the QM region for the GS, only the chromophore has been chosen. For the ES, the chromophore and the side chain of Ser146 have been considered because for many snapshot structures, a hydrogen bond between the phenolate group of the chromophore and the HG of Ser146 has been observed.

Figure 4. CAM-B3LYP/MM absorption and emission spectra of the WT_0, WT_1, and photoconverted (Phot. conv.) forms of LSSmOrange (in wavelengths). The vertical lines represent the experimental maximum absorption or emission wavelengths.

Figure 5. Overlay of the CAM-B3LYP/MM (denoted “CAM”) and ZIndo + PCM absorption (Abs.) and emission (Emiss.) spectra of the WT and photoconverted (Phot. conv.) LSSmOrange variants. Each CAM-B3LYP spectrum has been shifted to the maximum absorption or emission wavelength of the corresponding ZIndo spectrum.

width of the ZIndo spectrum of the photoconverted species agrees very well with the measured one, but the ZIndo WT_0 absorption spectrum has a slightly too low FWHM. With CAM-B3LYP, the line widths are both too small compared to the measured ones. The FWHMs in the ZIndo/CAM-B3LYP emission spectra are 68/41 nm, and for the WT_1 and photoconverted system, 76/45 nm, which is relatively similar to the experimental FWHMs of 48 and 55 nm. In the experimental absorption spectrum of the photoconverted form of LSSmOrange, the intensity drops down very sharply after the maximum absorption. This has not been observed in the calculated spectrum, which may indicate shortcomings in the accuracy of the MM FF used here for sampling structures.

As the ESPT most probably leads to a protonation of Asp161, we investigated the difference in the emission spectra of WT_1 with Asp161 protonated and deprotonated (see SI). We observed that the spectrum calculated for WT_1 Asp161 protonated is slightly shifted to a smaller wavelength compared to the one with deprotonated Asp161. Moreover, we performed a series of analyses to better understand the wavelength positions and shapes of the calculated spectra. Therefore, we investigated a possible correlation between several structural parameters of the snapshot structures and the calculated excitation energies (the results can be found in the SI). First, we performed a principal 12462


Article


an FF that describes the excited-state potential followed by vertical excitation-energy calculations, overall, appears to provide a reasonable model for the studied system (and most probably also for similar) FPs.

component analysis of the chromophore and plotted the excitation energies over the first two principal components of each studied system, but no correlation could be observed. We also investigated a possible correlation between the two dihedral angles ϕ and τ and the excitation energy of the chromophore, but we did not find any relationship. In the third step, we analyzed the chromophore bond lengths and angles and found that certain bond lengths correlate with the excitation wavelengths. We found that the sum of four distances over the ZIndo excitation wavelength yields a linear relationship with squared residues between 0.4 and 0.7 (see SI). Also, the CAM-B3LYP excitation energies of WT_0 correlate with the sum over the four distances. However, for the WT_1 and photoconverted forms, this correlation could not be observed (see the SI for details).

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b09815. MD trajectory analyses, structure−property analyses, additional calculated absorption and emission spectra, and the FF parameters (PDF)

■

■

AUTHOR INFORMATION

Corresponding Authors

CONCLUSIONS We first optimized the FF parameters for the GS and first excited singlet state of the LSSmOrange chromophore in its protonated and deprotonated forms. Then, we performed a careful analysis of the computational methodology employed for the QM part of the QM/MM calculation and studied the size of the QM region. For this specific chromophore, we found that TDDFT CAM-B3LYP with a rather small basis set (631G*) should already provide reasonable qualitative results. Moreover, we observed that for this specific system the QM region does not have to be chosen very large because the protein environment of the chromophore does not have a strong influence. The inclusion of the three closest charged amino acids with respect to the chromophore, Arg95, Glu148, and Glu215, and the inclusion of additional amino acids within 5 Å from the chromophore only slightly influence the absorption and emission wavelengths. Thus, we chose two small QM regions including only the chromophore or the chromophore and the Ser146 side chain, which has been observed to form hydrogen-bonding networks that might influence the excitation energies. From the calculated spectra, we observed that the ZIndo absorption and emission spectra are close to the experimental ones. Our calculations support the idea that the photoconverted species involves absorption by and fluorescence from an anionic chromophore, whereas WT involves absorption by the protonated and emission by the deprotonated chromophore. With CAM-B3LYP/6-31G*, an overestimation of the excitation energies by approximately 8−13 kcal/mol in comparison to the ZIndo and experimental results has been observed, which is in accordance with the previous results.70 Still, qualitatively, the CAM-B3LYP results reproduce the shifts in the different ZIndo absorption and emission spectra. In addition to this error in the absolute positions of the absorption and emission peaks with CAM-B3LYP, both the ZIndo and TDDFT approaches lead to other noticeable differences with respect to the experiment. For example, the line shape and width, although in fair agreement with the experiment, are not perfect. Also, as noted above, a small blue shift in the emission spectrum is observed experimentally upon photoconversion. Neither ZIndo nor CAM-B3LYP shows a blue shift. These discrepancies could be due to the shortcomings in the electronic structure method or the sampling protocol and warrant further investigation. Nevertheless, we conclude that the approach we have chosen here to obtain emission spectra, namely, the parametrization of

*E-mail: [email protected] (M.B.). *E-mail: [email protected] (J.N.H.). ORCID

Maike Bergeler: 0000-0003-0944-8605 Jeremy N. Harvey: 0000-0002-1728-1596 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work has been financially supported by the Schweizer Nationalfond (SNF Early PostDoc Scholarship). REFERENCES

(1) Lippincott-Schwartz, J.; Altan-Bonnet, N.; Patterson, G. H. Photobleaching and Photoactivation: Following Protein Dynamics in Living Cells. Nat. Cell Biol. 2003, 5, S7−S14. (2) Miyawaki, A. Fluorescent Proteins in a New Light. Nat. Biotechnol. 2004, 22, 1374−1376. (3) Lukyanov, K. A.; Chudakov, D. M.; Lukyanov, S.; Verkhusha, V. V. Photoactivatable Fluorescent Proteins. Nat. Rev. Mol. Cell Biol. 2005, 6, 885−891. (4) Tsien, R. Y. The Green Fluorescent Protein. Annu. Rev. Biochem. 1998, 67, 509−544. (5) Fron, E.; Keersmaecker, H. D.; Rocha, S.; Baeten, Y.; Lu, G.; Ujii, H.; der Auweraer, M. V.; Hofkens, J.; Mizuno, H. Mechanism Behind the Apparent Large Stokes Shift in LSSmOrange Investigated by Time-Resolved Spectroscopy. J. Phys. Chem. B 2015, 119, 14880− 14891. (6) Kogure, T.; Araki, S. K. T.; Saito, K.; Kinjo, M.; Miyawaki, A. A Fluorescent Variant of a Protein from the Stony Coral Montipora Facilitates Dual-Color Single-Laser Fluorescence Cross-Correlation Spectroscopy. Nat. Biotechnol. 2006, 24, 577−581. (7) Kogure, T.; Kawano, H.; Abe, Y.; Miyawaki, A. Fluorescence Imaging Using a Fluorescent Protein with a Large Stokes Shift. Methods 2008, 45, 223−226. (8) Piatkevich, K. D.; Hulit, J.; Subach, O. M.; Wu, B.; Abdulla, A.; Segall, J. E.; Verkhusha, V. V. Monomeric Red Fluorescent Proteins with a Large Stokes Shift. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 5369−5374. (9) Piatkevich, K. D.; Malashkevich, V. N.; Almo, S. C.; Verkhusha, V. V. Engineering ESPT Pathways Based on Structural Analysis of LSSmKate Red Fluorescent Proteins with Large Stokes Shift. J. Am. Chem. Soc. 2010, 132, 10762−10770. (10) Shcherbakova, D. M.; Hink, M. A.; Joosen, L.; Gadella, T. W. J.; Verkhusha, V. V. An Orange Fluorescent Protein with a Large Stokes Shift for Single-Excitation Multicolor FCCS and FRET Imaging. J. Am. Chem. Soc. 2012, 134, 7913−7923. (11) van Thor, J. J.; Gensch, T.; Hellingwerf, K. J.; Johnson, L. N. Phototransformation of Green Fluorescent Protein with UV and 12463


Article

The Journal of Physical Chemistry B Visible Light Leads to Decarboxylation of Glutamate 222. Nat. Struct. Biol. 2002, 9, 37−41. (12) Habuchi, S.; Cotlet, M.; Gensch, T.; Bednarz, T.; HaberPohlmeier, S.; Rozenski, J.; Dirix, G.; Michiels, J.; Vanderleyden, J.; Heberle, J.; et al. Evidence for the Isomerization and Decarboxylation in the Photoconversion of the Red Fluorescent Protein DsRed. J. Am. Chem. Soc. 2005, 127, 8977−8984. (13) Caricato, M.; Mennucci, B.; Tomasi, J. Solvent Effects on the Electronic Spectra: An Extension of the Polarizable Continuum Model to the ZINDO Method. J. Phys. Chem. A 2004, 108, 6248−6256. (14) Sanchez-Garcia, E.; Doerr, M.; Thiel, W. QM/MM Study of the Absorption Spectra of DsRed.M1 Chromophores. J. Comput. Chem. 2010, 31, 1603−1612. (15) Randino, C.; Nadal-Ferret, M.; Gelabert, R.; Lluch, M. M. J. M. A Time-Dependent DFT/Molecular Dynamics Study of the ProtonWire Responsible for the Red Fluorescence in the LSSmKate2 Protein. Theor. Chem. Acc. 2013, 132, 1327. (16) Schwabe, T.; Beerepoot, M. T. P.; Olsen, J. M. H.; Kongsted, J. Analysis of Computational Models for an Accurate Study of Electronic Excitations in GFP. Phys. Chem. Chem. Phys. 2015, 17, 2582−2588. (17) Isborn, C. M.; Götz, A. W.; Clark, M. A.; Walker, R. C.; Martínez, T. J. Electronic Absorption Spectra from MM and ab Initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein. J. Chem. Theory Comput. 2012, 8, 5092−5106. (18) Nadal-Ferret, M.; Gelabert, R.; Moreno, M.; Lluch, J. M. How Does the Environment Affect the Absorption Spectrum of the Fluorescent Protein mKeima? J. Chem. Theory Comput. 2013, 9, 1731− 1742. (19) Imhof, P. Computational Study of Absorption Spectra of the Photoconvertible Fluorescent Protein EosFP in Different Protonation States. J. Chem. Theory Comput. 2012, 8, 4828−4836. (20) Armengol, P.; Gelabert, R.; Moreno, M.; Lluch, J. M. Chromophore Interactions Leading to Different Absorption Spectra in mNeptune1 and mCardinal Red Fluorescent Proteins. Phys. Chem. Chem. Phys. 2016, 18, 16964−16976. (21) Chen, F.; Zeng, Q.; Zhuang, W.; Liang, W. Characterizing the Structures, Spectra, and Energy Landscapes Involved in the ExcitedState Proton Transfer Process of Red Fluorescent Protein LSSmKate1. J. Phys. Chem. B 2016, 120, 9833−9842. (22) Pedone, A.; Prampolini, G.; Monti, S.; Barone, V. Absorption and Emission Spectra of Fluorescent Silica Nanoparticles from TDDFT/MM/PCM Calculations. Phys. Chem. Chem. Phys. 2011, 13, 16689−16697. (23) Röhrig, U.; Frank, I.; Hutter, J.; Laio, A.; VandeVondele, J.; Röthlisberger, U. QM/MM Car-Parrinello Molecular Dynamics Study of the Solvent Effects on the Ground State and on the First Excited Singlet State of Acetone in Water. Chem. Phys. Chem. 2003, 4, 1177− 1182. (24) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671−690. (25) Vanommeslaeghe, K. mackerell.umaryland.edu/kenno/cgenff/ faq.php (accessed November 18, 2016). (26) Mayne, C. G.; Saam, J.; Schulten, K.; Tajkhorshid, E.; Gumbart, J. C. Rapid Parameterization of Small Molecules Using the Force Field Toolkit. J. Comput. Chem. 2013, 34, 2757−2770. (27) Humphrey, W.; Dalke, A.; Schulten, K. VMD − Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (29) Reuter, N.; Lin, H.; Thiel, W. Green Fluorescent Proteins: Empirical Force Field for the Neutral and Deprotonated Forms of the Chromophore. Molecular Dynamics Simulations of the Wild Type and S65T Mutant the Neutral and Deprotonated Forms of the

Chromophore. Molecular Dynamics Simulations of the Wild Type and S65T Mutant. J. Phys. Chem. B 2002, 106, 6310−6321. (30) 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, 361−373. (31) http://biophysics.cs.vt.edu/H++ (accessed March 4, 2016). (32) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: Automating pK Prediction and the Preparation of Biomolecular Structures for Atomistic Molecular Modeling and Simulations. Nucleic Acids Res. 2012, 40, W537−W541. (33) Myers, J.; Grothaus, G.; Narayanan, S.; Onufriev, A. A Simple Clustering Algorithm Can Be Accurate Enough for Use in Calculations of pKs in Macromolecules. Proteins 2006, 63, 928−938. (34) Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. H++: A Server for Estimating pKas and Adding Missing Hydrogens to Macromolecules. Nucleic Acids Res. 2005, 33, W368− W371. (35) http://www.charmm.org (accessed December 4, 2013). (36) 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, 187−217. (37) 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, 1781− 1802. (38) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (39) Harvey, J. N. Spin-Forbidden CO Ligand Recombination in Myoglobin. Faraday Discuss. 2004, 127, 165−177. (40) Ponder, J. W.; Richards, F. M. An Efficient Newton-Like Method for Molecular Mechanics Energy Minimization of Large Molecules. J. Comput. Chem. 1987, 8, 1016−1024. (41) Ridley, J. E.; Zerner, M. C. An Intermediate Neglect of Differential Overlap Technique for Spectroscopy: Pyrrole and the Azines. Theor. Chem. Acc. 1973, 32, 111−134. (42) Ridley, J. E.; Zerner, M. C. Triplet States via Intermediate Neglect of Differential Overlap: Benzene, Pyridine and the Diazines. Theor. Chim. Acta 1976, 42, 223−236. (43) Bacon, A. D.; Zerner, M. C. An Intermediate Neglect of Differential Overlap Theory for Transition Metal Complexes: Fe, Co and Cu Chlorides. Theor. Chim. Acta 1979, 53, 21−54. (44) Zerner, M. C.; Lowe, G. H.; Kirchner, R. F.; MuellerWesterhoff, U. T. An Intermediate Neglect of Differential Overlap Technique for Spectroscopy of Transition-Metal Complexes. Ferrocene. J. Am. Chem. Soc. 1980, 102, 589−599. (45) de Mello, P. C.; Zerner, M. C.; Hehenberger, M. Converging SCF calculations on excited states. Int. J. Quantum Chem. 1982, 21, 251−259. (46) Anderson, W. P.; Edwards, W. D.; Zerner, M. C. Calculated Spectra of Hydrated Ions of the First Transition-Metal Series. Inorg. Chem. 1986, 25, 2728−2732. (47) Hanson, L. K.; Fajer, J.; Thompson, M. A.; Zerner, M. C. Electrochromic Effects of Charge Separation in Bacterial Photosynthesis: Theoretical Models. J. Am. Chem. Soc. 1987, 109, 4728−4730. (48) Thompson, M. A.; Zerner, M. C. A Theoretical Examination of the Electronic Structure and Spectroscopy of the Photosynthetic Reaction Center from Rhodopseudomonas viridis. J. Am. Chem. Soc. 1991, 113, 8210−8215. (49) Zerner, M. C. Semiempirical Molecular Orbital Methods. Rev. Comput. Chem. 1991, 2, 313−366. (50) Yanai, T.; Tew, D.; Handy, N. A New Hybrid Exchange− Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (51) Neese, F. The ORCA Program System. WIREs Comput. Mol. Sci. 2012, 2, 73−78. 12464


Article

The Journal of Physical Chemistry B (52) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (53) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; et al. MOLPRO, version 2012.1, a package of ab initio programs, 2012. http://www.molpro.net. (54) Grant, B. J.; Rodrigues, A. P.; ElSawy, K. M.; McCammon, J. A.; Caves, L. S. Bio3d: An R Package for the Comparative Analysis of Protein Structures. Bioinformatics 2006, 22, 2695−2696. (55) Schrödinger, L. The PyMOL Molecular Graphics System, version 1.8; Schrödinger, LLC, 2015. (56) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. Benchmarks for Electronically Excited States: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, No. 134110. (57) Jacquemin, D.; Wathelet, V.; Perpete, E. A.; Adamo, C. Extensive TD-DFT Benchmark: Singlet-Excited States of Organic Molecules. J. Chem. Theory Comput. 2009, 5, 2420−2435. (58) Silva-Junior, M. R.; Sauer, S. P.; Schreiber, M.; Thiel, W. Basis Set Effects on Coupled Cluster Benchmarks of Electronically Excited States: CC3, CCSDR(3) and CC2. Mol. Phys. 2010, 108, 453−465. (59) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks of Electronically Excited States: Basis Set Effects on CASPT2 Results. J. Chem. Phys. 2010, 133, No. 174318. (60) Jacquemin, D.; Perpete, E. A.; Ciofini, I.; Adamo, C. Assessment of the ωB97 Family for Excited-State Calculations. Theor. Chem. Acc. 2011, 128, 127−136. (61) Jacquemin, D.; Mennucci, B.; Adamo, C. Excited-State Calculations with TD-DFT: From Benchmarks to Simulations in Complex Environments. Phys. Chem. Chem. Phys. 2011, 13, 16987− 16998. (62) Vydrov, O. A.; Scuseria, G. E. Assessment of a Long-Range Corrected Hybrid Functional. J. Chem. Phys. 2006, 125, No. 234109. (63) Vydrov, O. A.; Heyd, J.; Krukau, A.; Scuseria, G. E. Importance of Short-Range versus Long-Range Hartree-Fock Exchange for the Performance of Hybrid Density Functionals. J. Chem. Phys. 2006, 125, No. 074106. (64) Vydrov, O. A.; Scuseria, G. E.; Perdew, J. P. Tests of Functionals for Systems with Fractional Electron Number. J. Chem. Phys. 2007, 126, No. 154109. (65) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (66) Helms, V.; Winstead, C.; Langhoff, P. Low-Lying Electronic Excitations of the Green Fluorescent Protein Chromophore. J. Mol. Struct.: THEOCHEM 2000, 506, 179−189. (67) Epifanovsky, E.; Polyakov, I.; Grigorenko, B.; Nemukhin, A.; Krylov, A. I. Quantum Chemical Benchmark Studies of the Electronic Properties of the Green Fluorescent Protein Chromophore. 1. Electronically Excited and Ionized States of the Anionic Chromophore in the Gas Phase. J. Chem. Theory Comput. 2009, 5, 1895−1906. (68) Filippi, C.; Zaccheddu, M.; Buda, F. Absorption Spectrum of the Green Fluorescent Protein Chromophore: A Difficult Case for ab Initio Methods? J. Chem. Theory Comput. 2009, 5, 2074−2087. (69) Send, R.; Kaila, V. R. I.; Sundholm, D. Benchmarking the Approximate Second-Order Coupled-Cluster Method on Biochromophores. J. Chem. Theory Comput. 2011, 7, 2473−2484. (70) List, N. H.; Olsen, J. M.; Rocha-Rinza, T.; Christiansen, O.; Kongsted, J. Performance of Popular XC-Functionals for the Description of Excitation Energies in GFP-Like Chromophore Models. Int. J. Quantum Chem. 2012, 112, 789−800. (71) Daday, C.; Curutchet, C.; Sinicropi, A.; Mennucci, B.; Filippi, C. J. Chem. Theory Comput. 2015, 11, 4825. (72) De Keersmaecker, H.; Fron, E.; Rocha, S.; Kogure, T.; Miyawaki, A.; Hofkens, J.; Mizuno, H. Photoconvertible Behavior of LSSmOrange Applicable for Single Emission Band Optical Highlighting. Biophys. J. 2016, 111, 1014−1025.

12465


MM-Based Calculations of Absorption and Emission Spectra of

Recommend Documents