# MATLAB: Estimating viscoelastic spectrum using Mason's method

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
20.309: Biological Instrumentation and Measurement

```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;
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);
```