ADD appendix for matlab stuff

This commit is contained in:
Nathan Dwarshuis 2021-08-03 19:01:31 -04:00
parent e61817f192
commit fe9a57d519
6 changed files with 500 additions and 2 deletions

3
.gitignore vendored
View File

@ -13,3 +13,6 @@ figures/*
!tables
tables/*
!tables/*.tex
!code
!code/*

61
code/diffusion_mab.m Normal file
View File

@ -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

80
code/diffusion_stp.m Normal file
View File

@ -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

View File

@ -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

127
code/washing_out.tex Normal file
View File

@ -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}

View File

@ -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