fitFCS
A Matlab program to fit FCS autocorrelation curves.
Written by Peter Nagy (peter.v.nagy@gmail.com)
https://peternagyweb.hu
https://peternagygroup.com
Synopsis
The program fits FCS autocorrelation curves using a model involving 2-D and 3-D diffusion components and a triplet component. The program requires the Optimization Toolbox to run. Two advanced fitting algoritms, global search and simulated annealing, are only available if the Global Optimization Toolbox is installed in Matlab.
The program can be executed in two modes: command-prompt or graphical user interface (GUI)-mode.
GUI-mode
Type 'fitFCS' at the Matlab command prompt without any input or output arguments. The following GUI will be dsiplayed.
Variable names: You have to enter the names of variables containing time, the autocorrelation curves (one or more curves columnwise) and the weights for the fitting (one or more, columnwise). Variables can be Matlab arrays or tables. If a table variable type is used as an input for the autocorrelation curves, the table headers will be used for identifying the curve in the output. If the struct variable output of 'extractCurvesFromZeissFCS' or 'extractCurvesFromNikonFCS' is used, and multiple fields are to be concatenated, type the following in the variable field to import all curves in the 'ac' field of all conditions (assuming 'curvesNikon' is the name of the structure variable): [curvesNikon.ac]
Providing weights is optional. The squared difference will be multiplied by the weight during fitting. If weights are not specified, every value will have the same weight. Once the fields for variable names are filled in, press the "Load" button to import them. The curve (or the first curve if multiple ones are present) will be plotted. If multiple curves are present, you can select which one is plotted in the "Which curve" drop-down menu.
Batch processing: If this option is selected, all curves in the variable specified for the autocorrelation curves will be fitted. Only the curve specified in the "Which curve" drop-down menu will be displayed. After the fitting, not only the original curve, but the fitted one is shown as well.
Feed results to initial box: If this option is selected, the fit results will overwrite the initial values in the Parameters window.
Diffusion components: You can specify if the diffusion correlation time or the diffusion coefficient is present in the model. If the diffusion coefficient will be fitted, the radius of the point-spread function in the X-Y plane, parameter ω, will have to be included in the model as well.
Load/Save model: The model specified in the program can be saved to a structure variable in the Matlab workspace. The structure variable will have two fields, modelType and parameterList. This structure variable can be loaded with the Load model button, and it can also be used when running the program in command-prompt mode (see later).
Fitting algorithm: Three fitting algorithms are implemented.
- Local minimization: Once reasonable initial values for the fitted parameters have been found, this algorithm is both good and fast.
- Global search: The slowest, but most reliable algorithm that fits the measured autocorrelation curve even with relatively bad initial values.
- Simulated annealing: A fast global optimization algorithm that is not as good as the "Global search" option for finding perfect fits.
Model: You can add 2-D or 3-D diffusion components (up to 3) and one triplet component to the model by ticking the "Add" box and specifying the number of such components in the corresponding drop-down menu. The model is based on the following equations:
where
- G is the autocorrelation function
- I(t) and I(t+Δt) are the intensities measured at time t and Δt, respectively
- N is the number of particles in the confocal volume
- Ti the fraction of the corresponding triplet / non-fluorescent component
- τtriplet, τ2D and τ3D are the lifetime of the triplet / non-fluorescent component, the 2D diffusion correlation time and the 3D diffusion correlation time, respectively
- α is the anomalous exponent characterizing the deviation from normal diffusion. α=1 corresponds to normal diffusion.
- fi and fj are the fractions of the corresponding diffusion components
- κ is the beam shape parameter, i.e. the ratio of the Z-dimension and the horizontal dimension of the confocal volume
- ω is horizontal radius of the confocal volume
- D is the diffusion coefficient
As it can be seen, the offset of the autocorrelation function at t=∞ is 1. If the autocorrelation function to be fitted approaches zero at t=∞, you have to add one to it.
The parameters window on the right will be populated according to the model. For every parameter, you can specify whether it is a constant or a fitted parameter. If it is a constant, its value must be specified. If if is a fitted parameter, its initial value as well as lower and upper bounds must be given.
The fitted curve will be displayed in the same graph as the original data. By sliding the markers you can specify which part of the original data set will be used for fitting.
Output: The program generates three different kinds of ouput.
- A table, shown below, containing all the fitted parameters. The name, the value and the SD of the parameters are shown in the first three columns. If 'IsConstant' column contains 1, the correspoinding parameter wasn't fitted. If multiple autocorrelation curves were fitted, the DatasetID shows which curve the parameter corresponds to. If multiple diffusion components are present, the fraction of mobile entities diffusing according to the corresponding component is labeled by Fd2d_i or Fd3d_i. The last diffusion component does not have an associated fraction, since its fraction is equal to 1-the sum of all other fractions. The last two rows contain the squared deviation between the fitted and the orignal data (SS) and the degree of freedom (df).
- Figure structures: all the figures containing the original data and the fitted autocorrelation function, generated during the fit, are stored in a structure variable. In order to regenerate the plot, enter the following command at the Matlab command prompt:
struct2handle(fs(1),0);
where fs is the structure variable generated during batch fitting, and the command above specifies that the result of the first fit is to be shown. '0' specifies that the generated figure window will be the daughter of the root object.
- Plot data in table(s) is a cell array with a 6-column table for each dataset (if there are more datasets) or a single 6-column table (if there is only one dataset).
- 1st column: time (original resolution, full range)
- 2nd colmn: measured autocorrelation function plotted for time points in column 1
- 3rd column: time (original resolution, fitted range)
- 4th column: residuals (original resolution, fitted range)
- 5th column:
time for the fitted autocorrelation function (full range)
- 6th column: fitted autocorrelation function
Fitting is performed by clicking on the "Fit" button.
Command-prompt mode
Syntax:
8 input arguments:
[fittedParams,parameterNames,treatAsConstant,resultTable,figureStructures,plotData] = fitFCS(time,ACcurve,weights,range,modelType,parameterList,typeOfFittingAlgorithm,calculateD);
7 input arguments (calculateD is assumed to be false in this case):
[...]=fitFCS(time,ACcurve,weights,range,modelType,parameterList,typeOfFittingAlgorithm);
7 input arguments:
[...]=fitFCS(time,ACcurve,weights,range,modelStructure,typeOfFittingAlgorithm,calculateD);
6 input arguments (calculateD is assumed to be false in this case):
[...]=fitFCS(time,ACcurve,weights,range,modelStructure,typeOfFittingAlgorithm);
Input parameters:
- time: time column vector (array or table)
- ACcurve - autocorrelation curves culomnwise (array or table)
- weights - weights columnwise (array or table). The squared difference will be multiplied by this weight during fitting. If no weighing is to be
performed, provide an emtpy variable for weights, e.g. fitFCS(time,ACcurve,[],range,...)
- range - [lowestIndex, highestIndex] to fit
- modelType - a cell array in which the 1st, 3rd, 5th, etc. cell is a string defining a model component followed by a number defining how many such components are to be added to the model. Allowed model components: 'diff2d', 'diff3d', 'triplet'. E.g. {'diff3d',2,'diff2d',2,'triplet',1}
- parameterList - a cell array in which the 1st, 3rd, 5th, etc. cell is a string defining a parameter followed by an N x 2 or N x 4 matrix where N is the number of the components in which the parameter is present. The first number defines if it is a constant parameter (1) or to be fitted (0), the second number is the value of the constant parameter or the initial value of the fitted parameter. The 3rd and 4th numbers define the lower and upper bounds of the parameter, respectively, if it is a fitted parameter. If the parameter is a fitted one, and only two numbers are provided, then the upper and lower bounds be filled in with default values.
Allowed parameters for
- 2D diffusion: 'tauD2d' OR 'D2d', 'w' (omega)
- 3D diffusion: 'tauD3d', 'kappa' OR 'D3d', 'kappa', 'w'
- for 2D or 3D diffusion: Fd2d and Fd3d (fractions; these two matrices must have N-1 rows altogether, where N is the number of diffusion components (since Σ Fd2d + Σ Fd2d = 1), e.g. if there are two 2D diffusion components and two 3D diffusion components, Fd2d has 2 rows, and Fd3d has 1 row.
- triplet: 'tauTriplet', 'tripletFraction'
'N' is always a parameter.
E.g. 'tauD2d',[0 50;0 20],'tauTriplet',[0 0.1 0 0.2]
to add two 2D diffusion correlation times (both of them to be fitted, with initial values of 50 and 20, without user-defined lower and upper bounds) and a triplet lifetime to be fitted with an initial value of 0.1, and user defined lower and upper bounds of 0 and 0.2, respectively.
- modelStructure: The model and the corresponding parameters can also be defined with a structure that can be generated by the GUI and saved to the Matlab workspace. The structure has two fields: modelType and parameterList
modelType = struct with fields:
diff2d: 2 (these numbers show the number of the corresponding component type)
diff3d: 2
triplet: 1
parameterList = struct with fields:
N: [0 0.5000 0 Inf] the 1st number indicates whether the parameter is constant (1) or to be fitted (0)
the 2nd number indicates the constant value or the initial value of the parameter
the 3rd and 4th numbers are the lower and upper bounds, respectively
tauD2d: [2×4 double]
kappa: [1 6 0 Inf]
tauD3d: [2×4 double]
tauTriplet: [0 0.1000 0.001 50]
tripletFraction: [0 0.5000 0 1]
Fd2d: [2x4 double]
Fd3d: [1x4 double] the total number of rows in Fd2d and Fd3d altogether is (number of diffusion components)-1
If the model is specified using the modelType and parameterList cell arrays, the parameterList is filled with default values based on the modelType. Therefore, this function call, in which the parameterList cell array is empty, is also OK:
fitFCS(t,ac,error,[1 205],{'diff2d',1','triplet',1},{},1,1);
If, on the other hand, the model is specified with the modelStructure, the parameterList field must perfectly correspond to the modelType. Otherwise, the program will report an error.
- typeOfFittingAlgorithm - 1: fminunc, 2: global fitting, 3 - simulated annealing
- calculateD - 0: tauD is calculated, ω is NOT required in the parameter list; 1: D is calculated, ω is required in the parameter list
E.g. fitFCS(t,ac(:,1),error(:,1),[5 205],{'diff3d',2,'diff2d',2,'triplet',1},{'tauD3d',[0 100;0 200],'tauD2d',[0 50;0 20],'tauTriplet',[0 0.1],'tripletFraction',[0 0.5],'N',[0 0.5]},2,0);
to fit a model consisting of two 3D diffusion components, two 2D diffusion components and one triplet component using the correlation times (and not the diffusion coefficients) using global fitting
The same result is achieved with the following command:
fitFCS(t,ac(:,1),error(:,1),[5 205],modelStruct,2,0);
Output parameters:
- fittedParams - the fitted parameters. Parameters for different datasets are stored in different rows
- parameterNames - the names of ALL parameters (fitted and constant ones)
- treatAsConstant - 0 for constant, 1 for fitted parameter.
- resultTable - results of different datasets are displayed below each other. The last column contains the dataset IDs.The structure of this table is the same as that produced in the GUI-mode (see above).
- figureStructures - structures containing the figures that were generated (same as in GUI-mode, see above). During the execution of the program, all plots are displayed in the same figure; therefore, only the last one will be visible (in command prompt mode). But all the figures correspoinding to every dataset are stored in the figureStructures output variable. In order to restore the figure corresponding to the 2nd dataset, type the following command at the Matlab command prompt:
struct2handle(figStruct(2),0);
figStruct is the variable containing the figure structures
- plotData - a cell array containing a 6-column table for each dataset (if there are more datasets) or a single 6-column table (if there is only one dataset).
- 1st column: time (original resolution, full range)
- 2nd colmn: measured autocorrelation function plotted for time points in column 1
- 3rd column: time (original resolution, fitted range)
- 4th column: residuals (original resolution, fitted range)
- 5th column:
time for the fitted autocorrelation function (full range)
- 6th column: fitted autocorrelation function