image thumbnail

Monte Carlo Markov Chain for inferring parameters for an Ordinary Differential Equation model

by

 

This function uses a Monte Carlo Markov Chain algorithm to infer parameters for an ODE model

mcmc_ode

Contents

function mcmc_ode(N,parallel,diagnostics,nofigure)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name - mcmc_ode
% Creation Date - 2nd March 2012
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
%
% Description -
%    This uses a hierarchical Bayesian model to estimate ODE parameters for
%    a target cell limited model with eclipse phase.
%    This has narrower proposal distribution and has narrower mu, capsigma,
%    sigma. The initial values for mu, capsigma, and sigma are also
%    narrower and based on the paper equations
%	 Has blockwise updating of parameters.
%    beta fixed
%    3 moderately noisy mice
%
%
% Parameters -
%               N - Number of MCMC Iterations
%               parallel -
%                           If 1, then run ODE solving in parallel cores
%                           If 0, then run ODE solving serially
%               diagnostics -
%                           If 1, then print diagnostics
%                           If 0, then no diagnostics
%               nofigure -
%                           If 1, then do not display figures
%                           If 0, then display figures
%
% Example -
%               mcmc_ode(10,0,0,0)
%
%
%
% License - BSD
%
% Change History -
%                   2nd Mar 2012 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Load the data

fid_fileptr1 = importdata('knockout_RNA_div_500_1.tex');
passed_time_end1 = fid_fileptr1.data(end,1);
v_data1 = fid_fileptr1.data(:,2)';
time_data1 = fid_fileptr1.data(:,1)';

fid_fileptr2 = importdata('knockout_RNA_div_500_2.tex');
passed_time_end2 = fid_fileptr2.data(end,1);
v_data2 = fid_fileptr2.data(:,2)';
time_data2 = fid_fileptr2.data(:,1)';

fid_fileptr3 = importdata('knockout_RNA_div_500_3.tex');
passed_time_end3 = fid_fileptr3.data(end,1);
v_data3 = fid_fileptr3.data(:,2)';
time_data3 = fid_fileptr3.data(:,1)';

Initialize ODE parameters

target = 2.3e5; infected1 = 0; infected2 = 0; gamma = 44.43; k = 4;
beta = log10(1.388e-4); % fixed

Gibbs sampling

tic;

if matlabpool('size') == 0 && parallel == 1
    matlabpool open
end
Error using mcmc_ode (line 68)
Not enough input arguments.

Initialize the Gibbs sampler

% Set hyper-parameters
a = 23; b = 0.01; % a and b set using test_gamma.m optimal value of sigma_powerminus2
% is around 0.3 with tight distribution (see histogram generated)
% CHANGE HERE if data changed
% Order V0, p, delta
eta = [ log10(10^3.5) log10(27) log10(4.1) ]';
psi = diag([ 0.01^2 0.01^2 0.01^2 ]); % variance
omega = diag([ 10 10 10 ]); % adjusted after testing values of cap_sigma_powerminus1
nu = 2;

% Set other parameters
% CHANGE HERE if data changed
m1 = 5; m2 = 5; m3 = 5; % number of measurements for each subject
n = 3; % number of subjects

Set initial values for the Metropolis-Hastings algorithm

propsigma_V0 = 0.001; % standard dev. of proposal distribution for V0

propsigma_p = 0.001;  % standard dev. of proposal distribution for p

propsigma_delta = 0.001; % standard dev. of proposal distribution for delta

% Set initial values for Gibbs sampler
seed = 1; rand('state',seed); randn('state',seed); % set the random seed

Initialize hyperparameters

% sigma_powerminus2 = gamrnd(a,b); % from the equation in paper
sigma_powerminus2 = gamrnd(a + m1 + m2 + m3, ...
                           1/( 1/b + (1/2 * 0) ) ...
                          )

% mu = mvnrnd(eta,psi); % from the equation in paper
% cap_sigma_powerminus1 = wishrnd(omega,nu); % from the equation in paper

% CHANGE HERE if data changed
% Initialize theta_1, theta_2 and theta_3
V0_1 = log10(100);  p_1 = log10(20); delta_1 = log10(4);
V0_2 = log10(100);  p_2 = log10(20); delta_2 = log10(4);
V0_3 = log10(100);  p_3 = log10(20); delta_3 = log10(4);

