Particle Diffusion

From Bradwiki
Jump to navigation Jump to search

Monte Carlo Simulation of Basic Particle Diffusion

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

Matlab Script

  1 function [] = BasicParticleDiffusion()
  2 clc; close all; clear all;
  3 
  4 %=========================================%
  5 %	USER-ENTERED STARTING PARAMETERS
  6 %-----------------------------------------%
  7 
  8 Ndots = 200;		% number of particles in ES
  9 NSteps = 500;		% number of steps per trial
 10 NLoops = 10;		% number of loops to average over
 11 HiLoLoop = 5;		% number of diffusion rates to test between hi and low
 12 
 13 
 14 Sc = 1;			% scale of model
 15 dT = .1;		% time step (s)
 16 
 17 
 18 % Diffusion Rates [Choquet:		ES= .13		SynBSL= .05		SynLTP= .01]
 19 Des = (Sc*dT)* .1;			% D coef ES
 20 Dsb = (Sc*dT)* .05;			% D coef Syn Basal
 21 Dsp = (Sc*dT)* .01;			% D coef Syn LTP
 22 
 23 % Store a Diffusion Rate Copy
 24 D0 = Des;			% Diffusion Rate of Surf (D = L² / 2d*dT)
 25 D1 = Dsb;			% Diffusion Rate SYN1: BAS
 26 D2 = Dsp;			% Diffusion Rate SYN2: LTP
 27 
 28 DpDb = linspace(Dsp,Dsb,HiLoLoop);
 29 %=========================================%
 30 
 31 
 32 
 33 %=========================================%
 34 %			DO-DONT SWITCHES
 35 %-----------------------------------------%
 36 doESreg = 0;
 37 doAddDots = 1;
 38 doSubDots = 0;
 39 
 40 doTrace=0;
 41 doLiveScatter = 0;
 42 doLiveBar = 0;
 43 
 44 % Extra Stuff
 45 NdotsO = Ndots;
 46 %=========================================%
 47 
 48 
 49 
 50 
 51 %#############################################################
 52 for OLoop = 1:HiLoLoop
 53 %#############################################################
 54 
 55 
 56 dSb = DpDb(OLoop);
 57 dSp = Dsp;
 58 
 59 
 60 %=========================================%
 61 %	BASE DIFFUSION RATES EQUATIONS
 62 %-----------------------------------------%
 63 %{
 64 % Sc = Sc;					% scale of model (life:model)
 65 % dT = dT;					% time step (s)
 66 % dm = 2;                     % dimensions
 67 % Da = dSb*dT/Sc;			% Diffusion Rate A (D = L² / 2d*dT)
 68 % Db = dSp*dT/Sc;			% Diffusion Rate B
 69 % Dr = Da/Db;				% Ratio of Da:Des (1/Ls)^2;
 70 % Dn = Da/Dr;				% new D after scaling L
 71 % k = sqrt(dm*Da);			% stdev of D's step size distribution
 72 % L = sqrt(2*dm*Da);		% average diagonal (2D) step size
 73 % Lx = L/sqrt(2);           % average linear (1D) step size
 74 % Ls = 1/sqrt(Dr);			% scales Lx values for Dn
 75 % MSD = 2*dm*Da;			% mean squared displacement
 76 %}
 77 
 78 dm = 2;                     % dimensions
 79 k = sqrt(dm*Des);			% stdev of Des step size distribution
 80 
 81 Lsb1 = 1/sqrt(Des/dSb);		% scales Lx values for Dn
 82 Lsb2 = 1/sqrt(Des/dSp);		% scales Lx values for Dn
 83 
 84 
 85 XYL = zeros(2,Ndots);		% XY particle locations
 86 XYS = zeros(2,Ndots);		% XY particle step sizes
 87 %=========================================%
 88 
 89 
 90 
 91 
 92 
 93 
 94 %=========================================%
 95 %		SYNAPTIC LOCATION SETUP
 96 %-----------------------------------------%
 97 POLYSz = [3 6];		% size of ES (XY in µm)
 98 BOXwh = 0.8;		% size of ES (XY in µm)
 99 
