% Tumor Growth SIMIODE Activity Demonstration

% Demonstrating the one- and two-compartment models for tumor evolution.
% This code is set up to estimate the parameters for both the one- and
% two-compartment models and plot the final fits vs the given data set. It
% will also print out the SSE for both models for easy comparison.

% Estimated parameters can also be overridden to let students experiment
% with their own choices of parameters.

% NOTE: When choosing which data to load, consider the following:
% 'Part1_DataSet.mat' contains the data from Part 1 of the activity.
    % This data set should be used to demonstrate that both models can perform
    % well in the case where there is very little necrotic tissue in the tumor.
% 'Part2_DataSet.mat' contains the data from Part 2 of the activity with a
    % large necrotic core. This data set should be used to demonstrate that the
    % double ODE model performs far better than the single ODE model in
    % scenarios where the necrotic core is large.

clear all
close all

global data

% Load the data set
load('Part2_DataSet.mat')


% FIT SINGLE ODE MODEL TO DATA
initGuess = [.5 2];
lb = [0 0];
ub = [2 5];
opt = optimset('Display','off');
[params, mse] = fmincon(@sse_ocm, initGuess, [], [], [], [], lb, ub, [], opt); %this is a frequentist optimizer
fprintf('The SSE for the Single ODE Model is %f\n', mse);


% Simulate trajectory of single ODE model given optimized parameter values. 
% NOTE: Alternatively, you could change the values of the "params" vector
% here in order to let students see how changing parameter values affects
% the model fit.
v0 = .02; %Initial condition
[t1,vol1] = ode23(@(t,v)oneCompModel(t,v, params), 0:1:45, v0);




% FIT DOUBLE ODE MODEL TO DATA
clear initGuess lb ub params ss0
initGuess = [.5 2 .1 .5];
lb = [0 0 0 0];
ub = [2 5 1 1];
[params, ss0] = fmincon(@sse_tcm, initGuess, [], [], [], [], lb, ub, [], opt);
fprintf('The SSE for the Double ODE Model is %f\n', ss0);


% Simulate trajectory of double ODE model given optimized parameter values. 
% NOTE: Alternatively, you could change the values of the "params" vector
% here in order to let students see how changing parameter values affects
% the model fit.
v0 = [.02 0]; %All of tumor is viable to start
[t2,vol2] = ode23(@(t,v)twoCompModel(t,v, params), 0:1:45, v0);



% PLOT TOTAL TUMOR VOLUME DATA AGAINST BOTH MODEL FITS TO ASSESS WHICH IS BETTER
figure(1)
plot(data.day, data.totalVolume, 'ok',  'MarkerSize', 10, 'MarkerFaceColor','k')
hold on
plot(t1,vol1,'--r','Linewidth',2)
plot(t2,vol2(:,1)+vol2(:,2), '-b','Linewidth',2)
xlabel('Day')
ylabel('Volume')
legend('Tumor Volume Data', 'Single ODE Model, V(t)', 'Double ODE Model, V(t)+N(t)','Location','SouthEast')
title('Comparing Model Fits')
set(gca,'FontSize',18)



% HELPER FUNCTIONS FOR SOLVING ODE MODELS AND COMPUTING SSE
function dval = oneCompModel(t,V,params)
    
    dval(1,1) =  params(1)*V*(1-V/params(2));%-params(3)*V; %Include "-params(3)*V" piece if you want a death mechanism for viable cells

end


function dval = twoCompModel(t,val,params)

    V = val(1);
    N = val(2);
    
    dval(1,1) =  params(1)*V*(1-V/params(2))-params(3)*V;
    dval(2,1) =  params(3)*V-params(4)*N;  %The "-params(4)*N" could be deleted if you don't want a decay mechanism for necrotic cells
end



function sse = sse_ocm(params)
global data

v0 = .02; 
[t,vol] = ode23(@(t,v)oneCompModel(t,v, params), 0:1:45, v0);

sse = (sum((data.totalVolume - interp1(t, vol(:,1),  data.day)).^2))/length(data.totalVolume);

end


function sse = sse_tcm(params)
global data
v0 = [.02 0]; %initial conditions (assume all of tumor is viable to begin)
[t,vol] = ode23(@(t,v)twoCompModel(t,v, params), 0:1:45, v0);

totalVol = vol(:,1)+vol(:,2);

sse = sum((data.totalVolume - interp1(t, totalVol,  data.day)).^2);

end