Calculating MSD and Diffusion Coefficients
From Course Wiki
Revision as of 23:55, 16 October 2016 by Juliesutton (Talk | contribs)
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; 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; 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 % 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 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 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;