Assignment 4 part 3: Track microspheres over time

From Course Wiki
Jump to: navigation, search
20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg

This is Part 3 of Assignment 4.

Particle tracking

In this part of the assignment, you will write code to follow microscopic objects throughout a series of movie frames. Then, you'll record movies of small, fluorescent microspheres first diffusing in purely viscous solutions of glycerol-water, and (next week) moving in fibroblast cells after endocytosis. Calculating the mean squared displacement of their motion as a function of time interval will allow you to characterize their physical environment and behavior, first in terms of diffusivity and viscosity coefficients of the glycerol-water mixtures, next recognizing other material or transport properties in fibroblast cells.

Contextual background

Brownian motion

This section was adapted from

If you have ever looked at an aqueous sample through a microscope, you have probably noticed that every small particle you see wiggles about continuously. Robert Brown, a British botanist, was not the first person to observe these motions, but perhaps the first person to recognize the significance of this observation. Experiments quickly established the basic features of these movements. Among other things, the magnitude of the fluctuations depended on the size of the particle, and there was no difference between "live" objects, such as plant pollen, and things such as rock dust. Apparently, finely crushed pieces of an Egyptian mummy also displayed these fluctuations.

Brown noted: [The movements] arose neither from currents in the fluid, nor from its gradual evaporation, but belonged to the particle itself.

This effect may have remained a curiosity had it not been for A. Einstein and M. Smoluchowski. They realized that these particle movements made perfect sense in the context of the then developing kinetic theory of fluids. If matter is composed of atoms that collide frequently with other atoms, they reasoned, then even relatively large objects such as pollen grains would exhibit random movements. This last sentence contains the ingredients for several Nobel prizes!

Indeed, Einstein's interpretation of Brownian motion as the outcome of continuous bombardment by atoms immediately suggested a direct test of the atomic theory of matter. Perrin received the 1926 Nobel Prize for validating Einstein's predictions, thus confirming the atomic theory of matter.

Since then, the field has exploded, and a thorough understanding of Brownian motion is essential for everything from polymer physics to biophysics, aerodynamics, and statistical mechanics. One of the aims of this lab is to directly reproduce the experiments of J. Perrin that lead to his Nobel Prize. A translation of the key work is included in the reprints folder. Have a look – he used latex spheres, and we will use polystyrene spheres, but otherwise the experiments will be identical. In addition to reproducing Perrin's results, you will probe further by looking at the effect of varying solvent molecule size.

Diffusion coefficient of microspheres in suspension

According to theory,[1][2][3][4] the mean squared displacement of a suspended particle is proportional to the time interval as: $ \left \langle {\left | \vec r(t+\tau)-\vec r(t) \right \vert}^2 \right \rangle=2Dd\tau $, where r(t) = position, d = number of dimensions, D = diffusion coefficient, and $ \tau $= time interval.

Make a synthetic movie of diffusing particles

Generate a synthetic movie of diffusing microspheres using the function SimulatedBrownianMotionMovie below. (Hint: you may need to make some small modifications to the code to get it to work. Perhaps the code you wrote in Parts 1 and 2 will be helpful...)

function SimulatedMovie = SimulatedBrownianMotionMovie( varargin )
    % this is a fancy way to take care of input arguments
    % all of the arguments have default values
    % if you want to change a default value, call this function with name/value pairs like this:
    % foo = SimulatedBrownianMotionMovie( 'NumberOfParticles', 10, 'ParticleDiameter' 3.2E-6 );
    % it's safe to not understand this, if you like.
    p = inputParser();
    p.addParameter( 'NumberOfParticles', 5 );
    p.addParameter( 'ParticleDiameter', 1E-6 );
    p.addParameter( 'PixelSize', 6.9E-6 / 25 );
    p.addParameter( 'ImageSize', [ 300 300 ] );
    p.addParameter( 'FrameInterval', 1/10 );
    p.addParameter( 'TemperatureKelvin', 293 );
    p.addParameter( 'Viscosity', 1.0E-3 );
    p.addParameter( 'BoltzmannConstant', 1.38e-23 );
    p.addParameter( 'NumberOfFrames', 100 );
    p.addParameter( 'ShowImages', true );
    p.parse( varargin{:} );
    assignLocalVariable = @( name, value ) assignin( 'caller', name, value );
    allParameters = fields( p.Results );
    for ii = 1:numel(allParameters)
        assignLocalVariable( allParameters{ii}, p.Results.(allParameters{ii}) ); 
    % at this point, there will be local scope variables with the names of
    % each of the parameters specified above
    close all
    diffusionCoefficient = BoltzmannConstant * TemperatureKelvin / ( 3 * pi * Viscosity * ParticleDiameter );
    displacementStandardDeviation = sqrt( diffusionCoefficient * 2 * FrameInterval );
    % create initial particle table with random positions
    particleRadiusInPixels = ParticleDiameter / 2 / PixelSize;
    Position = bsxfun( @times, rand( NumberOfParticles, 2 ), ImageSize );
    Radius = particleRadiusInPixels * ones( NumberOfParticles, 1 );
    TotalPhotonNumber = sqrt(2) * particleRadiusInPixels^2 * 4095 * ones( NumberOfParticles, 1 );
    particleTable = table( Position, Radius, TotalPhotonNumber );
    SimulatedMovie = zeros( ImageSize(1), ImageSize(2), NumberOfFrames );
    for ii = 1:NumberOfFrames
        SimulatedMovie(:,:,ii) = SyntheticFluorescentMicrosphereImageWithNoise( particleTable, ImageSize ) / 4095;
        if( ShowImages )
            imshow( SimulatedMovie(:,:,ii) );
        % update the particle position according to the diffusion
        oldPosition = table2array( particleTable(:, 'Position') );
        Position = oldPosition + displacementStandardDeviation / PixelSize * randn( size( oldPosition ) );
        % coefficient
        particleTable(:, 'Position') = table( Position );

