Local Electric Field Controls Fluorescence Quantum Yield of Red and Far-Red Fluorescent Proteins

Genetically encoded probes with red-shifted absorption and fluorescence are highly desirable for imaging applications because they can report from deeper tissue layers with lower background and because they provide additional colors for multicolor imaging. Unfortunately, red and especially far-red fluorescent proteins have very low quantum yields, which undermines their other advantages. Elucidating the mechanism of nonradiative relaxation in red fluorescent proteins (RFPs) could help developing ones with higher quantum yields. Here we consider two possible mechanisms of fast nonradiative relaxation of electronic excitation in RFPs. The first, known as the energy gap law, predicts a steep exponential drop of fluorescence quantum yield with a systematic red shift of fluorescence frequency. In this case the relaxation of excitation occurs in the chromophore without any significant changes of its geometry. The second mechanism is related to a twisted intramolecular charge transfer in the excited state, followed by an ultrafast internal conversion. The chromophore twisting can strongly depend on the local electric field because the field can affect the activation energy. We present a spectroscopic method of evaluating local electric fields experienced by the chromophore in the protein environment. The method is based on linear and two-photon absorption spectroscopy, as well as on quantum-mechanically calculated parameters of the isolated chromophore. Using this method, which is substantiated by our molecular dynamics simulations, we obtain the components of electric field in the chromophore plane for seven different RFPs with the same chromophore structure. We find that in five of these RFPs, the nonradiative relaxation rate increases with the strength of the field along the chromophore axis directed from the center of imidazolinone ring to the center of phenolate ring. Furthermore, this rate depends on the corresponding electrostatic energy change (calculated from the known fields and charge displacements), in quantitative agreement with the Marcus theory of charge transfer. This result supports the dominant role of the twisted intramolecular charge transfer mechanism over the energy gap law for most of the studied RFPs. It provides important guidelines of how to shift the absorption wavelength of an RFP to the red, while keeping its brightness reasonably high.


