%This program simulates the diffusion of antibodies from a capillary into
%the surrounding tissue in a tumor.  Cylindrical geometry is used to
%simulate diffusion outward from a capillary.  A mixed (Robin) boundary
%condition is used at the capillary wall (see below), a no flux Neumann boundary
%condition is used at the edge of the Krogh cylinder.  This is in effect
%similar to a method of mirrors where an equal diffusion from surrounding
%capillaries prevents further diffusion from the original.  A biexponential
%can be used to simulate decaying plasma concentrations, or a
%fixed value can be used to simulate a constant infusion.  The method of
%lines is used to simulate the free antibody, free antigen, and bound
%antibody-antigen complex concentrations due to binding, release,
%diffusion, catabolism, synthesis, and degradation.  Effective diffusion
%constants are used to treat the pseudo-homogenous tissue.
%The total time and number of points can be adjusted to achieve total
%saturation and accurate modeling in the most efficient amount of
%computation time.

%flux = P(epsilon*R*Cp - Ce) where Ce is extracellular concentration, Cp is
%the plasma concentration, R is the partition coefficient, and epsilon is 
%the void fraction.

%This parent code (for diffusion in spheroids) was revised 10/4/04 to 
%define all concentrations in terms of total tumor volume.  The void 
%fraction epsilon is included to convert to effective concentrations where 
%needed.
%
%Greg Thurber, MIT, 2005
%
%A menu-driven front end was added to allow facile parameter inputs.
%Also, characteristic times are calculated and displayed.
%KDW 11/6/06


function [iflag] = Vascularized_Tumor_Penetration()

clear all;
close all;

tic;

iflag = 0;

%__________________________________________________________________________
%Between the above and next lines are the input parameters for the model


%k_off is calculated using these parameters:
Param.k_on = 1e5;    %per M*s
Param.k_off= 2.e-6;  %per s
abparamgo=99;
while abparamgo == 99;
    disp(['1) scFv in mouse']);
    disp(['2) scFv in human']);
    disp(['3) PEGylated scFv in mouse']);
    disp(['4) human IgG in mouse']);
    disp(['5) mouse IgG in mouse']);
    disp(['6) human IgG in human']);
    abparaminput=input('Choose one of the above:');
    if (abparaminput ==1)
        Param.Ab_diffusion = 80e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.3;    %tumor void fraction
        Param.P = 5e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 0.05;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 3;
        fraction_alpha = 0.8;
        abparamgo=0;
    end
    if (abparaminput ==2)
        Param.Ab_diffusion = 80e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.3;    %tumor void fraction
        Param.P = 5e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 0.42;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 5.32;
        fraction_alpha = 0.9;
        abparamgo=0;
    end
    if (abparaminput ==3)
        Param.Ab_diffusion = 16e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.1;    %tumor void fraction
        Param.P = 5e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 1.5;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 5.32;
        fraction_alpha = 0.9;
        abparamgo=0;
    end
    if (abparaminput ==4)
        Param.Ab_diffusion = 14e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.1;    %tumor void fraction
        Param.P = 3.e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 1.2;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 6.81;
        fraction_alpha = 0.66;
        abparamgo=0;
    end
    if (abparaminput ==5)
        Param.Ab_diffusion = 14e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.1;    %tumor void fraction
        Param.P = 3.e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 1.;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 7;
        fraction_alpha = 0.7;
        abparamgo=0;
    end
        if (abparaminput ==6)
        Param.Ab_diffusion = 14e-12;   %m^2/s antibody diffusion constant
        Param.epsilon = 0.1;    %tumor void fraction
        Param.P = 3.e-9; %permeability of capillary in m/s
        t_half_blood_alpha = 12.74;   %set in HOURS, converted to k_pharmacokinetic below
        t_half_blood_beta = 86.92;
        fraction_alpha = 0.43;
        abparamgo=0;
    end
end


Param.R_tumor = 100e-6;   %radius of tumor in meters
Param.R_capillary = 10e-6;   %Radius of capillary in meters

Param.Agen_initial = 3e-7;  %mol/tumor volume
Param.Abody_initial = 2e-6; %M (plasma conc!)  converted to tumor conc in pk function

%Endocytic trafficking parameters
metabhalflife=12;
Param.ke = log(2)/(metabhalflife*3600);    %1.6e-5 is approximately a 12 hour half-life
metabhalflifebound=12;
Param.ke2 = Param.ke;  %different rate if bound?
Param.Rs = Param.ke.*Param.Agen_initial;   %if steady state, Rs should equal ke*Agen_initial


Param.infusion_t = 3600*0;  %infusion time in SECONDS
%if performing integration for specified time period
time_interval = Param.infusion_t/3600 + 5.*max([t_half_blood_alpha, t_half_blood_beta]); %set total time in HOURS for simulation
      

