Particle Diffusion

Monte Carlo Simulation of Basic Particle Diffusion

A downloadable version of this code can be found here on Github gist.

Matlab Script

function [] = BasicParticleDiffusion()
clc; close all; clear all;

%=========================================%
%	USER-ENTERED STARTING PARAMETERS
%-----------------------------------------%

Ndots = 200;		% number of particles in ES
NSteps = 500;		% number of steps per trial
NLoops = 10;		% number of loops to average over
HiLoLoop = 5;		% number of diffusion rates to test between hi and low


Sc = 1;			% scale of model
dT = .1;		% time step (s)


% Diffusion Rates [Choquet:		ES= .13		SynBSL= .05		SynLTP= .01]
Des = (Sc*dT)* .1;			% D coef ES
Dsb = (Sc*dT)* .05;			% D coef Syn Basal
Dsp = (Sc*dT)* .01;			% D coef Syn LTP

% Store a Diffusion Rate Copy
D0 = Des;			% Diffusion Rate of Surf (D = L² / 2d*dT)
D1 = Dsb;			% Diffusion Rate SYN1: BAS
D2 = Dsp;			% Diffusion Rate SYN2: LTP

DpDb = linspace(Dsp,Dsb,HiLoLoop);
%=========================================%



%=========================================%
%			DO-DONT SWITCHES
%-----------------------------------------%
doESreg = 0;
doAddDots = 1;
doSubDots = 0;

doTrace=0;
doLiveScatter = 0;
doLiveBar = 0;

% Extra Stuff
NdotsO = Ndots;
%=========================================%




%#############################################################
for OLoop = 1:HiLoLoop
%#############################################################


dSb = DpDb(OLoop);
dSp = Dsp;


%=========================================%
%	BASE DIFFUSION RATES EQUATIONS
%-----------------------------------------%
%{
% Sc = Sc;					% scale of model (life:model)
% dT = dT;					% time step (s)
% dm = 2;                     % dimensions
% Da = dSb*dT/Sc;			% Diffusion Rate A (D = L² / 2d*dT)
% Db = dSp*dT/Sc;			% Diffusion Rate B
% Dr = Da/Db;				% Ratio of Da:Des (1/Ls)^2;
% Dn = Da/Dr;				% new D after scaling L
% k = sqrt(dm*Da);			% stdev of D's step size distribution
% L = sqrt(2*dm*Da);		% average diagonal (2D) step size
% Lx = L/sqrt(2);           % average linear (1D) step size
% Ls = 1/sqrt(Dr);			% scales Lx values for Dn
% MSD = 2*dm*Da;			% mean squared displacement
%}

dm = 2;                     % dimensions
k = sqrt(dm*Des);			% stdev of Des step size distribution

Lsb1 = 1/sqrt(Des/dSb);		% scales Lx values for Dn
Lsb2 = 1/sqrt(Des/dSp);		% scales Lx values for Dn


XYL = zeros(2,Ndots);		% XY particle locations
XYS = zeros(2,Ndots);		% XY particle step sizes
%=========================================%






%=========================================%
%		SYNAPTIC LOCATION SETUP
%-----------------------------------------%
POLYSz = [3 6];		% size of ES (XY in µm)
BOXwh = 0.8;		% size of ES (XY in µm)

% ES Size
POLYSz = POLYSz./Sc;		% scale enclosures 
XWIDE = POLYSz(1)/2;        % half X enclosure size (will double below)
YHIGH = POLYSz(2)/2;        % half X enclosure size (will double below)
BOARDER = [-POLYSz(1)/2 -POLYSz(2)/2 POLYSz(1) POLYSz(2)];


S1rad = .4;
S2rad = .4;

S1C = [0; YHIGH/2];
S2C = [0; -YHIGH/2];

S1BL = S1C - S1rad;
S2BL = S2C - S2rad;

% Radius from synapse counters
S1G = repmat(S1C,1,numel(XYL(1,:)));
S2G = repmat(S2C,1,numel(XYL(1,:)));

Box1 = sqrt(sum((XYL-S1G).^2)) < S1rad;
Box2 = sqrt(sum((XYL-S2G).^2)) < S2rad;


DenArea = POLYSz(1) * POLYSz(2);
S1Area = (pi*S1rad^2);
S2Area = (pi*S2rad^2);
SynAreas = S1Area + S2Area;
EssArea = DenArea - SynAreas;
%=========================================%




