Qual

From Bradwiki
Jump to navigation Jump to search

Qual Paper SI Material

Qual Project Source Code

Github2.png
All source code used in my qualifying exam project can be found here in my github repo for the project. An online version of my qual with enhanced features like citation data popup graphics and animations can be found on this page


SI Fig. 1 - Actin Filament Surface Hull


SI Fig. 2 - Brownian Motion Matlab Code

Brownian Motion Matlab Code



  1 function [boom] = BradsBrownianMotionScript()
  2 format compact;
  3 format short;
  4 close all;
  5 
  6 % Note: This matlab script requires 'msdanalyzer' toolbox
  7 %-------------############################------------------%
  8 %                 STARTING PARAMETERS
  9 %-------------############################------------------%
 10 D = .3;                     % Diffusion Rate
 11 Ds = .1;                    % Diffusion Scalar (Ds = Dn)
 12 Dn = D/(D/Ds);              % new D after scaling L
 13 d = 2;                      % dimensions
 14 dT = 1;                     % time step
 15 k = sqrt(d*D);	            % stdev of D's step size distribution
 16 MSD = 2*d*D;                % mean squared displacement
 17 L = sqrt(2*d*D);            % average diagonal (2D) step size
 18 Lx = L/sqrt(2);             % average linear (1D) step size
 19 Ls = 1/sqrt(D/Ds);          % scales Lx values for Dn
 20 
 21 
 22 MSDtest = [1 0 0];			% test: D, Dn, or L
 23 Scale = 1/10;				% scale of model
 24 Ndots = 100;
 25 Nsteps = Ndots;
 26 
 27 
 28 xyl = ones(2,Ndots);
 29 xyds = ones(2,Ndots);
 30 lims = ((D+1)^2)*10;
 31 
 32 
 33 %===========================================================%
 34 %               LIVE PARTICLE DIFFUSION
 35 %-----------------------------------------------------------%
 36 for t = 1:Nsteps
 37 
 38     xyds = STEPxyds(Ndots, k);
 39 	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
 40 	MAINPLOT(xyl, lims);
 41 
 42 end
 43 %===========================================================%
 44 
 45 
 46 %===========================================================%
 47 %               MSD RANDOM STEPS ANALYSIS
 48 %-----------------------------------------------------------%
 49 tracks = cell(Ndots, 1);
 50 
 51 stepN = 1;
 52 for t = 1:Nsteps 
 53 
 54 	xyds = STEPxyds(Ndots, k);
 55 	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
 56     [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);
 57 
 58 stepN = stepN+1;
 59 end
 60 MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);
 61 
 62 
 63 
 64 
 65 
 66 %===========================================================%
 67 %               MSD UNIFORM STEPS ANALYSIS
 68 %-----------------------------------------------------------%
 69 stepN = 1;
 70 for t = 1:Nsteps 
 71 
 72 xyds = stepsize(Ndots, Lx);
 73 		
 74 	[xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls);
 75     [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);
 76 
 77 stepN = stepN+1;
 78 end
 79 MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);
 80 
 81 
 82 boom = 'done';
 83 %-------------############################------------------%
 84 end %         ##    END MAIN FUNCTION   ##
 85 %-------------############################------------------% 
 86 %%
 87 
 88 
 89 
 90 
 91 
 92 %%
 93 %-------------############################------------------%
 94 %             ##      SUBFUNCTIONS      ##
 95 %-------------############################------------------% 
 96 
 97 
 98 %-------------------------------------------%
 99 % STEP SIZE GENERATOR