%converting time interval to seconds
tfinal = 3600*time_interval;   %seconds

%Input parameters in structure
Param.num_pts = 60;   %number of mesh points for simulation, ~60 gives decent (unvariable) results

paraminput=99;
while paraminput ~= 0;
    Param.K_equil=Param.k_off/Param.k_on;  %calculate K_equil
    disp([' ']);
    disp(['Binding parameters:']);
    disp(['(1) k_on = ',num2str(Param.k_on),' /M/s.']);
    disp(['(2) k_off = ',num2str(Param.k_off),' /M/s.']);
    disp(['    K_d = ',num2str(Param.K_equil),' M.']);
    disp([' ']);
    disp(['Tumor parameters:']);
    disp(['(3) [Antigen] =',num2str(Param.Agen_initial),' mol/tumor volume.']);
    disp(['(4) epsilon (antibody-accessible fraction) =',num2str(Param.epsilon)]);
    disp(['(5) Capillary radius = ',num2str(Param.R_capillary),' meters.']);
    disp(['(6) Free antigen metabolic half life =',num2str(metabhalflife),' hr.']);
    disp(['(7) Antibody-bound antigen metabolic halflife =',num2str(metabhalflifebound),' hr.']);
    disp(['(8) Radius of tumor cylinder considered = ',num2str(Param.R_tumor),' meters.']);
    disp([' ']);
    disp(['Transport parameters:']);
    disp(['(9) Antibody diffusivity = ',num2str(Param.Ab_diffusion),' m^2/s.']);
    disp(['(10) Capillary permeability coefficient = ',num2str(Param.P),' m/s.']);
    disp([' ']);
    disp(['Pharmacokinetic parameters:']);
    disp(['(11) Alpha half life = ',num2str(t_half_blood_alpha),' hr.']);
    disp(['(12) Beta half life = ',num2str(t_half_blood_beta),' hr.']);
    disp(['(13) Alpha fraction = ',num2str(fraction_alpha)]);
    disp([' ']);
    disp(['Antibody dosage:']);
    disp(['(14) [Antibody]initial = ',num2str(Param.Abody_initial),' mol/L.']);
    disp([' ']);
    disp(['Simulation parameters:']);
    disp(['(15) Final time = ',num2str(tfinal),' seconds.'])
    disp(['(16) Number of mesh points = ',num2str(Param.num_pts),' .'])
    disp([' ']);
    paraminput=input('Type 0 to execute simulation or parameter number to change:');
    if (paraminput == 1)
        Param.k_on=input('k_on (/M/s)?');
    end
    if (paraminput == 2)
        Param.k_off=input('k_off (/M/s)?');
    end
    if (paraminput == 3)
        Param.Agen_initial=input('[Antigen] (mol/tumor volume)?');
    end
    if (paraminput == 4)
        Param.epsilon=input('epsilon?');
    end
    if (paraminput == 5)
        Param.R_capillary=input('Capillary radius? (meters)?');
    end
    if (paraminput == 6)
        metabhalflife=input('Free antigen metabolic half life (hr)?');
        Param.ke = log(2)/(metabhalflife*3600);
        Param.Rs = Param.ke.*Param.Agen_initial;
    end
    if (paraminput == 7)
        metabhalflifebound=input('Antibody-antigen metabolic half life (hr)?');
        Param.ke2 = log(2)/(metabhalflifebound*3600);
    end
    if (paraminput == 8)
        Param.R_tumor=input('Radius of tumor cylinder Rcap (meters)?');
    end
    if (paraminput == 9)
        Param.Ab_diffusion=input('Antibody diffusivity D (m^2/s)?');
    end
    if (paraminput == 10)
        Param.P=input('Capillary permeability coefficient P (m/s)?');
    end
    if (paraminput == 11)
        t_half_blood_alpha=input('Alpha half life (hr)?');
    end
    if (paraminput == 12)
        t_half_blood_beta=input('Beta half life (hr)?');
    end
    if (paraminput == 13)
        fraction_alpha=input('Alpha fraction?');mi
    end
    if (paraminput == 14)
        Param.Abody_initial=input('[Antibody]initial (M)?');
    end
    if (paraminput == 15)
        tfinal=input('Final time (seconds)?');
    end
    if (paraminput == 16)
        Param.num_pts=input('Number of mesh points?');
    end
end
                                  

                    

%common parameters are changed above this line
%__________________________________________________________________________



%convert half life in minutes to k in s^-1
Param.pharma_k_alpha = log(2)/(t_half_blood_alpha*3600); 
Param.pharma_k_beta = log(2)/(t_half_blood_beta*3600); 
Param.fraction_alpha = fraction_alpha;

