CODEWITHSUNDEEP
matlab
BioMSD
BioMSD

Reputation: 43

Calculating mean-squared displacement (msd) with MATLAB

This is part of a matlab intro course for biologists. I have my data points (for a single particle!) in a matrix with 4 columns (time, x, y, z) and a few thousand rows. What I want to do is to calculate the mean-squared displacement for the particle using the xyz coordinates for all time steps. MSD is defined as MSD=average(r(t)-r(0))^2 where r(t) is the position of the particle at time t and r(0) is the initial position, so in a sense it is the distance traveled by the particle over time interval t. This is what I have so far. I'm obviously a MATLAB newbie, and I'd really appreciate som help/input.

dx=zeros(length(data),1);   %Create space for displacements.
dy=zeros(length(data),1);
dz=zeros(length(data),1);

for i=1:length(pos) 
    dx(i)=data(i+1,2)-data(1,2);   %Calculate the distance at each time step
    dy(i)=data(i+1,3)-data(1,3);   %back to the origin.
    dz(i)=data(i+1,4)-data(1,4);
end

What I need to do next is calculate mean square of these values. But one thing that worries me is how much of the information in the data-matrix is simply ignored if the initial position of the particle stays fixed. I figured an algorithm along these lines would make a better use (more accurate, better statistics) of the data, but I'm not an experienced enough programmer to implement it...

note: dt = delta t = time step, Pos = position

Now compare every position with Pos 2 instead.

Now compare every position with Pos 3

Position 4...

This way there would be many more data points for each dt. Does it make sense to do like this to anyone else? Would this make the calculation less 'noisy' (ie make the plot of msd vs time smoother?) Perhaps there are some handy matlab functions to make this problem easier? I'd really appreciate some input.

Upvotes: 4

Views: 36299

Answers (2)

Jan Eglinger
Jan Eglinger

Reputation: 4090

Note that there's a MATLAB class available to do MSD analysis: MSD analyzer, including extensive documentation and tutorials.

Basic usage:

ma = msdanalyzer(2, SPACE_UNITS, TIME_UNITS);
ma = ma.addAll(tracks);
ma = ma.computeMSD;

Upvotes: 0

Jonas
Jonas

Reputation: 74940

If you're calculating the displacement only with respect to the first position, then you don't actually average anything, since MSD(1) is the average distance your object travels within one time step. So your reasoning is spot on.

However, I would calculate MSD the other way around: Get all displacements at dt=1 (1-2,2-3,3-4,...) and average. This is MSD(1). Then, you get all displacements at dt=2 (1-3,2-4,...) and average. This is MSD(2). And so on.

A useful property of Matlab is that you can vectorize some computations, i.e. do the calculation on an entire array, rather than doing them element-by-element. So if you have an array a of 100-by-1 coordinates, the difference between each coordinate and the next one is b=a(2:100)-a(1:99), or more generally, b=a(2:end)-a(1:end-1), so b(1) is a(2)-a(1), b(2) is a(3)-a(2) etc.

To calculate MSD from your array data (where I assume that time is in equal steps!), you'd write

nData = size(data,1); %# number of data points
numberOfdeltaT = floor(nData/4); %# for MSD, dt should be up to 1/4 of number of data points

msd = zeros(numberOfDeltaT,3); %# We'll store [mean, std, n]

%# calculate msd for all deltaT's

for dt = 1:numberOfDeltaT
   deltaCoords = data(1+dt:end,2:4) - data(1:end-dt,2:4);
   squaredDisplacement = sum(deltaCoords.^2,2); %# dx^2+dy^2+dz^2

   msd(dt,1) = mean(squaredDisplacement); %# average
   msd(dt,2) = std(squaredDisplacement); %# std
   msd(dt,3) = length(squaredDisplacement); %# n
end

Upvotes: 5

Related Questions