Code covered by the BSD License  

Highlights from
Graphical User Interface for Solving Ordinary Differential Equations

image thumbnail

Graphical User Interface for Solving Ordinary Differential Equations

by

 

A Graphical User Interface for Solving Ordinary Differential Equations (Immune Response Model)

odecall_eclipse_tcl_jv_local(target,infected,...
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;
    

Contact us