From fe9a57d519bb90c276b6d7056b6be0ed9a29556f Mon Sep 17 00:00:00 2001 From: ndwarshuis Date: Tue, 3 Aug 2021 19:01:31 -0400 Subject: [PATCH] ADD appendix for matlab stuff --- .gitignore | 5 +- code/diffusion_mab.m | 61 +++++++++ code/diffusion_stp.m | 80 ++++++++++++ code/microcarrier_diffusion_washing.m | 179 ++++++++++++++++++++++++++ code/washing_out.tex | 127 ++++++++++++++++++ tex/thesis.tex | 50 ++++++- 6 files changed, 500 insertions(+), 2 deletions(-) create mode 100644 code/diffusion_mab.m create mode 100644 code/diffusion_stp.m create mode 100644 code/microcarrier_diffusion_washing.m create mode 100644 code/washing_out.tex diff --git a/.gitignore b/.gitignore index 776e9ea..0386601 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,7 @@ figures/* !tables tables/* -!tables/*.tex \ No newline at end of file +!tables/*.tex + +!code +!code/* \ No newline at end of file diff --git a/code/diffusion_mab.m b/code/diffusion_mab.m new file mode 100644 index 0000000..87aae2c --- /dev/null +++ b/code/diffusion_mab.m @@ -0,0 +1,61 @@ +function diffusion_mab + +% bulk properties +V_b = 1.075; % volume of bulk solution (cm^3) +n = 16000; % number of microcarriers in bulk volume + +% ligand properties (IgG) +m_L0 = 4; % mass of ligand (ug) +MW_L = 150000; % molecular weight of ligand (g/mol) +D_LW = 4.8e-7; % ligand diffusion coefficient in water (cm^2/s) +a = 1; % rxn stoic coeff (L per R) + +% fixed carrier properties +geometry = 0.190; % geometric diffusion factor of microcarriers +R = 0.01275; % microcarrier radius (cm) +C_R0 = 144; % receptor density in carriers (pM/cm^3, nM) +D_L_app = D_LW*geometry*60; % apparent ligand diffusion coeff. (cm^2/min) + +% init +t_0 = 0; % initial time (min) +t_f = 90; % final time (min) + +alpha = 0.999999999999; % fudge factor to prevent MATLAB implosion +r_i0 = R * alpha; % init interfatial radius in microcarriers (cm) + +C_Lb0 = m_L0 / MW_L * 1e6 / V_b; % init conc. of ligand in bulk (pmol/cm^3, nM) + +% y1 = r_i, y2 = C_Lb +dydt = @(t,y) [ + ((R * D_L_app * y(2) / a))/(C_R0 * y(1) * (y(1) - R)); + (4 * pi * n * D_L_app * y(2) * R * y(1)) / (V_b * (y(1) - R))]; + +[t,Y] = ode45(dydt,[t_0 t_f],[r_i0 C_Lb0]); + +r_i = Y(:,1); % interfacial radius (cm) +C_Lb = Y(:,2); % bulk ligand concentration (pmol/cm^3, nM) + +% fluxes (pmol/cm^2/min) +N_L = (D_L_app*C_Lb.*r_i)./(R^2*(r_i/R+(-1))); + +% NOTE: flow rates (pmol/min) will be the same at interface and at r=R due +% to pseudo-steady state assumption +w_L = -(4*pi*R^2*N_L); + +% concentration profiles (nM) as function of interfacial radius +C_r = diff(cumtrapz(t,-w_L))./diff(4/3*pi*r_i.^3); + +% total bound mass (ug) +m_bt_L = trapz(t,w_L)*MW_L/1e6*n; + +disp(['total bound mass at ',num2str(t(end)),' min']); +disp([num2str(m_bt_L),' ug, ', num2str(m_bt_L/m_L0*100), '% (ligand)']); + +figure; +plot(t,C_Lb*MW_L/1e6); +xlabel('t (min)'); +ylabel('concentration (ug/ml)'); +legend('ligand'); +title ('bulk concentration'); + +end diff --git a/code/diffusion_stp.m b/code/diffusion_stp.m new file mode 100644 index 0000000..6346e73 --- /dev/null +++ b/code/diffusion_stp.m @@ -0,0 +1,80 @@ +function diffusion_stp + +% ADJUST THIS TO MINIMIZE SSE +geometry = 0.190; % geometric diffusion factor of microcarriers + +% bulk parameters +V_b = 1.008; % volume of bulk solution (cm^3) +n = 16000; % number of microcarriers in bulk volume + +% ligand properties (streptavidin) +m_L0 = 40; % mass of ligand (ug) +MW_L = 55000; % molecular weight of ligand (g/mol) +D_LW = 6.2e-7; % ligand diffusion coeff. in water (cm^2/s) +a = 1; % rxn stoic coeff (L per R) + +% fixed carrier properties +R = 0.01275; % microcarrier radius (cm) +C_R0 = 3403; % receptor density in carriers (pM/cm^3, nM) +D_L_app = D_LW * geometry * 60; % effective diffusion coeff (cm^2/min) + +% init +t_0 = 0; % initial time (min) +t_f = 65; % final time (min) + +alpha = 0.999999999999; % fudge factor to prevent MATLAB implosion +r_i0 = R * alpha; % init. interfatial radius in microcarrier (cm) + +C_Lb0 = m_L0 / MW_L * 1e6 / V_b; % init. conc. of ligand in bulk (pmol/cm^3, nM) + +% y1 = r_i, y2 = C_Lb +dydt = @(t,y) [ + ((R * D_L_app * y(2) / a)) / (C_R0 * y(1) * (y(1) - R)); + (4 * pi * n * D_L_app * y(2) * R * y(1)) / (V_b * (y(1) - R))]; + +[t,Y] = ode45(dydt, [t_0 t_f], [r_i0 C_Lb0]); + +r_i = Y(:, 1); % interfacial radius (cm) +C_Lb = Y(:, 2); % bulk ligand concentration (pmol/cm^3, nM) + +% fluxes (pmol/cm^2/min) +N_L = (D_L_app * C_Lb .* r_i) ./ (R^2 * (r_i / R + (-1))); + +% NOTE: flow rates (pmol/min) will be the same at interface and at r=R due +% to pseudo-steady state assumption +w_L = -(4 * pi * R^2 * N_L); + +% concentration profiles (nM) as function of interfacial radius +C_r = diff(cumtrapz(t, -w_L)) ./ diff(4 / 3 * pi * r_i.^3); + +% total bound mass (ug) +m_bt_L = trapz(t, w_L) * MW_L / 1e6 * n; + +disp(['total bound mass at ', num2str(t(end)),' min']); +disp([num2str(m_bt_L), ' ug, ', num2str(m_bt_L / m_L0 * 100), '% (ligand)']); + +% experimental data +exp_time1 = [7.5 15 30 60]; +exp_time2 = [15 20 25 30]; +exp_conc1 = [20.20 17.15 13.78 13.99]; % bulk concentration in ug/ml +exp_conc2 = [18.73 16.74 14.91 14.75]; + +figure; +plot(exp_time1,exp_conc1, '.'); +hold on; +plot(exp_time2,exp_conc2, '.'); +plot(t, C_Lb * MW_L / 1e6); +hold off; +xlabel('t (min)'); +ylabel('concentration (ug/ml)'); +legend('exp1', 'exp2', 'ligand'); +title ('bulk concentration'); + +% fitting to experimental data (least squares algorithm) +e1 = interp1(t, C_Lb, exp_time1) * MW_L / 1e6 - exp_conc1; +e2 = interp1(t, C_Lb, exp_time2) * MW_L / 1e6 - exp_conc2; + +% MINIMIZE THIS (SSE) +disp(['SSE: ', num2str(sum([e1 e2].^2))]); + +end diff --git a/code/microcarrier_diffusion_washing.m b/code/microcarrier_diffusion_washing.m new file mode 100644 index 0000000..a468378 --- /dev/null +++ b/code/microcarrier_diffusion_washing.m @@ -0,0 +1,179 @@ +function microcarrier_diffusion_washing() + +% initial/reaction volume for carriers (cm^3) +V_b0 = 1; + +disp(['###Biotin washing###', sprintf('\n')]) +n_biotin = 10; % nmol biotin remaining after attachment +D_biotin = 5.0e-6; % diffusion coeff (cm^2/s) + +carrier_stages(V_b0, 3, n_biotin, D_biotin) + +disp(['###Streptavidin washing###', sprintf('\n')]) +m_stp = 15; % ug stp remaining after attachment +n_stp = m_stp / 55000 * 1000; +D_stp = 6.2e-7; % diffusion coeff (cm^2/s) + +carrier_stages(V_b0, 2, n_stp, D_stp) + +disp(['###Antibody washing###', sprintf('\n')]) +m_Ab = 1; % ug Ab added +n_Ab = m_Ab / 150000 * 1000; +D_Ab = 4.8e-7; % diffusion coeff (cm^2/s) + +carrier_stages(V_b0, 2, n_Ab, D_Ab) + +end + +function carrier_stages(V_b0, stages, n_L, D_L) + +% initial concentration of L (nM) +C_L = n_L / (V_b0 / 1000); + +% max fill for 15 ml conical tube is about 10ml at vortex speed = 6 without +% splashing the medium (PBS) onto the cap/rim and possibly compromising +% sterility and losing carriers +% +% wash volume, fill up to this to let biotin out of carriers (cm^3) +V_wash = 10; + +% dilution volume, fill up to this to reduce total biotin in solution (cm^3) +V_dilution = 15; + +C_L_bulk = C_L; + +for i = 1:stages + disp([sprintf('\n'), '---Stage ', num2str(i), '---', sprintf('\n')]) + [C_L, n_L_carrier, n_L_bulk] = carrier_diffusion_out(V_b0, ... + V_wash, C_L, C_L_bulk, D_L); + + C_L_bulk = C_L*V_wash / V_dilution; + n_L_bulk = n_L_bulk * V_b0 / V_dilution; +end + +disp([sprintf('\n'), 'final bulk amount = ' , num2str(n_L_bulk), ' nmol']); +disp(['final total amount = ', num2str(n_L_carrier + n_L_bulk), ' nmol']); +disp(['overall reduction = ', ... + num2str((1 - (n_L_carrier + n_L_bulk) / n_L) * 100), ' %', ... + sprintf('\n')]); + +end + +function [C_carrier_ave_final, n_Lc_f, n_Lb_f] = ... + carrier_diffusion_out(V_b0, V_bf, C_Lc_pw, C_Lb_pw, D_LW) + +% V_b0: prewash volume of bulk solution (cm^3) +% V_bf: wash volume of bulk solution (cm^3) +% C_Lc_pw: prewash concentration of ligand in carriers (pmol/cm^3, nM) +% C_Lb_pw: prewash concentration of ligand in bulk (pmol/cm^3, nM) +% D_LW: ligand diffusion coefficient in water (cm^2/s) + +% number of microcarriers in bulk volume +n = 16000; + +% geometric diffusion factor of microcarriers +geometry = 0.190; + +% microcarrier radius (cm) +R = 0.01275; + +% apparent ligand diffusion coeff in microcarrier (cm^2/min) +D_L_app = D_LW * geometry * 60; + +% void fraction +void = 0.95; + +% volume occupied by carriers (cm^3) +V_c = void * 4 / 3 * pi * R^3 * n; + +% prewash amount in carriers (nmol) +n_Lc_pw = (V_c / 1000) * C_Lc_pw; + +% prewash amount outside carriers (nmol) +n_Lb_pw = ((V_b0 - V_c) / 1000) * C_Lb_pw; + +t_0 = 0; % initial time (min) +t_f = 30; % final time (min) +C_Lb0 = n_Lb_pw / ((V_bf - V_c) / 1e3); % initial conc. in bulk (nM) + +% final conditions (after long time) +% we use this calculate an average bulk concentration to create a constant +% boundary in the BVP, otherwise we have a to solve a free BVP, +% which involves sacrificing a kitten... +n_L_trans = (V_bf * n_Lc_pw - V_c * n_Lb_pw) / (V_c + V_bf); + +% final concentration of ligand in bulk solution (pmol/cm^3, nM) +C_Lbf = (n_Lb_pw + n_L_trans) / (V_bf / 1e3); + +disp(['initial bulk concentration = ', num2str(C_Lb0), ' nM']); +disp(['final bulk concentration = ', num2str(C_Lbf), ' nM']); +disp(['bulk percent change = ', ... + num2str((1 - C_Lb0 / C_Lbf) * 100), '%']); +disp(['initial carrier concentration = ', num2str(C_Lc_pw), ' nM']); +disp(['initial carrier amount = ', num2str(n_Lc_pw), ' nmol']); +disp(['initial bulk amount = ', ... + num2str(n_Lb_pw), ' nmol', sprintf('\n')]); + +m = 2; % spherical +r = linspace(0, R, 50); +t_f = 1; +tolerance = 0.1; % iterate until center and outside conc are within this +increment = 1; % length of time increments (min) +diff = 1000000; % init to some huge number + +while diff > tolerance + t = linspace(t_0, t_f, 50); + + Y = pdepe(m,@pde,@init,@bound,r,t); + C = Y(:,:,1); + + % final concentration in center of carrier at final time + C_f_center = C(end,1); + + % test to see how far off the center is from bulk + diff = C_f_center / C_Lbf - 1; + t_f = t_f + increment; +end + +% average final concentration of carriers +C_carrier_ave_final = (C_f_center + C_Lbf) / 2; + +% nmol still remaining in carriers +n_Lc_f = C_carrier_ave_final * (V_c / 1000); + +% final nmol in bulk +n_Lb_f = n_Lb_pw + (n_Lc_pw - n_Lc_f); + +% percent nmol of ligand removed from carriers +perc_removed = (1 - n_Lc_f / n_Lc_pw) * 100; + +% amount actually transferred over theoretical +efficiency = (n_Lc_pw - n_Lc_f) / n_L_trans * 100; + +disp(['final time to equilibrium = ', num2str(t_f), ' min']); +disp(['Total removed from carriers = ', num2str(n_Lc_pw - n_Lc_f), ' nmol']); +disp(['Total remaining in carriers = ', num2str(n_Lc_f), ' nmol']); +disp(['Total remaining in bulk = ', num2str(n_Lb_f), ' nmol']); +disp(['Percent removed from carriers = ', num2str(perc_removed), ' %']); + +function [c, f, s] = pde(r, t, c, DcDr) + c = 1; + f = D_L_app * DcDr; + s = 0; +end + +function u0 = init(r) + u0 = C_Lc_pw; +end + +function [pl,ql,pr,qr] = bound(rl,cl,rr,cr,t) + pl = 0; + ql = 1; + % assume that the concentration boundary + % is the average of the initial and + % theoretical final concentration in bulk + pr = cr - (C_Lbf + C_Lbf)/2; + qr = 0; +end + +end diff --git a/code/washing_out.tex b/code/washing_out.tex new file mode 100644 index 0000000..a79c822 --- /dev/null +++ b/code/washing_out.tex @@ -0,0 +1,127 @@ +\begin{verbatim} +###Biotin washing### + + +---Stage 1--- + +initial bulk concentration = 879.6421 nM +final bulk concentration = 986.9753 nM +bulk percent change = 10.875% +initial carrier concentration = 10000 nM +initial carrier amount = 1.3197 nmol +initial bulk amount = 8.6803 nmol + +final time to equilibrium = 3 min +Total removed from carriers = 1.1882 nmol +Total remaining in carriers = 0.13145 nmol +Total remaining in bulk = 9.8686 nmol +Percent removed from carriers = 90.0394 % + +---Stage 2--- + +initial bulk concentration = 58.4119 nM +final bulk concentration = 69.8638 nM +bulk percent change = 16.3917% +initial carrier concentration = 996.0629 nM +initial carrier amount = 0.13145 nmol +initial bulk amount = 0.57641 nmol + +final time to equilibrium = 3 min +Total removed from carriers = 0.1221 nmol +Total remaining in carriers = 0.0093421 nmol +Total remaining in bulk = 0.69852 nmol +Percent removed from carriers = 92.8928 % + +---Stage 3--- + +initial bulk concentration = 4.1514 nM +final bulk concentration = 4.9653 nM +bulk percent change = 16.3917% +initial carrier concentration = 70.7919 nM +initial carrier amount = 0.0093421 nmol +initial bulk amount = 0.040967 nmol + +final time to equilibrium = 3 min +Total removed from carriers = 0.0086782 nmol +Total remaining in carriers = 0.00066396 nmol +Total remaining in bulk = 0.049645 nmol +Percent removed from carriers = 92.8928 % + +final bulk amount = 0.0033096 nmol +final total amount = 0.0039736 nmol +overall reduction = 99.9603 % + +###Streptavidin washing### + + +---Stage 1--- + +initial bulk concentration = 23.9902 nM +final bulk concentration = 26.9175 nM +bulk percent change = 10.875% +initial carrier concentration = 272.7273 nM +initial carrier amount = 0.035991 nmol +initial bulk amount = 0.23674 nmol + +final time to equilibrium = 14 min +Total removed from carriers = 0.032315 nmol +Total remaining in carriers = 0.003676 nmol +Total remaining in bulk = 0.26905 nmol +Percent removed from carriers = 89.7862 % + +---Stage 2--- + +initial bulk concentration = 1.6335 nM +final bulk concentration = 1.9538 nM +bulk percent change = 16.3917% +initial carrier concentration = 27.8557 nM +initial carrier amount = 0.003676 nmol +initial bulk amount = 0.01612 nmol + +final time to equilibrium = 15 min +Total removed from carriers = 0.0034097 nmol +Total remaining in carriers = 0.00026636 nmol +Total remaining in bulk = 0.019529 nmol +Percent removed from carriers = 92.7541 % + +final bulk amount = 0.001302 nmol +final total amount = 0.0015683 nmol +overall reduction = 99.4249 % + +###Antibody washing### + + +---Stage 1--- + +initial bulk concentration = 0.58643 nM +final bulk concentration = 0.65798 nM +bulk percent change = 10.875% +initial carrier concentration = 6.6667 nM +initial carrier amount = 0.00087977 nmol +initial bulk amount = 0.0057869 nmol + +final time to equilibrium = 17 min +Total removed from carriers = 0.00078901 nmol +Total remaining in carriers = 9.0766e-05 nmol +Total remaining in bulk = 0.0065759 nmol +Percent removed from carriers = 89.683 % + +---Stage 2--- + +initial bulk concentration = 0.040335 nM +final bulk concentration = 0.048242 nM +bulk percent change = 16.3917% +initial carrier concentration = 0.6878 nM +initial carrier amount = 9.0766e-05 nmol +initial bulk amount = 0.00039802 nmol + +final time to equilibrium = 18 min +Total removed from carriers = 8.4099e-05 nmol +Total remaining in carriers = 6.6675e-06 nmol +Total remaining in bulk = 0.00048212 nmol +Percent removed from carriers = 92.6543 % + +final bulk amount = 3.2141e-05 nmol +final total amount = 3.8809e-05 nmol +overall reduction = 99.4179 % +\end{verbatim} \ No newline at end of file diff --git a/tex/thesis.tex b/tex/thesis.tex index fe958ad..1680040 100644 --- a/tex/thesis.tex +++ b/tex/thesis.tex @@ -19,6 +19,31 @@ \usepackage[version=4]{mhchem} \usepackage{pgfgantt} \usepackage{setspace} +\usepackage{listings} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% my attempt to make MATLAB code look pretty + +\definecolor{dkgreen}{rgb}{0,0.6,0} +\definecolor{gray}{rgb}{0.5,0.5,0.5} +\definecolor{mauve}{rgb}{0.58,0,0.82} + +\lstset{frame=tb, + language=Matlab, + aboveskip=3mm, + belowskip=3mm, + showstringspaces=false, + columns=flexible, + basicstyle={\small\ttfamily}, + numbers=none, + numberstyle=\tiny\color{gray}, + keywordstyle=\color{blue}, + commentstyle=\color{dkgreen}, + stringstyle=\color{mauve}, + breaklines=true, + breakatwhitespace=true, + tabsize=3 +} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % benevolently force figures stay in their own subsection @@ -4446,7 +4471,30 @@ hosted using \gls{aws} using their proprietary Aurora implementation. The code is available here: \url{https://github.gatech.edu/ndwarshuis3/mdma}. -\chapter{reaction kinetics code} +\chapter{binding kinetics code} + +The \gls{stp} binding kinetic profile was fit and calculated using the following +MATLAB code. Note that the \inlinecode{geometry} parameter was varied so as to +minimize the \inlinecode{SSE} output. + +\lstinputlisting{../code/diffusion_stp.m} + +The geometric diffusivity from above (the \inlinecode{geometry} variable) was +used in the below code to obtain the reaction profile for the \gls{mab} binding +step. The model is the same except for the parameters which were changes to +reflect the \gls{mab} coating process. + +\lstinputlisting{../code/diffusion_mab.m} + +\chapter{washing kinetics code} + +The wash steps for the \gls{dms} were modeled using the following code: + +\lstinputlisting{../code/microcarrier_diffusion_washing.m} + +Complete output from this code is shown below: + +\input{../code/washing_out.tex} \chapter{references} \renewcommand{\chapter}[2]{} % noop the original bib section header