Local redox conditions in cells imaged via non-fluorescent transient states of NAD(P)H

The autofluorescent coenzyme nicotinamide adenine dinucleotide (NADH) and its phosphorylated form (NADPH) are major determinants of cellular redox balance. Both their fluorescence intensities and lifetimes are extensively used as label-free readouts in cellular metabolic imaging studies. Here, we introduce fluorescence blinking of NAD(P)H, as an additional, orthogonal readout in such studies. Blinking of fluorophores and their underlying dark state transitions are specifically sensitive to redox conditions and oxygenation, parameters of particular relevance in cellular metabolic studies. We show that such dark state transitions in NAD(P)H can be quantified via the average fluorescence intensity recorded upon modulated one-photon excitation, so-called transient state (TRAST) monitoring. Thereby, transitions in NAD(P)H, previously only accessible from elaborate spectroscopic cuvette measurements, can be imaged at subcellular resolution in live cells. We then demonstrate that these transitions can be imaged with a standard laser-scanning confocal microscope and two-photon excitation, in parallel with regular fluorescence lifetime imaging (FLIM). TRAST imaging of NAD(P)H was found to provide additional, orthogonal information to FLIM and allows altered oxidative environments in cells treated with a mitochondrial un-coupler or cyanide to be clearly distinguished. We propose TRAST imaging as a straightforward and widely applicable modality, extending the range of information obtainable from cellular metabolic imaging of NAD(P)H fluorescence.


ELECTRONIC STATE MODEL
The local electronic state probability distribution of fluorophores in a sample can be described by the state vector ̅ ( ̅ , ). In the case of the NADH model in Figure The one-photon excitation (OPE) rate is calculated as 01 ( ̅ ) = • Φ ( ̅ ), where is the fluorophore OPE cross section and Φ ( ̅ ) is the local excitation photon flux (related to by Φ = (ℎ ) ⁄ , where ℎ ⁄ is the excitation photon energy). A measurement of the free NADH fluorescence lifetime, = 0.4 ns, is also used to fix the value for the sum of the first excited singlet state decay rates, 10 + + = 1⁄ .
Note that for rectangular excitation pulses of constant laser intensity, the model matrix ( ̅ ) is time independent during excitation. This allows for efficient calculation of ̅ ( ̅ , ) using the matrix exponent , which can be implemented numerically or worked out analytically by diagonalizing ( ̅ ).
For most fluorophores, including NADH, only the first excited singlet state contributes significantly to the overall luminescence. In the case of a homogeneous sample, the detected luminescence is therefore , where ( ̅ ) is the collection efficiency function of the detection system, the overall detection quantum yield, is the fluorophore concentration and the fluorescence quantum yield. For a sample with multiple emissive states, the above equation is easily modified to include additional components of the state vector. Two-photon experiments used a Ti:Sapphire excitation laser with a repetition rate of = 76 MHz and a pulse width measured to 150 fs FWHM. Since the pulse period ( = 1⁄ ≈ 13 ns) is orders of magnitude shorter than the shortest TRAST pulse duration to be measured (2 μs), the excitation source can be treated as a continuous laser with an average excitation rate of , where Φ ( ̅ , , ′) describes the local excitation photon flux and is the TPE cross section. The integration variable ′ refers to micro-time in relation the laser sync signal and is used to describe the excitation dependence within a single laser pulse period. Slower excitation dependencies introduced by TRAST modulation (such as laser scanning, see section 3 in SI) are described by the macro-time, .

SCANNED TRAST GAUSSIAN PULSE CORRECTION (MACRO-TIME AVERAGING)
Performing scanned TRAST with a Gaussian excitation laser adds a time dependence to the excitation rate, 01 ( ̅ , ), with the irradiance varying continuously as the beam moves past a fixed point in the sample. These variations occur on the same time scale as the intended TRAST pulse duration and must therefore be taken into account through a time-dependent model matrix, ( ̅ , ). As a result, determining ( ) as outlined in eqs. (S3-S5) becomes considerably more computationally expensive. In this work, we treated scanned TRAST as if the excitation pulses were rectangular in time, thereby allowing the same treatment as for the case of a stationary modulated laser (section 1 in SI).
We thus determined effective rectangular pulse amplitudes and durations that best match the experimental Gaussian pulses ( Figure S1C). Simulation of the complete, time-dependent scenario showed that errors introduced by assuming rectangular pulses are minimized for effective pulse amplitudes of 0.71 ± 0.012 (95 % conf.) times the Gaussian peak amplitude and pulse durations of 0.90 ± 0.038 (95 % conf.) times the Gaussian −2 pulse width ( Figure S1B). The pulse optimization was performed for two very different photophysical models, representing NADH and Rhodamine 110 1 respectively. Both models returned the same pulse correction factors, irrespective of what photophysical rate parameters were used within those models ( Figure S1A, only data for NADH shown). It can be noted, within the given error margins, that the effective pulse amplitude found through simulations corresponds to the fluorescence intensity weighted average amplitude of a Gaussian pulse far from saturation ( 01 ( ̅ , ) ≪ 10 ).
, where 01 ( ̅ , ) = 01, is the Gaussian excitation rate. Here ( ̅ ) = 0 ( ̅ )⁄ represents the duration of the Gaussian excitation pulse and is directly related to the spatial −2 beam width in the scanning direction, 0 ( ̅ ), and to the scanner speed, . For such a rectangular pulse to preserve the total laser irradiance applied, i.e. for the area under the excitation pulse to stay the same, the effective pulse duration is , again matching the simulated result.

