# TechNote/Metal induced energy transfer with wide-field microscopy

Here we describe the steps to reproduce the results of AppNote/Metal induced energy transfer.

## Online data

All the data and the complete MATLAB script are available for download:
photonscore-miet-tech-note.zip
(1,9 GB, sha1: `ba91f0e70624af037bd6692d3e8ca320e393ea53`

).

## Loading the data

There are three datafiles we have out of the measurements. Lets load them all
into memory with `read_photons`

function:

```
tt = photonscore.read_photons('tirf.photons');
wf = photonscore.read_photons('wf.photons');
rwf = photonscore.read_photons('irf_wf.photons');
% get time channel width in ps.
dt_channel = photonscore.file_info('tirf.photons').dt_channel;
```

For examle `tirf.photon`

dataset:

```
>> tt
tt =
struct with fields:
x: [79339363×1 uint16]
y: [79339363×1 uint16]
dt: [79339363×1 uint16]
```

Here we see that `tt`

is pretty large dataset that has about 80 millions
photons. Therefore we can take a brief look on intesity image with `1024`

pixels
binning:

```
img = photonscore.hist_2d(tt.x, 0, 4096, 1024, tt.y);
imagesc(img);
colormap gray
colorbar
```

## Take a look on timing

Here for `hist_2d`

call we use the binning range from `0`

to `4096`

that is the whole range of the recorded positions and we specify `1024`

as a
number of bins in this range.

As we can see on the figure there is enough photons to go with 1 “megapixel”
binning. Ok, lets sort the data with `flim.sort`

function fixing
the binning size:

```
fl_tt = photonscore.flim.sort(tt.x, 0, 4096, 1024, tt.y, tt.dt);
fl_wf = photonscore.flim.sort(wf.x, 0, 4096, 1024, wf.y, wf.dt);
fl_rwf = photonscore.flim.sort(rwf.x, 0, 4096, 1024, rwf.y, rwf.dt);
```

The sorting routine is a compute-intense operation. The duration may vary depending on CPU speed and number of cores in the system. Once the data sorting is done, one can compute mean and median lifetimes:

```
[md_tt, me_tt] = photonscore.flim.medimean(fl_tt);
[md_wf, me_wf] = photonscore.flim.medimean(fl_tt);
```

Here we use `flim.iw_tau`

to visualize intensity
weighted median (or mean `me_tt`

and `me_wf`

) lifetimes:

```
subplot(1,2,1)
img_tt = photonscore.flim.iw_tau(fl_tt.image, md_tt);
imagesc(img_tt);
pbaspect([1 1 1]);
subplot(1,2,2)
img_wf = photonscore.flim.iw_tau(fl_wf.image, md_wf);
imagesc(img_wf);
pbaspect([1 1 1]);
```

The resulting images with a default `preview.png`

palette applied.

## Fitting the decays

And now we can take a look on the decay:

```
r = [0 4095];
h_tt = photonscore.hist_1d(fl_tt.time, r(1), r(2), r(2)-r(1));
h_wf = photonscore.hist_1d(fl_wf.time, r(1), r(2), r(2)-r(1));
h_rwf = photonscore.hist_1d(fl_rwf.time, r(1), r(2), r(2)-r(1));
semilogy(h_tt);
hold on;
semilogy(h_wf);
semilogy(h_rwf);
hold off;
legend 'TIRF' 'Wide-field' 'IRF'
```

The result of the above code should look similar to figure below where we can
see the whole of our timing signal. Obviously one would not need everything from
the range acquired. So the first thing we redefine the range changing the first
line in the fragment above to `r = [550 3930];`

to focus only on the sub-range
of the time window.

The second observation one shall do is IRF background that should be removed by
subtracting `950`

from the histogram. Let us introduce variable `h_rwf_nbg`

that
holds background-free IRF and proceed with discrete fitting:

```
h_rwf_nbg = max(0, h_rwf - 950);
clear x0 xu xl;
xl.background = 50;
x0.background = 250;
xu.background = 550;
xl.tau_ref = 0.2;
x0.tau_ref = 2;
xu.tau_ref = 8;
xl.irf_shift = -20;
x0.irf_shift = 0;
xu.irf_shift = 20;
% 3-component fit
xl.tau = [10; 10; 10];
x0.tau = [200; 40; 280];
xu.tau = [800; 800; 800];
x0.a = [1 1 1];
[x1wf, m1wf] = photonscore.flim.fit_decay(h_rwf_nbg, h_wf, x0, xl, xu);
photonscore.flim.plot_fit_decay_model(m1wf, h_rwf_nbg, h_wf, dt_channel);
```

The resulting discrete fit look like this:

In general it is tricky to get an initial parameter guess `x0`

because of
mathematical instability of optimization routine. Therefore the above section
was evaluated multiple times with 1, 2 and 3 components.

The purpose of the discrete fit was to find `irf_shift`

and `tau_ref`

parameters to proceed with MEM analysis:

```
>> x1tt
x1tt =
struct with fields:
irf_shift: -0.1356
tau_ref: 2.4623
background: 237.9709
a: [3×1 double]
tau: [3×1 double]
likelihood: 9.3878e+03
```

In this particular dataset the shift parameter is close to zero (`-0.1356`

). But
the reference dye lifetime (`2.4623`

) will be very helpful for further
analysis:

```
% generate exponentially spaced tau
tau = photonscore.flim.log_tau_range(400, [10 800]);
% Create a model for admixture expectation maximization routine
m = photonscore.flim.convolve(...
h_rwf_nbg/sum(h_rwf_nbg), ... % with normalized to 1 IRF
x1tt.irf_shift,... % IRF shift we got with discrete fit
length(h_tt), ... % number of time channels
tau,...
x1tt.tau_ref); % use delta function reconvolution
% append background component
m = [m ones(size(m,1), 1)/size(m,1)];
[x2tt m2tt] = photonscore.flim.admixture_em(m, h_tt, ones(size(m,2), 1), ...
'iterations', 20000);
```

The resulting time spectra:

And the model:

Now we can run per pixel analysis:

```
%% Per pixel EM
% make a mask with more than 50 photons per pixel
mask = int32(fl_tt.image > 50);
% get a linear index of those pixels
fm = find(mask);
% mark active pixels with 1,2,3...
mask(fm) = 1:length(fm);
% extract decays per pixel
dd = photonscore.flim.decay_from_mask(fl_tt, mask, r(1), r(2), r(2)-r(1));
% extract IRF per pixel
rr = photonscore.flim.decay_from_mask(fl_rwf, mask, r(1), r(2), r(2)-r(1));
% remove the first column that correspons to all the pixels with mark 0
dd = dd(:, 2:end);
rr = rr(:, 2:end);
% plot the tau specta
plot(x2tt(1:end-1));
% put maximas here
px_tau = [21 52 111 164];
% form the model matrix
m = photonscore.flim.convolve(h_rwf_nbg/sum(h_rwf_nbg), x1tt.irf_shift,...
size(dd, 1), ...
px_tau,...
x1tt.tau_ref);
x5tt = ones(size(m,2), size(dd, 2));
tic;
x5tt = photonscore.flim.admixture_em(m, dd, x5tt, ...
'iterations', 50, 'method', 3);
toc;
```

The resulting pixel densities: