# Calculating MSD and Diffusion Coefficients

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
20.309: Biological Instrumentation and Measurement

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.

```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 (m^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;```