DIFFUSION MODELLING
Equation (S3) models electronic state transitions but does not take the effects of diffusion into account. A more complete description can be based on the reaction-diffusion equation , where D is the diffusion coefficient which in this case was assumed to be constant for all species. Diffusion can speed up dark state recovery by exchanging dark molecules for fresh ones from outside the excitation volume, thereby reducing the observed TRAST amplitude. Since dark states with relaxation times much faster than the diffusion time through the detection volume are not affected in this way 2 , it is often sufficient to use a simplified diffusion model based on eq. (S3) with effective diffusion recovery rates inserted in the model matrix ( ̅ ) instead.
In this work we used the model matrix in eq. (S10), where the faster states are not directly affected by diffusion, but the long-lived NAD + is allowed a set of two effective diffusion recovery pathways back to 0 NADH. The reason to include two recovery relaxations is that the extended confocal detection volumes are not spherical, but elongated along the optical axis, and thus a slower recovery along this axis might need to be considered. These recovery rates reflect that the depletion zones of 0 NADH go beyond the excitation volumes and depend on the different excitation geometries under OPE and TPE respectively. The diffusion recovery rates were fitted globally for the OPE and TPE measurement data and the values are presented below in Supplementary Table 1.
The second diffusion recovery state, based on 2 and 3 , was found to be necessary to fit the OPE data. The corresponding TPE measurements did not require this additional dark state, likely reflecting the more confined excitation volume of TPE. The inclusion of these rates ( 2 and 3 ) also in the TPE global fit resulted in a second diffusion recovery state of near zero lifetime (1 2 ⁄ ≈ 0), while offering no improvement in fit quality.

SPATIALLY AVERAGED EXCITATION RATE
Typically, a gradient based diffusion model, see eq. (S9), is not required 2, 3 . The calculation of ( ) in eqs.
(S1-S5) can then be simplified considerably by removing the spatial dependence altogether. The average sample brightness becomes , where ADH( ) 1 is calculated using a space-invariant model matrix, . The excitation rate 01 ( ̅ ) is replaced by the average observed excitation rate. For stationary TRAST we use a weighted volume average 4,5 , where ̅̅̅̅̅̅̅̅̅ 1 ( ̅ ) = 01 ( ̅ )/( 01 ( ̅ ) + 10 ) is the 1 population at the onset of excitation, after equilibrium of the singlet states, but before the build-up of dark states. For scanned TRAST the spatial averaging is performed only in the two dimensions perpendicular to the scan direction , since the scan direction now effectively acts as the temporal domain. The averaging along the scan direction was discussed further in SI section 3.

DUTY CYCLE CONSIDERATIONS
In the TRAST analysis we assume that all pulses in a pulse train produce an identical fluorescence signal, ( ). In other words, the duty cycle must be low enough to allow a complete return to the same initial distribution, ̅ ( ̅ , 0), between each pulse. In previous works based on TRAST, a fixed duty cycle of = 0.01 was proven low enough for the organic fluorophores used 2, 4, 5 .
The influence of duty cycle in this study was tested by repeatedly measuring NADH for a few fixed pulse durations while varying the duty cycle. For stationary TRAST ( Figure S2A), as the duty cycle approaches = 1 (continuous excitation without recovery) all pulse durations obviously produce the same fluorescence signal, while for low enough duty cycles there is no dependence on duty cycle at all. When measuring the OPE TRAST curves in this work (Figure 1), as a compromise between low duty cycle (less distortion) and long illumination time (less noise), we performed our measurements at = 0.01.
For scanned TRAST measurements there is no connection between measurement time and duty cycle. Instead, we are free to select arbitrarily low duty cycles by increasing the length of the path to be scanned. In this case we used a duty cycle of = 0.0046, below the region where we expect distortions in the TRAST curves. As shown in Figure S2B, there is no change in the data if the duty cycle is lowered further.

MODEL EVALUATION FOR OPE (AIC/BIC)
A set of similar photophysical models were compared by use of both the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) 6 . For least squares fitting of multiple models to the same data, AIC and BIC scores can be calculated by , where is the residual sum of squares from the fit, the number of data points (constant) and is the number of fitted parameters in the model. The model with the lowest AIC/BIC score is preferred. The two criteria differ only in how much they penalise models with a larger number of fit parameters, . In our case, with a relatively large , the BIC is much more restrictive towards additional model parameters.
If the lowest AIC score is and model deviates from this score by ∆ ( ) = − , the relative likelihoods of each tested model can be expressed by the Akaike weights , where is the total number of models tested. Weights for the BIC scores can be calculated in the same way. Figure S4 shows the calculated weights for both AIC and BIC for 9 different variations of the model in Figure 2, with this base model receiving by far the highest likelihood according to both criteria. The BIC strongly prefers models with fewer parameters and assigns also a small probability to a model without any triplet state at all (M1). The AIC is more permissive towards additional rates and suggests small probabilities also for electron ejection from the triplet state (M2) and photon driven electron ejection from 1 (M4).