%display characteristic moduli
disp([' ']);

Thiele_modulus = ((Param.R_tumor^2).*(Param.Agen_initial./Param.epsilon).*Param.ke./(Param.Ab_diffusion.*(Param.Abody_initial./...
    (1+(Param.Ab_diffusion./(2.*Param.R_capillary.*Param.P))))))^0.5;

disp(['Thiele modulus = ',num2str(Thiele_modulus),]);


Clearance_modulus = (Param.R_tumor^2).*(Param.Agen_initial./Param.epsilon)./...
    (2.*Param.Abody_initial.*Param.R_capillary.*Param.P.*(Param.infusion_t + ...
    (exp(-Param.pharma_k_alpha.*Param.infusion_t).*Param.fraction_alpha./Param.pharma_k_alpha)+...
    (exp(-Param.pharma_k_beta.*Param.infusion_t).*(1-Param.fraction_alpha)./Param.pharma_k_beta)));

disp(['Clearance modulus = ',num2str(Clearance_modulus),]);

biot=2*Param.P*Param.R_capillary/Param.Ab_diffusion;
disp(['Biot number = ',num2str(biot),]);

disp([' ']);

endocytosis_reach = 1.7*1e6.*(Param.P.*(log(2)./Param.ke).*Param.R_capillary.*(Param.epsilon.*Param.Abody_initial./...
    Param.Agen_initial))^0.5;
disp(['Penetration distance limited by endocytosis = ',num2str(endocytosis_reach),' microns.']);

clearance_reach = 1.4*1e6.*(Param.P.*(Param.infusion_t + ...
    (exp(-Param.pharma_k_alpha.*Param.infusion_t).*Param.fraction_alpha./Param.pharma_k_alpha)+...
    (exp(-Param.pharma_k_beta.*Param.infusion_t).*(1-Param.fraction_alpha)./Param.pharma_k_beta))...
    .*Param.R_capillary.*(Param.epsilon.*Param.Abody_initial./Param.Agen_initial))^0.5;
disp(['Penetration distance limited by clearance = ',num2str(clearance_reach),' microns.']);

disp([' ']);
disp(['Characteristic times (seconds):']);
tdiff=(Param.R_tumor^2)./Param.Ab_diffusion;
disp(['Nonbinding diffusive time =',num2str(tdiff)]);
tclear=fraction_alpha*t_half_blood_alpha*3600 + (1-fraction_alpha)*t_half_blood_beta*3600;
disp(['Clearance time =',num2str(tclear)]);
tmetab=1./Param.ke;
disp(['Metabolic time =',num2str(tmetab)]);
tbindag=1./(Param.k_on.*Param.Agen_initial./Param.epsilon);
disp(['Binding time (Ab perspective) =',num2str(tbindag)]);
tbindab=1./(Param.k_on.*Param.Abody_initial./Param.epsilon);
disp(['Binding time (Ag perspective) =',num2str(tbindab)]);

%set up rows for initial conditions
initial_antigen = Param.Agen_initial*ones(1,Param.num_pts);
initial_antibody = zeros(1,Param.num_pts);
initial_bound = zeros(1,Param.num_pts);

initial = [initial_antibody initial_antigen initial_bound];

opt = odeset('AbsTol',1e-8);    %decrease tolerance

%Invoke ODE solver
[t,Y] = ode23s(@differential_equations, [0 tfinal], initial, opt, Param);

t_minutes = t./60;
%Plot output

%make x axis values

%Add the capillary radius to the dimensions so will have capillary in
%between
x_minus = [linspace(-Param.R_tumor, -Param.R_capillary, Param.num_pts), -Param.R_capillary, 0];
x_plus = [0, Param.R_capillary, linspace(Param.R_capillary, Param.R_tumor,Param.num_pts)];
x_axis = 1.e6*[x_minus, x_plus]';



for counter = 1:length(t)
    plasma_conc(counter) = calc_conc(Param.Abody_initial,t(counter),Param);
end;
%The above concentrations will be in overall concentrations, not plasma
%concentrations!
%figure;
%plot(t_minutes,plasma_conc/Param.epsilon);


%Bound receptors
%figure;
Z = [zeros(1,length(plasma_conc))', zeros(1,length(plasma_conc))', Y(:,(1+2*Param.num_pts):3*Param.num_pts)];

%Create mirror image so can visualize both "sides" of cylinder
for q = 1:Param.num_pts+2
    Zb(:,q) = Z(:,Param.num_pts + 3 - q);
end;
mirror = [Zb Z];


