Table of Contents Author Guidelines Submit a Manuscript
Abstract and Applied Analysis
Volume 2015 (2015), Article ID 859849, 17 pages
Research Article

On Data Space Selection and Data Processing for Parameter Identification in a Reaction-Diffusion Model Based on FRAP Experiments

1Industrial Mathematics Institute, Johannes Kepler University of Linz, Altenbergerstr 69, 4040 Linz, Austria
2Institute of Complex Systems, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Faculty of Fisheries and Protection of Waters, University of South Bohemia-České Budějovice, Zámek 136, 373 33 Nové Hrady, Czech Republic

Received 2 January 2015; Revised 15 May 2015; Accepted 18 May 2015

Academic Editor: Benito M. Chen-Charpentier

Copyright © 2015 Stefan Kindermann and Štěpán Papáček. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Fluorescence recovery after photobleaching (FRAP) is a widely used measurement technique to determine the mobility of fluorescent molecules within living cells. While the experimental setup and protocol for FRAP experiments are usually fixed, data (pre)processing represents an important issue. The aim of this paper is twofold. First, we formulate and solve the problem of relevant FRAP data selection. The theoretical findings are illustrated by the comparison of the results of parameter identification when the full data set was used and the case when the irrelevant data set (data with negligible impact on the confidence interval of the estimated parameters) was removed from the data space. Second, we analyze and compare two approaches of FRAP data processing. Our proposition, surprisingly for the FRAP community, claims that the data set represented by the FRAP recovery curves in form of a time series (integrated data approach commonly used by the FRAP community) leads to a larger confidence interval compared to the full (spatiotemporal) data approach.

1. Introduction

