81 lines
2.7 KiB
Mathematica
81 lines
2.7 KiB
Mathematica
|
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
|