function virus_vector = ddecall_tcl(target,infected1,infected2,...
virus,beta,p,delta,gamma,tau,...
time_phase,time_vector,filename,plotfig,...
figorssr,mac,interp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to call DDE solver and plot results and calculate SSR
%
% Name - ddecall_tcl
% Creation Date - 19th Aug 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
%
% Description - Function to call DDE solver and plot results and
% calculate SSR
%
% Parameters -
% target - initial target cell density (/mL)
% infected1 - initial infected but not producing virus
% cell density (/mL)
% infected2 - initial infected and producing virus
% 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) in log 10
% gamma - immune system clearance rate (/day)
% tau - duration of fixed delay (day)
% 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
%
%
% Assumptions - 1) all parameters passed in numerical values i.e. not
% logged except virus, beta, p (which are logged)
% 2) TCL model with correction term with two infected
% cell compartments and fixed delay
% 3) File has two columns - first column has time in days
% and second column has virus load in PFU/mL in log10
% 4) Model parameters passed locally in intial conditons
% to ode solver
% 5) ddecore called inline (saving of 0.03 sec on 2.2GHz
% Intel Core i7, 4core, 8GB RAM, parallel matlab
%
% 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 -
% 19th Aug 2011 - Creation by Soumya Banerjee
% 21st Aug 2011 - time_vector parameter added
% 7th Dec 2011 - Made delta logged parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
high_virus_value = 20; % a very high value of virus (in log10) that is
% returned when there is a problem with integration
options = odeset('RelTol', .00001, ...
'NonNegative', [1 2 3 4]);
sol = dde23(@ddecore_tcl_inline,[tau],[target ...
infected1 ...
infected2 ...
double(10^virus) ...
double(10^beta) ...
double(10^p) ...
gamma ...
double(10^delta)],...
[0, time_phase],...
options);
t = sol.x;
% 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(sol.y(4,:));
% 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);
plot(t,V,'-b','linewidth',2)
xlabel('Time post infection (in days)','fontsize',18); ylabel('Virus concentration in blood (log_1_0 (PFU/mL))','fontsize',18);
% 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',10)
% axis([-0.1 8.5 0 3.5])
legend('simulated','actual data','Location','NorthEast')
end
% If interpolation needed, then do
if interp == 1
% Need some error handling here
% e.g. for log10virus log10beta log10p log10delta
% 2.5000 -3.9612 1.1276 0.6248, the returned virus vector is
% -1.0458 -2.7730 -4.4999 -6.2563 NaN
% T0 = 2.3e4, k = 3
virus_vector = interp1(t, V, time_vector);
if ~isempty(find((isnan(virus_vector) == 1))) ...
|| ~isempty(find(isfinite(virus_vector) == 0)) ...
|| ~isempty(find(isreal(virus_vector) == 0))
% there is a NaN in virus_vector or complex number or infinite
% due to integration problem return something that gives high SSR
virus_vector = high_virus_value*ones(1,size(time_vector,2));
end
else
virus_vector = [];
end
function dydt = ddecore_tcl_inline(~,y,Z)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to solve DDEs
%
% Name - ddecore_tcl_inline
% Creation Date - 19th Aug 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
% Acknowledgements - "Solving Delay Differential Equations with dde23"
% by L.F. Shampine and S. Thompson
%
% Description - core function to solve DDEs, to be called from dde23
%
% Parameters -
% t - time
% y - vector containing values for all compartments
% Z - vector containing values for all compartments at
% (current time - delay)
%
% Assumptions - 1) All parameters passed in numerical values i.e. not
% logged
% 2) Target cell limited model with correction term and two
% compartments of infected cells
% 3) First argument (t) replaced by ~. It might give some
% performance benefit
% 4) Inline version of "ddecore_tcl"
%
% Comments - Make this function inline to speed up computation as this
% is called repeatedly by dde23
%
% License - BSD
%
% Change History -
% 19th Aug 2011 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dydt = y;
% Initial conditions
T = y(1); % Target cells
I1 = y(2); % Infected but not producing virus compartment
I2 = y(3); % Infected and producing virus compartment
V = y(4); % Virus
% Values of model parameters passed in y
beta = y(5);
p = y(6);
gamma = y(7);
delta = y(8);
% Time delayed solutions stored in this column vector (for all
% compartments)
ylag = Z; % there is only one delay time, and so Z is a 1 column vector
% if two delay times, then reference appropriate column
% System of delay differential equations
dydt(1) = -beta*T*V;
dydt(2) = beta*T*V - beta*ylag(1)*ylag(4);
dydt(3) = beta*ylag(1)*ylag(4) - delta*I2;
dydt(4) = p*I2 - gamma*V - beta*T*V;
% These are the model parameters. Since they do not change with time,
% their rate of change is 0
dydt(5) = 0;
dydt(6) = 0;
dydt(7) = 0;
dydt(8) = 0;