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

ddecall_tcl(target,infected1,infected2,...
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;
    
    

Contact us