Calculating MSD and Diffusion Coefficients

Example MSD calculation of sum and difference particles using MATLAB.

This page has example code for calculating the MSD and diffusion coefficient from a movie of fluorescent microspheres.

Diffusion measurements from particle tracking data

To calculate diffusion coefficients from the raw movie data, you must first run the FindCentroids and track functions (provided in PSet 2). Using the output trajectories of track, the function CalculateDiffusionCoefficientFromTrajectories will plot the mean squared displacement (MSD) as a function of time, and calculate the average diffusion coefficient for the particles. Additional input arguments include the pixel size (in the sample plane), the frame rate, and an optional figure handle for plotting the MSD.

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;