INTRODUCTION
Red fluorescent proteins (RFPs) and biosensors derived from them present an important addition to a rich palette of genetically encoded fluorescent probes widely employed in bioimaging (Tsien et al., 1998;Shaner et al., 2005;Day and Davidson, 2009;Wiedenmann et al., 2009;Chudakov et al., 2010). Their red-shifted absorption and fluorescence make it possible to report from deeper layers of tissues with less background autofluorescence compared to green fluorescent proteins (GFPs). The spectral red shifts in RFPs are due to a longer π-conjugated system in the chromophore structure, that includes an additional N-acylimine group appended to a common GFP chromophore, p-hydroxybenzylidineimidazolinone. Although the fluorescence quantum yield (QY) of the first discovered tetrameric red FP, DsRed, was quite large (∼0.7), the most popular red-shifted monomeric mutant variants carrying the same chromophore, mPlum and mCherry, show much lower QYs: 0.1-0.2 (Shaner et al., 2005). RFPs that are further red-shifted, such as eqFP670 and mGrape3, are even dimmer, with quantum yields of ∼0.05 (Lin et al., 2009;Shcherbo et al., 2010); see Figure 1. A number of red genetically encoded calcium indicators (GECIs) also show low QYs (∼0.2) even in their active fluorescent state (Molina et al., 2019).
Fast nonradiative relaxation leading to a low QY in the reddest variants of RFPs can be due to different mechanisms. First, the shift of the emission frequency of a molecule to the red (due to chemical modifications or interactions with environment) often results in acceleration of the internal conversion, following the "energy gap law" (Englman and Jortner, 1970;Jung et al., 2012). This mechanism is quite general and simply reflects the dependence of vibrational states distribution and coupling between them (through Franck-Condon factors) on the energy gap between the two electronic terms. Another mechanism of fast relaxation involves twisted intramolecular charge transfer (TICT) that becomes possible due to structural flexibility of the chromophore in some FP variants. In this case, the rotation about one or both of the bridge methine bonds (i.e., phenolate, P, or imidazolinone, I) drives the molecular system to a twisted state that corresponds to a conical intersection of the excited and ground states potential energy surfaces. Once in this state, the molecule undergoes an ultrafast transition from the excited to the ground state. Twisting in the excited state occurs in concert with significant charge transfer (CT) across the bridge from one side of the chromophore to another (Olsen and Smith, 2007;Simine et al., 2018;Park and Rhee, 2016;Sun et al., 2012, Moron et al., 2019. In the isolated anionic GFP chromophore, such TICT states are quasi-stable intermediates on the excited state potential energy surfaces for either I-or P-rotations (Martin et al., 2004;Altoe et al., 2005). Those states are close to conical intersection seams, but still are separated from the ground state surface by small gaps (Martin et al., 2004;Altoe et al., 2005). In contrast, twisting of the isolated anionic RFP chromophore around the P-bond leads directly to a conical intersection seam at twisting angles of ∼75°-90° (Olsen and Smith, 2007). Like in the GFP chromophore, this rotation is virtually barrierless and exergonic. In contrast, twisting around the I-bond has a barrier and is endergonic. Olsen and Smith suggested that, in contrast to the GFP chromophore, in the RFP chromophore the electronegativity of the acylimine (A) substituent plays a decisive role in selecting the P-pathway of the TICT process.
As for the chromophore inside a protein, several factors can impede its ultrafast relaxation along the TICT pathway (Tolbert et al., 2012;Jung et al., 2012). Steric clashes with the surrounding bulky groups is the most obvious one. These interactions are usually much more efficient for rotation around the I-bond and less so for the P-bond (Chudakov et al., 2003;Stiel et al., 2008) because of the larger moving volume in the former case. A volume-conserving hula-twist motion involving concert rotation around both exocyclic bonds was put forward to explain low quantum yields in some GFP mutants (Jung et al., 2005) and the RFP mPlum (Moron et al., 2019), although this mechanism was questioned for GFP in (Altoe et al., 2005).
The effect of electrostatic interactions of the chromophore with the protein surrounding (including hydrogen bonding) can also contribute to the dynamics of TICT. In one scenario, a stronger electric field directed from P to I would shift the phenolate π-conjugation resonance from a quinonoid to a benzenoid form. This would facilitate rotation around the P-bond because it becomes closer to single in character. This so-called electronic effect of controlling nonradiative relaxation was experimentally observed in certain GFP mutants. Variants having less hydrogen bonding of the phenolate oxygen (i.e., more quinonoid structure) generally showed an increased quantum yield and fluorescence lifetime compared to mutants with more hydrogen bonds (benzenoid structure) (Jung et al., 2005;Jung et al., 2012). Alternatively, if the charge transfer corresponding to a shift of electronic density from P to I upon P-rotation is a more important factor than the bond order, FIGURE 1 | Fluorescence spectra of a set of RFPs studied in this work. The amplitude of the fluorescence intensity is normalized such that the integral under the curve is proportional to the fluorescence quantum yield. than applying the field from I to P would speed up the rotation through an electrostatic driving force. This particular mechanism was theoretically predicted for GFPs (Park and Rhee, 2016;Simine et al., 2018) and RFPs (Olsen and Smith, 2007;Sun et al., 2012;Moron et al., 2019), however, it was not experimentally confirmed yet. Thus, it is crucial to experimentally measure the protein internal field. Such measurements however present a serious challenge, especially in the case of the two-dimensional RFP chromophore, in which both components of the field E projected onto the molecular axes (E x and E y ) can be important.
In general, any of the above mechanisms, including vibrational relaxation described by the energy gap law and chromophore twisting (possibly sensitive to the electric field) can contribute to a fast nonradiative relaxation in RFPs (Figure 1). Our goal is to understand which of these processes is most important, if any. To this end, we need to evaluate the electric fields created by the protein environment at the chromophore site. A correlation (or the absence of one) between the nonradiative relaxation rate (k nR ) and the amplitude and direction of the local electric field for a set of RFP variants can reveal the mechanism of relaxation. Is a redder shift inextricably connected to a faster nonradiative relaxation? We aim to answer this important question, and, eventually, reveal clues as to how to find RFP variants with large quantum yields.

Identification and Cloning of XRFP
XRFP was identified as part of a large survey of published raw mRNA-Seq data sets in the NCBI Sequence Read Archive (SRA) database. Transcriptomes were assembled using Trinity (Grabherr et al, 2011;Haas et al., 2013) and using the public Galaxy bioinformatics server (Afgan et al., 2018). Candidate FPencoding transcripts were identified by BLAST homology searching using avGFP as the query against the assembled transcriptome database. For each FP homolog found, the coding region was identified and a synthetic gene was designed to produce the encoded polypeptide sequence using codons optimized for Escherichia coli expression using an in-house BioXP3200 instrument (SGI-DNA, La Jolla, CA) or ordered as a gBlock double-stranded gene fragment (Integrated DNA Technologies, San Diego, CA). Fragments encoding FPs were inserted using Gibson assembly (Gibson et al., 2009) into the vector pNCS for expression (see below). XRFP was identified along with green-and orange-emitting FPs from the animal (data not shown).

Linear Absorption and Fluorescence Properties
Linear absorption spectra were measured with a Lambda 950 spectrophotometer (Perkin Elmer). To obtain the peak extinction coefficient of an anionic chromophore, we performed a stepwise alkaline titration up to complete denaturation of the protein. At each step, 10 μl of 0.1 M NaOH was added to 0.5 ml of protein solution, held in a spectroscopic cuvette with 1-cm optical path, followed by recording of absorption spectrum. To correct for different protein concentrations at each step, absorption spectrum was multiplied by a factor equal to the ratio of total volume to initial volume of solution. Usually a dependence of optical density of the native anionic chromophore peak versus optical density of the denatured protein peak (at 452 nm) shows a linear region. In this region of pH, there are only two species present: the anionic chromophore in native protein and denatured protein (in the same region the series of absorption spectra show an isosbestic point). The slope of this linear dependence is equal to the ratio of extinction coefficients of the two species. Using known extinction value of denatured protein, 44,100 M −1 cm −1 at 452 nm (Ward, 2005;Gross et al., 2000), we obtain the extinction coefficient of anionic species. In the cases where the immature green chromophore was present initially, its contribution to the final denatured protein concentration was subtracted.
Corrected fluorescence spectra and fluorescence quantum yields were measured with an integrating sphere fluorometer (Quantaurus-QY, Hamamtsu) using 1-cm quartz cuvettes. The peak OD of the samples was <0.1. The reference (buffer-only) measurements were done in the same cuvette.
Fluorescence lifetime was measured with a Digital Frequency Domain system ChronosDFD (ISS) appended to a PC1 ISS spectrofluorometer. Peak optical density of the samples and reference solutions in 1-cm cuvettes was kept below 0.1. Fluorescence was excited with a 518-nm laser diode (ISS, model 73292). The excitation was modulated with multiple harmonics in the range 10-300 MHz for eqFP670 and mPlum, and 5-150 MHz for the rest of the proteins. In this method, a fluorescence lifetime standard is used to obtain the instrumental response function in each individual measurement. For cross checking, we employed Rose Bengal (Sigma-Aldrich) in ethanol with τ 0.78 ns (Fleming et al., 1977;Cramer and Spears, 1978) or methanol, τ 0.55 ns, (Fleming et al., 1977;Cramer and Spears, 1978;Reed et al., 1981;Rodgers, 1981, Lakowicz et al., 1986; Rhodamine 6G in deionized water, τ 4.0 ns (Reisfeld et al., 1988;Magde et al., 2002;O'Hagan et al., 2001); and Rhodamine B in deionized water, τ 1.74 ns (Boens et al., 2007). For each RFP, we used two or three different standards and the results were similar (deviations <6%). Fluorescence of all the samples and standards was collected at 90°through two identical 561LP Edge Basic TM filters (Semrock) to cut off all excitation light. In the case of eqFP670 and mCherry (pH7.4), we used an additional HQ 650/20 (Chroma) filter to selectively collect fluorescence from the red-shifted form. For mCherry at pH 11.4, we used an additional HQ 577/10 filter (Chroma) to selectively collect fluorescence of the blue-shifted form. All the modulation ratio and phase delay curves were fitted to model functions corresponding to a single exponential fluorescence decay. The corresponding χ 2 values were in the range from 0.53 to 1.16. Corrected fluorescence excitation spectra were obtained with a LS-55 spectrofluorimeter (Perkin Elmer). For mCherry at pH 11.4, the registration wavelength was 575 nm, where the redshifted form does not fluoresce. For all other proteins, the registration wavelength was selected at the red side of the fluorescence spectrum near its maximum.

Two-photon Absorption Spectra, Cross Sections, and Polarization Ratio
Two-photon excitation spectra and cross sections were measured as described previously (Drobizhev et al., 2020). Briefly, we used an Insight DeepSee femtosecond laser (Spectra Physics) tunable in the 680-1,300 nm range coupled with a PC1 photon counting spectrofluorometer (ISS). LDS 798 in CDCl 3 :CHCl 3 (2:1) solution was used as a spectral shape standard and Rhodamine 6G in MeOH was used as a 2PA cross section standard, with σ 2 10 GM at 1060 nm. The laser beam was focused into the sample with a NIR achromatic lens, f 45 mm (Edmund Optics). The sample solutions were held in 3-mm spectroscopic cuvettes (Starna) and fluorescence was collected from the first 0.7-mm layer of solution, to avoid laser absorption by the solvent. A 745SP filter (Semrock) was placed after the sample to cut the scattered light. For mCherry pH 11.4, an additional HQ 577/10 filter (Chroma) was used.
To evaluate the two-photon polarization ratio Ω, the same experimental system was used. Additionally, a quarter-wave plate (Thorlabs) was placed in front of the entrance diaphragm of the fluorometer and the Glan polarizer (ISS) was set after the sample in the detection optical path. The quarter-wave plate was mounted on a rotation stage (Thorlabs) which made it possible to rotate it around both the laser propagation direction and the vertical axis. This helped to adjust the rotation and tilt angles to make the polarization very close to circular at each laser wavelength. The ellipticity of polarization was <7% in the tuning range 900-1,300 nm. By definition, Ω is the ratio of 2P cross sections obtained under circular and linear polarization of excitation. Experimentally, we collect fluorescence at 90°to excitation, first using circularly polarized light and then linearly (vertically) polarized light. To get rid of the effect of fluorescence polarization, at each polarization of laser, two measurements were done: one with the detection polarizer set vertically, and another with the detection polarizer set horizontally. The two-photon polarization ratio was then calculated as follows (Wan and Johnson, 1994): where F is the fluorescence signal, indexes O and V describe circular and vertical polarizations of excitation, respectively, and indexes V and H describe vertical and horizontal positions of detection polarizer. The sensitivity of the detector (PMT R928) for vertical versus horizontal polarization of fluorescence was checked by exciting the sample with horizontally polarized laser light and the ratio of the two signals (G-factor) was G 1.0. Standard deviations of Ω were estimated from 10 measurements of Ω in the same conditions at a number of different wavelengths.

MD Parameters
Force field parameters for the mCherry, mPlum, and DsRed chromophores were those published in (Dmitrienko et. al., 2006) with the exception of excluding the phenolic hydroxyl hydrogen in order to form the anionic chromophore. The optimized potentials for liquid simulations (OPLS) amino acid atom types suggested by Dmitrienko et. al. were also used. Standard OPLS-AA atom types were used for the extended chain atoms. A variant of the mCherry chromophore residue force field was also constructed, wherein equilibrium bond lengths and angles of the force field were set to those found in the 2H5Q pdb crystal structure rather than those suggested by Dmitrienko et al. for the neutral DsRed chromophore.
For both the DsRed and mCherry/mPlum chromophores, hydrogen addition rules were added to the OPLS-AA force field in such a way as to complete typical valency requirements for the anionic chromophores.

Simulation Details
All MD simulations were done with Gromacs-5.0.2 (Pronk et al., 2013). Initial atomic charges were found by protonating the pdb structure of a given mFruit protein using Gromacs, embedding it in a box of transferable intermolecular potential with 3 points (TIP3P) solvent, adding sodium ions in sufficient number to balance the charge of the protein, and then subjecting the chromophore to a semi-empirical INDO/S2 (ZINDO) analysis in the presence of the electric field due to non-chromophore protein atoms, water molecules, and ions in the simulation cell. Initial charges were calculated in this way for mCherry, mPlum, and DsRed, with a separate initial charge set determined for mCherry with its Glu215 residue protonated.
Following determination of an initial set of atomic charges for each mFruit protein chromophore based on the electrostatic environment, simulations were carried out. Each protein/ Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 simulation was treated with 100 ps of equilibration under the NVT ensemble with the protein heavy atoms constrained by a large harmonic restoring force, followed by 400 ps of NPT equilibration using the Berendsen barostat and a further 500 ps using the Parrinello-Rahman barostat, also with heavy protein atoms constrained. Finally, protein constraints were removed and 100 ns of dynamics was carried out, with positions saved each picosecond for a total of 100001 frames. The simulation was split into individual coordinate files, which were used for all subsequent analysis. For each simulation frame a ZINDO calculation was performed, which yielded excitation energy, wavelength, oscillator strength, transition dipole, and Δμ for the 59 lowest energy transitions. All atoms in the simulation cell that were not part of the quantum chromophore were used to calculate the electric field vector and electric potential at each chromophore atom, so that the electrostatic environment could be completely accounted for in the ZINDO calculation. The output of the ZINDO calculations was also used to determine the S 0 and S 1 state charges for each simulation frame. The S 1 -S 0 charge differences for each simulation were averaged, and those average values were used along with the atomic coordinates for each simulation frame to calculate the shift in chromophore excitation energy due to each protein atom that was not part of the quantum chromophore as well as for each water molecule and ion. This straightforward Coulombic summation allowed for the unambiguous assignment of Stark shifts due to protein, water and ions for each simulation frame as well.
Following the first set of 100 ns simulations based on initial MD charges, S 0 charge values for each quantum chromophore atom were plotted to assess charge convergence. The initial charges based on the solvated pdb structures were found to be quite close to those calculated for the chromophore during dynamics but were sufficiently different to warrant a second charge iteration. Average charge values for each MD chromophore atom were extracted from the above described plots, with an attempt made to identify portions of the trajectory where a given charge appeared to be approximately constant. These new MD chromophore atom charge values were used to begin a second set of 100 ns simulations of each RFP. Data for each simulation was analyzed using the ZINDO technique as described above, and chromophore charges during dynamics were found to be in excellent agreement with the static, assigned, second generation MD charges. After the 100 ns second charge iteration simulations, dynamics were continued for another 100 ps with frames collected every 2 fs. This high resolution data was analyzed in the same fashion as the coarser data sets described above.

Role of Radiative and Nonradiative Relaxation Rates in Controlling Fluorescence Quantum Yield
We selected seven RFPs with a large variation of absorption/ fluorescence peak wavelengths and fluorescence quantum yields.
Since the extinction coefficient does not change very much in the series, the key parameter determining molecular brightness (ε max φ) is the fluorescence quantum yield (φ). The quantum yield decreases about an order of magnitude when going from DsRed2 and mScarlet to the most red-shifted variant eqFP670, Table 1. The quantum yield and fluorescence lifetime depend on both the radiative relaxation rate (k R ) and nonradiative relaxation rate (k nR ): To understand the role of k R and k nR in controlling quantum yield, we calculated them from the measured φ and τ, and present them in Table 1. The dependence of k R on the cube of fluorescence peak frequency, ] f (in cm −1 ), is shown in Figure 2.
The dashed line is a linear fit representing the Einstein coefficient for spontaneous emission as a function of ] 3 f , Here n is the refractive index of the medium and h is the Planck constant. The slope of the fit is proportional to the matrix element of the transition dipole moment squared μ em 2 . The slope provides an average transition dipole moment, |μ em | (7.1 ± 0.2) D. In general, the k R value does not change more than twofold in the series. On the other hand, k nR varies about 20 times; see Table 1. We thus conclude that k nR plays a decisive role in controlling quantum yield.

Checking the Energy Gap Law
The energy gap law predicts an exponential decrease of k nR with the increase of the energy gap ΔE for fluorescence transition between the energy levels of S 1 and S 0 states (Englman and Jortner, 1970): Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 where C is a constant factor, ξ is a parameter that only slowly varies with ΔE and ] M is the maximum molecular normal vibration frequency (in Hz). Figure 3 shows the dependence of nonradiative relaxation rate (in logarithmic scale) on the fluorescence peak frequency ] f ΔE/ch for the series of RFPs studied here (red circles). Although the RFPs qualitatively follow the predicted dependence Eq. 7, the correlation is not very strong (Pearson's R −0.858). In particular, although the peak fluorescence frequencies of mScarlet and mCherry pH 11.4 are very close, their k nR values differ more than twofold, a difference much larger than the experimental errors. This is similarly true for the XRFP and mPlum pair. These quantitative inconsistencies suggest that the internal conversion through vibrionic coupling between S 1 and S 0 states, reflected in the energy gap law, is probably not the dominant mechanism of relaxation. Another, indirect support for a failure of this mechanism is obtained by considering the behavior of k nR in a series of proteins with green anionic chromophore, shown in Figure 3 by green triangles. The mutants and homologues of GFP of this series were characterized previously (Drobizhev et al., 2011;Molina et al., 2017). They were selected such that they all have high quantum yields, 0.67-0.91, (to exclude other possible deactivation mechanisms, e.g., twisting around the bridge bonds), and their fluorescence spectra span a broad range, from 482 (Rosmarinus) to 537 nm (phiYFP). Due to a close similarity in the structures of the green and red chromophores, as well as closeness in shape of the absorption and fluorescence spectra in the two series, one would expect that the energy gap law, if dominant in controlling k nR , would be observed in both series. However, the "green" set does not show any correlation between k nR and ] fl . We therefore conclude that the energy gap law is not a main mechanism of relaxation in the studied The columns show, in order, maximum absorption wavelength, maximum extinction coefficient, maximum fluorescence wavelength, fluorescence quantum yield, fluorescence lifetime, radiative relaxation rate, and non-radiative relaxation rate. In cases where multiple forms of chromophore were present in solution (e.g. red and green immature forms) we provide the extinction coefficient for a major (red) form. a Obtained in (Drobizhev et al., 2011) using the Strickler-Berg formula.
FIGURE 2 | Dependence of radiative relaxation rate on cube of maximum fluorescence frequency. Dashed line represents the best fit to the Einstein equation for spontaneous emission rate.
FIGURE 3 | Dependence of nonradiative decay rate (log scale) on peak fluorescence frequency. Red circles correspond to RFPs studied here and green triangles correspond to a number of GFPs variants, studied previously (Drobizhev et al., 2011;Molina et al., 2017). The dashed lines represent linear fits to the two sets of data.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 set of RFPs (and GFPs) and will consider alternative mechanism in Section "Describing the Rate of TICT Process with Marcus Formalism". To this end, we must first investigate the internal electric field created by the protein at the chromophore site.
Physical Model for Determining the Protein Internal Electric Field at the Chromophore Site The protein surrounding the chromophore creates an internal electric field that may play an important role in the mechanism of nonradiative relaxation. To determine the components of this field in the plane of the chromophore, we developed a physical model based on the following assumptions: 1) The local electric field does not change neither upon electronic excitation of the chromophore (environment is polarizable only in the ground state) nor upon twisting around exocyclic P-bond in the excited state.
2) The potentially non-homogeneous field varying from one atom to another on the chromophore will produce the same effect as a homogeneous effective field E. The first assumption was recently proven to a first approximation for DsRed (List et al., 2012) and a series of GFP variants (Nifosi et al., 2019). The effective field obtained after averaging the fields on several atoms of the chromophore's π-conjugation system was also shown to describe optical properties reasonably well (Nifosi et al., 2019). In our approach of evaluating the local electric field, we use the permanent dipole moment of the chromophore μ (g,e) , where indices g and e correspond to the ground and excited states, respectively, as a linear metric of the field. In fact, if an electric field of a large magnitude (|E| ∼10 7 -10 8 V/cm) is applied to a dipolar molecule with nonzero polarizability α, its dipole moment μ (g,e) will acquire an additional, induced, part equal to μ (g,e) ind αE comparable to the vacuum dipole moment μ (g,e) 0 . (The underline denotes the tensorial nature of polarizability.) Experiments and calculations show that such strong fields are indeed present in proteins (Geissinger et al., 1995;Manas et al., 1999;Schweitzer-Stenner, 2008;Callis and Burgess, 1997;Vivian and Callis, 2001). Therefore, for the linear in the field approximation, the total dipole moment reads (Atkins and Friedman, 1997): For simplicity, we consider the chromophore to be planar in the ground electronic state, although small twist and tilt angles are present in some RFPs (Shu et al., 2006). We also assume here that μ (e) and α (e) correspond to an excited state with positions of the atomic nuclei unchanged compared to the ground state. Applying Eq. 8 separately to the ground and the excited state and then subtracting the former from the latter, we obtain: We now select x and y coordinate axes such that they correspond to the main axes of the 2x2 tensor of the polarizability change, Δα, i.e., the frame where it is diagonal with components Δα xx and Δα yy . In the same coordinate frame, the Δμ μ (e) − μ (g) vector has components Δμ x and Δμ y ; Δμ 0,x and Δμ 0,y are the components of Δμ 0 . Equation 9 can be projected onto the x and y coordinate axes resulting in Δμ y Δμ 0,y + Δα yy E y .
(Note that in our previous papers (Drobizhev et al., 2009;Drobizhev et al., 2012a) the coefficient ½ was used in the second term of the right-hand side of Eq. 9 following an erroneous presentation in some previous literature, see e.g. (Berlin et al., 2006). A correct formula (Atkins and Friedman, 1997) does not contain it.) If Δμ x and Δμ y were known, then it would be possible to calculate the field components E x and E y , by solving Eqs 10, 11: Now, we will demonstrate how to obtain Δμ x and Δμ y values for a chromophore within a protein, using one-and two-photon absorption spectroscopy. In our method we also rely on the quantum-mechanically calculated parameters of the isolated chromophore: transition frequency, ] 0 , and components of polarizability tensor, Δα xx and Δα yy , and dipole moment vector, Δμ 0,x and Δμ 0,y . This approach allows us to limit the size of the molecular system that has to be calculated quantummechanically to the chromophore group. We can then evaluate the parameters sensitive to the protein environment purely from experiment.
Two experimental values have to be measured in the region of the pure electronic S 0 -> S 1 transition: 1) the one-photon absorption frequency ] and 2) the change of the permanent dipole moment upon excitation, Δμ |Δμ|. (Here and throughout we call the pure electronic, or 0-0, transition, the one that does not involve excitation of high frequency vibrations, ] v > 1,000 cm −1 ). The first value can be obtained from the onephoton absorption (1PA) spectrum. The second requires the two-photon absorption (2PA) spectrum in absolute crosssection values, as well as the two-photon polarization ratio (i.e., the ratio of fluorescence signals obtained upon circular and linear two-photon excitation).
Due to the Stark effect, the chromophore experiences a spectral shift of its one-photon absorption transition frequency from ] 0 in vacuum to ] inside the protein. Since we are dealing with strong fields such that Δμ itself depends on the field, we should include the quadratic terms in the field dependence (quadratic Stark effect). The optical transition energy for the chromophore in protein will read (Atkins and Friedman, 1997): We use the e.s.u. system of units here and throughout, unless specified otherwise. Frequency ] is expressed in cm −1 , h is the Planck constant, and c is the speed of light in vacuum. Substituting Eq. 12 and Eq. 13 into Eq. 14, we obtain the relation between the transition energy and the Δμ 0,x , Δμ 0,y , Δμ x , and Δμ y components: To find Δμ x and Δμ y , it is convenient to re-group Eq. 15 as follows: Equation 16 represents a conic section in the (Δμ x , Δμ y ) coordinate plane, which, depending on the signs of Δα xx and Δα yy , could be either an ellipse or a hyperbola (if Δα xx ≠ Δα yy ).
The curve is fully determined if the value of ] is measured and Δμ 0,x , Δμ 0,y , Δα xx , Δα yy , and ] 0 are calculated quantum mechanically.
The second equation for finding the Δμ x and Δμ y components comes from the expression of the 2PA cross section, deduced from the second order perturbation theory in the two-level approximation of the S 0 → S 1 transition (Drobizhev et al., 2009) for linearly polarized degenerate excitation. The two-level approximation was justified by quantum mechanical calculations for the RFP chromophore (Drobizhev et al., 2009). The 2PA cross section corresponding to the pure electronic peak, σ 2 (0-0), depends on several factors. These are: Δμ 2 , the extinction coefficient ε(0-0) (in M −1 cm −1 ), transition frequency ] (we assume that the one-photon absorption peak coincides with the 0-0 transition), and c, the angle between the transition dipole moment μ and the change of the permanent dipole moment Δμ : where the factor A is equal to A 4 10 3 π ln 10f 2 opt 5hc 2 nN A .
Here, N A is the Avogadro's number, f opt is the local field factor at optical frequency ]/2, which is usually assumed to be of a Lorentz form (Bloebmergen, 1965): f opt (n 2 +2)/3. The angle c can be obtained experimentally by comparing the 2PA cross sections measured under circularly (σ M 2 ) and linearly (σ h 2 ) polarized light. Within the two-level approximation of the degenerate 2PA process, the ratio of the two measurements Ω depends on c as follows (Meath and Power, 1984): As one can see, the polarization ratio spans the range between 2/3 and 3/2 when c varies from 0°to 90°. Solving Eq. 19 for c results in c ± arccos 2(1 − Ω) 4Ω − 1 + πn; n 0, 1, 2, . . .
Solving Eq. 17 for Δμ 2 and noticing that Δμ 2 Δμ 2 x + Δμ 2 y , we get Equation 21 presents a circle in the (Δμ x , Δμ y ) coordinate plane with the radius equal to The curve Eq. 21 is fully determined if the values of ], σ 2 (0−0), ε(0−0), and c are measured and the constant A is calculated according to Eq. 18. Generally, the system of two conical sections Eq. 16 and Eq. 21 can have up to four possible solutions for the vector Δμ (Δμ x , Δμ y ). In Section "Electric Fields in RFPs" we will show how independent knowledge of the direction of the transition dipole moment in the molecular frame, as well as MD simulations can help to select the right solution.

