Brownian Motions

Brownian motion describes the stochastic diffusion of particles as they travel through n-dimensional spaces filled with other particles and physical barriers. Here the term particle is a generic term that can be generalized to describe the motion of molecule (e.g. H₂O) or proteins (e.g. NMDA receptors); note however that stochastic diffusion can also apply to things like the price index of a stock (see random walk) or the propagation heat energy across a surface. Brownian motion is among the simplest continuous-time stochastic processes, and a limit of various probabilistic processes (see random walk). As such, Brownian motion is highly generalizable to many applications, and is directly related to the universality of the normal distribution. In some sense, stochastic diffusion is a pure actuation of the basic statistical properties of probability distributions - it is distribution sampling translated into movements.

SPECIFICATION - Physical experiments suggest that Brownian motion has

• continuous increments
• increments of a particle over disjoint time intervals are independent events
• each increment is assumed to result from collisions with many molecules
• each increment is assumed to have a normal probability distribution
  • (Central Limit Theorem of probability theory: the sum of a large number of independent identically distributed random variables is approximately normal)
• the mean increment is zero as there is no preferred direction
• the position of a particle spreads out with time
• the variance of the increment is proportional to the length of time that Brownian motion has been observed

The probability density of a normally distributed random variable with mean μ and standard deviation σ is given by:

$f(x,\mu ,\sigma )={\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}({\frac {x-\mu }{\sigma }})^{2}}$

Mathematically, the random process called Brownian motion is denoted here as B(t) and defined for times t ≥ 0; the probability density of Brownian particles at the end of time period [0, t] is obtained by substituting μ = 0 and σ = √t, giving:

${\frac {1}{\sqrt {t2\pi }}}e^{-{\frac {1}{2}}({\frac {x}{\sqrt {t}}})^{2}}$

where x denotes the value of random variable B(t). The probability distribution of the increment B(t + u) − B(t) is:

$P[B(t+u)-B(t)\leq a]=\int _{x=-\infty }^{a}{\frac {1}{\sqrt {u2\pi }}}e^{-{\frac {1}{2}}({\frac {x}{\sqrt {u}}})^{2}}dx$

Common Diffusion Terms & Parameters

I've found it most helpful to learn about the terms and parameters used to characterize stochastic diffusion as they apply to diffusion phenomena along a flat 2D surface. These parameters include:

D  : diffusion rate coefficient

$D={L^{2} \over 2d\cdot t}$

L : step length

d : dimensions

t : time

D (in units: µm²/s) is the mean diffusion rate per unit time (velocity), often in µm²/s for biological motion on a molecular level (n.b. µm²/s are the units for particle diffusion on a surface/membrane; the surface can actually be curved or ruffled such that the surface fills a 3D space, however the diffusion equations for those instances are slightly more complex. Units for 3D space are naturally: µm³/s). This refers to how fast, on average, the particle moves along its surface trajectories. This value is often of ultimate interest, that is, the goal of single-particle tracking studies is often to define the average diffusion rate of a molecule and report it as a Diffusion Coefficient; once D is known, this Diffusion Coefficient can then be easily implemented in computer simulations of Brownian motion. However under most circumstances, the diffusion rate cannot be observed directly in empirical experiments - this would require the ability to visualize all the microscopic particles and collisions that dictate the particle's movement on a nanosecond timescale. In the animation on the right, the diffusion rate can actually be quantified directly; but what is often seen when observing particle diffusion through a microscope would more closely resemble this:

Instead, D is often calculated from the mean squared diffusion (MSD) path of the particle, defined below.

MSD  : mean squared displacement

$MSD={2d\cdot D}$

d : dimensions

D : diffusion coefficient

MSD (in units: µm²) is the mean square displacement of a particle over some time period. If D (the diffusion coefficient) is unknown, MSD can be computed using the following equation:

Error creating thumbnail: File missing

...where the quantity 'x ' refers to the displacement of a particle over 'n ' total steps, with 'k ' indicating the current step.

K  : standard deviation of the normal distribution for D

$k={\sqrt {d\cdot D}}$

d : dimensions

D : diffusion coefficient

k (in units: µm) is the standard deviation (σ) of the normal distribution of step lengths that, when randomly sampled, will give rise to a diffusion rate D. This value is useful for simulating Brownian motion for a particular diffusion rate.

L  : mean step length

$L={\sqrt {2d\cdot D}}$

d : dimensions

D : diffusion coefficient

L or Λ (in units: µm) is the average step length per interval of observation. In diffusion simulations, this is the step size per iteration (equivalent form: L = (2d * D).5).