if doLiveScatter
%=========================================%
%               FIGURE SETUP
%-----------------------------------------%
Flh = figure(1);
set(Flh,'Units','pixels','OuterPosition',[1e2  1e2  7e2  7e2],'Color',[.95,.95,.95]);  
xlim = [-XWIDE*1.01 XWIDE*1.01]; ylim = [-YHIGH*1.01 YHIGH*1.01];
%---
	subplot('Position',[.05 .05 .65 .9])
Ph1 = scatter(XYL(1,:),XYL(2,:),5,[0 0 1]);
	set(Ph1,'Marker','o','SizeData',60,'LineWidth',.5,...
		'MarkerFaceColor',[.95 .1 .1],'MarkerEdgeColor','none')
	axis([xlim, ylim]);
	% axis off
	set(gca,'XTick',[],'YTick',[],'Color',[.8 .8 .8])
	hold on


rectangle('Position',BOARDER); hold on;
rectangle('Position',[S1BL(1),S1BL(2),S1rad*2,S1rad*2],'Curvature',[1,1]); hold on;
rectangle('Position',[S2BL(1),S2BL(2),S2rad*2,S2rad*2],'Curvature',[1,1]); hold on;


pause(.5)

	Ah3 = subplot('Position',[.75 .05 .22 .9]);
Ph3 = bar([1 1]);
	set(Ah3,'Ylim',[0 100])

pause(.5)
%=========================================%
end %if doLiveScatter


%=========================================%
%		MEMORY PREALLOCATION
%-----------------------------------------%
%SMxE = zeros(NLoops);
% preallocate for trace dot
XYLp = zeros(2,NSteps);	
% preallocate for data saves
DatN = NSteps/100;