Parameters of the Isolated RFP Chromophore Calculated Quantum Mechanically
To obtain the Δμ 0,x , Δμ 0,y , Δα xx , and Δα yy values, we started from the available structure of the mCherry protein (Shu et al., 2006). We cut the chromophore group (residues 65 and 66 in 2H5Q pdb file) out of the protein, added the hydrogen atoms to it, and optimized its geometry using Gaussian 09 (Frisch et al., 2016) and the B3LYP/6-311++G(d,p) density functional and basis set (Figure 4, bottom). In the optimization process, we forced the I and P rings to lie in the same plane (i.e., flat chromophore). We performed geometry optimization for a chromophore in a set of electric field values applied separately along each axis, x and y. Molecular system of coordinates was selected such that the change of polarizability tensor Δα is diagonal in it, as described before (Drobizhev et al., 2012a), see Figure 4, bottom. Using a version of Zerner's INDO/S2-CIS ("ZINDO") method (Ridley and Zerner, 1973) modified to add oxygen parameters suggested by the Truhlar group (Li et al., 1999) we calculated Δμ x and Δμ y at each value of the field for the vertical optical transition. (Note that the TD-DFT often underestimates the Δμ values (Jacquemin, 2016) and ZINDO performs better than TD-DFT in calculating the RFP chromophore photophysical properties (Schaefer et al., 2007)). Figure 4 shows the dependences of Δμ x and Δμ y on E x (a) and E y (b) obtained for the RFP chromophore in vacuum. We find that at zero field, Δμ 0,x 5.53D and Δμ 0,y −3.17D. The signs of the two components suggest that at least for the isolated chromophore, the electronic density shifts from P to I and also from I to A part upon excitation (by definition, the electric dipole moment is directed from negative charge to positive charge). The dependences of Δμ components on the field were fitted to linear functions (Eqs 10,11) in the ranges of the field E x 0 ÷ 20 MV/cm and E y 0 ÷ 50 MV/cm, typical for the set of RFPs studied here (vide infra) with fixed zero-field values Δμ 0,x and Δμ 0,y . The slopes of these regressions constitute the elements of the Δα tensor: Δα xx −61.3 ± 0.2 Å 3 and Δα yy 6.9 ± 0.3 Å 3 . As expected, the main component of Δα is directed along the molecular axis encountering the greater number of π-electrons, i.e., the line connecting the centers of two π-conjugated rings. The minor component is due to the A part and/or parts of conjugated system of P and I extended in a direction perpendicular to x. The vacuum transition frequency ] 0 of the red FP chromophore was previously obtained in (Taguchi et al., 2009) using fragment molecular orbital (FMO) calculations and the method of configuration interaction singles with perturbative doubles, including higher-order corrections and partial renormalization, i.e., PR-CIS(Ds). This method provided an excellent agreement between the theoretical and experimental transition frequencies for the RFP chromophore within proteins. The corresponding differences between the two were 0 cm −1 for DsRed, −560 cm −1 for mCherry, and 320 cm −1 for mStrawberry. Transition frequencies of the isolated model chromophore with slightly different frozen conformations, corresponding to DsRed, mCherry, and mStrawberry geometry in low-temperature crystals, were calculated to be 15,890, 16,370, and 16,530 cm −1 , respectively. For our purposes we take the average value of these three: ] 0 16,260 cm −1 . We assume that it will represent well all possible conformational variations of the chromophore.