λ  : mean 1D step length component of L

$\lambda ={L \over {\sqrt {2}}}$

λ (in units: µm) is the average 1-dimensional step length for each component (X,Y,Z) dimension of L. For example, simulating 2D particle diffusion will require the generation of individual step lengths for both the X and Y dimension. The total step distance from the origin will be the length of the hypotenuse created by the individual X and Y component step lengths. In fact, the equation: λ = L / √2 is derived from the Pythagorean theorem for right triangles, such that 2λ² = λ² where 2λ² represents a² + b² and λ² represents c².

ε  : step length scalar coefficient

$\epsilon ={\sqrt {\delta \over {D}}}$

δ : new desired diffusion rate

ε (in units: au) is a coefficient value that, when multiplied by each λ component step length, will scale those lengths to achieve a new diffusion rate δ.

Brownian Motion Matlab Code

function [out] = BradsBrownianMotionScript()
format compact;
format short;
close all;

% Note: This matlab script requires 'msdanalyzer' toolbox
%-------------############################------------------%
%                 STARTING PARAMETERS
%-------------############################------------------%
D = .3;                     % Diffusion Rate
Ds = .1;                    % Diffusion Scalar (Ds = Dn)
Dn = D/(D/Ds);              % new D after scaling L
d = 2;                      % dimensions
dT = 1;                     % time step
k = sqrt(d*D);	            % stdev of D's step size distribution
MSD = 2*d*D;                % mean squared displacement
L = sqrt(2*d*D);            % average diagonal (2D) step size
Lx = L/sqrt(2);             % average linear (1D) step size
Ls = 1/sqrt(D/Ds);          % scales Lx values for Dn


MSDtest = [1 0 0];			% test: D, Dn, or L
Scale = 1/10;				% scale of model
Ndots = 100;
Nsteps = Ndots;


xyl = ones(2,Ndots);
xyds = ones(2,Ndots);
lims = ((D+1)^2)*10;


%===========================================================%
%               LIVE PARTICLE DIFFUSION
%-----------------------------------------------------------%
for t = 1:Nsteps

    xyds = STEPxyds(Ndots, k);
	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
	MAINPLOT(xyl, lims);

end
%===========================================================%


%===========================================================%
%               MSD RANDOM STEPS ANALYSIS
%-----------------------------------------------------------%
tracks = cell(Ndots, 1);

stepN = 1;
for t = 1:Nsteps 

	xyds = STEPxyds(Ndots, k);
	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
    [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);

stepN = stepN+1;
end
MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);





%===========================================================%
%               MSD UNIFORM STEPS ANALYSIS
%-----------------------------------------------------------%
stepN = 1;
for t = 1:Nsteps 

xyds = stepsize(Ndots, Lx);
		
	[xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls);
    [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);

stepN = stepN+1;
end
MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);


boom = 'done';
%-------------############################------------------%
end %         ##    END MAIN FUNCTION   ##
%-------------############################------------------% 
%%





%%
%-------------############################------------------%
%             ##      SUBFUNCTIONS      ##
%-------------############################------------------% 


%-------------------------------------------%
% STEP SIZE GENERATOR
%-------------------------------------------%
function xyds = STEPxyds(Ndots, k)

    xyds = (k * randn(2,Ndots));

end


%-------------------------------------------%
% MOVE PARTICLES MAIN FUNCTION
%-------------------------------------------%
function [xyl] = AMPARSTEP(Ndots, xyds, xyl)
	
	for j = 1:Ndots
        xyl(:,j) = xyl(:,j)+xyds(:,j);
	end
	
end


%-------------------------------------------%
% LIVE DIFFUSION PLOT
%-------------------------------------------%
function [] = MAINPLOT(xyl, lims)
%-------------------------------------------%
xlim = [-lims lims];
ylim = [-lims lims];
zlim = [-5 5];

%=================================%
%       MAIN 2D PLOT
%---------------------------------%
figure(1)
subplot(2,1,1), 
AMPARPlot = gscatter(xyl(1,:),xyl(2,:));
axis([xlim, ylim]);
set(AMPARPlot,'marker','.','markersize',[6],'color',[1 0 0])


%=================================%
%           3D PLOT
%---------------------------------%
figure(1);
subplot(2,1,2), 
gscatter(xyl(1,:),xyl(2,:)); view(20, 30);
axis normal;
grid off
axis([xlim, ylim, zlim]);
set(gca, 'Box', 'on');

end


