MATLAB: Estimating viscoelastic spectrum using Mason's method

From Course Wiki
Revision as of 01:24, 5 October 2013 by Steven Wasserman (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg

function [Gstar, GPrime, GDoublePrime, Tau ] = CalculateViscoelasticSpectrumMason(MsdValues, SamplingFrequency, ParticleRadius, Temperature, Smoothing)
%
    Tau = (1:length(MsdValues))' / SamplingFrequency;

    kB   = 1.38e-23;            % Boltzmann constant

    if(nargin < 5)
        Smoothing = 3;
    end

    if(nargin < 4)
        Temperature = 293;
    end

    smoothLogMsdValues =  log(smooth(MsdValues, Smoothing));
    smoothTimeInterval = Tau;

    % create a time interval axis to match diff(MSD)
    Tau = (smoothTimeInterval(1:(end-1)) + smoothTimeInterval(2:end)) / 2;
    shiftedMsd = (smoothLogMsdValues(1:(end-1)) + smoothLogMsdValues(2:end)) / 2;

    % compute G*, G' and G'' by Mason method
    alpha = diff(smoothLogMsdValues) ./ diff(log(smoothTimeInterval));
    gammaOfAlpha = 0.457 * (1 + alpha).^2 - 1.36 * (1 + alpha) +1.9; % approimation value from Mason
    %gammaOfAlpha = gamma(1 + alpha);

    Gstar = abs(2 * kB * Temperature ./ (3 * pi * ParticleRadius * shiftedMsd .* gammaOfAlpha));
    GPrime = Gstar .* cos(pi * alpha / 2);
    GDoublePrime = Gstar .* sin(pi * alpha / 2);

end

% this code tests the CalculateGStarGPrimeGDoublePrime function by
% simulating MSDs for two extreme cases: a particle on a spring and a
% freely diffusing particle

close all

numberOfFrames = 1800;
numberOfMsds = 1799;
sampleRate = 10;
particleRadius = 0.5e-6;
tau = (1:numberOfMsds) / sampleRate;

% generate synthetic MSD versus tau curve for a freely diffusing particle
diffusionCoefficient = 1e-15;
msdFree = diffusionCoefficient * tau;

% compute MSD versus time interval vector and G*, G', G''
[gStarFree, gPrimeFree, gDoublePrimeFree, outputTau] = CalculateGStarGPrimeGDoublePrime(msdFree, sampleRate, particleRadius);

% now simulate a confined microsphere (see Mason equation 4)
saturationMsd = sqrt(.01e-15);
diffusionTimeConstant = 1 / 3;
msdConfined = saturationMsd ^ 2 * (1 - exp(-tau / diffusionTimeConstant));

% compute MSD versus time interval vector and G*, G', G''
[gStarConfined, gPrimeConfined, gDoublePrimeConfined, outputTau] = CalculateGStarGPrimeGDoublePrime(msdConfined, sampleRate, particleRadius);

figure;
loglog( outputTau, gPrimeFree, outputTau, gDoublePrimeFree, outputTau, gPrimeConfined, outputTau, gDoublePrimeConfined );
legend( 'G'' Free', 'G'''' Free', 'G'' Confined', 'G'''' Confined');

% 
% k = sqrt(expectedDiffusionCoefficient / 2);
% 
% dx = k * randn(numberOfFrames,1);
% dy = k * randn(numberOfFrames,1);
% 
% x = cumsum(dx);
% y = cumsum(dy);
% 
% % format centroids to match track.m output -- 3rd column is frame number,
% % fourth column is particle number
% centroids = [x y (1:numberOfFrames)' ones(numberOfFrames,1)];
% 
% % compute MSD versus time interval vector and G*, G', G''
% results = CalculateDiffusionCoefficientFromTrajectories(centroids, 1, sampleRate);
% [GStarConfined, GPrimeConfined, GDoublePrimeConfined] = CalculateGStarGPrimeGDoublePrime(results.particles{1}.msd(1:numberOfMsds), sampleRate, particleRadius);