100 %-------------------------------------------%
101 function xyds = STEPxyds(Ndots, k)
102 
103     xyds = (k * randn(2,Ndots));
104 
105 end
106 
107 
108 %-------------------------------------------%
109 % MOVE PARTICLES MAIN FUNCTION
110 %-------------------------------------------%
111 function [xyl] = AMPARSTEP(Ndots, xyds, xyl)
112 	
113 	for j = 1:Ndots
114         xyl(:,j) = xyl(:,j)+xyds(:,j);
115 	end
116 	
117 end
118 
119 
120 %-------------------------------------------%
121 % LIVE DIFFUSION PLOT
122 %-------------------------------------------%
123 function [] = MAINPLOT(xyl, lims)
124 %-------------------------------------------%
125 xlim = [-lims lims];
126 ylim = [-lims lims];
127 zlim = [-5 5];
128 
129 %=================================%
130 %       MAIN 2D PLOT
131 %---------------------------------%
132 figure(1)
133 subplot(2,1,1), 
134 AMPARPlot = gscatter(xyl(1,:),xyl(2,:));
135 axis([xlim, ylim]);
136 set(AMPARPlot,'marker','.','markersize',[6],'color',[1 0 0])
137 
138 
139 %=================================%
140 %           3D PLOT
141 %---------------------------------%
142 figure(1);
143 subplot(2,1,2), 
144 gscatter(xyl(1,:),xyl(2,:)); view(20, 30);
145 axis normal;
146 grid off
147 axis([xlim, ylim, zlim]);
148 set(gca, 'Box', 'on');
149 
150 end
151 
152 
153 %-------------------------------------------%
154 % MANUAL STEP SIZE FUNCTION
155 %-------------------------------------------%
156 function xyds = stepsize(Ndots, Lx)
157 %-------------------------------------------%
158 
159    Lx(1:2,1:Ndots) = Lx;
160    xyd = randi([0 1],Ndots,2)';
161    xyd(xyd == 0) = -1;
162    xyds = (Lx.*xyd);
163    
164 end
165 
166 
167 %-------------------------------------------%
168 % MSD SCALED STEPS FUNCTION
169 %-------------------------------------------%
170 function [xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls)
171 	
172 	for j = 1:Ndots
173         xyds(:,j) = xyds(:,j)*Ls;
174         xyl(:,j) = xyl(:,j)+xyds(:,j);
175 	end	
176 	
177 end
178 
179 
180 %-------------------------------------------%
181 % MSD TRACKS GENERATOR
182 %-------------------------------------------%
183 function [tracks] = MSDfun(stepN, Nsteps, tracks, xyds)
184     time = (0:Nsteps-1)';
185     xymsd = xyds';
186     xymsd = cumsum(xymsd,1);
187     tracks{stepN} = [time xymsd];
188 
189 end
190 
191 
192 %-------------------------------------------%
193 % MSD TRACKS ANALYSIS
194 %-------------------------------------------%
195 function [] = MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest)
196 
197 SPACE_UNITS = 'µm';
198 TIME_UNITS = 's';
199 N_PARTICLES = Ndots;
200 N_TIME_STEPS = Nsteps;
201 N_DIM = 2;
202 
203 oD = D;				% raw		µm^2/s
204 D  = D*Scale;       % to-scale	µm^2/s
205 
206 oDn = Dn;			% raw		µm^2/s
207 Dn = Dn*Scale;		% to-scale	µm^2/s
208 
209 oL = L;				% raw		µm
210 L = L*Scale;		% to-scale	µm
211 
212 dTbase = dT;		% raw time-step 
213 dT = dT*Scale;		% to-scale time-step
214 k = k;				% stdv of step distribution
215 
216 ma = msdanalyzer(2, SPACE_UNITS, TIME_UNITS);
217 ma = ma.addAll(tracks);
218 disp(ma)
219 
220 figure
221 ma.plotTracks;
222 ma.labelPlotTracks;
223 
224 ma = ma.computeMSD;
225 ma.msd;
226 
227 t = (0 : N_TIME_STEPS)' * dT;
228 [T1, T2] = meshgrid(t, t);
229 all_delays = unique( abs(T1 - T2) );
230 
231 figure
232 ma.plotMSD;
233 
234 
235 cla
236 ma.plotMeanMSD(gca, true)
237 
238 mmsd = ma.getMeanMSD;
239 t = mmsd(:,1);
240 x = mmsd(:,2);
241 dx = mmsd(:,3) ./ sqrt(mmsd(:,4));
242 errorbar(t, x, dx, 'k')
243 
244 [fo, gof] = ma.fitMeanMSD;
245 plot(fo)
246 ma.labelPlotMSD;
247 legend off
248 
249 
250 ma = ma.fitMSD;
251 
252 good_enough_fit = ma.lfit.r2fit > 0.8;
253 Dmean = mean( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
254 Dstd  =  std( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
255 
256 Dheader1 = ['Raw Unscaled Values'];
257 Dhead1 = ['    D        Dn        L'];
258 Ddat1 = [oD oDn oL];
259 disp(' ')
260 disp(Dheader1)
261 disp(Dhead1)
262 disp(Ddat1)
263 
264 
265 yourtesthead = ['YOU ARE TESTING DIFFUSION FOR:'];
266 if MSDtest(1)
267 	yourtest = ['   D:   original diffusion rate'];
268 elseif MSDtest(2)
269 	yourtest = ['   Dn:  new diffusion rate'];
270 elseif MSDtest(3)
271 	yourtest = ['   L:  step length'];
272 else
273 	yourtest = ['   generic diffusion rate'];
274 end
275 disp(yourtesthead)
276 disp(yourtest)
277 
278 disp(' ')
279 fprintf('Estimation of raw D coefficient from MSD:\n')
280 fprintf('D = %.3g ± %.3g (mean ± std, N = %d)\n', ...
281     Dmean, Dstd, sum(good_enough_fit));
282 
283 
284 
285 
286 % Retrieve instantaneous velocities, per track
287  trackV = ma.getVelocities;
288 
289  % Pool track data together
290  TV = vertcat( trackV{:} );
291 
292  % Velocities are returned in a N x (nDim+1) array: [ T Vx Vy ...]. So the
293  % velocity vector in 2D is:
294  V = TV(:, 2:3);
295 
296  % Compute diffusion coefficient
297 varV = var(V);
298 mVarV = mean(varV); % Take the mean of the two estimates
299 Dest = mVarV / 2 * dT;
300 
301 
302 
303 Dheader2 = ['Scaling to model...'];
304 Dhead2 = ['    D        Dn        L'];
305 Ddat2 = [D Dn L];
306 
307 disp(' ')
308 disp(Dheader2)
309 disp(Dhead2)
310 disp(Ddat2)
311 fprintf('Estimation from velocities histogram:\n')
312 fprintf('Tested D = %.3g %s, compare to scaled Des value of %.3g %s\n', ...
313     Dest, [SPACE_UNITS '²/' TIME_UNITS], D, [SPACE_UNITS '²/' TIME_UNITS]);
314 
315 % printf('D.psd target value was %.3g %s\n', ...
316 %     Dest, msdDpsd, [SPACE_UNITS '²/' TIME_UNITS]);
317 
318 end



SI Fig. 3 - Title

Text Text Text


SI Fig. 4 - Title

Text Text Text


SI Fig. 5 - Title

Text Text Text