% CAUTION - mu initially initialized to theta_1
mu = [V0_1 p_1 delta_1]' % later to mean of theta_1, 2 and 3
D = inv(omega) + ...
    ([V0_1 p_1 delta_1]' - mu)*([V0_1 p_1 delta_1]' - mu)' + ...
    ([V0_2 p_2 delta_2]' - mu)*([V0_2 p_2 delta_2]' - mu)' + ...
    ([V0_3 p_3 delta_3]' - mu)*([V0_3 p_3 delta_3]' - mu)';

cap_sigma_powerminus1 = wishrnd(inv(D), ...
                                n + nu ...
                                   )

B = n*cap_sigma_powerminus1 + inv(psi);
C = cap_sigma_powerminus1 * ...
        ([V0_1 p_1 delta_1]' + [V0_2 p_2 delta_2]' + [V0_3 p_3 delta_3]') ...
            + psi\eta;
mu = mvnrnd(B\C, inv(B))'


% initialize state vectors and cell arrays that will hold the values
% of updated hyperparameters at each MCMC step
state_sigma_powerminus2 = zeros(N,1);
state_mu = cell(N,1); state_cap_sigma_powerminus1 = cell(N,1);
% CAUTION - change here if number of parameters inferred decreases
state_theta1 = zeros(N,3); state_theta2 = zeros(N,3); state_theta3 = zeros(N,3);

t = 1; % initialize iteration at 1
state_sigma_powerminus2(t,1) = sigma_powerminus2; % save the current state
state_mu{t} = mu;
state_cap_sigma_powerminus1{t} = cap_sigma_powerminus1;
state_theta1(t,:) = [V0_1 p_1 delta_1];
state_theta2(t,:) = [V0_2 p_2 delta_2];
state_theta3(t,:) = [V0_3 p_3 delta_3];

% for diagnostics
ssr_new_array = zeros(1,N); ssr_old_array = zeros(1,N);
ssr_new_array(1,1) = 0; ssr_old_array(1,1) = 0;

if diagnostics == 1
    temp_V0 = zeros(N); temp_newV0 = zeros(N); temp_pratio = zeros(N);
    temp_alpha = zeros(N); temp_u = zeros(N);
    retval_old_array = zeros(1,N); retval_new_array = zeros(1,N);
    temp_p = zeros(1,N); temp_newp = zeros(1,N);
    temp_delta = zeros(1,N); temp_newdelta = zeros(1,N);

    temp_V0(1) = 0; temp_newV0(1) = 0; temp_pratio(1) = 0;
    temp_alpha(1) = 0; temp_u(1) = 0;
    retval_old_array(1,1) = 0; retval_new_array(1,1) = 0;
    temp_p(1,1) = 0; temp_newp(1,1) = 0;
    temp_delta(1,1) = 0; temp_newdelta(1,1) = 0;
end

Start sampling

while t < N % Iterate until we have N-1 samples
    t = t + 1;

    % display progress
    if mod(t,100) == 0
        t
    end

    if diagnostics == 1
        t
    end

Gibbs sampling algorithm

Get new value of sigma^(-2)

    totalssr = 0;
    parfor iInnerCount = 1:n
        if iInnerCount == 1
            totalssr = totalssr + ...
                odecall_eclipse_tcl_jv_local_fileptr_mcmc(...
                        target,infected1,infected2,beta,...
                        [ V0_1 p_1 delta_1 ]',gamma,k,passed_time_end1,...
                        time_data1,fid_fileptr1,0,0,1,1,v_data1);
        elseif iInnerCount == 2
                totalssr = totalssr + ...
                odecall_eclipse_tcl_jv_local_fileptr_mcmc(...
                        target,infected1,infected2,beta,...
                        [ V0_2 p_2 delta_2 ]',gamma,k,passed_time_end2,...
                        time_data2,fid_fileptr2,0,0,1,1,v_data2);
        elseif iInnerCount == 3
                totalssr = totalssr + ...
                odecall_eclipse_tcl_jv_local_fileptr_mcmc(...
                        target,infected1,infected2,beta,...
                        [ V0_3 p_3 delta_3 ]',gamma,k,passed_time_end3,...
                        time_data3,fid_fileptr3,0,0,1,1,v_data3);
        end
    end
    sigma_powerminus2 = gamrnd(a + m1 + m2 + m3, ...
                               1/( 1/b + (1/2 * totalssr) ) ...
                              )

    % Get new value of mu

    B = n*cap_sigma_powerminus1 + inv(psi);
    C = cap_sigma_powerminus1 * ...
        ([V0_1 p_1 delta_1]' + [V0_2 p_2 delta_2]' + [V0_3 p_3 delta_3]') ...
            + psi\eta;

    mu = mvnrnd(B\C, inv(B))'

    % Get new value of cap_sigma_powerminus1
    D = inv(omega) + ...
    ([V0_1 p_1 delta_1]' - mu)*([V0_1 p_1 delta_1]' - mu)' + ...
    ([V0_2 p_2 delta_2]' - mu)*([V0_2 p_2 delta_2]' - mu)' + ...
    ([V0_3 p_3 delta_3]' - mu)*([V0_3 p_3 delta_3]' - mu)';

    cap_sigma_powerminus1 = wishrnd(inv(D), ...
                                    n + nu ...
                                   )

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Metropolis-Hastings algorithm

Get new value of theta_1, theta_2 & theta_3

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Get new value for theta_1

Propose a new value for V0 (for theta_1)

    prop_array = mvnrnd( [V0_1 p_1 delta_1], diag([propsigma_V0 propsigma_p propsigma_delta]) );
    new_V0_1    = prop_array(1,1);
    new_p_1     = prop_array(1,2);
    new_delta_1 = prop_array(1,3);
    % for diagnostics
    if diagnostics == 1
        temp_V0(t) = V0_1;
        temp_newV0(t) = new_V0_1;
        temp_p(t) = p_1;
        temp_newp(t) = new_p_1;
        temp_delta(t) = delta_1;
        temp_newdelta(t) = new_delta_1;
    end

    [retval_new ssr_new] = ...
         postdist([ new_V0_1 new_p_1 new_delta_1]',sigma_powerminus2,mu,...
                        cap_sigma_powerminus1,fid_fileptr1,...
                        target,infected1,infected2,beta,...
                        gamma,k,passed_time_end1,time_data1,v_data1,diagnostics...
                        );
    [retval_old ssr_old] = ...
         postdist([ V0_1     p_1     delta_1]',sigma_powerminus2,mu,...
                        cap_sigma_powerminus1,fid_fileptr1,...
                        target,infected1,infected2,beta,...
                        gamma,k,passed_time_end1,time_data1,v_data1,diagnostics...
                        );

    pratio = retval_new/retval_old;
    alpha = min( [ 1 pratio ] ); % Calculate the acceptance ratio
    u = rand; % Draw a uniform deviate from [ 0 1 ]

    % for diagnostics
    if diagnostics == 1
        temp_pratio(t) = pratio;
        temp_alpha(t) = alpha;
        temp_u(t) = u;
        retval_old_array(1,t) = retval_old;
        retval_new_array(1,t) = retval_new;
    end

    ssr_new_array(1,t) = ssr_new;
    ssr_old_array(1,t) = ssr_old;

    if u < alpha % Do we accept this proposal?
       V0_1    = new_V0_1; % proposal becomes new value
	   p_1	= new_p_1;
	   delta_1 = new_delta_1;
    end

Get new value for theta_2

Propose a new value for V0 (for theta_2)

    prop_array  = mvnrnd( [V0_2 p_2 delta_2], diag([propsigma_V0 propsigma_p propsigma_delta]) );
    new_V0_2    = prop_array(1,1);
    new_p_2     = prop_array(1,2);
    new_delta_2 = prop_array(1,3);
    pratio = ...
         postdist([ new_V0_2 new_p_2 new_delta_2]',sigma_powerminus2,mu,...
                        cap_sigma_powerminus1,fid_fileptr2,...
                        target,infected1,infected2,beta,...
                        gamma,k,passed_time_end2,time_data2,v_data2,diagnostics...
                        ) / ...
         postdist([     V0_2 p_2 delta_2]',sigma_powerminus2,mu,...
                        cap_sigma_powerminus1,fid_fileptr2,...
                        target,infected1,infected2,beta,...
                        gamma,k,passed_time_end2,time_data2,v_data2,diagnostics...
                        );
    alpha = min( [ 1 pratio ] ); % Calculate the acceptance ratio
    u = rand; % Draw a uniform deviate from [ 0 1 ]
    if u < alpha % Do we accept this proposal?
       V0_2 = new_V0_2; % proposal becomes new value
	   p_2	= new_p_2;
	   delta_2 = new_delta_2;
    end

Get new value for theta_3

Propose a new value for V0 (for theta_3)

    prop_array  = mvnrnd( [V0_3 p_3 delta_3], diag([propsigma_V0 propsigma_p propsigma_delta]) );
    new_V0_3    = prop_array(1,1);
    new_p_3     = prop_array(1,2);
    new_delta_3 = prop_array(1,3);


    [retval_new ssr_new] = ...
    postdist([ new_V0_3 new_p_3 new_delta_3]',sigma_powerminus2,mu,...
                    cap_sigma_powerminus1,fid_fileptr3,...
                    target,infected1,infected2,beta,...
                    gamma,k,passed_time_end3,time_data3,v_data3,diagnostics...
                    );
    [retval_old ssr_old] = ...
    postdist([     V0_3 p_3 delta_3]',sigma_powerminus2,mu,...
                    cap_sigma_powerminus1,fid_fileptr3,...
                    target,infected1,infected2,beta,...
                    gamma,k,passed_time_end3,time_data3,v_data3,diagnostics...
                    );
    pratio = retval_new/retval_old;
    alpha = min( [ 1 pratio ] ); % Calculate the acceptance ratio
    u = rand; % Draw a uniform deviate from [ 0 1 ]
    if u < alpha % Do we accept this proposal?
       V0_3 = new_V0_3; % proposal becomes new value
	   p_3	= new_p_3;
	   delta_3 = new_delta_3;
    end

Save state

    state_sigma_powerminus2(t,1) = sigma_powerminus2;
    state_mu{t} = mu;
    state_cap_sigma_powerminus1{t} = cap_sigma_powerminus1;
    state_theta1(t,:) = [V0_1 p_1 delta_1];
    state_theta2(t,:) = [V0_2 p_2 delta_2];
    state_theta3(t,:) = [V0_3 p_3 delta_3];
end

disp('Initial values:')
[state_theta1(1,1) state_theta1(1,2) state_theta1(1,3)]

disp('Final values:')
[state_theta1(end,1) state_theta1(end,2) state_theta1(end,3)]

if diagnostics == 1
    disp('state_cap_sigma_powerminus1')
    state_cap_sigma_powerminus1{1}
    state_cap_sigma_powerminus1{2}
    state_cap_sigma_powerminus1{3}
    state_cap_sigma_powerminus1{4}
    state_cap_sigma_powerminus1{5}
    state_cap_sigma_powerminus1{6}
    state_cap_sigma_powerminus1{7}
    state_cap_sigma_powerminus1{8}
    state_cap_sigma_powerminus1{9}
    state_cap_sigma_powerminus1{10}

    disp('state_sigma_powerminus2')
    state_sigma_powerminus2(:,1)

    disp('state_mu')
    state_mu{1}
    state_mu{2}
    state_mu{3}
    state_mu{4}
    state_mu{5}
    state_mu{6}
    state_mu{7}
    state_mu{8}
    state_mu{9}
    state_mu{10}
end

Plot parameters

[ST,~] = dbstack();
function_name = ST.name;

if nofigure ~= 1
    figID = figure;
    plot([1:N],state_theta1(:,1),'-k','linewidth',2)
    ylabel('log_1_0 V_0');
    xlabel('Number of MCMC Iterations');
    print(figID, '-djpeg', sprintf('V0_iter_%s_%d_%s.jpg', function_name, N, date));

    figID = figure;
    plot([1:N],state_theta1(:,2),'-r','linewidth',2)
    ylabel('log_1_0 p');
    xlabel('Number of MCMC Iterations');
    print(figID, '-djpeg', sprintf('p_iter_%s_%d_%s.jpg', function_name, N, date));

    figID = figure;
    plot([1:N],state_theta1(:,3),'-k','linewidth',2)
    ylabel('log_1_0 delta');
    xlabel('Number of MCMC Iterations');
    print(figID, '-djpeg', sprintf('delta_iter_%s_%d_%s.jpg', function_name, N, date));

    figID = figure;
    hist(state_theta1(:,1),1000)
    xlabel('log_1_0 V0');
    title('Histogram');
    print(figID, '-djpeg', sprintf('V0_hist_%s_%d_%s.jpg', function_name, N, date));

    figID = figure;
    hist(state_theta1(:,2),1000)
    xlabel('log_1_0 p');
    title('Histogram');
    print(figID, '-djpeg', sprintf('p_hist_%s_%d_%s.jpg', function_name, N, date));

    figID = figure;
    hist(state_theta1(:,3),1000)
    xlabel('log_1_0 delta');
    title('Histogram');
    print(figID, '-djpeg', sprintf('delta_hist_%s_%d_%s.jpg', function_name, N, date));
end

Save evolution of parameters in file

Fields - log10(V0),log10(p) ,log10(delta) all in log10 scale show and write results to file

fid_output = fopen(sprintf('results_evolution_%s_%d_%s.csv', function_name, N, date),'w');
if diagnostics == 1
    fprintf(fid_output,...
        'log10V0,log10p,log10delta,tempoldV0,tempnewV0,temp_pratio,temp_alpha,temp_u,ssr_new,ssr_old,retval_new_array,retval_old_array,tempoldp,tempnewp,tempolddelta,tempnewdelta\r\n');
else
    fprintf(fid_output,...
        'log10V0_1,log10p_1,log10delta_1,log10V0_2,log10p_2,log10delta_2,log10V0_3,log10p_3,log10delta_3,ssr_new,ssr_old,sigma_powerminus2,log10V0_pop,log10p_pop,log10delta_pop\r\n');
end

% output result array to file
for iCount=1:N
% for iCount = 1:4 for testing loop
    if diagnostics == 1
        fprintf(fid_output,'%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f\r\n',...
                        state_theta1(iCount,1),state_theta1(iCount,2),...
                        state_theta1(iCount,3),...
                        temp_V0(iCount),temp_newV0(iCount),temp_pratio(iCount),...
                        temp_alpha(iCount),temp_u(iCount),...
                        ssr_new_array(1,iCount),ssr_old_array(1,iCount),...
                        retval_new_array(1,iCount),retval_old_array(1,iCount),...
                        temp_p(1,iCount),temp_newp(1,iCount),...
                        temp_delta(1,iCount),temp_newdelta(1,iCount)...
                        );
    else
        fprintf(fid_output,'%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f,%12.10f\r\n',...
                        state_theta1(iCount,1),state_theta1(iCount,2),...
                        state_theta1(iCount,3),...
                        state_theta2(iCount,1),state_theta2(iCount,2),...
                        state_theta2(iCount,3),...
                        state_theta3(iCount,1),state_theta3(iCount,2),...
                        state_theta3(iCount,3),...
                        ssr_new_array(1,iCount),ssr_old_array(1,iCount),...
                        state_sigma_powerminus2(iCount,1),...
                        state_mu{iCount}(1),state_mu{iCount}(2),state_mu{iCount}(3)...
                        );
    end
end


% close output file
fclose(fid_output);

toc;

% Yeah, clean up after yourself
if matlabpool('size') > 0 && parallel == 1
    matlabpool close
end

Auxiliary functions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Posterior distribution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [retval ssr] = postdist(param_vector,sigma_powerminus2,mu,...
                        cap_sigma_powerminus1,fileptr,...
                        target,infected1,infected2,beta,...
                        gamma,k,time_phase,time_data,v_data,diagnostics...
                        )
ssr = odecall_eclipse_tcl_jv_local_fileptr_mcmc(...
                        target,infected1,infected2,beta,...
                        param_vector,gamma,k,time_phase,...
                        time_data,fileptr,0,0,1,1,v_data);
if diagnostics == 1
    disp('ssr')
    ssr
    disp('whole term')
    (-sigma_powerminus2/2)*ssr ...
        - ( (param_vector - mu)' * cap_sigma_powerminus1 * (param_vector - mu) )/2
    disp('first term')
    (-sigma_powerminus2/2)*ssr
    disp('second term')
    - ( (param_vector - mu)' * cap_sigma_powerminus1 * (param_vector - mu) )/2
end

retval = exp( (-sigma_powerminus2/2)*ssr ...
                - ( (param_vector - mu)' * cap_sigma_powerminus1 * (param_vector - mu) )/2 ...
            );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ODE solvers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ssr = odecall_eclipse_tcl_jv_local_fileptr_mcmc(target,infected1,...
                            infected2,beta,param_vector,gamma,k,...
                            time_phase,time_vector,fileptr,plotfig,...
                            figorssr,mac,interp,emp_virus_vector)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to call ODE solver and plot results and calculate SSR
%
% Name - odecall_eclipse_tcl_jv_local_fileptr_mcmc
% Creation Date - 9th Dec 2011
% Author - Soumya Banerjee
% Website - www.cs.unm.edu/~soumya
%
%
% Description - function to call ODE solver with eclipse phase, target cell
%                   limited model, and plot results and
%                   calculate SSR, for MCMC
%
% Parameters -
%               target - initial target cell density (/mL)
%               infected1 - initial latently infected cell density (/mL)
%               infected2 - initial productively infected cell density
%               (/mL)
%               beta - infectivity log10
%               param_vector(1,4) has
%                   virus - initial virus density (/mL) in log 10
%
%                   p - replication rate (PFU/day) in log 10
%                   delta - infected cell death rate (/day) in log 10
%               gamma - immune system clearance rate (/day)
%               k - eclipse phase (rate of transition from I1 to I2) (/day)
%               time_phase - duration of simulation(days)
%               time_vector - vector of measured times (days)
%               fileptr - handle to 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
%               emp_virus_vector - empirically measured virus data (1,n)
%                                   row vector
%
%
%
%
% 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) Target cell limited model with eclipse phase
%               6) File pointer handle provided
%               7) Returns SSR on logged data and simulation
%
% 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 -
%                   29th Aug 2011 - Creation by Soumya Banerjee
%                    7th Sep 2011 - Modification by Soumya Banerjee
%                                   plots now saved and parameters
%                                   printed in title
%                    7th Nov 2011 - Added testing lines to print
%                                   parameters
%                    9th Nov 2011 - Added code to handle infinite, complex
%                                   numbers and NaN
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

high_virus_value = 20; % a very high value of virus (in log10) that is
% returned when there is a problem with integration

virus = param_vector(1,1);
p = param_vector(2,1);     delta = param_vector(3,1);

options = odeset('RelTol', .00001, ...
                 'NonNegative', [1 2 3 4]);

[t1 z1] = ode15s(@odecore_jv_tcl_eclipse_inline,[0, time_phase], ...
                                            [target ...
                                            infected1 ...
                                            infected2 ...
                                            double(10^virus) ...
                                            double(10^beta) ...
                                            double(10^p) ...
                                            gamma ...
                                            double(10^delta) ...
                                            k],...
                                            options);
% FOR TESTING
% [virus beta p delta]

% 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(:,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);
    figID = figure
    plot(t1,V,'-b','linewidth',2)
    xlabel('Time post infection (in days)','fontsize',25); 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])
    %     title(strcat(strcat(strcat(strcat(strcat(strcat(strcat('V0 = ',num2str(10^virus)),...
    %         strcat('beta = ',num2str(10^beta))), 'p = '),...
    %         num2str(10^p))),'delta = '),num2str(10^delta)),'fontsize',18);

    title(strcat(strcat('T0 =',num2str(target)),strcat(strcat(strcat(strcat(strcat('k ='...
        ,num2str(k)),strcat(strcat(strcat(strcat('V0 = ',num2str(10^virus)),...
       strcat('beta = ',num2str(10^beta))), 'p = '),...
       num2str(10^p)))),'delta = '),num2str(10^delta))),'fontsize',18);

    legend('simulated','actual data','Location','NorthEast')
    print(figID, '-djpeg', sprintf('tcl_jv_plot_T0%s_k%s.jpg', num2str(target),num2str(k)));
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

    sim_virus_vector = interp1(t1, V, time_vector);