The method of fluorescence recovery after photobleaching (FRAP) is based on the measurement of the change of fluorescence emitted by autofluorescent molecules or fluorescently tagged compounds (e.g., green fluorescence proteins (GFP) in a region of interest (ROI), usually a 2D Euclidean bounded domain, in response to a high-intensity laser pulse provided by confocal laser scanning microscopy (CLSM). The initial steady state is thus perturbed by an external stimulus, the so-called bleach. The bleach (or bleaching) causes an irreversible loss of fluorescence in the bleached area, apparently without any significant damage to intracellular structures. After bleaching, the observed recovery in fluorescence reflects the mobility of fluorescence compounds from the area outside the bleach [1].

To quantify the mobility of photosynthetic proteins of different microbial species is the main research interest of biologists. In Figure 1, we observe an example of FRAP data in form of a time series of both 2D and 1D fluorescence profiles reflecting the photosynthetic proteins mobility on a membrane [2, 3]. In the observed specimen (thylakoid membrane of the red algae Porphyridium cruentum, cf. Figure 1(a)), the experimentalists usually select a 2D rectangular ROI with the larger side perpendicular to the bleach strip in order to perform the averaging along the shorter axis. The resulting 1D fluorescence profiles (Gaussian-shaped signals in the central bleached region) are plotted on Figure 1(b). Based on spatiotemporal FRAP data, the diffusion constant is usually estimated by fitting a closed form model to the preprocessed FRAP data; see, for example, [49]. Since all closed form equation models need some simplified conditions to be fulfilled, for example, specific assumptions about initial and boundary conditions, they have limited accuracy. The simulation based models are more general although computationally more expensive, and they are using the raw data; the underlying process is modeled numerically; see, for example, PDE based model [10] or Monte Carlo simulations based model [11, 12], and model parameters are inferred within an optimization procedure.

Figure 1: Experimental data from FRAP experiments [2, 3]. (a) Representative FRAP image sequence for a single cell of red algae Porphyridium cruentum for phycobilisome fluorescence. First, a fluorescence image before bleaching was detected ((A) before bleach), and then the phycobilisome fluorescence was bleached out across the middle of the cell in the vertical direction (red rectangle) and the postbleach images were recorded every 8 seconds (5 plots are displayed only). The total measurement time was 200 seconds; the length of the scale bar is 3 μm. (b) One-dimensional fluorescence profiles. The abscissa represents the position along the axis perpendicular to the bleach strip. The ordinate quantifies the corresponding values of fluorescence (in arbitrary units) averaged along the axis parallel to the bleach strip. In the central region we see the step-wise recovery of the Gaussian-shaped signals, from the lowest values of the first postbleach profile to the highest prebleach (steady-state) values on the top.

Let us mention that CLSM allows the selection of a thin cross section of the sample by rejecting the information coming from the out-of-focus planes. However, the small energy level emitted by the fluorophore and the amplification performed by the photon detector is the cause of significant measurement noise. Therefore, FRAP images are in general very noisy with a small signal-to-noise ratio (SNR); that is, in order to get reliable results for the diffusion coefficient (often modeled as time and space invariant), an adequate technique residing in regularization is mandatory [13]. The analysis of the ill-posedness of the parameter identification of reaction-diffusion models based on spatiotemporal FRAP images was treated in [14, 15] and previously also in [16]. Currently, our approach is implemented into the special software CA-FRAP based on the UFO (Universal Functional Optimization) system [17]. For the real FRAP data with the red algae Porphyridium cruentum we obtained a good agreement with reference values (the range of result 10−14 m2s−1, i.e., 10−2μm2s−1) [3, 15, 18].

In this paper, we focus on another rather theoretical issue with great practical relevance residing in FRAP data (pre)processing. The key concept relies on sensitivity analysis; confer [19]. Sensitivity analysis, performed to test the sensitivity of the parameter inference of particular modelling conditions, is well established. However, we are using sensitivity analysis to reduce data set size prior to conducting parameter inference; that is, we calculate the sensitivities and afterwards select the “relevant” data space where the sensitivity is “sufficiently high.” The importance of this point is notably raised when the cost of data acquisition (besides considering the computational issues related to the “irrelevant” data processing) cannot be neglected. As far as we know, this approach is novel in both the FRAP-related and statistic literature. Furthermore, based on the same sensitivity analysis, we analyze two approaches of FRAP data processing: (i) the integrated data approach commonly used by the FRAP community and (ii) the full (spatiotemporal) data case. The inertia of the common FRAP data processing methods almost prohibits promoting new approaches with few exceptions [10, 20, 21].

Our paper is organized as follows. After the introductory section we describe the problem and introduce the sensitivity analysis in Section 2. Then in Section 3, we develop a new theoretical approach allowing a data space reduction. Furthermore, in Section 4, we investigate a common way of treating the data (in the FRAP community), that is, using the so-called FRAP recovery curve in form of a time series (integrated data approach), and we state a proposition exhibiting larger sensitivities (and hence smaller confidence intervals) of the parameter identification problem with full data compared to that with integrated data. The novelty and advantages of our approach as well as outlooks for further research are discussed in Section 5.

2. Preliminaries

2.1. Reaction-Diffusion System

Consider a model of reaction-diffusion system describing the fluorescent solute concentration in time and space. Right now we assume that two model parameters, (i) a diffusion coefficient and (ii) a reaction factor , are time-dependent and constant in space. The introduced reaction term reflects the continuous photobleaching during data acquisition and is proportional to the fluorescent particle concentration (which is proportional to the fluorescence signal measured by CLSM); that is, it is described as a first-order reaction [10]. The governing equation for the spatiotemporal fluorescence signal is a diffusion equation with reaction term as follows: Boundary conditions could be, for example, We also consider the simplest case of unbounded domains, , in which we set appropriate decay conditions at , . Equations (1) (and variants) are the basis for all the further analysis.

In the case of constant coefficients, the solution to this problem can be expressed by means of the Green function for the heat equation. Consider It is well-known that solves the problem for and and that solves (1) for constant parameters .

Some frequently used cases are that of a diffusion in free space, that is, on the whole domain without boundary conditions. In that case, the Green function (in 1 and 2 spatial dimensions, resp.) is the well-known heat kernel. Consider

In FRAP experiments, the initial condition, that is, the first postbleach profile (with the background or prebleach signal subtracted), is often modeled as a Gaussian [4, 7, 8] (cf. Figure 1(b)) which leads in the one- or two-dimensional case to initial conditions of the form where is the maximum depth at time (i.e., for ) and is the half-width of the bleach at height (depth) ; see, for example, [7]. The explicit solution for in the free space case is then given by

2.2. Parameter Identification Problem

We now discuss the parameter identification problem, where we try to infer about the parameters , by using direct measurements of in some space-time domain. That is, we assume that one of the following kind of data are observed: The first formulation (DD) is a discrete one, and the second one (CD) is an idealized continuous version of it. In the second one, the set is some observation (monitored) region, which does not have to be the full region . Of course, if the data grid is small enough, the discrete data case can be seen as an approximation to the continuous case. We will use this at several instances, for example, for an approximation of sums in the discrete case by the corresponding integrals.

Remark 1. A FRAP data structure usually consists of a time sequence of rectangular matrices (2D fluorescence profiles, cf. Figure 1(a)), where each entry quantifies the fluorescence intensity at a particular point in a 2D domain (e.g., by a number between 0 and 255). Consider where is the spatial index uniquely identifying the position where the signal is measured and is the time index (the initial condition corresponds to ).
Moreover, the data points are uniformly distributed in both time (the time interval between two consecutive measurements is constant) and space (on an equidistant 1D or 2D mesh). Nevertheless, in the following we adopt the simplified notation consisting in using only one index () for all data in the space-time domain.

We now define the forward map (also called parameter-to-data map) where is as in (5). Our regression model is now where the data are modeled as contaminated with additive white noise as follows: Here denotes the true coefficients and is a data error vector which we assume to be normally distributed with variance . Consider

Given some data, the aim of the parameter identification problem is to find such that (12) is satisfied in some appropriate sense. Since (12) usually consists of an overdetermined system (there are more data points than unknowns), it cannot be expected that (12) holds with equality, but instead an appropriate notion of solution (which we adopt for the rest of the paper) is that of a least-squares solution (with denoting the Euclidean norm on ). One has

The above defined parameter identification problem is usually ill-posed for nonconstant coefficients such that regularization has to be employed; see, for example, [13]. The solution of a practical example based on FRAP data was presented in [15]. If we a priori restrict the coefficients to be constant, then the identification problem becomes well-posed.

Further, in Section 3, we solve and afterward discuss the problem of data space selection based on sensitivity analysis. The prerequisite for this study is the confidence intervals estimation which we introduce in the next subsection.

2.3. Sensitivity Analysis and Confidence Intervals

For the sensitivity analysis we require the Fréchet-derivative of the forward map ; that is, A corresponding quantity is the Fisher information matrixAccording to Bates and Watts [22], we can estimate confidence intervals. Suppose we have computed as least-squares solutions in the sense of (15). Let us define the residual as Then according to [22], it is possible to quantify the error between the computed parameters and the true parameters . In fact, we have an approximate confidence interval as follows:where corresponds to the upper quantile of the Fisher distribution with and degrees of freedom.

The parameter can be either disregarded for its low value [4] or measured by another set of FRAP experiments. These new experiments allow the identification of parameters related to the fluorescence loss due to the bleaching during scanning; see, for example, [10]. Consequently, in case when is either fixed (known) or disregarded (), we only have as unknown. Then a similar result holds and the Fisher information matrix collapses into the scalar quantity , and the confidence interval for full observations is described as follows:with

Remark 2. In (20), several simplifications are possible. Note that according to our noise model, the residual term is an estimator of the error variance [22] such that the approximation holds for being large. This approximation will be used in several instances in the following, and the term in (20) can be viewed as rather independent of or .
Moreover, we remind the reader that the Fisher distribution with and degrees of freedom converges to the -distribution as . Hence, the term can approximately be viewed as independent of as well and of moderate size.

3. Theoretical Results on Data Space Selection

In this section we present the theoretical results concerning an (in further defined sense) optimal data space selection for the problem of parameter identification based on FRAP data. For all the results following, we always consider only the constant diffusivity as unknown; set . Hence, (20) is the central estimate which we use in the sequel. We mention that most of the analysis can be extended to the case of two unknowns as well. However, the case of an unknown nonconstant diffusivity is out of the scope of this paper.

3.1. Data Reduction: Relevant versus Irrelevant Data

The confidence interval estimate (20) gives us useful information on the quality of a least-squares estimate. A central ingredient in this estimate is the sensitivity (or the Fisher information matrix), which is the main factor controlling the error . As a next step, starting from (20), we might change the data observation and consider the question if we can use less data without increasing the error too much.

Let us assume now that we have different kind of data (i.e., less data points, also called reduced data in the following). One has where and we compute a least-squares estimate using these data. Without loss of generality let us impose the following ordering of the data , that is, that the new data are just the first values of the original data. We generally impose the condition that is large enough such that certain estimators (like ) are close to their expected value, that is, that a law of large numbers is reasonable. As a rule of thumb, we impose in all of the following that . The case of extreme data reduction to a few points requires further more detailed analysis, which we omit here.

The corresponding confidence interval estimate (20) still holds with a (usually larger) confidence interval Note that for the residual we now have to take the Euclidean norm in in the corresponding expression (21) (i.e., a sum of the form ).

Suppose we would like to have a similar confidence interval as before with the new data set instead. We can thus estimate using (25) Let us further introduce an upper bound such that then we finally find a confidence interval using the reduced data. It will be announced in the following lemma.

Lemma 3. Let be the number of reduced data points and suppose that is chosen according to (27). Then using reduced data, one finds a confidence interval of the form

Remark 4. The relevance of this observation is that we have almost the same type of estimates as in the full data case (20), but we might ignore irrelevant data points. The Fisher distribution converges to the -distribution as and thus it holds that if both and . For reasonable values of and (e.g., and ), it can be verified that . This situation corresponds to our initial assumption of sufficiently large . Note that a similar analysis as below could be done if this assumption is not valid by incorporating the factor into the analysis.
Thus, in practice and for the later analysis we replace the estimate (28) by the simpler one (forgetting the first factor). Consider

Remark 5. The previous Lemma 3 can be seen as a criterion for selecting . If this is chosen as an upper bound, then (27) is relating to . Then this will be the main condition for choosing (as small as possible such that (27) holds). With this choice and if is small, we basically obtain almost the same confidence interval than when using the full data. Hence, out of the original data we can pick a subset which contains almost the same information. If (27) holds with small , we thus call the new data the relevant data and the complement the irrelevant data.

Let us mention, following Remark 2, that the expressions and are both similar to for and are large. More precisely, the estimator has variance , and hence the ratio of these two terms is roughly of the order of With our assumption of large data sizes, , and for reasonable , the two expressions can hence be considered equal in our analysis. This explains, why we considered the right-hand sides of the corresponding estimates (20) and (28) (resp., (29)) almost equal. In fact, our analysis is rather independent of the and (as long as we impose the condition of being large as above) because we compare two situations (full data and reduced data) with similar and the corresponding factors cancel out. The dependence of the confidence intervals itself on the noise variance is, of course, well known and (in the well-posed case) linear; see, for example, [6, 9] and Figure 4.

3.2. Data Reduction Based on the Data Noise Level

There is also the possibility of using other selection criteria as well for the data reduction step. The following one is based on comparing the sensitivity with the noise level. In general, for deterministic errors, the data noise level refers to the norm of the noise. In our case, that is, that of stochastic errors, we identify the noise level with the noise variance .

Lemma 6. Let be the number of reduced data points and let hold. (This means that the relative difference with respect to is less than 100%.) Furthermore, we assume the following selection criterion with some freely selected parameter : Then using reduced data, we find a confidence interval of the form

Proof. The proof follows the steps of that preceding Lemma 3. We obtain the estimate in Lemma 6 by using assumptions (30) and in the second part of (26) as which immediately yields the result.

Remark 7. If the approximation of residual (22) is valid and if the ratio of as in Lemma 3, then Lemma 6 leads to a confidence interval of the order which is again of the same order as the full data case; the length of the confidence interval in the reduced data case is multiplied by the factor compared to the full data case.

Remark 8. The previous analysis can be extended to the case of an unknown reaction rate as well. The selection criterion corresponding to Lemmas 3 and 6 should be formulated as respectively.

Example 9. In this example we underline the previous theoretical results by numerical computations. To do so, we compute a least-squares estimate using different kind of data. In a first experiment we consider a rectangular spatiotemporal data grid with space interval and time interval . The grid size in both space and time direction was set to . We simulated data by assuming an exact diffusion coefficient with bleach radius and computed the data for the 1D case by (8). These data were perturbed by normally distributed additive noise (white Gaussian noise) with standard deviation . Based on these data we computed a least-squares estimate of the diffusion coefficient using a Gauss-Newton procedure. We refer to this setup as the full data case.
Next, we computed a similar estimate using less data. For this, we considered the data point selection criterion (27) for various values of . The data regions were selected by first ordering all data points according to their sensitivities Then we selected the reduced data from the region of points , using this ordering and the smallest such that (27) holds for the choice of . In a sense, we take the most sensitive data points as relevant data. In Figure 2(a), the gray region corresponds to that region where the data points were selected and the complementary region, the white one, represents irrelevant data that can be neglected without much loss of quality. A similar process (using the same sensitivity ordering) was performed using selection criterion (30) with values of . The corresponding data regions are displayed in Figure 2(b). Below of each relevant data region, the ratio of is shown; that is, represents the percentage of data reduction. For instance, in the first row, rightmost plot, , only 50% of the original data are used.
Having identified several regions of relevant data using different selection criteria, we considered the quality of the computed estimate using the reduced data. For the seven cases (full data and six cases of reduced data), we sampled 1000 times a normally distributed noise vector (with fixed standard deviation ), and for each of the noise realizations, we computed a least-squares fit by a Gauss-Newton iteration and computed the squared error This yields 1000 samples out of the error distribution of . For the seven cases, a boxplot of this error can be seen in Figure 2(c). Here the box interior line corresponds to the mean of the (over 1000 realizations) and the edges of the box in each column represents the 25th and 75th percentiles. The cross marks + are considered outliers. The size of the box (more precisely the top of the box) can be interpreted as an indicator of the size of a confidence interval.
This plot corresponds quite well with our theoretical findings. For instance, the case , in the second column of the boxplot, corresponds to the case of using only 10% of the original data, that is, 90% of data reduction. According to Lemma 3, the theoretical confidence interval bound is roughly times larger compared to the full data case. And indeed, we observe the box in the second column in Figure 2(c) being roughly 2 times larger compared to the full data case in column one. The other cases with and behave similarly. Although the boxes in the last three columns for are within the theoretical bound in Lemma 6, it seems that this bound is not very precise, since the actual error is quite smaller than the bound in (31). In that respect, the data selection method of Lemma 3 is more accurate and preferable.
Note that for the case (fourth column), we have virtually the same error distribution as for the full data case even though we have only used 50% of the original data.
There is a rather esthetic issue of the lotus flower-like plots in Figure 2(a); what is the meaning of these lotus-like grey regions? The answer is simple: those regions exhibit the sensitivity to the fluorescence loss in photobleaching (FLIP); see [18] for review and [23] for a more creative exploration of the FRAP and FLIP synergic effect.

Figure 2: Examples of data reduction. Gray color indicates the data space, where the relevant data are taken. White color represents the data regions of irrelevant data points. (a) Relevant data region using (27) with . (b) Relevant data region using (30) with . (c) Boxplot of the distribution of the squared error based on 1000 noisy signal samples using full and reduced data; cross marks (+) indicates the outliers. The data were simulated using (8), with an exact diffusion coefficient , bleach radius , and additive white Gaussian noise with standard deviation .

Remark 10. In Example 9, we have—in a sense—selected the optimal region by ordering the sensitivities. However, the data selection criteria of Lemmas 3 and 6 do not depend on this ordering but can be performed in different situations. In fact, in practice, one probably wants to select a rectangular region as a relevant data region instead of the complicated one in Figure 2. This is of course possible in Lemmas 3 and 6. One can compare all subregions of rectangular form and accept those which satisfy the criterion in the above lemmas. Such a concept will be pursued in the following subsections.

Remark 11. Another way of reducing the amount of data could be to use a coarser data grid; that is, instead of , one might perform a subsampling and just use every second grid point (in space or time direction). However, such an idea does not fit into our context because a subsampling will in general increase the confidence interval significantly. Looking at (25), one can observe that a subsampling by taking every second data point will lead to Fréchet-derivative which is about half the size of the original one and thus leading to a doubling of the confidence intervals. On the contrary, in our framework we consider cases where the confidence intervals stay almost the same. Thus subsampling will lead to an which is not of the order of . Still, such a procedure makes sense, if one has an idea of the wanted size of the confidence intervals. Then the grid size can be selected such that estimate (25) fits to one’s needs. For a rectangular grid on say a space-time region , this requires an estimate of , which can roughly be obtained by the asymtotics as . Thus from estimates of the integral on the right-hand side, can be adapted to reach a desired confidence.

3.3. Data Reduction for the Free Space Case

We now study the data reduction procedure in view of practical relevance. That is, we want to find simple formulas for selecting a rectangular spatiotemporal region which yield approximately similar error bounds than when using all the spatial points as data.

Let us consider the problem on free space with and unknown and a Gaussian as in (7) as initial condition. From the solutions (8) in the one- and two-dimensional case, we obtain by some calculations for the derivativeBy an appropriate scaling based on the experimental design factor (yielding the new spatial dimensionless variable ) and on the diffusion coefficient (yielding the second dimensionless variable ) we may simplify the expressions above.

Lemma 12. In the setting as before one can express the sensitivity as with the functions

In the one-dimensional case let us further assume the usual situation that full data are given on a space-time cylinder , where is the space interval of observations and is the time interval as follows: Here and are the spacings of the grid in space and time direction, respectively. The sensitivity can be approximated by an integral Forgetting the factor (which cancels out anyway later), we find withand the error function .

Analogously in the two-dimensional case, we can consider data on a space-time cylinder , yielding that with

A particular case happens if in the above integrals. Then respectively,

The data reduction process involves estimates of the sensitivities and hence of as in (27). The previous calculations lead to the following lemma, which allows for data reduction in the spatial variable compared to full spatial information; see Example 14 below.

Lemma 13. Let and be defined as in (43) and (45), respectively. Let and . If thenwhere

Proof. With (48) we have that for all . With , (43) reads where we used the fact that is monotonically increasing for (This can easily be seen from its derivative.) The same proof holds with in place of . Note that is again monotonically increasing for .

This lemma leads to simple rules for selecting the length of data intervals with respect to the spatial dimension, as the following example explains.

Example 14. Let us now again exemplify the data reduction procedure. Suppose we want to find the relevant data on a space-time cylinder which yields a comparable confidence interval to full data given on the whole spatial domain with being large and with equal grid spacings . Assuming being large, we can approximately view the full data interval as and in the following treat this as the case when the full data are given on the whole space interval.
The condition we impose now is (27), with some chosen small . We have the situation that the ideal corresponds to data on , while the actual data are given on . Condition (27), assuming that the data points are dense such that the sums can be approximated by integrals, can be stated as Rearranging terms and with the notation as before we come to the condition or equivalently We are now in the situation, where we can apply Lemma 13. Using (48) with and and given , we have to find a such that , respectively, This inequality can be solved by numerical means; for instance, if we allow for a 1% increase of the confidence interval, we take and find numerically a value of . For (0.1% increase) we find a value . Thus with the above setting of variables, we have the following recommendation for the length of data interval :which leads to the confidence intervals that are almost identical to the case of full spacial data.
A similar procedure holds in the two-dimensional case. If instead of full data on we want to take only relevant data on a cylinder , we have to find a such that (55) holds with instead of . A numerical examination of this function reveals the recommendation for the spatial radius of the data-cylinder. One hasNote that in these examples we only restrict the spatial coordinates while the time interval for the full data case and the restricted data case is the same. The question comes up if it is possible to reduce the time interval as well and to find a smaller interval such that (27) holds with a small and hence yielding similar confidence intervals. Unfortunately it is not possible to do so without some loss in precision of the resulting parameter estimates. In fact, suppose that the spatial interval is chosen as in (56), then we can approximate the sensitivity on the interval well by . Hence (27) holds using this approximation if According to Lemma 3, the ratio of the confidence interval using data on and is then approximately given by The last ratio can be easily evaluated using (46). It turns out that this ratio is of the order of only if . Thus it is always better to have a larger time interval. Still, if one can quantify the “costs” of using a large time interval and adds these to the ratio in (59), this can lead to an optimization problem in which case an optimal time interval of data measurements can be easily calculated.

Example 15. Let us illustrate the findings of this subsection by numerical examples. The main point in the previous analysis was that a very large spatial data domain can be replaced without loss of quality by a smaller one using formula (56). For the results in Figure 3, we first considered data given on a rectangular grid with large spatial size. The full data case corresponds to data given on the spatial interval and on the time interval . The spacing of the grid was in space direction and in time direction. For all the computations in this example we used an exact diffusivity and normally distributed additive random noise with . As in Example 9, we computed an estimate for the diffusivity by using different kind of data. For each data region we used 1000 samples of noise yielding the same number of squared errors .
We considered seven different kinds of rectangular data regions. The corresponding data regions (with labels (I)–(VII)) are indicated in Figure 3(a) by different shadings. The corresponding boxplot of the error is given in Figure 3(b).
The first boxplot in Figure 3(b) (labeled as (I)) corresponds to the full data case with data given on the full grid on . For the second boxplot (labeled as (II)) we used reduced data on an interval with computed by (56) which yields . Note that the time interval here is the same as that for the full data case. According to the above analysis such a procedure with (56) should not change the confidence interval significantly. Indeed, we observe in Figure 3(b)that the cases (I) and (II) are almost equal. Note that in case (II), we only use 54% of the full data with similar results, which again indicates the big amount of irrelevant data in the full data case.
The plots on the right, labeled as (III)–(VII), correspond to a choice of a smaller time interval with the corresponding length of the spatial interval given by (56). We have already indicated that it is always better to have a larger time interval, because a smaller time interval will yield a larger confidence interval according to (59).
In case (III) we considered a time interval with which yields a very small data region indicated by the (hardly visible) darkest color in the right plot close to the origin of space-time coordinate system. Case (IV) corresponds to , case (V) corresponds to , and cases (VI) and (VII) correspond to and , respectively. The according spatial interval is set by (56). The boxplots on the right-hand side correspond to the theoretical findings above. Note that case (III) has a too large error to fit to the scaling of the plot (hence, there is no corresponding box in Figure 3(b)). We observe that all the boxes in cases (III)–(VII) have larger confidence intervals than those in the full data case (I). Only the last case (VII) is similar to the one in the full data case (I). This shows that a reduction of the time interval (Another issue to be considered is the enlargement of the time interval between two consecutive measurements as the FRAP experiment proceeds.) has to be paid by an enlargement of the confidence intervals.
In Figure 4 we performed the same experiments as in Figure 3 for some selected regions, but with different noise variances . The first 5 plots correspond to the confidence intervals of Figure 3(b) for regions (I), (II), (IV), (V), and (VII), but now with a different for each plot. From left to right and top to bottom we have It can be observed that the relative size of two confidence intervals corresponding to the same region, that is, (I), (II), (IV), (V), and (VII), is almost the same in all cases and hence rather independent of . Of course, the absolute size depends on , as can be seen from the rightmost bottom plot, where we indicate the intervals for region (II) versus different . The roughly linear dependence of the mean of the on is in agreement with theory; see (20) with (22).

Figure 3: Comparison of the squared error for 7 rectangular spatiotemporal regions. (a) Data regions for the cases (I)–(VII); full data are represented by the white region , case (I). Reduced data, in cases (II)–(VII), are indicated by different shadings. (b) Boxplot of the error corresponding to the use of data from the respective regions.
Figure 4: Dependence of on the data noise level. Boxplot of error for regions (I), (II), (IV) (V), and (VII) for 5 different . (f) shows the confidence intervals for region (II) plotted against the noise variance . Note: be aware of different scaling along the vertical axes of each of 6 plots.

4. FRAP Recovery Curves: Integrated Data Approach

In this section we investigate the different way of treating the data; in the following we call it the integrated data approach, which is a common procedure in the evaluation of FRAP experiments [1, 12]. Almost without exception, the experimental biologists are using for their purposes the so-called FRAP recovery curve, the space averaged signal (; see below), instead of the spatiotemporal data . Our aim is to compare both approaches in view of sensitivity and confidence intervals.

By integrating the data, we have to take into account that the data error is integrated as well. This will in general reduce the variance and this fact causes the false opinion about the higher precision of the integrated data approach.

Proposition 16. Let us consider that the full data , , are given on a rectangular spatiotemporal grid where is the number of points in the ROI (either 2D or 1D) and is the number of time points. The (space-)integrated data are given on a temporal grid. Assume that the data error is normally distributed and that the number of data points is large enough such that is almost equal to the data error variance in either case and that the Fisher-quantile satisfies .
Then the full data approach has equal or smaller confidence interval bounds for the parameter estimates of .

Proof. Let us introduce for the integrated data case the parameter-to-data mapping As the spatial errors are distributed according to , then for the integrated errors it holds that and we have the regression model Note that the error variance is now reduced by a factor
After a regression using this model, we may employ a similar estimate for the confidence interval as before (20), Note that the sensitivity for the integrated data approach is (after changing the order of sum and derivative)Compare this to the corresponding expression in the full (nonintegrated) data case, which reads With the assumption that the residual error corresponds to the error variance, we find We may thus compare the sensitivities for the two approaches: Now using the Cauchy-Schwarz inequality in the last line and the assumption , we find that the upper bound for the integrated data case is always larger than that of the full data case. Consider which finishes the proof.

An explanation for this proposition is that although the data error variance is smaller for the integrated data case, the sensitivity is smaller as well, and this fact causes the greater or equal confidence interval bounds.

Remark 17. The full data case always leads to a smaller confidence interval and is thus preferable. For the one-dimensional ROI, with the uniform spatiotemporal grid where is the number of points on the space axis, is the number of time points, and is the time step and for the the reasonable case that , we found the empirical estimate This inequality can be partially verified by analytical considerations. Indeed, the ratio on the left of (72) is roughly the ratio between the corresponding sensitivities in (69). Using a large data-size approximation, that is, approximating the sums by integrals, we can evaluate the corresponding integrals as in (43) and find withwhere . A plot of the function reveals that it is close to if , close to if , and monotonically decaying from to for . Thus if , then . In particular if , then in the integrand is close to and the intervals are almost the same. This corresponds to the right-hand side of (72) being , hence verifying (72). Conversely if , then is almost , and the ratio in (72) is close to . The case corresponds to . Within our assumptions of this leads to the condition and hence the right-hand side of (72) being small () as well.

Remark 18. As the proof of Proposition 16 is based on the Cauchy-Schwarz inequality, we can consider the cases in which it is sharp, that is, when the left-hand side of (69) is close to the right-hand side. This happens when the two vectors in the inequality are collinear; in our case this means that for all . Thus, by this reasoning, if the sensitivity does not vary much in the spacial direction, then the confidence intervals in the full data and integrated data case should be approximately equal. We can roughly interpret this case as that where all data points are highly relevant and no data reduction is necessary. Conversely, if the sensitivities are not constant in spacial direction, then the full data case yields significantly smaller confidence intervals.

Example 19. Again, we underpin the theoretical study in this section by numerical experiments. For the plots in Figure 5, we compared the squared errors for the full data case and the integrated data case for four kinds of data regions (labeled as (I)–(IV)). As data we used again a rectangular region with different values of . The grid spacing , the value of the exact diffusivity , and the noise variance were the same as in Examples 9 and 15. As before we used 1000 samples for each of cases (I)–(IV) and for both the full and the integrated data we computed an error distribution of .
The length of the spatial interval was as follows: for case (I): , for case (II): , for case (III): , and for case (IV): . The time interval was the same for all four cases.
In the plot in Figure 5, we display four pairs of boxplots of the squared error corresponding to the four cases of data regions. In each of the pairs, the left boxplot corresponds to using full data as before, while the right boxplot corresponds to using spatially integrated data (60) on the same region. The error plots are in full accordance with Proposition 16, because in each pair, the left box is smaller or equal to the right box, which means that the error using full data (left) is always smaller or equal to that of using integrated data (right). The results also fit quantitatively to the estimate in (72). Note that the factor is equal to , for cases (I), (II), (III), and (IV). Thus, in case (I), we have such that (72) predicts a comparable error for the full and integrated data approach. Indeed, the first pair in the boxplot has comparable sizes of the boxes. In the other three case we have ; hence (72) predicts a larger error for the integrated data case. This can be observed again for cases (II)–(IV) where the integrated data case has a larger error with increasing value of . The explanation of such a behavior resides in the inclusion of more and more larger region of irrelevant data into the spatial integration.

Figure 5: Comparison of the squared errors for the full data approach and the integrated data approach for four data regions with different length of the spatial interval (the factor is equal to , resp.) labeled as (I)–(IV); time interval was the same for all cases . In the pairs, the left and right boxplot correspond to using full data (Full) and spatially integrated data (Int), respectively.

Remark 20. Unexpectedly, Figure 5 shows another interesting phenomenon, from left to right, with the growing ROI (more precisely for increasing ratio ); while the squared error for the full data approach monotonically decreases, there is a minimum for the integrated data approach; that is, there is an optimum (maximum sensitivity) for some value . Indeed, when the ROI is too small we disregard a huge amount of the relevant data, and conversely, when the ROI is too big, the integrated approach kills the valuable spatial resolution causing the loss of sensitivity of the output data to the parameter value. Hence, this is an explanation of the existence of an optimal size of ROI (i.e., an optimal value of ) for the integrated data approach.

5. Conclusion and Future Directions

In this paper, we have introduced and solved the problem of data space selection-reduction and data processing for parameter identification in a reaction-diffusion model based on FRAP experiments. First, either for a known data noise level (determined by parameter ) or for previously chosen coefficient of enlargement of the confidence interval for the parameter estimate (parameter ), we propose a new method to separate the relevant and irrelevant data sets.

Afterwards, we show that the data set represented by the FRAP recovery curves (the integrated data approach) leads to a larger confidence interval compared to the spatiotemporal data. This rigorous statement is worth to be promoted in the FRAP community because it goes against the common knowledge and practice.

Four numerical examples support our conclusions. The corresponding graphical visualizations, despite their simplicity, help us in discovering new features of FRAP and FLIP techniques and inspire new investigation related to the optimal design of photobleaching experiments. For instance, the lotus flower-like plot (Figure 2) led us to the investigation related to the combination of FRAP and FLIP experiments [23]. The next step in the same direction was done in our conference paper [24], where we use the sensitivity analysis developed here to find specific conditions for optimizing the radius of bleach spot. Again, the inspiration to optimize a design factor was induced by Figure 5; that is, the existence of an optimal size of monitored region for the integrated data approach(cf. Remark 20) led us to look for an optimal value of bleach size.

Once proving the concept, our future goal is even more ambitious. Based on the same sensitivity analysis we shall suggest the optimal values for the experimental design variables, that is, the experimental protocol modifications.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This work was supported by the OeAD (Austrian Agency for International Mobility and Cooperation in Education, Science and Research) within the programme “Aktion Oesterreich-Tschechien (AOeCZ-Universitaetslehrerstipendien)’’ and by the Ministry of Education, Youth and Sports of the Czech Republic, Projects CENAKVA (no. CZ.1.05/2.1.00/01.0024) and CENAKVA II (no. LO1205 under the NPU I program).


  1. F. Mueller, D. Mazza, T. J. Stasevich, and J. G. McNally, “FRAP and kinetic modeling in the analysis of nuclear protein dynamics: what do we really know?” Current Opinion in Cell Biology, vol. 22, no. 3, pp. 403–411, 2010. View at Publisher · View at Google Scholar · View at Scopus
  2. R. Kaňa, O. Prášil, and C. W. Mullineaux, “Immobility of phycobilins in the thylakoid lumen of a cryptophyte suggests that protein diffusion in the lumen is very restricted,” FEBS Letters, vol. 583, no. 4, pp. 670–674, 2009. View at Publisher · View at Google Scholar · View at Scopus
  3. R. Kaňa, E. Kotabová, M. Lukeš, and et al, “Phycobilisome mobility and its role in the regulation of light harvesting in red algae,” Plant Physiology, vol. 165, no. 4, pp. 1618–1631, 2014. View at Publisher · View at Google Scholar
  4. D. Axelrod, D. E. Koppel, J. Schlessinger, E. Elson, and W. W. Webb, “Mobility measurement by analysis of fluorescence photobleaching recovery kinetics,” Biophysical Journal, vol. 16, no. 9, pp. 1055–1069, 1976. View at Publisher · View at Google Scholar · View at Scopus
  5. J. Ellenberg, E. D. Siggia, J. E. Moreira et al., “Nuclear membrane dynamics and reassembly in living cells: targeting of an inner nuclear membrane protein in interphase and mitosis,” The Journal of Cell Biology, vol. 138, no. 6, pp. 1193–1206, 1997. View at Publisher · View at Google Scholar · View at Scopus
  6. P. Jönsson, M. P. Jonsson, J. O. Tegenfeldt, and F. Höök, “A method improving the accuracy of fluorescence recovery after photobleaching analysis,” Biophysical Journal, vol. 95, no. 11, pp. 5334–5348, 2008. View at Publisher · View at Google Scholar · View at Scopus
  7. C. W. Mullineaux, M. J. Tobin, and G. R. Jones, “Mobility of photosynthetic complexes in thylakoid membranes,” Nature, vol. 390, no. 6658, pp. 421–424, 1997. View at Publisher · View at Google Scholar · View at Scopus
  8. S. Seiffert and W. Oppermann, “Systematic evaluation of FRAP experiments performed in a confocal laser scanning microscope,” Journal of Microscopy, vol. 220, no. 1, pp. 20–30, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  9. T.-T. Tsay and K. A. Jacobson, “Spatial fourier analysis of video photobleaching measurements: principles and optimization,” Biophysical Journal, vol. 60, no. 2, pp. 360–368, 1991. View at Publisher · View at Google Scholar · View at Scopus
  10. O. N. Irrechukwu and M. E. Levenston, “Improved estimation of solute diffusivity through numerical analysis of FRAP experiments,” Cellular and Molecular Bioengineering, vol. 2, no. 1, pp. 104–117, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. J. Mai, S. Trump, R. Ali et al., “Are assumptions about the model type necessary in reaction-diffusion modeling? A FRAP application,” Biophysical Journal, vol. 100, no. 5, pp. 1178–1188, 2011. View at Publisher · View at Google Scholar · View at Scopus
  12. I. F. Sbalzarini, Analysis, Modeling and Simulation of Diffusion Processes in Cell Biology, VDM Dr. Muller, 2009.
  13. H. Engl, M. Hanke, and A. Neubauer, Regularization of Ill-Posed Problems, Kluwer, Dortrecht, The Netherlands, 1996.
  14. R. Kaňa, C. Matonoha, Š. Papáček, and J. Soukup, “On estimation of diffusion coefficient based on spatio-temporal FRAP images: an inverse ill-posed problem,” in Programs and Algorithms of Numerical Mathematics 16, J. Chleboun, K. Segeth, J. Šístek, and T. Vejchodský, Eds., pp. 100–111, Academy of Sciences of the Czech Republic, 2013. View at Google Scholar · View at MathSciNet
  15. Š. Papáček, R. Kaňa, and C. Matonoha, “Estimation of diffusivity of phycobilisomes on thylakoid membrane based on spatio-temporal FRAP images,” Mathematical and Computer Modelling, vol. 57, no. 7-8, pp. 1907–1912, 2013. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  16. K. S. Zadeh, H. J. Montas, and A. Shirmohammadi, “Identification of biomolecule mass transport and binding rate parameters in living cells by inverse modeling,” Theoretical Biology and Medical Modelling, vol. 3, article 36, 2006. View at Publisher · View at Google Scholar · View at Scopus
  17. L. Lukšan, M. Tuma, C. Matonoha et al., “UFO 2014—interactive system for universal functional optimization,” Tech. Rep. V-1191, Institute of Computer Science, Academy of Sciences of the Czech Republic, Prague, Czech Republic, 2014, View at Google Scholar
  18. R. Kaňa, “Mobility of photosynthetic proteins,” Photosynthesis Research, vol. 116, no. 2-3, pp. 465–479, 2013. View at Publisher · View at Google Scholar · View at Scopus
  19. A. Cintrón-Arias, H. T. Banks, A. Capaldi, and A. L. Lloyd, “A sensitivity matrix based methodology for inverse problem formulation,” Journal of Inverse and Ill-posed Problems, vol. 17, no. 6, pp. 545–564, 2009. View at Publisher · View at Google Scholar · View at MathSciNet
  20. K. Sadegh Zadeh and H. J. Montas, “A class of exact solutions for biomacromolecular diffusion-reaction in live cells,” Journal of Theoretical Biology, vol. 264, no. 3, pp. 914–933, 2010. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  21. K. Sadegh Zadeh, “A synergic simulation-optimization approach for analyzing biomolecular dynamics in living organisms,” Computers in Biology and Medicine, vol. 41, no. 1, pp. 24–36, 2011. View at Publisher · View at Google Scholar · View at Scopus
  22. D. M. Bates and D. G. Watts, Nonlinear Regression Analysis and Its Applications, John Wiley & Sons, New York, NY, USA, 1988. View at Publisher · View at Google Scholar · View at MathSciNet
  23. Š. Papáček, J. Jablonský, C. Matonoha, R. Kaňa, and S. Kindermann, “FRAP & FLIP: two sides of the same coin?” in Proceedings of the 3rd International Conference on Bioinformatics and Biomedical Engineering (IWBBIO '15), Granada, Spain, April 2015, F. Ortuño and I. Rojas, Eds., vol. 9044, pp. 444–455, Springer, 2015. View at Google Scholar
  24. Š. Papáček, J. Jablonský, R. Kaňa, C. Matonoha, and S. Kindermann, “From data processing to experimental design and back again: a parameter identification problem based on FRAP images,” in Proceedings of the 17th International Conference on Digital Image Processing (ICDIP '15), Dubai, UAE, January 2015.