SMxE = zeros(fliplr(POLYSz .* 100));
%=========================================%


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for Ln = 1:NLoops
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



	%==================================================%
	for Nt = 1:NSteps 
	%==================================================%

	% generates step sizes
		XYS = (k * randn(2,Ndots));	

		
	% Sticky Functions	
		Box1 = sqrt(sum((XYL-S1G).^2)) < S1rad;
		Box2 = sqrt(sum((XYL-S2G).^2)) < S2rad;
		XYS(:,Box1) = XYS(:,Box1)*(Lsb1);
		XYS(:,Box2) = XYS(:,Box2)*(Lsb2);
		

	% adds step to location
		XYL = XYL+XYS;	


	% rebound at enclosure walls
		[XYL] = ENCLOSE(Nt,XYL,XWIDE,YHIGH,Ndots); 


	% counts particles in each roi
		if mod(Nt,100)==0	
			SYN1n(Nt/100) = sum(Box1>0);
			SYN2n(Nt/100) = sum(Box2>0);
		end



		%------------
	% PLOTS - Scatter, Bar, Trace, iSc
	
		if doLiveScatter
			set(Ph1,'XData',XYL(1,:),'YData',XYL(2,:));
			drawnow
		end

		if doLiveBar
			if mod(Nt,100)==0
			set(Ph3,'YData',[SYN1n(Nt/100) SYN2n(Nt/100)]);
			drawnow
			end
		end

		if doTrace
			XYLp(:,Nt) = XYL(:,1);		% save step of first dot (for trace)
			DUALPLOT(Nt,XYL,Ph1,XYLp,Ph2)
			XL(Nt,:) = XYL(1,:);
			YL(Nt,:) = XYL(2,:);
		end

		
		%------------
		if mod(Nt,400)==0
		%if Nt >= 400; keyboard; end;
		%%
		
		pXYL = XYL;
		GR1x = XYL(1,:)+XWIDE;
		GR1y = XYL(2,:)+YHIGH;
		GR1c = round(GR1x*100);
		GR1r = round(GR1y*100);
		GR1c(GR1c < 1) = 1;
		GR1c(GR1c > (XWIDE*200)) = (XWIDE*200);
		GR1r(GR1r < 1) = 1;
		GR1r(GR1r > (YHIGH*200)) = (YHIGH*200);
		pXYL(1,:) = GR1c;
		pXYL(2,:) = GR1r;

		G1xy = SMxE;
		for xy = 1:numel(GR1c)
		G1xy(GR1r(xy),GR1c(xy)) = 1;
		end
		
		%%
		end
		%------------
	
	

	%if Nt >= 1; keyboard; end;
	%==================================================%
	end % for Nt = 1:Nsteps 
	%==================================================%

	
	S1nPD = (SYN1n)./ S1Area;
	S2nPD = (SYN2n)./ S2Area;

  S1Pn{Ln} = S1nPD;
  S2Pn{Ln} = S2nPD;
  S1S2r{Ln} = ((SYN1n+1) ./ (SYN2n+1))';
  EsPn{Ln} = ((Ndots - (SYN1n+SYN2n))')./ EssArea;

  MeanESSn(Ln) = mean(EsPn{Ln});
  MeanSYNbN(Ln) = mean(S1Pn{Ln});
  MeanSYNpN(Ln) = mean(S2Pn{Ln});
  MeanRatio(Ln) = mean(S1S2r{Ln});



  if mod(Ln,10)==0;disp(Ln); end;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


mESSn(OLoop) = mean(MeanESSn);
mSYNbN(OLoop) = mean(MeanSYNbN);
mSYNpN(OLoop) = mean(MeanSYNpN);
mRATbp(OLoop) = mean(MeanRatio);

DboxD(:,OLoop) = [Des; dSb; dSp];
DboxR(OLoop) = dSp/dSb;



%#############################################################
end % OLoop
%#############################################################


DPLOT(DboxD,DboxR,mRATbp,mESSn,mSYNbN,mSYNpN,NSteps,Des,dSb,dSp)

FINALPLOT(EsPn,S1Pn,S2Pn,S1S2r,Des,dSb,dSp)

%-------------------------------------------%
end %varargout={S1Pn,S2Pn};
%======================================================================%
%% MAIN FUNCTION END
%======================================================================%







%======================================================================%
% ENCLOSE: keep particles inside polygon
%-------------------------------------------%
function [XYL] = ENCLOSE(Nt,XYL,XWIDE,YHIGH,Ndots)

	for j = 1:Ndots 
		if XYL(1,j)>(XWIDE) || XYL(1,j)<(-XWIDE)
			XYL(1,j) = sign(XYL(1,j))*(XWIDE);
		elseif XYL(2,j)>(YHIGH) || XYL(2,j)<(-YHIGH)
			XYL(2,j) = sign(XYL(2,j))*(YHIGH);
		end
	end

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


%======================================================================%
% DUALPLOT: plots diffusion
%-------------------------------------------%
function [] = DUALPLOT(Nt,XYL,Ph1,XYLp,SPh2)
%-------------------------------------------%

LLg = 3;
if Nt > LLg
xp = (XYLp(1,(Nt-LLg):Nt));
yp = (XYLp(2,(Nt-LLg):Nt));
else
xp = (XYLp(1,Nt));
yp = (XYLp(2,Nt));
end

set(Ph1,'XData',XYL(1,:),'YData',XYL(2,:));
drawnow

plot(xp,yp,'Parent',SPh2);
drawnow

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



%======================================================================%
% FINALPLOT: plots steady-state averages
%-------------------------------------------%
function [] = DPLOT(DboxD,DboxR,mRATbp,mESSn,mSYNbN,mSYNpN,NSteps,Des,Db1,Db2)


%%
%------------------------------------------%
Fh = figure;
set(Fh,'OuterPosition',[200 200 1400 500])
% Color
neongreen = [.1 .9 .1];
liteblue = [.2 .9 .9];
hotpink=[.9 .1 .9];
%------------------------------------------%
axes('Position',[.07 .14 .4 .78]);
%------------------------------------------%
% [ph1] = line(1:numel(DboxR),DboxR);
[ph1] = plot(DboxR);
legend(ph1, 'D Ratio');
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
[ph2] = plot(mRATbp);
legend([OUTH;ph2],OUTM{:},'SS Ratio');
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
%------------------------------------------%
% xt = (get(gca,'XTick'));
% xt = linspace(0,NSteps,numel(xt)); %.* dT./60
% set(gca,'XTickLabel', sprintf('%.0f|',xt))
set(legend,'Location','NorthWest');
%------------------------------------------%
MS1 = 9; c1=neongreen; c2=hotpink; c3=liteblue;
set(ph1,'LineStyle','-','Color',c1,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
set(ph2,'LineStyle',':','Color',c2,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
hTitle  = title('Diffusion Coefficient Ratio vs Particle Steady-State Ratio');
hXLabel = xlabel('Trial');
hYLabel = ylabel('D Ratio or Particle-Density Ratio (+/- SEM)');
set(gca,'FontName','Helvetica');
set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
set([hXLabel, hYLabel],'FontSize',10);
set( hTitle,'FontSize',12);
set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
'XMinorTick','on','YMinorTick','on','YGrid','on', ...
'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
haxes=axis;
ylim([0 haxes(4)*1.2 ]);
% xlim([0 (haxes(2)*.9)]);
haxes=axis;
hold on

estxt = horzcat(' ES: ', num2str(Des), '  ');
s1txt = horzcat(' S1: ', num2str(Db1), '  ');
s2txt = horzcat(' S2: ', num2str(Db2), '  ');
text(haxes(2)/2.5,haxes(4)/1.05,strcat(s1txt,' : ',s2txt, ' \bullet (um²/s; ES:',estxt,')' ),...
	'FontSize',10,'HorizontalAlignment','left');

%------------------------------------------%
axes('Position',[.55 .14 .4 .78]);
%------------------------------------------%
% [ph1] = line(1:numel(DboxR),DboxR);
[ph1] = plot(mSYNbN);
legend(ph1, 'S1');
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
[ph2] = plot(mSYNpN);
legend([OUTH;ph2],OUTM{:},'S2');
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
[ph3] = plot(mESSn);
legend([OUTH;ph3],OUTM{:},'ES');
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
%------------------------------------------%
xt = (get(gca,'XTick')) .* 0;
xt(1:2:end) = DboxR;

clear xtt
xtt{1,numel(xt)} = [];
xtt = num2cell(xt);
xtt(2:2:end) = {[]};

set(gca,'XTickLabel', sprintf('%.2f|',xtt{1,:}))
set(legend,'Location','NorthEast');
%------------------------------------------%
MS1 = 9; c1=neongreen; c2=hotpink; c3=liteblue;
set(ph1,'LineStyle','-','Color',c1,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
set(ph2,'LineStyle',':','Color',c2,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
set(ph3,'LineStyle','-.','Color',c3,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c3,'MarkerFaceColor',c3);
hTitle  = title('Diffusion Rate Effects on Particle Steady-State');
hXLabel = xlabel('S1:S2 Diffusion Coefficient Ratio');
hYLabel = ylabel('Particle Per Square Micron (+/- SEM)');
set(gca,'FontName','Helvetica');
set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
set([hXLabel, hYLabel],'FontSize',10);
set( hTitle,'FontSize',12);
set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
'XMinorTick','on','YMinorTick','on','YGrid','on', ...
'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
haxes=axis;
ylim([0 haxes(4)*1.2 ]);
% xlim([0 (haxes(2)*.9)]);

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


%======================================================================%
% FINALPLOT: plots steady-state averages
%-------------------------------------------%
function [] = FINALPLOT(EsPn,S1Pn,S2Pn,S1S2r,Des,dSb,dSp)

%%
Fh = figure;
set(Fh,'OuterPosition',[200 200 650 700])

% Color
neongreen = [.08 .88 .08];
liteblue = [.2 .9 .9];
hotpink=[.9 .1 .9];
%------------------------------------------%
[ph1 hax1 mu1] = CIplot(EsPn,liteblue,5);
legend(ph1, horzcat('ES _{', num2str(Des), ' um²/s}'));
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
ph11 = plot(mu1);
hold on
[ph2 hax2 mu2] = CIplot(S1Pn,neongreen,5);
legend([OUTH;ph2],OUTM{:},strcat('S1 _{', num2str(dSb), ' um²/s}'));
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
ph22 = plot(mu2);
hold on
pause(.5)
[ph3 hax3 mu3] = CIplot(S2Pn,hotpink,5);
legend([OUTH;ph3],OUTM{:},strcat('S2 _{', num2str(dSp), ' um²/s}'));
[LEGH,OBJH,OUTH,OUTM] = legend;
hold on
ph33 = plot(mu3);
hold on
%------------------------------------------%
hleg = legend;
set(hleg,'FontSize',14)

xt = (get(gca,'XTick'))*100;
set(gca,'XTickLabel', sprintf('%.0f|',xt))
%------------------------------------------%
MS1 = 9; c1=liteblue; c2=neongreen; c3=hotpink;
% set(ph1,'LineStyle','-','Color',c1,'LineWidth',1,...
% 'Marker','o','MarkerSize',6,'MarkerEdgeColor',c1,'MarkerFaceColor','none');
set(ph11,'LineStyle','none','Color',c1,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
set(ph22,'LineStyle','none','Color',c2,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
set(ph33,'LineStyle','none','Color',c3,'LineWidth',2,...
'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c3,'MarkerFaceColor',c3);
hTitle  = title('Particle Steady-State Density in Dendritic ROIs');
hXLabel = xlabel('Time');
hYLabel = ylabel('Particle Per Square Micron (+/- SEM)');
set(gca,'FontName','Helvetica');
set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
set([hXLabel, hYLabel],'FontSize',10);
set( hTitle,'FontSize',12);
set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
'XMinorTick','on','YMinorTick','on','YGrid','on', ...
'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
haxes=axis;
ylim([0 haxes(4)*1.2 ]);
% xlim([0 (haxes(2)*.9)]);
%%
end
%======================================================================%