Finding Dipole Moment Difference Δμ and the Angle γ between the Vectors μ and Δμ
In this section we will find the absolute values of Δμ and the angle between μ and Δμ using one-and two-photon absorption spectroscopy. This information is necessary for evaluation of local electric fields, performed in Section "Electric Fields in FIGURE 4 | Top: Calculated permanent dipole moment differences, Δμ x (black squares) and Δμ y (red circles), between the excited and ground states of the RFP chromophore as function of applied electric field E x (A) and E y (B). In the calculations, the chromophore structure was optimized for each value of the field with the phenolate and imidazolinone rings constrained to lie in the same plane (flat chromophore). Bottom: Chromophore structure optimized at zero field with the directions of axes x and y shown.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 RFPs". The top panels of Figure 5 show one-photon absorption, fluorescence, and two-photon absorption spectra in the region of the S 0 ->S 1 transition of RFPs. We found that all the proteins studied here can be divided in two groups, according to their Stokes shifts values. DsRed2, mCherry, and mScarlet have smaller Stokes shifts (<1,200 cm −1 ), and mPlum, XRFP, and FIGURE 5 | Top panels: one-photon absorption (dark blue line, right y-axis), fluorescence normalized to one-photon absorption peak (cyan line, right y-axis), and two-photon absorption (purple symbols, left y-axis) spectra in the region of the S 0 -> S 1 transition. Green triangles represent the calculated Δμ 2 values (left y-axis) as a function of excitation wavelength. Bottom panels show the two-photon polarization ratio Ω (red symbols) as a function of excitation wavelength. Horizontal dashed line represents a limiting case of Ω 2/3, c 0. The top x-axis corresponds to transition wavelength, that is the same for both 1P and 2P excitation, and the bottom x-axis corresponds to laser wavelength, used for 2P excitation.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 10 eqFP670 have larger Stokes shifts (≥1,500 cm −1 ). We assume that in proteins with larger Stokes shifts there are additional mechanisms of relaxation in the excited state.
It is known that the absorption spectra of mCherry at pH 7.4 and pH 11.4 (Shu et al., 2006), mPlum (Wang et al., 2004), and DsRed (Matz et al, 1999) contain contributions from at least two different forms. In particular, the mCherry spectrum always contains two contributions, one dominating at acidic pH and another at alkaline pH. In addition to a major red peak, mPlum and DsRed display a shorter wavelength minor peak, corresponding to the immature green chromophore (Gross et al., 2000;Shu et al., 2006). Therefore, in Figure 5 we present absorption spectra in the form of corrected fluorescence excitation spectra scaled to extinction coefficient of the major form. These excitation spectra were measured with the fluorescence monochromator set at a wavelength where the contributions of the minor forms were negligible. The 2PA spectra were all measured in the form of fluorescence excitation spectra and scaled to the two-photon absorption cross section. In these measurements, we again selected the excitation of the main form by choosing the appropriate observation wavelength.
At the long-wavelength side of the pure electronic transition, the shapes of the 1PA and 2PA spectra perfectly match (when plotted against transition wavelength). This demonstrates that the pure electronic transition is simultaneously allowed for both oneand two-photon absorption and its broadening does not depend on the mode of excitation. This behavior is expected for dipolar chromophores, such as the one considered here. Differences in the 2PA and 1PA spectral shapes at shorter wavelengths, i.e., in the region of vibronic transitions, can be explained by the Herzberg-Teller contribution to the Δμ factor, present only in 2PA (Drobizhev et al., 2012b). Green downward triangles in the top panels represent the wavelength dependence of |Δμ| 2 calculated according to Eqs 18, 21 in the region of pure electronic transition from independently measured σ 2 (0-0), ε(0-0), and c.
The bottom panels of Figure 5 show the two-photon polarization ratio Ω as a function of wavelength. This value allows to calculate the angle c. Analysis of the previous literature shows that the direction of transition dipole moment μ generally does not depend on the local environment in proteins with the same chromophore structure (List et al., 2012;Ansbacher et al., 2012;Nifosi et al, 2019, Myskova et al, 2020. On the other hand, the direction of Δμ is more sensitive to local electrostatics (see below). Therefore, we can consider Ω as a qualitative metric for the direction of Δμ within the chromophore coordinates. Furthermore, Ω can resolve spectrally overlapping transitions with different directions of Δμ. These transitions can belong to conformers with different structures of the chromophore environment creating different local electric fields.
For the proteins with a small Stokes shift (four upper graphs in Figure 5), Ω is virtually constant in the region 950-1,250 nm. This points to the presence of a single electronic S 0 -> S 1 electronic transition in this region, cf. (Masters et al., 2018). Ω values between 0.67 and 0.71 correspond to small c angles in these proteins, Table 2. mCherry at pH 11.4 presents an exception, where Ω slightly increases from 0.67 to 0.71 when going from the main peak toward the very red tail of the spectrum at λ > 1200 nm. This can be due to a minor, red-shifted, form contributing to the fluorescence signal at the long excitation wavelengths (despite selective detection of fluorescence). This is probably the form that dominates the spectrum at neutral pH with a different electrostatic environment of the chromophore (Shu et al., 2006). In fact, mCherry at pH 7.4 shows Ω 0.71 across most of the spectrum and its absorption is shifted to the red, peaking near 1,200 nm ( Figure 5).
For the proteins with a large Stokes shift (three bottom graphs in Figure 5), only eqFP670 displays a constant value of Ω across the spectrum. In mPlum and XRFP, Ω first increases and then reaches a plateau in the red tail of the 0-0 transition. Correspondingly, c changes from 19°to 26°in mPlum and from 0°to 22°in XRFP. This suggests a presence of at least two different conformations of the chromophore environment in these proteins.
It is known that in mPlum a large Stokes shift of fluorescence is due to a fast reorganization of the chromophore environment in the excited state (Konold et al., 2014;Faraji and Krylov, 2015;Yoon et al., 2016). Specifically, a direct hydrogen bond between the acylimine oxygen of the chromophore and the E16 carboxyl oxygen of the protein reorganizes to a water-mediated hydrogen bond at the same site. This causes the chromophore to switch from a red emitting form to a far-red emitting form. Since the switch between the two forms occurs on a picosecond time scale, the far-red form dominates in steady-state fluorescence. On the other hand, a small fraction of the far-red form is present in the ground state and can be directly excited (Faraji and Krylov, 2015;Yoon et al., 2016). Our experiment supports this observation because Ω, being independent of emission transition dipole moment by definition (Wan and Johnson, 1994), points to the presence of a minor, far-red form in the ground state with a different direction of Δμ. The latter could be assigned to watermediated hydrogen bonding at the acylimine oxygen. The spectral variation of Ω can be used to resolve the contributions of the two forms (we label the red and far-red forms a and b, respectively), similarly to the case of 1P fluorescence anisotropy measured as a function of excitation wavelength (Lakowicz, 2006). Using the property of additivity of Ω and assuming that both forms fluoresce with the same quantum yield φ b because the red form rapidly converts to the far-red one in the excited state, we can write (Rehms and Callis, 1993): where f a (λ) and f b (λ) are the fractional contributions to the 2PA of the two forms. Specifically, f a (λ) ρ a σ 2,a (λ) ρ a σ 2,a (λ)+ρ b σ 2,b (λ) , and f b (λ) 1 − f a (λ), where ρ a and ρ b are the fractional concentrations of forms a and b. Combining the last relation with Eq. 23, we can write: where Ω a 0.698 and Ω b 0.730 are the corresponding pure Ω values found experimentally. The spectral contributions of the two forms to the total two-photon absorption spectrum, σ 2,total (λ) ρ a σ 2,a (λ) + ρ b σ 2,b (λ), are calculated as follows: ρ a σ 2,a (λ) f a (λ)σ 2,total (λ), and are shown in Figure 6, top, together with the total 2PA spectrum.
Applying the same method to XRFP, we resolved its 2PA spectrum into two forms, Figure 6, bottom. In the case of eqFP670, the independence of Ω on excitation wavelength suggests that there is a single ground state conformation in the range of 1,200-1,300 nm. The large Stokes shift can still be tentatively explained by a conformational change of the chromophore's environment in the excited state. Table 2 summarizes the magnitude of vector Δμ, and the angles it makes with the vector μ (c) and the x-axis (δ), as well as other spectroscopic parameters that we used to obtain this information. For comparison, the table also includes the calculated data for the chromophore in vacuum and relevant literature data. For mPlum and XRFP, where two forms were found in the absorption spectra, we provide a set of parameters for each form. The |Δμ| values were previously measured using Stark spectroscopy for DsRed and mPlum in frozen water/glycerol solution, up to a constant local field factor f (Lounis et al., 2001;Abbyad et al., 2007). Both of them are consistent with our data if we assume f 1.75, which is in the range predicted in (Fried et al., 2013) for the frozen solutions, such as those used in Stark measurements. The angle c 13 0 found for DsRed in (Lounis et al., 2001) using Stark spectroscopy is close to our value (9 0 ). Our experimental value of |Δμ| for DsRed2 also agrees well with the quantum mechanical calculations (List et al., 2012) where different models were explored to describe protein environment. The |Δμ| values obtained for the red anionic chromophore in vacuum (Nifosí et al., 2007;List et al., 2012) with DFT calculations are lower than our semi-empirical calculation probably because the DFT often underestimates |Δμ| (Jacquemin, 2016).

Electric Fields in RFPs
A graphical representation of the curves described by Eqs 16, 21 in the Δμ x , Δμ y plane provides the solutions of this system of equations as the intersections of these two curves. Figure 7 plots the equations of the hyperbola Eq. 15 and the circle Eq. 21 for the DsRed2 protein, using the isolated chromophore parameters calculated in Section "Parameters of the Isolated RFP Chromophore Calculated Quantum Mechanically" and the experimental values of ] and |Δμ| taken from Table 2. There are four possible solutions for Δμ found at the intersections of the red and blue curves in each of the quadrants I-IV. To resolve this ambiguity, we take into consideration the following two points. 1) The Δμ vector in DsRed protein was calculated (List et al., 2012) for a polarizable environment (PE-QM), a frozen polarizable environment (FPE) and a non-polarizable environment (NPE). The purple dashed arrows in Figure 7 correspond to these three cases. The angle between Δμ and the x-axis, calculated  Column 2 shows 1PA 0-0 transition frequency, column 3 -two-photon cross section of the 0-0 transition, column 4two-photon polarizion ratio, (Ω) of the pure electronic transition, column 5 -angle γ between μ and Δμ, column 8 − angle β between μ and x-axis. The relative random errors of measurements of σ 2 (0-0) and Δμ are shown in %. Literature data are presented in italic. *Calculated in this work; (a) Calculated in (Taguchi et al., 2009); (b) Calculated (TD-DFT CAM-B3LYP) for different conformations of red chromophore in vacuum (List et al., 2012); (c) Calculated (DFT BLYP) for red chromophore in vacuum (Nifosí et al., 2007). (d) Calculated in vacuum, (Ansbacher et al., 2012); (e) Stark spectroscopy measurement (Lounis et al., 2001); (f) Calculated (TD-DFT PE-CAM-B3LYP) for DsRed with different physical models of protein environment (List et al., 2012); (g) Stark spectroscopy measurement (Abbyad et al., 2007); (h) Quantum mechanical calculation (Moron et al., 2019); (i) Measured in crystal using absorption of polarized light (Myskova et al., 2020 clockwise, increases in the order: PE-QM, FPE, and NPE, respectively. All three vectors are found in quadrant IV and are close to our experimental result obtained for quadrant IV (purple solid arrow).
2) The direction of the oscillating transition dipole moment μ was also calculated in (List et al., 2012), both for chromophore in vacuum and in protein with different models describing the environment. In all cases, the μ vector fell within the quadrants II and IV and made an angle β −(3-9)°with the x-axis, Figure 7, showing that the protein environment did not affect it much. The direction of μ in an isolated RFP chromophore, optimized in vacuum, was also calculated in (Ansbacher et al., 2012) with a result similar to (List et al., 2012), β 10°, Figure 7.
Knowing the direction of Δμ in a protein and angle c, we can predict the direction of μ from our own experiments. Since Ω measures only an absolute value of c, there will be two possible directions of μ for each pair of quadrants: I-III and II-IV. Suppose that Δμ lies within either quadrants II or IV, then two solutions for μ, μ + and μ − would be possible (Figure 7). The direction of μ + agrees better with the calculations of (List et al., 2012;Ansbacher, et al., 2012). Solutions corresponding to quadrants I and III result in worse agreement with theoretically predicted β. Therefore, all the evidence presented above suggests that the Δμ vector most probably falls within quadrant IV, as shown in Figure 7. Its numerical components (Δμ x and Δμ y ) are presented in Table 3.
Note that the direction of Δμ corresponds to a flow of electronic density globally from P to I and further to A.
The dipole vector diagrams for other proteins (similar to Figure 7) are shown in Supplementary Figure S1. It is interesting to note that the direction of μ in mCherry (β −8.70°) was recently measured independently in crystal (Myskova et al., 2020) and it is very close to what was predicted theoretically for the isolated chromophore (β -(6÷10°) (List et al., 2012;Ansbacher, et al., 2012) and for the DsRed protein (β −(3-6)°) (List et al., 2012). The β angles obtained in liquid solution with our model for μ + (β −19°for DsRed2 and −28°for mCherry pH 7.4) are close, but systematically larger in absolute value than those either calculated or measured in crystals. We think that this discrepancy is because the chromophore can visit a larger set of conformations in solution, allowing the acylimine tail to come into better conjugation with the rest of the molecule, thus turning μ closer to the direction of the tail. Note that the direction of the μ vector is quite conservative for the whole set of proteins studied FIGURE 7 | Graphical solution of Eq. 16 (blue curves) and Eq. 21 (red curve), to find the Δμ x and Δμ y components of Δμ for DsRed2. The structure of the DsRed2 chromophore is displayed in the background, illustrating the molecular system of coordinates. Green dashed bi-directional arrows depict the direction of the oscillating transition dipole moment μ, calculated in (List et al., 2012), for the chromophore in vacuum (light green) and in DsRed protein (dark green). The dotted light green arrow corresponds to a direction of μ calculated in vacuum in (Ansbacher et al. 2012). Purple dashed arrows correspond to the quantum mechanically calculated Δμ vector in DsRed protein environment described by three different models (List et al., 2012); see text. The purple solid arrow is the selected Δμ solution of Eqs 16,21. It was selected based on the known Δμ and μ vectors calculated in (List et al., 2012) and our MD simulations of the electric field; see text. Green solid bi-directional arrows represent the two possible directions of the μ vector based on our measurement of c and the selected direction of Δμ.
FIGURE 6 | Resolution of two-photon transitions corresponding to two different conformations of the chromophore's environment in mPlum (top) and XRFP (bottom). Top panels: Purple circles represent the total effective 2PA cross section σ 2,total (λ) a , blue rhombs-contribution from the red from (ρ a σ 2,a (λ)) and orange squaresfrom the far-red form (ρ b σ 2,b (λ)).  Figure S1). Our next step is to calculate the E x and E y components of the electric field from the experimentally defined Δμ x and Δμ y , using Eqs 12, 13. The results, that we call "quasi-empirical", are presented in Table 3 and Supplementary Table S1. To validate our "quasi-empirical" electric fields and to obtain an additional support for the direction of Δμ vector not only in DsRed but in other RFPs, we performed a series of MD simulations on DsRed, mCherry pH 7.4, mCherry pH 11.4, and mPlum. In this case, electric fields were obtained as negative derivatives of the time-averaged potentials on chromophore atoms along the π-conjugation pathway and projected onto the x-or y-direction, as was done before for GFPs (Drobizhev et al., 2015). We used the chains of atoms: CD2-CG2-CB2-CA2-C2 for E x and O2-C2-CA2 for E y evaluations (see Figure 4 for atoms labeling). The results of these calculations are presented in Supplementary Figure S2 and Supplementary  Table S1. Although in some cases the agreement between the experimental and calculated values of the fields is almost quantitative, i.e., within the standard deviations (mCherry pH 7.4 and pH 11.4, if Δμ is in quadrant IV), in others they agree qualitatively (DsRed and mPlum). Considering the direction of Δμ in DsRed and mCherry (at pH 7.4 and 11.4), the best match between our "quasi-empirical" model and MD simulations is obtained for Δμ falling in quadrant IV. For mPlum, the best match corresponds to quadrant III, with quadrant IV being the next closest. If Δμ of mPlum in fact would lie in quadrant III, the μ vector would adopt a quite unexpected direction, through quadrants I and III. We consider this improbable because of the "conservation" law for the direction of μ, see (List et al., 2012;Ansbacher, et al., 2012), Table 2, and Supplementary Figure S1. Therefore, we assume that Δμ of mPlum is within quadrant IV, like in the other three proteins. The components of the Δμ vector and electric field in all the proteins and in different conformational states, are summarized in Table 3.

Describing the Rate of TICT Process With Marcus Formalism
We now consider another possible mechanism of fluorescence quenching-ultrafast nonradiative transition (jump) from the excited to the ground state, when the RFP chromophore adopts a twisted conformation at a conical intersection seam. Prior to this event, the chromophore in the excited state proceeds through a slower process of twisting of its phenolate group around the bridging CG2-CB2 bond (Olsen and Smith, 2007;Sun et al, 2012;Moron et al., 2019). According to calculations, the twisting occurs in concert with significant CT across the bridge from the P to the I and A groups. When the angle between P and I rings becomes ∼75°-90°, a whole electronic density of the excited state orbital localizes on the I and A parts (Olsen and Smith, 2007). The charge q 0.32e (e is the electron charge) is transferred from P to A and I in the TICT process (Olsen and Smith, 2007).
TICT is common for many molecules featuring electrondonating and accepting groups capable of rotating relative to each other around a bridging bond (Grabowski et al., 2003). If the final CT state is stabilized by a polar solvent, the fluorescence quantum yield drops dramatically (Grabowski et al., 2003), pointing to the dependence of the CT rate on the thermodynamic free energy difference. We hypothesize that in RFPs, the protein electric field can contribute to stabilization (or destabilization) of the TICT state in RFPs. If, for example, the field directed from I to P (E x ) increases (corresponding to concentration of more positive charges on the imidazolinone site and/or more negative charges on the phenyl side in a certain mutant), the CT is expected to accelerate and the conical intersection seam will be reached faster.
In our model, we assume that the potential energy barrier for rotation from a planar to a strongly twisted chromophore is due to steric interactions of the P ring with nearby residues above and below the ring plane. In the RFPs studied here, the most important residues in these positions making direct van der Waals contacts with P are Pro63 (substituted to Thr in eqFP670), Met163 (Lys in DsRed2), and Ile197 (Ala in DsRed2, Arg in mScarlet and eqFP670), as well as Lys70 in DsRed2; see Supplementary Table S3. For our model we assume a priori that the barrier is equal for all proteins, i.e., two parabolic potentials corresponding to the initial and final states in the Marcus model always have same stiffness. The electric field, however, changes between proteins, as we have seen.
To correlate the barrier height with a change of free energy, we use Marcus theory formalism. We assume that the rate of TICT in the excited state is a limiting stage for nonradiative relaxation. Therefore, according to the Marcus theory in the high temperature limit, we can write (Sutin, 1999;Bixon and Jortner, 1999): where ΔG°is the standard Gibbs free energy change in the CT reaction, λ is the reorganization energy, k B is the Boltzmann factor, and T is the absolute temperature. The high temperature limit (h] V < k B T, where ] V is the frequency of vibration coupled with CT) can be justified for different low-frequency modes, such as collective vibrations of protein beta-barrel, ] V (1.5-9) × 10 11 s −1 , (Martin and Matyushov, 2012;Hu et al., 2012;Perticaroli et al., 2014), or twisting modes of the chromophore itself, (3-6) × 10 12 s −1 (Martin et al., 2004). The pre-exponential factor B for a single low-frequency mode ] V can be presented as B κ] V , where κ is the Landau-Zener dimensionless transmission coefficient (Bixon and Jortner, 1999). If κ << 1, the process is considered non-adiabatic, and if κ 1, it is adiabatic. In the non-adiabatic limit, where V is the electronic coupling parameter. Note that Eqs 28, 29 can be derived independently (of Marcus theory) on the basis of the more general Fermi golden rule for nonradiative electronvibrational transitions when only low-frequency modes are involved (Bixon and Jortner, 1999;Callis and Liu, 2004).
Taking the logarithm of both sides in Eq. 28, we obtain Suppose that a negative point charge of an absolute value q is transferred by a distance Δx opposite to the x-axis direction and a distance Δy opposite to the y-axis direction, see Figure 4 for definition of x and y directions. For simplicity, we consider a point charge moving from the geometric center of the initial charge distribution on P (corresponding to a locally excited state) to the geometrical center of the final charge distribution on I and A (corresponding to a TICT state). We assume that the parameters q, Δx, and Δy are the same for all the proteins. If an electric field E (E x , E y ) is present in course of the charge transfer, the electrostatic potential energy of the system will change by We assume that the field does not change during the process of charge transfer and is the same as that found in previous section. The total free energy change, ΔG 0 , consists of a pure molecular (vacuum) part, ΔG 0 vac , and an electrostatic part (imposed by the protein) ΔU: Substituting this into Eq. 30, we obtain which is a second-order polynomial in terms of E x + ηE y . In Eq 33, the constant parameters are defined as follows: Using our experimental results, we plotted lnk nR vs. (E x + ηE y ), varying the parameter η from −0.2 to 0.1 (Supplementary Figure  S3). Out of seven proteins, five (DsRed2, mCherry pH 7.4, mCherry pH 11.4, and red-shifted forms of XRFP and mPlum) fit satisfactorily to a second order polynomial in this range of η. mScarlet and eqFP670 drop significantly out of the correlation and therefore we do not include them in the fit (see discussion below). For each η value, the residual sum of squares of the fit was calculated and the minimum was achieved at η −0.1 (Supplementary Figure S4). The final plot of lnk nR vs (E x + ηE y ) corresponding to η −0.1 is shown in Figure 8.
In their consideration of TICT, Olsen and Smith calculated the parameters G 0 vac and q for the isolated RFP chromophore; see Figure 7 of (Olsen and Smith, 2007). We use these values as well as three experimental coefficients, b 0 , b 1 , and b 2 , corresponding to the best fit of our data to the Marcus model ( Figure 8) to find the final three parameters: λ, Δx, and B from the system of Eqs 34-36. The transversal displacement Δy was then calculated according to Δy η Δx. All the parameters are reported in  ] V mentioned above. This suggests non-adiabatic character of CT (κ << 1) and justifies applicability of Eq 28. It is interesting that the obtained charge displacement Δx 6.7 Å is close to the actual chromophore size and to the displacements predicted theoretically in (Olsen and Smith, 2007;Moron et al., 2019). The positive sign of Δx supports the prediction that the electronic density is moving from P to I, and the negative sign and magnitude of Δy (−0.67 Å) suggests a further slight shift from I to A, in qualitative agreement with predictions (Olsen and Smith, 2007).
To estimate more accurately the theoretically predicted magnitude of charge displacement (Olsen and Smith, 2007), we first notice that a part of the LUMO density, localized on the phenolate group before the transfer, is centered close to the middle point between the CE1 and CE2 atoms (Gross et al., 2000;Mochizuki et al., 2007;Hasegawa et al., 2010;List et al., 2012). Furthermore, the electronic density in a 90 0 P-twisted conformation is delocalized between I and A groups (Olsen and Smith, 2007), with a much higher weight on I than on A (Moron et al., 2019). Therefore, we can assume that the charge density moves from the point between the CE1 and CE2 atoms (whose position does not move during the P rotation) to the C1 atom of imidazolinone. The corresponding crystallographic values (see, e.g., 2H5Q pdb structure for mCherry) in this case are Δx 6.1 Å and Δy −0.9 Å, both in good agreement with our findings (Table 4). These results provide strong independent support both for our approach of finding protein electric fields and for describing the field effect in the framework of TICT. Although five proteins of the series fit quite well to the theoretical model, mScarlet and eqFP670 show significant discrepancies.