Develop particle tracking code

Use the code you developed in parts 1 and 2 of this assignment to write a function that finds the weighted centroids of the particles in each frame of a movie. Make sure to use the WeightedCentroid property instead of Centroid in the regionprops command (and be sure you understand the difference between the two!). Your function should return an N x 3 matrix with the following three columns (in order):

  1. the Y centroid coordinate,
  2. the X centroid coordinate, and
  3. the frame number.

Test out your function using your synthetic movie.

Once you have your N x 3 matrix, use the track function (an announcement with a link to this function was sent out to the class) to add a fourth column to the matrix containing a number that identifies each different particle. The N x 4 matrix will be the input to CalculateDiffusionCoefficientsFromTrajectories (below).


Make an x-y plot of two or more example bead trajectories from your synthetic movie. In other words, plot the y-centroid position vs. the x-centroid position for two or more particles. Remember, in MATLAB, you can isolate the coordinates of particle number 5, for example, by using a conditional statement:

particle5 = Nby4Matrix(:,4) == 5;

(the above statement is true when the fourth column of the matrix particle number is equal to 5). Then you can extract the x and y coordinates of particle 5 by indexing your matrix using this condition:

xyCoordinatesOfParticle5 = Nby4Matrix(particle5,1:2);

(Hint: If you subtract the initial position from each trajectory, then you can plot multiple trajectories on a single set of axes.)

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;
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];
        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 ));
        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;
% 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;
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;


Use CalculateDiffusionCoefficientsFromTrajectories to plot the MSD of each bead trajectory from your synthetic movie. Report the measured average diffusion coefficient and its uncertainty (with appropriate significant figures). Does this agree with the diffusion coefficient that you used to make the movie? (Hint: if it's orders of magnitude off, you may have a bug in your code.)

Estimating the diffusion coefficient by tracking suspended microspheres

Imaging chamber for fluorescent microspheres diffusing in water:glycerol mixtures

1. Track some 0.84μm Pink Spherotech polystyrene beads in water-glycerin mixtures (Samples A and B contain 80% and 85% glycerin, respectively).

Notes: Fluorescent microspheres have been mixed for you by the instructors into water-glycerin solutions A and B. (a) Vortex the stock Falcon tube, and then (b) transfer the bead suspension into its imaging chamber (consisting of a microscope slide, double-sided tape delimiting a 2-mm channel, and a 22mm x 40mm No. 1.5 coverslip, and sealed at both ends nail polish).
Tip 1: Ensure that the focal plane you choose to image is not near the coverslip or the slide. If some particles don't move or hardly move at all, it is likely that they are stuck to the coverslip. Adjust the focus so that you are viewing a plane near the middle of the sample. (A good way to do this is to focus on the top and bottom of the sample chamber and then split the difference.)
Tip 2: Each full frame of the full camera field of view takes up almost 2.5 MB of memory, so movies can get large very fast. Try to limit image data variables to a reasonable size by keeping the length of the movie short or limiting the Region of Interest (ROI) to a fraction of the full field of view.( A full-field, three minute move takes up about 4 GB, which is certain to push MATLAB over the edge.)

2. Record movies of beads diffusing in the two glycerol solutions and use your newly developed code to estimate the diffusion coefficient of each sample. Consider how many particles you should track and for how long. What factors determine the uncertainty and accuracy of your estimate?

  1. Include a snapshot of the 0.84 μm fluorescent beads monitored.
  2. Discuss any important considerations you made when preparing your samples and capturing your images. For example, how did you choose your exposure time, frame rate, number of particles in the region of interest, choice of sample plane, etc?
  3. Make an x-y plot of two or more example bead trajectories for each of the glycerin samples. (Hint: Subtract the initial position from each trajectory to plot multiple trajectories on a single set of axes.)
  4. Plot the average MSD vs τ results for the two glycerin samples (A and B); use log-log axes. Use the minimum number of axes that can convey your results clearly.
  5. Include a table of the diffusion coefficient, viscosity and measured glycerin/water ratio for each of the samples (A and B)
  6. Is the viscosity you measured close to the theoretical value predicted by this website?
  7. List several error sources present in your measurement of the MSD of diffusing particles. Classify each error as random or systematic, and technical or fundamental. Which sources do you think most significantly affect your measured diffusion coefficient? Note that we will quantify the magnitude of these errors in the next assignment.


Back to 20.309 Main Page


  1. A. Einstein, On the Motion of Small Particles Suspended in Liquids at Rest Required by the Molecular-Kinetic Theory of Heat, Annalen der Physik (1905).
  2. E. Frey and K. Kroy, Brownian motion: a paradigm of soft matter and biological physics, Ann. Phys. (2005). Published on the 100th anniversary of Einstein’s paper, this reference chronicles the history of Brownian motion from 1905 to the present.
  3. R. Newburgh, Einstein, Perrin, and the reality of atoms: 1905 revisited, Am. J. Phys. (2006). A modern replication of Perrin's experiment. Has a good, concise appendix with both the Einstein and Langevin derivations.
  4. M. Haw, Colloidal suspensions, Brownian motion, molecular reality: a short history, J. Phys. Condens. Matter (2002).