Code covered by the BSD License  

Highlights from
A Graphical User Interface for Solving Delay Differential Equations and doing Local Search

image thumbnail

A Graphical User Interface for Solving Delay Differential Equations and doing Local Search

by

 

A GUI for solving a set of delay differential equations (DDEs) and local search for best solutions

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

Contact us