function virus_vector = odecall_eclipse_tcl_jv_local(target,infected,...
virus,beta,p,delta,gamma,...
rho,eta,t_i,...
time_phase,time_vector,filename,plotfig,...
figorssr,mac,interp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to call ODE solver and plot results and calculate SSR
%
% Name - odecall_eclipse_tcl_jv_local
% Creation Date - 24th Dec 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
%
% Description - function to call ODE solver with eclipse phase, immune
% limited model, and plot results and
% calculate SSR
%
% Parameters -
% target - initial target cell density (/mL)
% infected1 - initial latently infected cell density (/mL)
% infected2 - initial productively infected cell density
% (/mL)
% virus - initial virus density (/mL) in log 10
% beta - infectivity log10
% p - replication rate (PFU/day) in log 10
% delta - infected cell death rate )/day)
% gamma - immune system clearance rate (/day)
% rho - efficacy of antibody neutralization (/PRNT50 /day)
% eta - rate of antibody production (/day)
% t_i - time at which adaptive antibody kicks in (days)
% time_phase - duration of simulation (days)
% time_vector - vector of measured times (days)
% filename - name of data file
% plotfig - 1 if data and model simulation needs to be
% plotted,
% 0 if no plot needed
% figorssr - 1 if there is a need to access the data file
% (needed when need to calculate SSR or
% need to plot the data in the figure),
% 0 if access to data file is not needed
% mac - 1 if this is my mac,
% 0 if my linux desktop
% interp - 1 if interpolation of simulation needed,
% 0 if not needed
%
%
% Example usage - virus_vector = odecall_eclipse_tcl_jv_local(2.3e5,0,...
% 0,4.50638,0.00219,log10(5.13),log10(2.36),44.43,4,...
% 10,[2 4 6 8 10],'knockout_RNA_div_500.tex',1,1,...
% 1,1)
%
%
% Assumptions - 1) all parameters passed in numerical values i.e. not
% logged except virus, beta, p and delta
% (which are logged to base 10)
% 2) Phase 1 model with correction term with one infected
% cell compartment
% 3) File has two columns - first column has time in days
% and second column has virus load in PFU/mL in log10
% (see attached file)
% 4) Model parameters passed locally in intial conditons
% to ode solver
% 5) Immune system limited model without eclipse phase
%
% Comments - 1) Make this function inline to speed up computation
% 2) Consider tinkering with ode options (using odeset)
% to speed up for your specific case
%
% License - BSD
%
% Change History -
% 24th Dec 2011 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = odeset('RelTol', .00001, ...
'NonNegative', [1 2 3]);
[t1 z1] = ode15s(@odecore_jv_tcl_noeclipse_inline,[0, time_phase], ...
[target ...
infected ...
double(10^virus) ...
double(10^beta) ...
double(10^p) ...
gamma ...
double(10^delta) ...
rho ...
eta ...
t_i],...
options);
% V is the fourth row of the output array
% CAUTION - change this if not true, e.g. in an ODE with one infected cell
% compartment
V = log10(z1(:,3)');
% If figure being plotted or ssr required, then import data
if figorssr == 1 || plotfig == 1
% If run on my mac then data file initial pathname is different
if mac == 1
fileptr = importdata(strcat('',filename), '\t');
else
% if run on my linux desktop data file initial pathname is
% different
fileptr = importdata(strcat('~/matlab/linux_stage98/',filename), '\t');
end
end
% Plot the results
if plotfig == 1
% subplot (2,2,1); plot(t,log10(zT))
% xlabel('Time (in days)'); ylabel('Count in log base 10'); title('Target cells progression');
%
% subplot (2,2,2); plot(t,zI)
% xlabel('Time (in days)'); ylabel('Count in log base 10'); title('Infected cells progression');
%
% subplot (2,2,3);
% figure
plot(t1,V,'-b','linewidth',2)
xlabel('Days post infection','fontsize',15); ylabel('Virus concentration(log_1_0 (PFU/mL))','fontsize',15);
% axis([-0.1 7])% 0 13])
hold on
% a = [7.8 9.8 10.3 10.3 8.4 1.8 0.85];
% subplot (2,2,3);
plot(fileptr.data(:,1),fileptr.data(:,2),'ro','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',8)
% axis([-0.1 8.5 0 3.5])
legend('simulated','data','Location','NorthEast')
end
% If interpolation needed, then do
if interp == 1
virus_vector = interp1(t1, V, time_vector);
else
virus_vector = [];
end
function dydt = odecore_jv_tcl_noeclipse_inline(t,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to solve ODEs
%
% Name - odecore_jv_tcl_noeclipse_inline
% Creation Date - 24th Dec 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
% Acknowledgements - Drew Levin
%
% Description - core function to solve ODEs, to be called from ode45, etc
% with eclipse phase for JV paper
%
% Assumptions - 1) All parameters passed in numerical values i.e. not
% logged
% 2) Immune limited model with correction term
% 3) No Eclipse phase
% 4) Inline version
% 5) Intermediate variables removed
%
% Comments - Make this function inline to speed up computation as this
% is called repeatedly by ode45 or ode15s
% - Input argument not used and replaced by ~
%
% License - BSD License
%
% Change History -
% 24th Dec 2011 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dydt = y;
% Initial conditions
T = y(1); % Target cells
I = y(2); % Infected cells
V = y(3); % Virus
% Values of model parameters passed in y
beta = y(4);
p = y(5);
gamma = y(6);
delta = y(7);
rho = y(8); % Ab efficacy
eta = y(9); % Ab production rate
t_i = y(10); % time at which Ab production kicks in
% System of differential equations
dydt(1) = -beta*T*V;
dydt(2) = beta*T*V - delta*I;
if t <= t_i
dydt(3) = p*I - gamma*V - beta*T*V;
else
dydt(3) = p*I - gamma*V - beta*T*V - rho*eta*(t - t_i)*V;
end
% These are the model parameters. Since they do not change with time,
% their rate of change is 0
dydt(4) = 0;
dydt(5) = 0;
dydt(6) = 0;
dydt(7) = 0;
dydt(8) = 0;
dydt(9) = 0;
dydt(10) = 0;