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: