Regression Analysis of Confocal FRAP and its Application to Diffusion in Membranes

In most biological processes, diffusion plays a critical role in transferring various bio-molecules to transfer desirable locations in an effective and energy-efficient manner. How fast molecules are transferred is measured by diffusion coefficients. Since each bio-molecules, in particular, signaling molecules have their unique diffusion coefficients and quantifying the diffusion coefficients help us to understand various time scales of both physiological and pathological processes in biological systems. Moreover, since diffusion profiles of a diffusant vary in different micro-environments of cell membranes, accurate diffusion coefficient also can provide a good picture of membrane landscapes as well as interactions of different membrane constituents. Currently, only a few experimental methods are available to assess the diffusion coefficient of a biomolecule of interest in live cells including Fluorescence Recovery After Photobleaching (FRAP). FRAP was developed to study diffusion processes of biomolecules in the cell membranes in the 1970s. Albeit its long history, the main principle of FRAP analysis has remained unchanged since its inception: fitting FRAP data to a theoretical diffusion model for the best fitting diffusion coefficient or using the relation between the half time of recovery and ROI size. In this study, we developed a flexible yet versatile confocal FRAP data analysis framework based on linear regression analysis which allows FRAP users to determine the diffusion from either single or multiple FRAP data points without data fitting. We also validated this approach for a series of fluorescently labeled soluble and membrane-bound proteins and lipids.


Introduction
Diffusion is a major transporting tool living organisms utilize without the consumption of energy. By diffusion process as part of larger cell signaling process, cells regulate the traffic of proteins, various critical gases, biochemical substances, and ions not only across the membranes but also between intra-and extra-cellular spaces through osmosis, ion channel, and other mechanisms [1,2]. On the plasma membrane of cells, many membrane proteins are in constant translational motion which allows them to move to various locations in less time and with significantly less effort [3]. The diffusion rate of spherical particles through a liquid with a low Reynolds number can be affected by protein's intrinsic properties such as protein sizes as well as surrounding environments [4,5]. How rapidly diffusion occurs is characterized by the diffusion coefficient D, a parameter that is dependent on temperature (T), dynamic viscosity ( ) and the radius of a diffusant ( ) with the Boltzmann's constant, B as a proportionality constant: which is referred to as the Stokes-Einstein equation [6,7]. D can be understood as a measure of the mean squared displacement per unit time of the diffusant. From a diffusion medium point of view, cells also have mechanisms to immobilize specific membrane proteins, as well as ways of confining both membrane protein and lipid molecules to particular domains in a continuous lipid bilayer which may affect protein's diffusion significantly [4,5]. Therefore, to calibrate a protein's diffusion rate is critically important to understand protein's biochemical properties, interactions with other cellular components, and cell membrane structures.
To measure the diffusion coefficient of a protein in a live cell is a challenging task, and fluorescence-based techniques such as single particle tracking (SPT), fluorescence correlation microscopy (FCS), photoactivation, and fluorescence recovery after photobleaching (FRAP) are currently used in the lab. Among these, FRAP is one of the most accessible methods to calibrate biological diffusion, since FRAP can be done by most modern laser scanning confocal microscopes (LSCMs) without any additional preparations.
FRAP takes advantage of the photo-illuminating and photobleaching properties of fluorescence molecules. Fluorescence molecules can be illuminated by proper strength of laster intensity within absorption spectrum wavelengths. In general, the fluorescence intensity from the fluorescence molecules is proportional to illumination laser intensity up to a certain threshold level of laser intensity. Beyond this threshold level of laser intensity, fluorescence molecules irreversibly lose fluorescent property, which is called photobleaching. For FRAP experiments, subject bio-molecules labeled by fluorescent tags are introduced to a cell and a user-defined region of interest (ROI) within the cell is photobleached using a strong pulse of laser for a short period of time. By reversing the laser intensity to illumination level, recovery of fluorescence in the ROI can be observed due to the mixing of photobleached molecules in the ROI and fluorescently intact molecules from outside of the ROI by random movements, or diffusion. The time course of fluorescence recovery in the ROI is called a FRAP curve, or FRAP data. By fitting this FRAP data by a series of theoretical FRAP curves generated using different diffusion coefficients, the best fitting diffusion coefficient is determined by minimizing the sum of squared error (SSE) between the FRAP data and theoretical FRAP curve [8][9][10].
For some FRAP users, the data fitting step can be technically difficult. As an alternative, for a circular ROI, many FRAP users used the ratio between the half time of recovery ( 1∕2 ) and the radius of ROI, or nominal radius ( r n ) to estimate the diffusion coefficient where 1∕2 is the time required for a bleach spot to recover halfway between the initial fluorescence intensity and steady-state fluorescence intensity [11]; Recently, we also improve this relation taking the diffusion during the photobleaching into consideration in terms of effective radius ( r e ) [12][13][14][15]; where r e is the effective radius measured from the postbleach profile ( Fig. 1A) modeled by Here, K is the bleaching depth parameter and r is the radial distance from the center of an ROI. Therefore, for a normalized postbleach profile, the maximum fluorescence intensity, F max = 1 , the minimum fluorescence intensity, F min = f (0) = 1 − K , and the effective radius, r e is the halfwidth of a post-bleach profile at 1 − Ke −2 height (Fig. 1A).
Different from data fitting, since this r e ∕ 1∕2 approach (Eq. 2) uses only one data point of a FRAP curve, the calculation of D is simple enough to be done by hand. On the contrary, since Eq. (2) uses only one data point, it may yield erroneous results when FRAP data are noisy . F i , F 0 and F ∞ are the prebleach steadystate fluorescence intensity, postbleach initial fluorescence intensity, and postbleach steady-state fluorescence intensity, respectively. t 2∕3 is the time when the recovery curve reaches 2/3 of the total recovery from F 0 to F ∞ . (C) Transformed normalized FRAP curve, g(t) (Eq. 7) for linear regression analysis. As the FRAP curve is getting closer to the steady-state, g(t) shows large fluctuations especially near t = 1∕2 . Still, even for this disadvantage, analytic r e ∕ 1∕2 expressions in Eq. (2) allow us to predict what will happen to D for different halftime of recovery ( 1∕2 ) and/or bleaching spot radius ( r n ), and/or diffusion during the photobleaching ( r e ) in intuitive ways. In this study, by applying linear regression analysis on confocal FRAP data, we will derive a generalized analytic expression similar to Eq. (2), which allows FRAP users to calculate diffusion coefficients from any number of FRAP data points without data fitting.

Theory
In our previous study [12], we derived a simplified FRAP equation for lateral diffusion confocal FRAP data from a circular ROI as: where D = r 2 e ∕(4D) , = r n ∕r e , and K is a bleaching depth parameter. The mobile fraction, M f is defined as where F i , F 0 and F ∞ are the prebleach steady-state fluorescence intensity, postbleach initial fluorescence intensity, and postbleach steady-state fluorescence intensity, respectively (Fig. 1B). For normalized FRAP data, F i = 1.
To apply linear regression analysis, we first transform Eq. (4) to a linear model. First, notice that since F(0) = F 0 , by letting t = 0 in Eq. (4),

By solving for K,
On the other hand, by solving Eq. (4) for K∕ ( If we use Eq. (6) for K, where we took the reciprocal of both sides in the last step. Therefore, for a given FRAP data F Data (t) , if we transform the FRAP data by then g(t) satisfies a linear relation in t with a slope m = 2∕[ D (1 + 2 )] and with 0 as the y-intercept; where the diffusion time, D = r 2 e ∕(4D). As we can expect from the definition of g(t) in Eq. (7) when F Data (t) is getting closer to F ∞ the denominator becomes smaller and, therefore g(t) will show larger fluctuations (Fig. 1B, C). Largely, a FRAP curve can be divided into two parts: The recovery phase and the steady steady phase. In the recovery phase, a FRAP curve shows monotone increasing from F 0 to F ∞ . In the steady steady phase, the FRAP curve fluctuates up and down near the postbleach steady-state fluorescence intensity, F ∞ . Mostly, FRAP data points in the steady steady phase form outliers in the transform data, g(t) (Fig. 1C). To avoid this, we choose the data points from the recovery phase, 0 ≤ t ≤ t 2∕3 , where t 2∕3 is the time the recovery curve reaches 2/3 of the total recovery ( Fig. 1B, C) to stay in the recovery phase of FRAP data. Now, suppose that g(t k ) is transformed FRAP data at t = t k and let y = mt k be the best fitting linear regression line for g(t k ) for k = 1, 2, 3, ⋯ , n . Then, the Sum of Squared Error (SSE) is defined as which is quadratic in m (i.e. a parabola). Therefore, SSE will have the minimum when the derivative, d dm SSE = 0 (the vertex of a parabola).
By solving for m, Once the minimizing slope m for SSE is determined, by using Eq. (8), we can calculate the best fitting D from, By using D = r 2 e ∕(4D) and 2 = r 2 n ∕r 2 e , we finally obtain This equation can be used to find the relationship between D and 1∕2 by using Therefore, by using only one data point, 1∕2 , F( 1∕2 ) in Eq. (10), which is identical to Eq. (2). In this sense, Eq. (10) can be regarded as a generalization of Eq. (2). In another study [13], we also showed that r e can be calculated from r 1∕2 by where r 1∕2 is the half-width at the half depth in the postbleach profile (Fig. 1A).
Combining Eqs. (7) and (8), Furthermore, by replacing g(t k ) with the original FRAP data ( F Data (t k ) ) using Eq. (7), Equations (9) and (10) indicate that a diffusion coefficient from FRAP data can be thought of as a weighted mean of transformed FRAP data over different time points. Therefore, by choosing a different set of FRAP data points, we can find a diffusion coefficient specific to that portion of the FRAP data. In our analysis, since we chose the FRAP data points from the recovery phase of a FRAP curve up to t 2∕3 , the diffusion coefficients found here will be the weighted mean of diffusion coefficients from FRAP data points up to t 2∕3 (Fig. 1).

Results
To validate our new FRAP analysis framework, we reanalyzed the FRAP data of several membrane proteins and lipid probes in COS-7 cell membranes from our previous publications [13,16]. The data analyzed here include FRAP data of Alexa488-conjugated Cholera toxin B subunit (Alexa488-CTxB), YFP-GL-GPI, YFP-GT46, and DiI16, which are standard model proteins for lateral diffusion in cell membranes (Fig. 2). CTxB is a homopentameric protein, containing binding sites for 5 GM1 gangliosides and is commonly used as a potential marker of lipid rafts. CTxB is localized to the outer leaflet of the plasma membrane and is well known for exceptionally slow diffusion for its molecular weight when attached to the gangliosides using its five B subunit. This is presumably due to either direct or indirect interactions with the actin cytoskeleton [16,17]. YFP-GL-GPI consists of yellow fluorescent protein (YFP) attached to a single GPI lipid anchor. Due to the lack of a transmembrane domain nor cytoplasmic tail, YFP-GL-GPI is located exclusively outer leaflet of the plasma membrane similar to CTxB. Even though YFP-GL-GPI is associated with cortical actin and possibly lipid-raft as well, YFP-GL-GPI has much faster diffusivity than CTxB even though YFP-GL-GPI has a larger molecular weight (28.4 kDa) than Alexa488-CTxB (12.2 kDa) [16][17][18]. While Alexa488-CTxB and YFP-GL-GPI are peripheral membrane proteins, YFP-GT46 is an artificial transmembrane protein consisting of an extracellular YFP and a short cytoplasmic tail. YFP-GT46 is not raft-associated [16,17,19].
Lastly, DiIC16 is a lipophilic fluorescent stain for labeling membranes. Once incorporated into membranes, DiIC16 diffuse laterally within the cellular plasma membranes with lipid molecules. Therefore, diffusion of DiIC16 is not restricted by the cytoskeleton and shows extremely high diffusivity [16].
To calculate D from FRAP data, we first measure the effective radius ( r e ) using Eq. (11) from the half depth of the postbleach profile, where F max and F min are maximum and minimum fluorescence intensities in the postbleach profiles (Fig. 1C), and then r 1∕2 was measured from the half-width at f 1∕2 in the postbleach profile. The effective radius, r e calculated from r 1∕2 were 3.0 ± 0.4 m for CTxB, 3.1 ± 0.3 m for GLGPI, 3.0 ± 0.3 m for GT46, and 3.7 ± 0.6 m for DiIC16 which show no significant differences from r e determined from data fitting (two-sample student t-test, p < 0.05 , n r 1∕2 = 12 , n fit = 12 ). Figure 3 show the postbleach profiles with r e from f 1∕2 versus data fitting in comparison with the postbleach data for Alexa488-CTxB, YFP-GL-GPI, YFP-GT46, and DiIC16.   We next transformed the data for linear regression using the relation Eq. (7) and choose data points up to t 2∕3 for regression analysis (Fig. 4). The transformed FRAP data for Alexa488-CTxB, YFP-GL-GPI, YFP-GT46, and DiIC16 showed a linear relation as expected and were well described by linear regression lines up to t 2∕3 (Fig. 4A1-D1). To confirm that the diffusion coefficients found from the linear regression analysis on truncated FRAP data up to t 2∕3 describe the whole FRAP data well, heoretical FRAP curves were plotted using the diffusion coefficients found from the linear regression analysis and compared with entire FRAP data for Alexa488-CTxB, YFP-GL-GPI, YFP-GT46, and DiIC16 (Fig. 4A2-D2). The results indicated that diffusion coefficients found from linear regression analysis on partial FRAP data up to t 2∕3 successfully describe the whole FRAP data. The diffusion coefficients found from the linear regression analysis were 0.19 ± 0.11 for Alexa488-CTxB ( n = 12 ), 1.21 ± 0.40 for YFP-GL-GPI ( n = 12 ), 0.49 ± 0.23 for YFP-GT46 ( n = 12 ), and 2.39 ± 0.74 for DiIC16 ( n = 12 ), which were in excellent agreement with what we reported in our previous study using data fitting methods (two-sample t-test, p < 0.05 ) (Fig. 4).

Discussion
Since the calculated D was obtained from multiple data points, the linear regression based FRAP analysis (Eq. 12) requires more steps compared to r e ∕ 1∕2 approach using Eq. (2). However, it should provide more reliable results. Moreover, this approach provides a comprehensive way to calculate the diffusion coefficient not only from multiple FRAP data points but also from a single data point. If we choose the data point (t i , F Data (t i )) form the FRAP data, then the corresponding diffusion coefficient will be which coincides with the result from our previous study even for the different origin of the equation [13]. In particular, at the half time of recovery, we have and therefore we obtain , as in Eq. (2). Similarly at t = t 2∕3 , and therefore we obtain which provides another layer of safety net to double-check the result from r e ∕ 1∕2 approach using Eq. (2). Importantly, if we use two data points ( 1∕2 , F Data ( 1∕2 )) and (t 2∕3 , F Data (t 2∕3 )) then the averaged diffusion coefficient from two data points should not be but the weighted average; from Eq. (12).
In summary, we developed a novel linear regression-based non-fitting FRAP analysis framework to calculate the diffusion coefficient from user-selected single or multiple confocal FRAP data points. This approach has a few major advantages over the other FRAP analysis approaches. First of all, since Eq. (12) works for any number of FRAP data points, there is no need to measure 1∕2 from the FRAP curve as in Eq. (2). Secondly, our new approach can be used to analyze separate parts of FRAP data in the different time scales of interest by choosing FRAP data points differently. In this way, timedependent diffusion coefficients can also be obtained from the FRAP data [13]. Thirdly, for a small number of FRAP data points, the explicit formula Eq. 12 becomes simple enough to handle by hand without any computer software. Lastly, the calculated diffusion coefficient is the true minimizer of the sum of squared errors (SSE) from the FRAP data and does not experience any pitfalls from data fitting algorithms. We believe that the flexibility and versatility of our new approach can be beneficial to many FRAP users.

r e -measurements
The postbleach images were normalized by dividing by the prebleach images, and then the postbleach profiles along the center of the normalized images were averaged over multiple data sets (n = 5-6). The mean postbleach profiles were then fitted to Eq. (3) for the bleaching depth parameter (K) and an effective radius ( r e ) as fitting parameters using a nonlinear least-squares fitting routine (nlinfit.m) available in MATLAB (version 9.2.0, R2017a; The Mathworks, Inc, Natick, MA). r e was also measured from r 1∕2 using Eq. (11), where r 1∕2 is the half-width at the half depth of a postbleach profile ( f 1∕2 ) where lower and upper 5 percent of data points in a postbleach profile were used to determine F max and F min .
Funding No funding was received to assist with the preparation of this manuscript.

Data Availability
The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Declarations
Ethics Approval Not applicable.

Conflicts of Interest
The author declares that there is no conflict of interest.