%     V

    if ~isempty(find((isnan(sim_virus_vector) == 1))) ...
            || ~isempty(find(isfinite(sim_virus_vector) == 0)) ...
            || ~isempty(find(isreal(sim_virus_vector) == 0))
        % there is a NaN in sim_virus_vector or complex number or infinite
        % due to integration problem return something that gives high SSR
        sim_virus_vector = high_virus_value*ones(1,size(time_vector,2));
    end
else
    sim_virus_vector = [];
end

ssr = sum((sim_virus_vector - emp_virus_vector).^2);


function dydt = odecore_jv_tcl_eclipse_inline(~,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to solve ODEs
%
% Name - odecore_jv_tcl_eclipse_inline
% Creation Date - 27th Aug 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) Target cell limited model with correction term
%               3) Eclipse phase included
%               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
%
% Change History -
%                   27th Aug 2011 - Creation by Soumya Banerjee
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    dydt = y;

%     % Initial conditions
%     T  = y(1); % Target cells
%     I1 = y(2); % Infected cells not producing virus (latently infected)
%     I2 = y(3); % Infected cells producing viru (productively infected)
%     V  = y(4); % Virus

%     % Values of model parameters passed in y
%     beta   = y(5);
%     p      = y(6);
%     gamma  = y(7);
%     delta  = y(8);
%     k      = y(9); % eclipse phase (rate of transition from I1 to I2)

    % System of differential equations
	dydt(1) = -y(5)*y(1)*y(4);
	dydt(2) = y(5)*y(1)*y(4) - y(9)*y(2);
    dydt(3) = y(9)*y(2) - y(8)*y(3);
	dydt(4) = y(6)*y(3) - y(7)*y(4) - y(5)*y(1)*y(4);

    % 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;
    dydt(9) = 0;

Contact us