Difference between revisions of "Calculating MSD and Diffusion Coefficients"

Example MSD calculation of sum and difference particles using MATLAB.

This page contains example MATLAB code for calculating the mean squared displacement (MSD) and diffusion coefficient from a movie of fluorescent microspheres.

Diffusion measurements from particle tracking data

The function CalculateDiffusionCoefficientFromTrajectories can be used to plot the MSD as a function of time for fluorescent particles and calculate their average diffusion coefficient. The input argument Centroids can be generated by the function track, which creates a 4 column matrix listing the x- and y- coordinates, the frame number and the particle number. Additional input arguments to CalculateDiffusionCoefficientFromTrajectories include the pixel size, the movie frame rate, and an optional figure handle for plotting the MSD.

Note that this is an updated version of CalculateDiffusionCoefficientFromTrajectories from the one distributed in PSet 2 (Fall 2016). In particular, the MSD algorithm has been improved for speed, a bug has been fixed to avoid unnecessarily excluding data points, and the code can now handle particles that appear in only one frame.

function out = CalculateDiffusionCoefficientFromTrajectories(Centroids, PixelSize, FrameRate, FigureHandle)
 
Centroids( :, 1:2 ) = Centroids( :, 1:2 ) * PixelSize;
 
numberOfParticles = max( Centroids(:,4) );
particles = cell( 1, numberOfParticles );
longestTrack = 0;
 
particlesToSkip = [];

if( nargin < 4 || isempty( FigureHandle ))
    FigureHandle = figure;
end
 
figure( FigureHandle )

for ii = 1:numberOfParticles
    % get the centroids for a single particle
    particle = Centroids(Centroids(:,4) == ii, :);
 
    numberOfSamples = size(particle,1); % number of samples in trajectory
    particles{ii} = struct();
     
    numberOfMsds = round(numberOfSamples-1);% / 8);
    particles{ii}.msd = zeros(1, numberOfMsds);
    particles{ii}.msdUncertainty = zeros(1, numberOfMsds);
    particles{ii}.number = numberOfMsds;
     
    if numberOfMsds==0
        particlesToSkip = [particlesToSkip ii];
    else
        
        for tau = 1:numberOfMsds
            dx = particle(1:1:(end-tau), 1) - particle((tau+1):1:end, 1);
            dy = particle(1:1:(end-tau), 2) - particle((tau+1):1:end, 2);
            drSquared = dx.^2+dy.^2;
            msd = mean( drSquared );
            
            particles{ii}.msd(tau) = msd;
            particles{ii}.msdUncertainty(tau) = std( drSquared ) / sqrt(length( drSquared ));
        end
        
        particles{ii}.diffusionCoefficient = particles{ii}.msd(1) / (2 * 2 * FrameRate^-1);
        
        longestTrack = max( longestTrack, numberOfSamples );
        loglog( (1:numberOfMsds)/FrameRate, particles{ii}.msd(1:numberOfMsds), 'color', rand(1,3) );
        hold on;
    end
end
 
% compute the average value of D and the ensemble average of msd
d = zeros(1,numberOfParticles);
averageMsd = zeros( 1, longestTrack );
averageMsdCount = zeros( 1, longestTrack);
 
for ii = 1:numberOfParticles
    if sum(particlesToSkip==ii)==0 
        numberOfMsds = length(particles{ii}.msd);
        d(ii) = particles{ii}.diffusionCoefficient;
        averageMsd(1:numberOfMsds) = averageMsd(1:numberOfMsds) + particles{ii}.msd;
        averageMsdCount(1:numberOfMsds) = averageMsdCount(1:numberOfMsds) + 1;
    end
end
 
figure( FigureHandle );
averageMsd = averageMsd ./ averageMsdCount;
loglog((1:longestTrack)/FrameRate, averageMsd, 'k', 'linewidth', 3);
ylabel('MSD (nm^2)');
xlabel('Time (s)');
title('MSD versus Time Interval');
%hold off;
 
out = struct();
out.diffusionCoefficient = mean(d);
out.diffusionCoefficientUncertainty = std(d)/sqrt(length(d));
out.particles = particles;