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