%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 mirros 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
%


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





%Input parameters in structure
Param.num_pts = 60;   %number of mesh points for simulation, ~60 gives decent (unvariable) results
Param.Ab_diffusion = 80e-12;   %m^2/s antibody diffusion constant

%k_off is calculated using these parameters:
Param.k_on = 1e5;    %per M*s
Param.K_equil = 20e-12;  %Molar - sm3E is 2e-11


Param.R_tumor = 100e-6;   %radius of tumor 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
Param.epsilon = 0.3;    %tumor void fraction

Param.R_capillary = 10e-6;   %Radius of capillary in meters
Param.P = 5e-9; %permeability of capillary in m/s

%Trafficking parameters
Param.ke = 1.6e-5;    %3.5e-5 is approximately a 5 hour half-life
Param.ke2 = Param.ke;  %different rate if bound?
Param.Rs = Param.ke.*Param.Agen_initial;   %if steady state, Rs should equal ke*Agen_initial

%set half life for antibody in blood for global pharmacokinetic model


t_half_blood_alpha = 0.05;   %set in HOURS, converted to k_pharmacokinetic below
t_half_blood_beta = 3;
fraction_alpha = 0.8;

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
                                        

                    

%common parameters are changed above this line
%__________________________________________________________________________


%calculate off rate
Param.k_off = Param.k_on*Param.K_equil;  %per s

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

%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;

%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 = [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;
meshc(x_axis, t_minutes, mirror./Param.Agen_initial);
xlabel('radius (0 at center)');
ylabel('Time (min)');
zlabel('Fraction of bound receptor');
title('Fraction of Bound Receptors vs. Time and Radius');

%Free antigen
figure;
Z1 = [zeros(1,length(plasma_conc))', zeros(1,length(plasma_conc))', Y(:,(1+Param.num_pts):2*Param.num_pts)];

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

meshc(x_axis, t_minutes, mirror2);
xlabel('radius (0 at center)');
ylabel('Time (min)');
zlabel('Concentration of free antigen (M)');
title('Concentration of Free Antigen vs. Time and Radius');

%Free antibody
figure;
Z2 = [plasma_conc'./1000, plasma_conc'./1000, Y(:,1:Param.num_pts)];

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

meshc(x_axis, t_minutes, mirror3);
xlabel('radius (0 at center)');
ylabel('Time (min)');
zlabel('Concentration of free antibody(M)');
title('Concentration of Free Antibody vs. Time and Radius');

%Total antigen material balance
figure;
total_antigen = Y(:,(1+2*Param.num_pts):3*Param.num_pts) + Y(:,(1+Param.num_pts):2*Param.num_pts);


meshc(total_antigen);
xlabel('radius (0 at center)');
ylabel('Time (min)');
zlabel('Concentration of free antibody(M)');
title('Concentration of Free Antibody vs. Time and Radius');

pretargeting_number = 1/(Param.ke.*(Param.fraction_alpha./Param.pharma_k_alpha+(1-Param.fraction_alpha)./Param.pharma_k_beta));
disp(['The pretargeting number for this simulation is: ',num2str(pretargeting_number)]);

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(['The endocytosis reach for this simulation is: ',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(['The clearance reach for this simulation is: ',num2str(clearance_reach),' microns.']);

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

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

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;