Difference between revisions of "Assignment 4 part 3: Track microspheres over time"

This is Part 3 of Assignment 4.

Overview

In this part of the assignment, you will develop code to track fluorescent microspheres in a sequence of frames. You'll check your code with a synthetic movie of diffusing beads, before recoding actual movies of fluorescent 0.84 μm microspheres diffusing in glycerol solutions. From these movies, you'll plot the mean squared displacement and measure the diffusion coefficient - and thus the viscosity of the glycerol.

Next week, you will use the same code to characterize the diffusive motion of particles in different rheological environments (including live cells).

Develop particle tracking code

Generate a synthetic movie of diffusing microspheres using the function SimulatedBrownianMotionMovie below (you may need to make some small modifications to the code to get it to work...).

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( 'NumerOfParticles', 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', 7.4E-6 / 40 );
    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}) ); 
    end
    % 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 );
    figure
    
    for ii = 1:NumberOfFrames
        ii
        SimulatedMovie(:,:,ii) = SyntheticFluorescentMicrosphereImageWithNoise( particleTable, ImageSize ) / 4095;
        if( ShowImages )
            imshow( SimulatedMovie(:,:,ii) );
            drawnow
        end
        
        % update the particle position according to the diffusion
        
        oldPosition = table2array( particleTable(:, 'Position') );
        Position = oldPosition + displacementStandardDeviation / PixelSize * randn( size( oldPosition ) );
        
        % coefficient
        particleTable(:, 'Position') = table( Position );
    end
end

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 unique identifier for each particle. The N x 4 matrix will be the input to CalculateDiffusionCoefficientsFromTrajectories (below).

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;

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?

Include a snapshot of the 0.84 μm fluorescent beads monitored. 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? 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.) 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. Include a table of the diffusion coefficient, viscosity and measured glycerin/water ratio for each of the samples (A and B) Is the viscosity you measured close to the theoretical value predicted by this website? 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.

To characterize the performance of your microscope for particle tracking, you will make a 180 s movie of fixed, 0.84 μm microspheres with your microscope and plot the mean squared displacement of pairs of particles in the movie.

Stability of microscope for particle tracking

Your stability plot should look something like this.

Once you have developed and tested your code, use the following procedure to measure the location precision of your microscope:

1. Bring a slide with fixed beads into focus.

• • Choose a field of view in which you can see at least 3 beads using the 40× objective.
  • Limit the field of view to only those beads by choosing a region of interest (ROI) in UsefulImageAcquisition tool. (Otherwise, you will have a gigantic movie that is hard to work with.)

1. Acquire a movie of beads for 3 minutes at a frame rate of at least 10 fps.
2. Use the code you developed to track the particles and plot MSD versus time interval.

The MSD for the difference trajectory of two fixed particles (why the difference trajectory?) should start out less than 100 nm² and still be less than 1000 nm² for t = 180 s. If the MSD is larger than that, refine your apparatus and methodology until you achieve that goal.

Turn in the code you developed and a plot of MSD for sum and difference tracks versus time interval.

Navigation

• Overview
• Part 1: Make a fake image
• Part 2: Measure resolution
• Part 3: Track microspheres over time

Back to 20.309 Main Page