Special Case: mScarlet
In the case of mScarlet, the experimental nonradiative rate is much slower than what is expected from the values of the field. This can be explained by the very tight surrounding of the chromophore that strongly elevates the steric energy barrier for P rotation (Bindels et al., 2017). More specifically, instead of the flexible, neutral Ile197 in mFruits, mScarlet has an Arg residue that faces the phenolate ring and makes multiple van der Waals contacts with it, impeding the P-rotation. Moreover, the Arg side group in mScarlet makes a cobweb of numerous hydrogen bonds that holds it more strongly in place compared to Ile in mFruits and Ala in DsRed2. In the original publication (Bindels et al, 2017), these strong local interactions were invoked to explain the exceptionally planar structure of mScarlet chromophore observed in crystal.

Special Case: eqFP670
In the case of eqFP670, the fluorescing state most probably corresponds to a different conformation of the chromophore/ environment in the excited state compared to the ground state, see Section "Finding Dipole Moment Difference Δμ and the Angle c between the Vectors μ and Δμ". Unfortunately, this conformation cannot be interrogated with 2PA (in contrast to mPlum and XRFP), because it is not present in the ground state. Therefore, we can only speculate about the values of the electric field corresponding to this conformer. Suppose that Δμ and the field in the excited state are not very different from the ground state, and that the Δμ vector lies in quadrant IV (Supplementary Figure S1). Then the nonradiative rate turns out to be several times too fast compared to the expected value based on the model fit (left red square in Figure 8).
However, one cannot rule out the possibility that the Δμ of eqFP670 occupies another quadrant. For instance, the direction of Δμ x can flip if E x reaches some critical value. This flipping can theoretically occur at E p x −Δμ 0,x /Δα xx 27 MV/cm; see Eq. 10. The k nR value steadily increases in concert with E x in the series of DsRed2, mCherry pH 11.4, XRFP, mCherry pH 7.4, and mPlum (Table 3). Therefore, one would expect an even larger E x value for eqFP670 than for mPlum (i.e., >19 MV/cm), instead of the predicted E x 14 MV/cm (if Δμ is in quadrant IV). If this were true, it could lead to a negative Δμ x . This situation would correspond to a position of Δμ in quadrants II or III with E x 39.5 MV/cm for both cases. Choosing between quadrants II and III depends on the direction of Δμ y . If Δμ y > 0, i.e., Δμ is flipped relatively to its direction in IV and lies in quadrant II, then the predicted E y value becomes very large (∼250 MV/cm) compared to the rest of the proteins (30-60 MV/cm). This is difficult to explain considering that the electrostatic environment in eqFP670 is not so different from the others. (Estimations show that an additional positive Arg197 and a conserved, presumably negative Glu215, can explain some increase in E y value of eqFP670 versus mFruits, but not larger than up to a factor of 2.) Therefore, we assume that a position of Δμ in quadrant III is possible, and we show a data point for eqFP670 corresponding to this quadrant in Figure 8 (right red square). This point is reasonably close to the model fit and provides E x 39.5 MV/ cm and E y 28.3 MV/cm.
Mechanistically, a large value of E x in eqFP670 (relative to other RFPs) can tentatively be explained by the presence of the Asn143 residue close to the phenolate oxygen of the chromophore. As the crystal structure (4EDS pdb) shows, the carbonyl oxygen atom OD1 of the Asn143 side chain comes in unusually close contact to the phenolate oxygen (∼2.4 Å), which is shorter than the sum of their van der Waals radii. A hydrogen bonding network built around the chromophore phenolate and the Asn143 residue probably helps to hold the two electronegative atoms in such a close contact. Since the Asn carbonyl oxygen carries a partial negative charge, this can result in pushing the negative charge of the phenolate oxygen toward other parts of the chromophore (e.g., I or A). This effect is analogous to an increase of the effective  (Olsen and Smith, 2007). See text for definitions of parameters.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 electric field in the x-direction that potentially flips the sign of Δμ x .