100 % ES Size
101 POLYSz = POLYSz./Sc;		% scale enclosures 
102 XWIDE = POLYSz(1)/2;        % half X enclosure size (will double below)
103 YHIGH = POLYSz(2)/2;        % half X enclosure size (will double below)
104 BOARDER = [-POLYSz(1)/2 -POLYSz(2)/2 POLYSz(1) POLYSz(2)];
105 
106 
107 S1rad = .4;
108 S2rad = .4;
109 
110 S1C = [0; YHIGH/2];
111 S2C = [0; -YHIGH/2];
112 
113 S1BL = S1C - S1rad;
114 S2BL = S2C - S2rad;
115 
116 % Radius from synapse counters
117 S1G = repmat(S1C,1,numel(XYL(1,:)));
118 S2G = repmat(S2C,1,numel(XYL(1,:)));
119 
120 Box1 = sqrt(sum((XYL-S1G).^2)) < S1rad;
121 Box2 = sqrt(sum((XYL-S2G).^2)) < S2rad;
122 
123 
124 DenArea = POLYSz(1) * POLYSz(2);
125 S1Area = (pi*S1rad^2);
126 S2Area = (pi*S2rad^2);
127 SynAreas = S1Area + S2Area;
128 EssArea = DenArea - SynAreas;
129 %=========================================%
130 
131 
132 
133 
134 if doLiveScatter
135 %=========================================%
136 %               FIGURE SETUP
137 %-----------------------------------------%
138 Flh = figure(1);
139 set(Flh,'Units','pixels','OuterPosition',[1e2  1e2  7e2  7e2],'Color',[.95,.95,.95]);  
140 xlim = [-XWIDE*1.01 XWIDE*1.01]; ylim = [-YHIGH*1.01 YHIGH*1.01];
141 %---
142 	subplot('Position',[.05 .05 .65 .9])
143 Ph1 = scatter(XYL(1,:),XYL(2,:),5,[0 0 1]);
144 	set(Ph1,'Marker','o','SizeData',60,'LineWidth',.5,...
145 		'MarkerFaceColor',[.95 .1 .1],'MarkerEdgeColor','none')
146 	axis([xlim, ylim]);
147 	% axis off
148 	set(gca,'XTick',[],'YTick',[],'Color',[.8 .8 .8])
149 	hold on
150 
151 
152 rectangle('Position',BOARDER); hold on;
153 rectangle('Position',[S1BL(1),S1BL(2),S1rad*2,S1rad*2],'Curvature',[1,1]); hold on;
154 rectangle('Position',[S2BL(1),S2BL(2),S2rad*2,S2rad*2],'Curvature',[1,1]); hold on;
155 
156 
157 pause(.5)
158 
159 	Ah3 = subplot('Position',[.75 .05 .22 .9]);
160 Ph3 = bar([1 1]);
161 	set(Ah3,'Ylim',[0 100])
162 
163 pause(.5)
164 %=========================================%
165 end %if doLiveScatter
166 
167 
168 %=========================================%
169 %		MEMORY PREALLOCATION
170 %-----------------------------------------%
171 %SMxE = zeros(NLoops);
172 % preallocate for trace dot
173 XYLp = zeros(2,NSteps);	
174 % preallocate for data saves
175 DatN = NSteps/100;
176 
177 SMxE = zeros(fliplr(POLYSz .* 100));
178 %=========================================%
179 
180 
181   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182   for Ln = 1:NLoops
183   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184 
185 
186 
187 	%==================================================%
188 	for Nt = 1:NSteps 
189 	%==================================================%
190 
191 	% generates step sizes
192 		XYS = (k * randn(2,Ndots));	
193 
194 		
195 	% Sticky Functions	
196 		Box1 = sqrt(sum((XYL-S1G).^2)) < S1rad;
197 		Box2 = sqrt(sum((XYL-S2G).^2)) < S2rad;
198 		XYS(:,Box1) = XYS(:,Box1)*(Lsb1);
199 		XYS(:,Box2) = XYS(:,Box2)*(Lsb2);
200 		
201 
202 	% adds step to location
203 		XYL = XYL+XYS;	
204 
205 
206 	% rebound at enclosure walls
207 		[XYL] = ENCLOSE(Nt,XYL,XWIDE,YHIGH,Ndots); 
208 
209 
210 	% counts particles in each roi
211 		if mod(Nt,100)==0	
212 			SYN1n(Nt/100) = sum(Box1>0);
213 			SYN2n(Nt/100) = sum(Box2>0);
214 		end
215 
216 
217 
218 		%------------
219 	% PLOTS - Scatter, Bar, Trace, iSc
220 	
221 		if doLiveScatter
222 			set(Ph1,'XData',XYL(1,:),'YData',XYL(2,:));
223 			drawnow
224 		end
225 
226 		if doLiveBar
227 			if mod(Nt,100)==0
228 			set(Ph3,'YData',[SYN1n(Nt/100) SYN2n(Nt/100)]);
229 			drawnow
230 			end
231 		end
232 
233 		if doTrace
234 			XYLp(:,Nt) = XYL(:,1);		% save step of first dot (for trace)
235 			DUALPLOT(Nt,XYL,Ph1,XYLp,Ph2)
236 			XL(Nt,:) = XYL(1,:);
237 			YL(Nt,:) = XYL(2,:);
238 		end
239 
240 		
241 		%------------
242 		if mod(Nt,400)==0
243 		%if Nt >= 400; keyboard; end;
244 		%%
245 		
246 		pXYL = XYL;
247 		GR1x = XYL(1,:)+XWIDE;
248 		GR1y = XYL(2,:)+YHIGH;
249 		GR1c = round(GR1x*100);
250 		GR1r = round(GR1y*100);
251 		GR1c(GR1c < 1) = 1;
252 		GR1c(GR1c > (XWIDE*200)) = (XWIDE*200);
253 		GR1r(GR1r < 1) = 1;
254 		GR1r(GR1r > (YHIGH*200)) = (YHIGH*200);
255 		pXYL(1,:) = GR1c;
256 		pXYL(2,:) = GR1r;
257 
258 		G1xy = SMxE;
259 		for xy = 1:numel(GR1c)
260 		G1xy(GR1r(xy),GR1c(xy)) = 1;
261 		end
262 		
263 		%%
264 		end
265 		%------------
266 	
267 	
268 
269 	%if Nt >= 1; keyboard; end;
270 	%==================================================%
271 	end % for Nt = 1:Nsteps 
272 	%==================================================%
273 
274 	
275 	S1nPD = (SYN1n)./ S1Area;
276 	S2nPD = (SYN2n)./ S2Area;
277 
278   S1Pn{Ln} = S1nPD;
279   S2Pn{Ln} = S2nPD;
280   S1S2r{Ln} = ((SYN1n+1) ./ (SYN2n+1))';
281   EsPn{Ln} = ((Ndots - (SYN1n+SYN2n))')./ EssArea;
282 
283   MeanESSn(Ln) = mean(EsPn{Ln});
284   MeanSYNbN(Ln) = mean(S1Pn{Ln});
285   MeanSYNpN(Ln) = mean(S2Pn{Ln});
286   MeanRatio(Ln) = mean(S1S2r{Ln});
287 
288 
289 
290   if mod(Ln,10)==0;disp(Ln); end;
291   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292   end
293   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294 
295 
296 mESSn(OLoop) = mean(MeanESSn);
297 mSYNbN(OLoop) = mean(MeanSYNbN);
298 mSYNpN(OLoop) = mean(MeanSYNpN);
299 mRATbp(OLoop) = mean(MeanRatio);
300 
301 DboxD(:,OLoop) = [Des; dSb; dSp];
302 DboxR(OLoop) = dSp/dSb;
303 
304 
305 
306 %#############################################################
307 end % OLoop
308 %#############################################################
309 
310 
311 DPLOT(DboxD,DboxR,mRATbp,mESSn,mSYNbN,mSYNpN,NSteps,Des,dSb,dSp)
312 
313 FINALPLOT(EsPn,S1Pn,S2Pn,S1S2r,Des,dSb,dSp)
314 
315 %-------------------------------------------%
316 end %varargout={S1Pn,S2Pn};
317 %======================================================================%
318 %% MAIN FUNCTION END
319 %======================================================================%
320 
321 
322 
323 
324 
325 
326 
327 %======================================================================%
328 % ENCLOSE: keep particles inside polygon
329 %-------------------------------------------%
330 function [XYL] = ENCLOSE(Nt,XYL,XWIDE,YHIGH,Ndots)
331 
332 	for j = 1:Ndots 
333 		if XYL(1,j)>(XWIDE) || XYL(1,j)<(-XWIDE)
334 			XYL(1,j) = sign(XYL(1,j))*(XWIDE);
335 		elseif XYL(2,j)>(YHIGH) || XYL(2,j)<(-YHIGH)
336 			XYL(2,j) = sign(XYL(2,j))*(YHIGH);
337 		end
338 	end
339 
340 end
341 %======================================================================%
342 
343 
344 %======================================================================%
345 % DUALPLOT: plots diffusion
346 %-------------------------------------------%
347 function [] = DUALPLOT(Nt,XYL,Ph1,XYLp,SPh2)
348 %-------------------------------------------%
349 
350 LLg = 3;
351 if Nt > LLg
352 xp = (XYLp(1,(Nt-LLg):Nt));
353 yp = (XYLp(2,(Nt-LLg):Nt));
354 else
355 xp = (XYLp(1,Nt));
356 yp = (XYLp(2,Nt));
357 end
358 
359 set(Ph1,'XData',XYL(1,:),'YData',XYL(2,:));
360 drawnow
361 
362 plot(xp,yp,'Parent',SPh2);
363 drawnow
364 
365 end
366 %======================================================================%
367 
368 
369 
370 %======================================================================%
371 % FINALPLOT: plots steady-state averages
372 %-------------------------------------------%
373 function [] = DPLOT(DboxD,DboxR,mRATbp,mESSn,mSYNbN,mSYNpN,NSteps,Des,Db1,Db2)
374 
375 
376 %%
377 %------------------------------------------%
378 Fh = figure;
379 set(Fh,'OuterPosition',[200 200 1400 500])
380 % Color
381 neongreen = [.1 .9 .1];
382 liteblue = [.2 .9 .9];
383 hotpink=[.9 .1 .9];
384 %------------------------------------------%
385 axes('Position',[.07 .14 .4 .78]);
386 %------------------------------------------%
387 % [ph1] = line(1:numel(DboxR),DboxR);
388 [ph1] = plot(DboxR);
389 legend(ph1, 'D Ratio');
390 [LEGH,OBJH,OUTH,OUTM] = legend;
391 hold on
392 [ph2] = plot(mRATbp);
393 legend([OUTH;ph2],OUTM{:},'SS Ratio');
394 [LEGH,OBJH,OUTH,OUTM] = legend;
395 hold on
396 %------------------------------------------%
397 % xt = (get(gca,'XTick'));
398 % xt = linspace(0,NSteps,numel(xt)); %.* dT./60
399 % set(gca,'XTickLabel', sprintf('%.0f|',xt))
400 set(legend,'Location','NorthWest');
401 %------------------------------------------%
402 MS1 = 9; c1=neongreen; c2=hotpink; c3=liteblue;
403 set(ph1,'LineStyle','-','Color',c1,'LineWidth',2,...
404 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
405 set(ph2,'LineStyle',':','Color',c2,'LineWidth',2,...
406 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
407 hTitle  = title('Diffusion Coefficient Ratio vs Particle Steady-State Ratio');
408 hXLabel = xlabel('Trial');
409 hYLabel = ylabel('D Ratio or Particle-Density Ratio (+/- SEM)');
410 set(gca,'FontName','Helvetica');
411 set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
412 set([hXLabel, hYLabel],'FontSize',10);
413 set( hTitle,'FontSize',12);
414 set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
415 'XMinorTick','on','YMinorTick','on','YGrid','on', ...
416 'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
417 haxes=axis;
418 ylim([0 haxes(4)*1.2 ]);
419 % xlim([0 (haxes(2)*.9)]);
420 haxes=axis;
421 hold on
422 
423 estxt = horzcat(' ES: ', num2str(Des), '  ');
424 s1txt = horzcat(' S1: ', num2str(Db1), '  ');
425 s2txt = horzcat(' S2: ', num2str(Db2), '  ');
426 text(haxes(2)/2.5,haxes(4)/1.05,strcat(s1txt,' : ',s2txt, ' \bullet (um²/s; ES:',estxt,')' ),...
427 	'FontSize',10,'HorizontalAlignment','left');
428 
429 %------------------------------------------%
430 axes('Position',[.55 .14 .4 .78]);
431 %------------------------------------------%
432 % [ph1] = line(1:numel(DboxR),DboxR);
433 [ph1] = plot(mSYNbN);
434 legend(ph1, 'S1');
435 [LEGH,OBJH,OUTH,OUTM] = legend;
436 hold on
437 [ph2] = plot(mSYNpN);
438 legend([OUTH;ph2],OUTM{:},'S2');
439 [LEGH,OBJH,OUTH,OUTM] = legend;
440 hold on
441 [ph3] = plot(mESSn);
442 legend([OUTH;ph3],OUTM{:},'ES');
443 [LEGH,OBJH,OUTH,OUTM] = legend;
444 hold on
445 %------------------------------------------%
446 xt = (get(gca,'XTick')) .* 0;
447 xt(1:2:end) = DboxR;
448 
449 clear xtt
450 xtt{1,numel(xt)} = [];
451 xtt = num2cell(xt);
452 xtt(2:2:end) = {[]};
453 
454 set(gca,'XTickLabel', sprintf('%.2f|',xtt{1,:}))
455 set(legend,'Location','NorthEast');
456 %------------------------------------------%
457 MS1 = 9; c1=neongreen; c2=hotpink; c3=liteblue;
458 set(ph1,'LineStyle','-','Color',c1,'LineWidth',2,...
459 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
460 set(ph2,'LineStyle',':','Color',c2,'LineWidth',2,...
461 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
462 set(ph3,'LineStyle','-.','Color',c3,'LineWidth',2,...
463 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c3,'MarkerFaceColor',c3);
464 hTitle  = title('Diffusion Rate Effects on Particle Steady-State');
465 hXLabel = xlabel('S1:S2 Diffusion Coefficient Ratio');
466 hYLabel = ylabel('Particle Per Square Micron (+/- SEM)');
467 set(gca,'FontName','Helvetica');
468 set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
469 set([hXLabel, hYLabel],'FontSize',10);
470 set( hTitle,'FontSize',12);
471 set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
472 'XMinorTick','on','YMinorTick','on','YGrid','on', ...
473 'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
474 haxes=axis;
475 ylim([0 haxes(4)*1.2 ]);
476 % xlim([0 (haxes(2)*.9)]);
477 
478 %%
479 end
480 %======================================================================%
481 
482 
483 %======================================================================%
484 % FINALPLOT: plots steady-state averages
485 %-------------------------------------------%
486 function [] = FINALPLOT(EsPn,S1Pn,S2Pn,S1S2r,Des,dSb,dSp)
487 
488 %%
489 Fh = figure;
490 set(Fh,'OuterPosition',[200 200 650 700])
491 
492 % Color
493 neongreen = [.08 .88 .08];
494 liteblue = [.2 .9 .9];
495 hotpink=[.9 .1 .9];
496 %------------------------------------------%
497 [ph1 hax1 mu1] = CIplot(EsPn,liteblue,5);
498 legend(ph1, horzcat('ES _{', num2str(Des), ' um²/s}'));
499 [LEGH,OBJH,OUTH,OUTM] = legend;
500 hold on
501 ph11 = plot(mu1);
502 hold on
503 [ph2 hax2 mu2] = CIplot(S1Pn,neongreen,5);
504 legend([OUTH;ph2],OUTM{:},strcat('S1 _{', num2str(dSb), ' um²/s}'));
505 [LEGH,OBJH,OUTH,OUTM] = legend;
506 hold on
507 ph22 = plot(mu2);
508 hold on
509 pause(.5)
510 [ph3 hax3 mu3] = CIplot(S2Pn,hotpink,5);
511 legend([OUTH;ph3],OUTM{:},strcat('S2 _{', num2str(dSp), ' um²/s}'));
512 [LEGH,OBJH,OUTH,OUTM] = legend;
513 hold on
514 ph33 = plot(mu3);
515 hold on
516 %------------------------------------------%
517 hleg = legend;
518 set(hleg,'FontSize',14)
519 
520 xt = (get(gca,'XTick'))*100;
521 set(gca,'XTickLabel', sprintf('%.0f|',xt))
522 %------------------------------------------%
523 MS1 = 9; c1=liteblue; c2=neongreen; c3=hotpink;
524 % set(ph1,'LineStyle','-','Color',c1,'LineWidth',1,...
525 % 'Marker','o','MarkerSize',6,'MarkerEdgeColor',c1,'MarkerFaceColor','none');
526 set(ph11,'LineStyle','none','Color',c1,'LineWidth',2,...
527 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c1,'MarkerFaceColor',c1);
528 set(ph22,'LineStyle','none','Color',c2,'LineWidth',2,...
529 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c2,'MarkerFaceColor',c2);
530 set(ph33,'LineStyle','none','Color',c3,'LineWidth',2,...
531 'Marker','o','MarkerSize',MS1,'MarkerEdgeColor',c3,'MarkerFaceColor',c3);
532 hTitle  = title('Particle Steady-State Density in Dendritic ROIs');
533 hXLabel = xlabel('Time');
534 hYLabel = ylabel('Particle Per Square Micron (+/- SEM)');
535 set(gca,'FontName','Helvetica');
536 set([hTitle, hXLabel, hYLabel],'FontName','AvantGarde');
537 set([hXLabel, hYLabel],'FontSize',10);
538 set( hTitle,'FontSize',12);
539 set(gca,'Box','off','TickDir','out','TickLength',[.02 .02], ...
540 'XMinorTick','on','YMinorTick','on','YGrid','on', ...
541 'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'LineWidth',1);
542 haxes=axis;
543 ylim([0 haxes(4)*1.2 ]);
544 % xlim([0 (haxes(2)*.9)]);
545 %%
546 end
547 %======================================================================%