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_Lb0 + C_Lbf) / 2;
    qr = 0;
end

end