Qual

From Bradwiki
Jump to: navigation, 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