%-------------------------------------------%
% MANUAL STEP SIZE FUNCTION
%-------------------------------------------%
function xyds = stepsize(Ndots, Lx)
%-------------------------------------------%

   Lx(1:2,1:Ndots) = Lx;
   xyd = randi([0 1],Ndots,2)';
   xyd(xyd == 0) = -1;
   xyds = (Lx.*xyd);
   
end


%-------------------------------------------%
% MSD SCALED STEPS FUNCTION
%-------------------------------------------%
function [xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls)
	
	for j = 1:Ndots
        xyds(:,j) = xyds(:,j)*Ls;
        xyl(:,j) = xyl(:,j)+xyds(:,j);
	end	
	
end


%-------------------------------------------%
% MSD TRACKS GENERATOR
%-------------------------------------------%
function [tracks] = MSDfun(stepN, Nsteps, tracks, xyds)
    time = (0:Nsteps-1)';
    xymsd = xyds';
    xymsd = cumsum(xymsd,1);
    tracks{stepN} = [time xymsd];

end


%-------------------------------------------%
% MSD TRACKS ANALYSIS
%-------------------------------------------%
function [] = MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest)

SPACE_UNITS = 'µm';
TIME_UNITS = 's';
N_PARTICLES = Ndots;
N_TIME_STEPS = Nsteps;
N_DIM = 2;

oD = D;				% raw		µm^2/s
D  = D*Scale;       % to-scale	µm^2/s

oDn = Dn;			% raw		µm^2/s
Dn = Dn*Scale;		% to-scale	µm^2/s

oL = L;				% raw		µm
L = L*Scale;		% to-scale	µm

dTbase = dT;		% raw time-step 
dT = dT*Scale;		% to-scale time-step
k = k;				% stdv of step distribution

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

figure
ma.plotTracks;
ma.labelPlotTracks;

ma = ma.computeMSD;
ma.msd;

t = (0 : N_TIME_STEPS)' * dT;
[T1, T2] = meshgrid(t, t);
all_delays = unique( abs(T1 - T2) );

figure
ma.plotMSD;


cla
ma.plotMeanMSD(gca, true)

mmsd = ma.getMeanMSD;
t = mmsd(:,1);
x = mmsd(:,2);
dx = mmsd(:,3) ./ sqrt(mmsd(:,4));
errorbar(t, x, dx, 'k')

[fo, gof] = ma.fitMeanMSD;
plot(fo)
ma.labelPlotMSD;
legend off


ma = ma.fitMSD;

good_enough_fit = ma.lfit.r2fit > 0.8;
Dmean = mean( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
Dstd  =  std( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;

Dheader1 = ['Raw Unscaled Values'];
Dhead1 = ['    D        Dn        L'];
Ddat1 = [oD oDn oL];
disp(' ')
disp(Dheader1)
disp(Dhead1)
disp(Ddat1)


yourtesthead = ['YOU ARE TESTING DIFFUSION FOR:'];
if MSDtest(1)
	yourtest = ['   D:   original diffusion rate'];
elseif MSDtest(2)
	yourtest = ['   Dn:  new diffusion rate'];
elseif MSDtest(3)
	yourtest = ['   L:  step length'];
else
	yourtest = ['   generic diffusion rate'];
end
disp(yourtesthead)
disp(yourtest)

disp(' ')
fprintf('Estimation of raw D coefficient from MSD:\n')
fprintf('D = %.3g ± %.3g (mean ± std, N = %d)\n', ...
    Dmean, Dstd, sum(good_enough_fit));




% Retrieve instantaneous velocities, per track
 trackV = ma.getVelocities;

 % Pool track data together
 TV = vertcat( trackV{:} );

 % Velocities are returned in a N x (nDim+1) array: [ T Vx Vy ...]. So the
 % velocity vector in 2D is:
 V = TV(:, 2:3);

 % Compute diffusion coefficient
varV = var(V);
mVarV = mean(varV); % Take the mean of the two estimates
Dest = mVarV / 2 * dT;



Dheader2 = ['Scaling to model...'];
Dhead2 = ['    D        Dn        L'];
Ddat2 = [D Dn L];

disp(' ')
disp(Dheader2)
disp(Dhead2)
disp(Ddat2)
fprintf('Estimation from velocities histogram:\n')
fprintf('Tested D = %.3g %s, compare to scaled Des value of %.3g %s\n', ...
    Dest, [SPACE_UNITS '²/' TIME_UNITS], D, [SPACE_UNITS '²/' TIME_UNITS]);

% printf('D.psd target value was %.3g %s\n', ...
%     Dest, msdDpsd, [SPACE_UNITS '²/' TIME_UNITS]);

end