Difference between revisions of "Calculating MSD and Diffusion Coefficients"
From Course Wiki
Juliesutton (Talk | contribs) (Created page with "Category:20.309 Category:Matlab {{Template:20.309}} Example MSD calculation of sum and difference particles using MATLAB. T...") |
Juliesutton (Talk | contribs) m |
||
Line 12: | Line 12: | ||
− | <pre> | + | <pre>function out = CalculateDiffusionCoefficientFromTrajectories(Centroids, PixelSize, FrameRate,FigureHandle) |
− | function out = CalculateDiffusionCoefficientFromTrajectories(Centroids, PixelSize, FrameRate, FigureHandle) | + | |
Centroids( :, 1:2 ) = Centroids( :, 1:2 ) * PixelSize; | Centroids( :, 1:2 ) = Centroids( :, 1:2 ) * PixelSize; | ||
Line 21: | Line 20: | ||
longestTrack = 0; | longestTrack = 0; | ||
+ | particlesToSkip = []; | ||
+ | |||
if( nargin < 4 || isempty( FigureHandle )) | if( nargin < 4 || isempty( FigureHandle )) | ||
FigureHandle = figure; | FigureHandle = figure; | ||
Line 26: | Line 27: | ||
figure(FigureHandle) | figure(FigureHandle) | ||
+ | |||
for ii = 1:numberOfParticles | for ii = 1:numberOfParticles | ||
% get the centroids for a single particle | % get the centroids for a single particle | ||
Line 38: | Line 40: | ||
particles{ii}.number = numberOfMsds; | particles{ii}.number = numberOfMsds; | ||
− | for tau = 1: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 | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
end | end | ||
Line 61: | Line 68: | ||
for ii = 1:numberOfParticles | for ii = 1:numberOfParticles | ||
− | numberOfMsds = length(particles{ii}.msd); | + | 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 | end | ||
Line 78: | Line 87: | ||
out.diffusionCoefficient = mean(d); | out.diffusionCoefficient = mean(d); | ||
out.diffusionCoefficientUncertainty = std(d)/sqrt(length(d)); | out.diffusionCoefficientUncertainty = std(d)/sqrt(length(d)); | ||
− | out.particles = particles; | + | out.particles = particles;</pre> |
− | </pre> | + |
Revision as of 02:02, 17 October 2016
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;