DISCUSSION
It was suggested on theoretical grounds that the nonradiative relaxation of the red FP chromophore in vacuum (Olsen and Smith, 2007) or in proteins (Sun et al., 2012;Moron et al., 2019) generally proceeds via TICT with the phenolate ring turning from 0°to ∼75°-90°around the adjacent C-C bridging bond. This process leads to a conical intersection of the S 1 and S 0 potential energy surfaces that opens a doorway for almost instantaneous S 1 ∼> S 0 relaxation. The phenolate rotation is barrierless for chromophore in vacuum, but encounters a potential barrier in proteins due to steric interactions between the phenolate ring and few nearby amino acid side chains. Since the twisting is coupled to significant charge transfer in the excited state, an external (to the chromophore) electric field could either facilitate or hinder this process, depending on its strength and direction. We apply the Marcus formalism to correlate the barrier height and the free energy difference (ΔG 0 ) between the final and initial states of the charge transfer. In this model, the ΔG 0 contains a contribution equal to the change of electrostatic potential energy between the two states. We observe that for five out of seven proteins, including DsRed2, mCherry (pH 7.4 and pH 11.4), mPlum, and XRFP, the radiationless relaxation rate follows the Marcus theory.
For mScarlet, the measured nonradiative rate is much slower than predicted by the model. We explain this by a much higher potential barrier for rotation of phenolate because of steric clashes with the R198 and Pro64 side chains. The former occupies a unique position on top of the phenolate group and makes several van der Waals contacts with it. The second flanks the phenolate from the opposite side and due to its rigid nature provides an additional barrier for P-rotation. Therefore, we can conclude that a combination of these two amino acids in the chromophore pocket is probably a main factor securing high brightness of mScarlet, compared to, e.g., XRFP (experiencing similar electric field, Figure 8).
In contrast, the nonradiative rate of eqFP670 appears to be larger than predicted by our model, for both possible directions of Δμ, i.e. in quadrants III and IV, Figure 8. Although fieldfacilitated TICT can qualitatively explain the very fast nonradiative rate of eqFP670 (see previous section), other mechanisms cannot be ruled out. One is vibrational relaxation, described by the energy gap law ("Checking the Energy Gap Law" section). For a far-red emitting protein such as eqFP670, this mechanism can strongly contribute to an increase of k nR . Another possibility is that some other molecular parameters involved in TICT are different in eqFP670 compared to other RFPs. These include the charge q, and CT distances Δx or Δy, as well as the potential barrier to rotation. Differences in any of the first three are probable, especially if the acylimine tail adopts a different conformation during the process of excited state relaxation similar to what was suggested in mPlum (Moron et al., 2019). There are two different factors that could affect the height of potential barrier for P-rotation. First, similar to mScarlet, the phenolate of eqFP670 is flanked by Arg197 with four van der Waals contacts between them (although in mScarlet there are six). Also, Arg197 participates in a number of hydrogen bonds that hold its position fixed. On the other hand, the more flexible Thr60 found on the opposite side of the phenolate group in eqFP670, replacing the rigid cycle of the Pro63 residue, probably makes phenolate rotation easier.
If we only consider the RFPs with a "moderately-soft" internal pocket around phenolate group, i.e., similar to those containing FIGURE 9 | 3D plot representing quantum yield (A) and shift of peak absorption transition frequency relative to the vacuum frequency (B) as a function of E x and E y in a set of RFPs with a "moderately-soft" phenolate pocket (see text). Five representative RFPs are shown by red dots. Large quantum yields in (A) correspond to red and small quantum yields-to purple color of the surface. The calculation is based on model Eqs 2, 33-36 and model parameters presented in Table 4. Large negative shifts in (B) corresponds to purple and large positive shifts-to red surface color. The prediction is based on Eq. 15 and parameters from "Parameters of the Isolated RFP Chromophore Calculated Quantum Mechanically" section.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 Pro63 and Ile197 (or Ala197), we can establish some qualitative structure-property relationships for photophysical parameters that can potentially help to guide future RFP engineering. Figure 9A shows a 3D plot of the quantum yield as a function of the electric field components E x and E y . The quantum yield is calculated using Eq. 2 with k R 0.158 × 10 9 s −1 (mean of all seven k R values, Table 1), and k nR described by Eqs 33-36. As one can see, for each fixed E y value, the quantum yield first decreases as E x increases, then reaches a minimum, and then starts to increase. For all reasonable E y values, in the region of negative E x values, the potential barrier for TICT via P-rotation becomes very high due to the combined action of the "charge locking" effect of E x and steric clashes of phenolate with nearby residues. This results in φ 1. With increasing positive E x , the barrier becomes lower, thus accelerating nonradiative relaxation. The k nR reaches its maximum (minimum φ) when ΔG 0 −λ, in accord with the Marcus model. Even larger E x values correspond to a so-called inverted Marcus region (Bixon and Jortner, 1999), where the k nR rate starts to decrease again and quantum yield to increase. For a fixed E x within the normal Marcus region, φ increases with the increase of E y , but the absolute values of gradients are much smaller than for E x . This behavior follows from the fact that the charge is transferred counter to the E y field, but at a much shorter distance compared to the x-direction. For an illustration, we show the data points corresponding to DsRed2, mCherry pH 11.4, mCherry pH 7.4, mPlum, and XRFP. The transition from DsRed2 to mCherry pH 11.4 corresponds to the drop of quantum yield by 30%, following an increase of E x by 3 MV/cm (with E y unchanged), whereas transition from mCherry pH 7.4 to mPlum also corresponds to a drop of φ by 30%, but with the decrease of E y by 12 MV/cm (with almost constant E x value).
The absorption frequency shift (relatively to vacuum ] 0 16,260 cm −1 , λ abs,0 615 nm) as a function of E x and E y is presented in Figure 9B. It is calculated according to Eq. 14. Here, the further red-shifted variants correspond to small positive (or, better, even negative) E y values simultaneously with E x values around 20-30 MV/cm. To have both a red-shifted absorption and a reasonably high quantum yield, one has to keep both E x and E y in the range of small positive or negative values, i.e., from −10 to +10 MV/cm. A favorably low E x of the DsRed2 protein leading to a high quantum yield is most likely due to the positively charged Lys163 residue hydrogen-bonded to the phenolate oxygen. In mFruits and mScarlet, E x is larger, probably because of lacking of the charge in that position.
Possible mutagenesis strategy for optimizing both parameters (φ and λ abs ) using the DsRed scaffold could contain the following steps. First, one should keep a positively charged residue H-bonded to the phenolate oxygen (to have a small E x ). Second, one should concentrate more positive charges close to the acylimine oxygen (to decrease E y ) and/or more negative charge on the imidazolinone carbonyl oxygen (to further decrease E x and decrease E y ). The last change may be a difficult task however, because the Arg95 residue hydrogen bonding to the imidazolinone oxygen seems to be conservative in RFPs. It is interesting that similar recipes for improving quantum yield and shifting fluorescence to the red were put forward on the basis of quantum calculations of several different mPlum conformers (Moron et al., 2019).
We finally would like to note that in the case of two-photon excitation, peak absorption wavelengths of RFPs automatically fall in the tissue transparency window and therefore the 2P brightness, but not the position of 2PA peak becomes a most critical parameter. The relative 2P brightness can be estimated to a first approximation as a product Δμ 2 φ. Both Δμ 2 and φ depend on the electric field and our model, see Figure 10, predicts that for getting 2P brighter variants one has to decrease or even flip the sign of E x component, with E y playing not too critical role. The structural guidelines for increasing 2P brightness are therefore similar to what was described above for obtaining higher quantum yield.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
MD and TH conceived the idea of the paper and supervised the project, MD and RM wrote the paper, MD and RM made experimental measurements, RM expressed and purified all the FIGURE 10 | 3D plot representing the 2P brightness (roughly proportional to Δμ 2 φ) in the region of the 0-0 excitation frequency as a function of E x and E y . Large brightness corresponds to red color on the surface and small brightness-to purple. Four representative RFPs are shown by red dots.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 633217 proteins, MD developed the physical models, PC and JS accomplished quantum mechanical and molecular dynamics simulations, GL, AS, and NS engineered XP fluorescent protein.