function output_cell_array = ...
dde_auto_fit_param(target,infected1,infected2,virus,beta,p, ...
delta,gamma,tau,time_phase,filename,parallel,...
curvefit,virus_lb,beta_lb,p_lb,delta_lb, ...
virus_ub,beta_ub,p_ub,delta_ub ...
)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to call delay differential equation (DDE) solver, plot results
% and use lsqcurvefit to get best solution
%
% Name - dde_auto_fit_param
% Creation Date - 7th Dec 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
%
% Description - Function to call delay differential equation (DDE) solver,
% plot results and use lsqcurvefit to get best solution
%
% 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)
% beta - infectivity
% p - replication rate (PFU/day)
% delta - infected cell death rate (/day)
% gamma - immune system clearance rate (/day)
% tau - eclipse phase (delay) in DDE (day)
% time_phase - duration of simulation (days)
% time_vector - vector of measured times (days)
% filename - name of data file
% parallel - If 1, then use multiple cores else not
% curvefit - If 1, then do curvefit else only plot initial
% guess
% virus_lb - lower bound on initial virus density (/mL)
% beta_lb - lower bound on infectivity
% p_lb - lower bound on replication rate (PFU/day)
% delta_lb - lower bound on infected cell death rate (/day)
% virus_ub - upper bound on initial virus density (/mL)
% beta_ub - upper bound on infectivity
% p_ub - upper bound on replication rate (PFU/day)
% delta_ub - upper bound on infected cell death rate (/day)
%
% Output - output_cell_array
% 1st cell is virus_vector - model simulated viraemia
% at experimental data points
% 2nd cell is coeff - optimized coefficients that
% minimize SSR (from lsqcurvefit)
%
%
% Assumptions - 1) all parameters passed in numerical values and are
% positive real numbers
%
% Comments -
%
% License - BSD license
%
% Change History -
% 7th Dec 2011 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize
coeff_init = [log10(virus),log10(beta),log10(p),log10(delta)]; % V0 logged, beta logged, p logged, delta logged
fid_filename = importdata(filename);
v_data = fid_filename.data(:,2)';
jv_time_vector = fid_filename.data(:,1)';
if curvefit ~= 1
%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial guess %
%%%%%%%%%%%%%%%%%%%%%%%%%
V0 = coeff_init(1); beta = coeff_init(2); p = coeff_init(3);
delta = coeff_init(4);
virus_vector = ddecall_tcl(target,infected1,infected2,V0,beta,p,...
delta,gamma,tau,...
time_phase,jv_time_vector,filename,1,1,1,1);
else
%%%%%%%%%%%%%%%%%%%%%%%%%
% Call curvefit %
%%%%%%%%%%%%%%%%%%%%%%%%%
% check to see if pool of workers is open, if not open
if matlabpool('size') == 0 && parallel == 1
matlabpool open
end
tic;
% Get lower and upper bounds on parameters
coeff_lb = [log10(virus_lb),log10(beta_lb),log10(p_lb),log10(delta_lb)]; % V0 logged, beta logged, p logged, delta logged
coeff_ub = [log10(virus_ub),log10(beta_ub),log10(p_ub),log10(delta_ub)]; % V0 logged, beta logged, p logged, delta logged
options = optimset('UseParallel','always');
ode_call = @(coeff_array,t) ddecall_tcl(target,infected1,infected2,...
coeff_array(1),coeff_array(2),...
coeff_array(3),coeff_array(4),...
gamma,tau,time_phase,...
jv_time_vector,filename,0,1,1,1);
coeff = lsqcurvefit(ode_call,coeff_init,[1:time_phase],v_data,coeff_lb,...
coeff_ub,options);
% plot the best solution
V0 = coeff(1); beta = coeff(2); p = coeff(3); delta = coeff(4);
virus_vector = ddecall_tcl(target,infected1,infected2,V0,beta,p,...
delta,gamma,tau,...
time_phase,jv_time_vector,filename,1,1,1,1);
% check to see if pool of workers is open, if not open
if matlabpool('size') > 0 && parallel == 1
matlabpool close
end
toc;
end
output_cell_array{1} = virus_vector;
if curvefit == 1
output_cell_array{2} = coeff;
else
output_cell_array{2} = [];
end