%This program simulates the diffusion of antibodies into microspheroid %tumors. The method of lines simulates the free antibody, free %antigen, and bound antibody-antigen complex concentrations due to binding, %release, diffusion, catabolism, synthesis, and degradation. % %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. %This code 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, 2004 % function iflag = antibody_diffusion_pharmacokinetic_version() 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 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 = 130e-6; %radius of tumor in meters Param.Agen_initial = 1e-6; %mol/tumor volume Param.Abody_initial = 3e-9; %M (plasma conc!) converted to tumor conc in pk function Param.epsilon = 0.3; %tumor void fraction %trafficking parameters Param.ke = 1.6e-5; %1.6e-5 is approximately a 12 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 %if performing integration for specified time period time_interval = 10; %set total time in HOURS for simulation %set half life for antibody in blood for global pharmacokinetic model %(single compartment, no normal tissue) t_half_blood = 21; %set in MINUTES, converted to k_pharmacokinetic below %Note, if want constant concentration (i.e. for infusion or in vitro), set %in pk function at end of program %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 = log(2)/(t_half_blood*60); %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 x_axis = linspace(-Param.R_tumor, Param.R_tumor, Param.num_pts*2); %Bound receptors figure; Z = Y(:,(1+2*Param.num_pts):3*Param.num_pts); %Create mirror image so can visualize both "sides" of sphere for q = 1:Param.num_pts Zb(:,q) = Z(:,Param.num_pts + 1 - q); end; mirror = [Zb Z]; surfc(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; surfc(x_axis, t_minutes, (mirror./Param.Agen_initial)); xlabel('radius (0 at center)'); ylabel('Time (min)'); zlabel('Fraction of Ag bound'); title('Concentration of Bound Receptors vs. Time and Radius'); axis([-Param.R_tumor Param.R_tumor 0 t_minutes(end) 0 1]); %Free antigen figure; Z1 = Y(:,(1+Param.num_pts):2*Param.num_pts); %Create mirror image so can visualize both "sides" of sphere for q = 1:Param.num_pts Zc(:,q) = Z1(:,Param.num_pts + 1 - q); end; mirror2 = [Zc Z1]; surfc(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 = Y(:,1:Param.num_pts); %Create mirror image so can visualize both "sides" of sphere for q = 1:Param.num_pts Zd(:,q) = Z2(:,Param.num_pts + 1 - q); end; mirror3 = [Zd Z2]; surfc(x_axis, t_minutes, mirror3); xlabel('radius (0 at center)'); ylabel('Time (min)'); zlabel('Concentration of free antibody(uM)'); title('Concentration of Free Antibody 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; %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); if (i==1) %Neumann BC's at center of spheroid Abplus = Y(i+1); Abminus = (4/3)*Ab - (1/3)*Abplus; %quadratic approximation and no flux else if (i==Param.num_pts) %Dirichelt BC's for outer edge Abplus = Ab_exterior; %define outer concentration Abminus = Y(i-1); else %interior points Abplus = Y(i+1); Abminus = Y(i-1); end end %compute derivatives %Antibody dY_dt(i) = D.*(((Abminus - 2*Ab + Abplus)./(dR^2)) + (Abplus - Abminus)./(dR*i*dR)) - 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. function Ab_conc = calc_conc(Ab_initial,time,Param) Co_plasma = Ab_initial; %rename parameters t = time; k = Param.pharma_k; Co = Param.epsilon.*Co_plasma; %convert plasma conc to tumor conc Ab_conc = Co*exp(-k*t); %calculates exponential decay %DO NOT USE Co_plasma or Ab_initial (have not been corrected for void fraction! Ab_conc = Co; %use if constant concentration return;