%meshc(x_axis, t_minutes, mirror);
%xlabel('radius (0 at center)');
%ylabel('Time (min)');
%zlabel('Concentration of bound receptor (M)');
%title('Concentration of Bound Receptors vs. Time and Radius');

figure;
plot(x_axis,mirror(length(t),:)/Param.Agen_initial);
xlabel('Radius (microns)');
ylabel('Fraction of bound receptor');
title('Bound antibody profile at final time point');


figure;
meshc(x_axis, t_minutes, mirror./Param.Agen_initial);
xlabel('Radius (microns; 0 at center)');
ylabel('Time (min)');
zlabel('Fraction of bound receptor');
title('Fraction of Bound Receptors vs. Time and Radius');



toc;
time = toc;

iflag = 1;

return;



%__________________________________________________________________________
%__________________________________________________________________________
%__________________________________________________________________________




%This function defines the derivatives that will be solved by a stiff
%solver.  The origin of these equations comes from material balances on the
%antibody, antigen, and bound complex concentrations.

function dY_dt = differential_equations(t,Y, Param)


dY_dt = zeros(3*Param.num_pts,1);
dR = Param.R_tumor/Param.num_pts;   %define step distance


%unpack parameters
D = Param.Ab_diffusion;
k_on = Param.k_on;
k_off = Param.k_off;
ke = Param.ke;
ke2 = Param.ke2;
Rs = Param.Rs;
Ab_exterior = calc_conc(Param.Abody_initial,t,Param);  %call function
epsilon = Param.epsilon;
R_capillary = Param.R_capillary;
P = Param.P;

%Define boundary conditions
for i = 1:Param.num_pts
    
    Ab = Y(i);
    Ag = Y(i+Param.num_pts);
    B = Y(i+2*Param.num_pts);
    
    %Because there is a void area for the capillary, the current radius
    %must be computed and used in the second derivative.  This value
    %calculates the correct current radius for each step.
    R_current = R_capillary + i*dR;
        
    if (i==1) %Robin BC at capillary wall.
              
        Abplus = Y(i+1);
        Abminus = (2*P*dR*Ab_exterior + D*(4*Ab - Abplus))./(2*P*dR + 3*D);
        
        %Ab_exterior is overall concentration, not plasma conc.
        %For higher plasma concentrations, the second term should be
        %subtracted from the first, so using the minus sign means the
        %derivative should be positive
        
    else
        if (i==Param.num_pts)  %Neumann no flux BC at outer edge
        
        Abplus = (4/3)*Y(i) - (1/3)*Y(i-1);
        Abminus = Y(i-1);
        
        else    %interior points

            Abplus = Y(i+1);  
            Abminus = Y(i-1);                       
            
        end
    end
    
%compute derivatives
%All concentrations are overall, not effective.

%Antibody
%Only change from spherical coordinates to cylindrical is the addition of a
%2 in the denominator of the second part of the 2nd derivative.  Normally, 
%a 2 in the numerator that is the shape factor for spheres is cancelled out
%by a 2 in the denominator from the central difference method.  The shape 
%factor for a cylinder is 1 (and 0 for planar).
dY_dt(i) = D.*(((Abminus - 2*Ab + Abplus)./(dR^2)) + (Abplus - Abminus)./(2*dR*R_current)) - k_on.*Ab.*Ag./epsilon + k_off.*B;

%Antigen
dY_dt(i+Param.num_pts) = Rs - k_on*Ab*Ag./epsilon + k_off*B - ke*Ag;

%Bound complex
dY_dt(i+2*Param.num_pts) = k_on*Ab*Ag./epsilon - k_off*B - ke2*B;

end

return;



%__________________________________________________________________________
%__________________________________________________________________________
%__________________________________________________________________________




%This function provides a single compartment pharmacokinetic model for the
%external concentration of antibody.  It provides the antibody
%concentration in the tumor, NOT the effective concentration


%This function provides a single compartment pharmacokinetic model for the
%external concentration of antibody.  It is given plasma concentrations but
%returns the overall concentration.


function Ab_conc = calc_conc(Ab_initial,time,Param)

Co_plasma = Ab_initial;    %rename parameters
t = time;
k_alpha = Param.pharma_k_alpha;
k_beta = Param.pharma_k_beta;
A = Param.fraction_alpha;
B = (1-A);
Co = Param.epsilon.*Co_plasma;  %convert plasma conc to tumor conc

if time < Param.infusion_t
    
    Ab_conc = Co;   %use if constant concentration

else
    adjusted_t = time - Param.infusion_t;
    
    Ab_conc = Co*(A.*exp(-k_alpha*adjusted_t)+B.*exp(-k_beta*adjusted_t)); %calculates exponential decay
    
end;


return;