% Model of prostate cancer dynamics under intermittent 
% androgen suppression therapy (adapted from [Liu2015])
% Includes IPA estimators
% Includes normalized cost function

% Author: Julia L. Fleck
% Date: April 9, 2015


clear all;
clc;

% Reset the random number generator so that rand produces the same random
% numbers as if we restarted Matlab
rng('default');


% Model parameters
% Patient 15
% alpha1 = 0.0213;
% alpha2 = 0.0242;
% beta1 = 0.01;
% beta2 = 0.0168;
% k1 = 7*10^-9;
% k2 = 1.0;
% k3 = 7.4*10^-9;
% k4 = 2.0;
% m1 = 0.00005;
% x30 = 8.0;
% tau = 12.5;
% lambda1_ON = 0.03;
% lambda1_OFF = 0.01;
% mu1 = 0.02;
% mu3 = 0.02;
% muTau = mu3 + (x30/tau);
% d = 1.0;

% Patient 1
alpha1 = 0.0204;
alpha2 = 0.0242;
beta1 = 0.0201;
beta2 = 0.0168;
k1 = 10.0*10^-9;
k2 = 1.0;
k3 = 10.0*10^-9;
k4 = 2.0;
m1_ON = 0.00005;
m1_OFF = 0.00003;
x30 = 12.0;
tau = 12.5;
lambda1_ON = 0.03;
lambda1_OFF = 0.01;
mu1_ON = 0.05;
mu1_OFF = 0.2;
mu3 = 0.02;
d = 1.0;

% Number of iterations
iterations = 10;

% Initial values
theta1 = 7.0; % Lower PSA threshold
theta2 = 15.0; % Upper PSA threshold

history = zeros(6,iterations); % matrix to store optimization results
                               % each column corresponds to one iteration
                               % first row is iteration
                               % second row is J
                               % third row is dJdtheta1
                               % fourth row is dJdtheta2
                               % fifth row new theta1
                               % sixth row new theta2


% Gradient-based optimization scheme
for iter = 1:iterations
    
    % Simulation parameters
    T = 2500.0; % Total simulation time in days
    clock = 0.0; % Global clock
    mode = 1; % Operational mode of system
              % mode = 1 means therapy in ON
              % mode = 0 means therapy is OFF
    flagLow = 0; % Flag to indicate if lower PSA threshold was reached
    flagHigh = 0; % Flag to indicate if upper PSA threshold was reached

    % Data structure
    HSCpop = 0.0; % Vector with HSC population size along simulation time
    CRCpop = 0.0; % Vector with CRC population size along simulation time
    androgen = 0.0; % Vector with androgen level along simulation time
    timeODE = 0.0; % Vector with simulation time of each run of ODEs
    dx1dtheta1_ON = 0.0; % Vector with x'1,1(t) discretized over all ON periods
    dx1dtheta2_ON = 0.0; % Vector with x'1,2(t) discretized over all ON periods
    dx2dtheta1_ON = 0.0; % Vector with x'2,1(t) discretized over all ON periods
    dx2dtheta2_ON = 0.0; % Vector with x'2,2(t) discretized over all ON periods
    dx1dtheta1_OFF = 0.0; % Vector with x'1,1(t) discretized over all OFF periods
    dx1dtheta2_OFF = 0.0; % Vector with x'1,2(t) discretized over all OFF periods
    dx2dtheta1_OFF = 0.0; % Vector with x'2,1(t) discretized over all OFF periods
    dx2dtheta2_OFF = 0.0; % Vector with x'2,2(t) discretized over all OFF periods
    ON_period = zeros(4,1000); % Matrix with start and end times of ON periods
                                % Each column is one ON period
                                % First row is start time
                                % Second row is end time
                                % Third row is dtau_dtheta1 at end of ON period
                                % Fourth row is dtau_dtheta2 at end of ON period
    OFF_period = zeros(4,1000); % Matrix with start and end times of OFF periods
                                % Each column is one OFF period
                                % First row is start time
                                % Second row is end time
                                % Third row is dtau_dtheta1 at end of OFF period
                                % Fourth row is dtau_dtheta2 at end of OFF period
    eventDynamics = zeros(11,1000); % Matrix to store IPA info
                                    % First row is time of event occurrence
                                    % Second row is type of event
                                    % Third row is dtaudtheta1
                                    % Fourth row is dtaudtheta2
                                    % Fifth row is dx1dtheta1
                                    % Sixth row is dx1dtheta2
                                    % Seventh row is dx2dtheta1
                                    % Eighth row is dx2dtheta2
                                    % Ninth row is size of HSC population
                                    % Tenth row is size of CRC population
                                    % Eleventh row is androgen level

%     % Initial conditions for Patient 15
%     HSCinit = 32.3; % Initial population size of HSC cells
%     CRCinit = 0.1; % Initial population size of CRC cells
%     PSAinit = HSCinit + CRCinit; % Initial PSA level
%     androgenInit = 16.0; % Initial androgen level

    % Initial conditions for Patient 1
    PSAinit = 29.0; % Initial PSA level
    HSCinit = 28.9; % Initial population size of HSC cells
    CRCinit = 0.1; % Initial population size of CRC cells
    androgenInit = 16.0; % Initial androgen level
    
    index = 1;
    eventCounter = 1;
    eventDynamics(9,1) = HSCinit;
    eventDynamics(10,1) = CRCinit;
    eventDynamics(11,1) = androgenInit;
    
    % Discretized time interval for solving ODEs
%     tInterval = 0:0.5:500; % Patient 15
    tInterval = 0:0.5:350; % Patient 1
    
    % Discretized time interval for estimating population dynamics
    tEstimate = 0:0.5:10;
    
    
    % Simulation
    while clock(end) < T

        % ON state
        if mode == 1
            % Solve system of ODEs
            [tx,x] = ode45(@prostateON, tInterval, [HSCinit;CRCinit;androgenInit]);

            % Store start time of this ON period
            ON_period(1,eventCounter) = clock(end)+tx(1);

            % Determine if and when lower PSA threshold was crossed
            for i = 1:size(tx,1)

                if x(i,1) + x(i,2) <= theta1

                    % Transition from q_ON to q_OFF
                    % Event e1 occurred
                    eventCounter = eventCounter+1;

                    % Duplicates existing column of eventDynamics
                    eventDynamics(:,eventCounter) = eventDynamics(:,eventCounter-1);

                    % Store ODE results
                    HSCpop = [HSCpop; x(1:i,1)];
                    CRCpop = [CRCpop; x(1:i,2)];
                    androgen = [androgen; x(1:i,3)];
                    timeODE = [timeODE; tx(1:i)];
                    eventDynamics(9,eventCounter) = x(i,1);
                    eventDynamics(10,eventCounter) = x(i,2);
                    eventDynamics(11,eventCounter) = x(i,3);

                    % Store event time information
                    eventDynamics(1,eventCounter) = clock(end)+tx(i);
                    eventDynamics(2,eventCounter) = 1;

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Evolution of state derivatives over [clock(end), clock(end)+tx(i))
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                    tInit = 0; % Start of integration interval
                    tFinal = tx(i); % End of integration interval
                    tStep = 0.0001; % Time step for discretization
                    tVec = tInit:tStep:tFinal; % Discretized time interval
                    x3_init = max(x(1,3),0); % Value of x3 at start of integration interval
                    dx11_init = eventDynamics(5,eventCounter); % Value of dx1dtheta1 at start of integration interval
                    dx12_init = eventDynamics(6,eventCounter); % Value of dx1dtheta2 at start of integration interval
                    dx21_init = eventDynamics(7,eventCounter); % Value of dx2dtheta1 at start of integration interval
                    dx22_init = eventDynamics(8,eventCounter); % Value of dx2dtheta2 at start of integration interval
                    dx11_ON = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx1dtheta1 in ON state
                    dx12_ON = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx1dtheta2 in ON state
                    dx21_ON = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx2dtheta1 in ON state
                    dx22_ON = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx2dtheta2 in ON state
                    x1_event = max(x(i,1),0); % Value of x1 at time of event occurrence (end of interval)
                    x2_event = max(x(i,2),0); % Value of x2 at time of event occurrence (end of interval)
                    x3_event = max(x(i,3),0); % Value of x3 at time of event occurrence (end of interval)

                    % State variable x1
                    % First term in A(t)
                    f1 = @(t) alpha1./(1+(exp(-(x3_init*exp(-(t-tInit)/tau)+tau*mu3*((1-exp(-(t-tInit)/tau)))-k1)*k2)));
                    Q1 = integral(f1,tInit,tFinal);

                    % Second term in A(t)
                    f2 = @(t) beta1./(1+(exp(-(x3_init*exp(-(t-tInit)/tau)+tau*mu3*((1-exp(-(t-tInit)/tau)))-k3)*k4)));
                    Q2 = integral(f2,tInit,tFinal);

                    % Third term in A(t)
                    f3 = @(t) (m1_ON/x30)*(x3_init*exp(-(t-tInit)/tau)+(tau*mu3*(1-exp(-(t-tInit)/tau))));
                    % f3 = @(t) (m1/x30)*(x3_init*exp(-(t-tInit)/tau)+(tau*mu3*(1-exp(-(t-tInit)/tau))));
                    Q3 = integral(f3,tInit,tFinal);

                    % Fourth term in A(t)
                    Q4 = (m1_ON+lambda1_ON)*(tFinal-tInit);
                    % Q4 = (m1+lambda1_ON)*(tFinal-tInit);

                    % A(t)
                    A = Q1 - Q2 - Q3 - Q4;

                    % Value of dx1dtheta1 at end of integration interval
                    dx11_final_ON = dx11_init*exp(A);

                    % Value of dx1dtheta2 at end of integration interval
                    dx12_final_ON = dx12_init*exp(A);

                    % Trajectory of x'1(t) over [tInit, tFinal)
                    Aloop = 0;
                    for j = 1:length(tVec)-1
                        delta_t = tVec(j+1) - tVec(j);
                        hON = x3_init*exp(-(tVec(j+1) - tVec(1))/tau) + mu3*tau*(1-exp(-((tVec(j+1) - tVec(1))/tau)));
                        R1 = alpha1./(1+exp(-(hON-k1)*k2));
                        R2 = beta1./(1+exp(-(hON-k3)*k4));
                        R3 = (m1_ON/x30)*hON;
                        % R3 = (m1/x30)*hON;
                        % R4 = (m1+lambda1_ON);
                        R4 = (m1_ON+lambda1_ON);
                        Aloop = Aloop + R1*delta_t - R2*delta_t - R3*delta_t - R4*delta_t;
                        dx11_ON(j) = dx11_init*exp(Aloop);
                        dx12_ON(j) = dx12_init*exp(Aloop);
                    end

                    % Stores discretized state derivative
                    dx1dtheta1_ON = [dx1dtheta1_ON; dx11_ON];
                    dx1dtheta2_ON = [dx1dtheta2_ON; dx12_ON];

    %                 % Plots trajectory of x'1(t) over [tInit, tFinal)
    %                 plot(tVec(1:end-1),dx11_ON);
    %                 plot(tVec(1:end-1),dx12_ON);

                    % State variable x2
                    % B1(t)
                    h1 = @(t) (alpha2*(1-((d/x30)*(x3_init*exp(-(t-tInit)/tau)+(tau*mu3*(1-exp(-(t-tInit)/tau)))))))-beta2;
                    B1 = integral(h1,tInit,tFinal);

                    % B2(t)
                    % Must use discretized approximation because of the nested integral
                    B1loop = 0;
                    B21 = 0;
                    B22 = 0;
                    for j = 1:length(tVec)-1
                        delta_t = tVec(j+1) - tVec(j);
                        hON = x3_init*exp(-(tVec(j+1) - tVec(1))/tau) + mu3*tau*(1-exp(-((tVec(j+1) - tVec(1))/tau)));
                        B1loop = B1loop + ((alpha2*(1-((d/x30)*hON))) - beta2)*delta_t;
                        % Y1 = m1*(1-(hON/x30));
                        Y1 = m1_ON*(1-(hON/x30));
                        Y21 = dx11_ON(j);
                        Y22 = dx12_ON(j);
                        Y3 = exp(-B1loop);
                        Y4 = exp(B1loop);
                        B21 = B21 + Y4*Y1*Y21*Y3*delta_t;
                        B22 = B22 + Y4*Y1*Y22*Y3*delta_t;
                        dx21_ON(j) = (dx21_init*exp(B1loop)) + B21;
                        dx22_ON(j) = (dx22_init*exp(B1loop)) + B22;
                    end

                    % Value of dx2dtheta1 at end of integration interval
                    dx21_final_ON = (dx21_init*exp(B1)) + B21;

                    % Value of dx2dtheta2 at end of integration interval
                    dx22_final_ON = (dx22_init*exp(B1)) + B22;

                    % Stores discretized state derivative
                    dx2dtheta1_ON = [dx2dtheta1_ON; dx21_ON];
                    dx2dtheta2_ON = [dx2dtheta2_ON; dx22_ON];

    %                 % Plots trajectory of x'2(t) over [tInit, tFinal)
    %                 plot(tVec(1:end-1),dx21_ON);
    %                 plot(tVec(1:end-1),dx22_ON);

                    %%%%%%%%%%%%%%%%%%%%%%%%
                    % Event time derivative
                    %%%%%%%%%%%%%%%%%%%%%%%%

                    % Event e1
                    % Transition from q_ON to q_OFF

                    % Estimate random processes at time of event occurrence
                    myrand_x1 = (0.001*randn(1)+0.0); % Sample from Gaussian white noise distribution
                    myrand_x2 = (0.001*randn(1)+0.0); % Sample from Gaussian white noise distribution

                    % Denominator of event time derivative
                    % mysum_e1 = (R1 - R2 - lambda1_ON)*x1_event + mu1 + ((alpha2*(1-((d/x30)*hON))) - beta2)*x2_event + myrand_x1 + myrand_x2;
                    mysum_e1 = (R1 - R2 - lambda1_ON)*x1_event + mu1_ON + ((alpha2*(1-((d/x30)*hON))) - beta2)*x2_event + myrand_x1 + myrand_x2;
                    
                    %mysum_e1 = abs((x(i,1)-x(i-2,1))/(tx(i)-tx(i-2))) + abs((x(i,2)-x(i-2,2))/(tx(i)-tx(i-2))); % Estimated using time interval immediately before event occurrence
                    %mysum_e1 = (beta1 + m1) + (alpha2 + m1);

                    % Event time derivatives
                    dtau_dtheta1 = (1 - dx11_final_ON - dx21_final_ON)/mysum_e1;
                    dtau_dtheta2 = (- dx12_final_ON - dx22_final_ON)/mysum_e1;

                    % Store event time derivative
                    eventDynamics(3,eventCounter) = dtau_dtheta1;
                    eventDynamics(4,eventCounter) = dtau_dtheta2;
                    ON_period(3,eventCounter-1) = dtau_dtheta1;
                    ON_period(4,eventCounter-1) = dtau_dtheta2;

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Update state derivatives
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%
                    
                    % Update initial conditions
                    HSCinit = x(i,1);
                    CRCinit = x(i,2);
                    androgenInit = x(i,3);
                    
                    % State variable x1
%                     firstTerm = R1 - (alpha1/(1+exp(-(x3_event-k1)*k2)));
%                     secondTerm = R2 - (beta1/(1+exp(-(x3_event-k3)*k4)));
%                     thirdTerm = R3 - ((m1/x30)*x3_event);
%                     delta_f1_e1 = (firstTerm - secondTerm + thirdTerm)*x1_event;
                    
                    % Estimate delta_f1_e1 within time interval around event occurrence
                    % Before event: cell population is decreasing
                    % After event: cell population is increasing
                    % Solve system of ODEs for OFF period to estimate population dynamics at event time
                    %[ty,y] = ode45(@prostateOFF, tEstimate, [HSCinit;CRCinit;androgenInit]);
%                     before_e1 = abs((x(i,1)-x(i-2,1))/(tx(i)-tx(i-2)));
%                     after_e1 = abs((y(3,1)-y(1,1))/(ty(3)-ty(1)));
%                     delta_f1_e1 = before_e1 - after_e1;

                    delta_f1_e1 = (beta1 + m1_ON) - (alpha1 + mu1_ON);
                    % delta_f1_e1 = (beta1 + m1) - (alpha1 + mu1);
                    dx11_ON_OFF = dx11_final_ON + (delta_f1_e1*dtau_dtheta1);
                    dx12_ON_OFF = dx12_final_ON + (delta_f1_e1*dtau_dtheta2);

                    
                    % State variable x2
%                     firstTerm = ((alpha2*d)/x30)*(x3_event - hON);
%                     secondTerm = (m1/x30)*(hON - x3_event);
%                     delta_f2_e1 = (firstTerm*x2_event) - (secondTerm*x1_event);

                    % Estimate delta_f2_e1 within time interval around event occurrence
                    % Before event: cell population is increasing
                    % After event: cell population is decreasing
%                     before_e1 = abs((x(i,2)-x(i-2,2))/(tx(i)-tx(i-2)));
%                     after_e1 = abs((y(3,2)-y(1,2))/(ty(3)-ty(1)));
%                     delta_f2_e1 = before_e1 - after_e1;

                    delta_f2_e1 = (alpha2 + m1_ON) - beta2;
                    % delta_f2_e1 = (alpha2 + m1) - beta2;
                    dx21_ON_OFF = dx21_final_ON + (delta_f2_e1*dtau_dtheta1);
                    dx22_ON_OFF = dx22_final_ON + (delta_f2_e1*dtau_dtheta2);

                    % Store state derivatives
                    eventDynamics(5,eventCounter) = dx11_ON_OFF;
                    eventDynamics(6,eventCounter) = dx12_ON_OFF;
                    eventDynamics(7,eventCounter) = dx21_ON_OFF;
                    eventDynamics(8,eventCounter) = dx22_ON_OFF;

                    % Perform state transition
                    mode = 0;

                    % Store end time of this ON period
                    ON_period(2,eventCounter-1) = clock(end)+tx(i);

                    % Update global clock
                    clock = [clock; clock(end)+tx(1:i)];

                    % Flag to indicate that lower threshold was reached
                    flagLow = 1;

                    % Exit the for loop
                    break;
                end
            end

            if flagLow == 0
                % Store ODE results
                HSCpop = [HSCpop; x(:,1)];
                CRCpop = [CRCpop; x(:,2)];
                androgen = [androgen; x(:,3)];
                timeODE = [timeODE; tx(1:end)];

                % Update initial conditions
                HSCinit = x(end,1);
                CRCinit = x(end,2);
                androgenInit = x(end,3);

                % No state transition (patient remains on therapy)
                % Repeat ON therapy cycle
                mode = 1;

                % Store end time of this ON period
                ON_period(2,eventCounter) = clock(end)+tx(end);
                eventCounter = eventCounter+1;

                % Updates global clock
                clock = [clock; clock(end)+tx(1:end)];
            end
        end

        % Simulates system in OFF state
        if mode == 0
            % Solves system of ODEs
            [ty,y] = ode45(@prostateOFF, tInterval, [HSCinit;CRCinit;androgenInit]);

            % Store start time of this OFF period
            OFF_period(1,eventCounter) = clock(end)+ty(1);

            % Determines if and when upper PSA threshold was crossed
            for i = 1:size(ty,1)

                if y(i,1) + y(i,2) >= theta2

                    % Transition from q_OFF to q_ON
                    % Event e2 occurred
                    eventCounter = eventCounter+1;

                    % Duplicates existing column of eventDynamics
                    eventDynamics(:,eventCounter) = eventDynamics(:,eventCounter-1);

                    % Store ODE results
                    HSCpop = [HSCpop; y(1:i,1)];
                    CRCpop = [CRCpop; y(1:i,2)];
                    androgen = [androgen; y(1:i,3)];
                    timeODE = [timeODE; ty(1:i)];
                    eventDynamics(9,eventCounter) = y(i,1);
                    eventDynamics(10,eventCounter) = y(i,2);
                    eventDynamics(11,eventCounter) = y(i,3);

                    % Store event time information
                    eventDynamics(1,eventCounter) = clock(end)+ty(i);
                    eventDynamics(2,eventCounter) = 2;

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Evolution of state derivatives over [clock(end), clock(end)+ty(i))
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                    tInit = 0; % Start of integration interval
                    tFinal = ty(i); % End of integration interval
                    tStep = 0.0001; % Time step for discretization
                    tVec = tInit:tStep:tFinal; % Discretized time interval
                    x3_init = max(y(1,3),0); % Value of x3 at start of integration interval
                    dx11_init = eventDynamics(5,eventCounter); % Value of dx1dtheta1 at start of integration interval
                    dx12_init = eventDynamics(6,eventCounter); % Value of dx1dtheta2 at start of integration interval
                    dx21_init = eventDynamics(7,eventCounter); % Value of dx2dtheta1 at start of integration interval
                    dx22_init = eventDynamics(8,eventCounter); % Value of dx2dtheta2 at start of integration interval
                    dx11_OFF = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx1dtheta1 in OFF state
                    dx12_OFF = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx1dtheta2 in OFF state
                    dx21_OFF = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx2dtheta1 in OFF state
                    dx22_OFF = zeros(length(tVec)-1, 1); % Vector to store trajectory of dx2dtheta2 in OFF state
                    x1_event = max(y(i,1),0); % Value of x1 at time of event occurrence (end of interval)
                    x2_event = max(y(i,2),0); % Value of x2 at time of event occurrence (end of interval)
                    x3_event = max(y(i,3),0); % Value of x3 at time of event occurrence (end of interval)

                    % State variable x1
                    % First term in C(t)
                    g1 = @(t) alpha1./(1+(exp(-(x3_init*exp(-(t-tInit)/tau)+((tau*mu3)+x30)*((1-exp(-(t-tInit)/tau)))-k1)*k2)));
                    U1 = integral(g1,tInit,tFinal);

                    % Second term in C(t)
                    g2 = @(t) beta1./(1+(exp(-(x3_init*exp(-(t-tInit)/tau)+((tau*mu3)+x30)*((1-exp(-(t-tInit)/tau)))-k3)*k4)));
                    U2 = integral(g2,tInit,tFinal);

                    % Third term in C(t)
                    g3 = @(t) (m1_OFF/x30)*(x3_init*exp(-(t-tInit)/tau)+(((tau*mu3)+x30)*(1-exp(-(t-tInit)/tau))));
                    % g3 = @(t) (m1/x30)*(x3_init*exp(-(t-tInit)/tau)+(((tau*mu3)+x30)*(1-exp(-(t-tInit)/tau))));
                    U3 = integral(g3,tInit,tFinal);

                    % Fourth term in C(t)
                    U4 = (m1_OFF+lambda1_OFF)*(tFinal-tInit);
                    % U4 = (m1+lambda1_OFF)*(tFinal-tInit);

                    % C(t)
                    C = U1 - U2 - U3 - U4;

                    % Value of dx1dtheta1 at end of integration interval
                    dx11_final_OFF = dx11_init*exp(C);

                    % Value of dx1dtheta2 at end of integration interval
                    dx12_final_OFF = dx12_init*exp(C);

                    % Trajectory of x'1(t) over [tInit, tFinal)
                    Cloop = 0;
                    for j = 1:length(tVec)-1
                        delta_t = tVec(j+1) - tVec(j);
                        hOFF = x3_init*exp(-(tVec(j+1) - tVec(1))/tau) + ((mu3*tau)+x30)*(1-exp(-((tVec(j+1) - tVec(1))/tau)));
                        S1 = alpha1./(1+exp(-(hOFF-k1)*k2));
                        S2 = beta1./(1+exp(-(hOFF-k3)*k4));
                        S3 = (m1_OFF/x30)*hOFF;
                        % S3 = (m1/x30)*hOFF;
                        % S4 = (m1+lambda1_OFF);
                        S4 = (m1_OFF+lambda1_OFF);
                        Cloop = Cloop + S1*delta_t - S2*delta_t - S3*delta_t - S4*delta_t;
                        dx11_OFF(j) = dx11_init*exp(Cloop);
                        dx12_OFF(j) = dx12_init*exp(Cloop);
                    end

                    % Stores discretized state derivative
                    dx1dtheta1_OFF = [dx1dtheta1_OFF; dx11_OFF];
                    dx1dtheta2_OFF = [dx1dtheta2_OFF; dx12_OFF];

    %                 % Plots trajectory of x'1(t) over [tInit, tFinal)
    %                 plot(tVec(1:end-1),dx11_OFF);
    %                 plot(tVec(1:end-1),dx12_OFF);

                    % State variable x2
                    % D1(t)
                    r1 = @(t) (alpha2*(1-((d/x30)*(x3_init*exp(-(t-tInit)/tau)+(((tau*mu3)+x30)*(1-exp(-(t-tInit)/tau)))))))-beta2;
                    D1 = integral(r1,tInit,tFinal);

                    % D2(t)
                    % Must use discretized approximation because of the nested integral
                    D1loop = 0;
                    D21 = 0;
                    D22 = 0;
                    for j = 1:length(tVec)-1
                        delta_t = tVec(j+1) - tVec(j);
                        hOFF = x3_init*exp(-(tVec(j+1) - tVec(1))/tau) + ((mu3*tau)+x30)*(1-exp(-((tVec(j+1) - tVec(1))/tau)));
                        D1loop = D1loop + ((alpha2*(1-((d/x30)*hOFF))) - beta2)*delta_t;
                        P1 = m1_OFF*(1-(hOFF/x30));
                        % P1 = m1*(1-(hOFF/x30));
                        P21 = dx11_OFF(j);
                        P22 = dx12_OFF(j);
                        P3 = exp(-D1loop);
                        P4 = exp(D1loop);
                        D21 = D21 + P4*P1*P21*P3*delta_t;
                        D22 = D22 + P4*P1*P22*P3*delta_t;
                        dx21_OFF(j) = (dx21_init*exp(D1loop)) + D21;
                        dx22_OFF(j) = (dx22_init*exp(D1loop)) + D22;
                    end

                    % Value of dx2dtheta1 at end of integration interval
                    dx21_final_OFF = (dx21_init*exp(D1)) + D21;

                    % Value of dx2dtheta2 at end of integration interval
                    dx22_final_OFF = (dx22_init*exp(D1)) + D22;

                    % Stores discretized state derivative
                    dx2dtheta1_OFF = [dx2dtheta1_OFF; dx21_OFF];
                    dx2dtheta2_OFF = [dx2dtheta2_OFF; dx22_OFF];

    %                 % Plots trajectory of x'2(t) over [tInit, tFinal)
    %                 plot(tVec(1:end-1),dx21_OFF);
    %                 plot(tVec(1:end-1),dx22_OFF);

                    %%%%%%%%%%%%%%%%%%%%%%%%
                    % Event time derivative
                    %%%%%%%%%%%%%%%%%%%%%%%%

                    % Event e2
                    % Transition from q_OFF to q_ON

                    % Estimate random processes at time of event occurrence
                    myrand_x1 = (0.001*randn(1)+0.0); % Sample from Gaussian white noise distribution
                    myrand_x2 = (0.001*randn(1)+0.0); % Sample from Gaussian white noise distribution

                    % Denominator of event time derivative
                    % mysum_e2 = (S1 - S2 - lambda1_OFF)*x1_event + mu1 + ((alpha2*(1-((d/x30)*hOFF))) - beta2)*x2_event + myrand_x1 + myrand_x2;
                    mysum_e2 = (S1 - S2 - lambda1_OFF)*x1_event + mu1_OFF + ((alpha2*(1-((d/x30)*hOFF))) - beta2)*x2_event + myrand_x1 + myrand_x2;
                    
                    %mysum_e2 = abs((y(i,1)-y(i-2,1))/(ty(i)-ty(i-2))) + abs((y(i,2)-y(i-2,2))/(ty(i)-ty(i-2))); % Estimated using time interval immediately before event occurrence
                    %mysum_e2 = (alpha1 + mu1) + beta2;

                    % Event time derivatives
                    dtau_dtheta1 = (- dx11_final_OFF - dx21_final_OFF)/mysum_e2;
                    dtau_dtheta2 = (1 - dx12_final_OFF - dx22_final_OFF)/mysum_e2;

                    % Store event time derivatives
                    eventDynamics(3,eventCounter) = dtau_dtheta1;
                    eventDynamics(4,eventCounter) = dtau_dtheta2;
                    OFF_period(3,eventCounter-1) = dtau_dtheta1;
                    OFF_period(4,eventCounter-1) = dtau_dtheta2;


                    %%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Update state derivatives
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%
                    
                    % Update initial conditions
                    HSCinit = y(i,1);
                    CRCinit = y(i,2);
                    androgenInit = y(i,3);

                    % State variable x1
%                     firstTerm = S1 - (alpha1/(1+exp(-(x3_event-k1)*k2)));
%                     secondTerm = S2 - (beta1/(1+exp(-(x3_event-k3)*k4)));
%                     thirdTerm = S3 - ((m1/x30)*x3_event);
%                     delta_f1_e2 = (firstTerm - secondTerm + thirdTerm)*x1_event;

                    % Solve system of ODEs to estimate population dynamics at event time
                    %[tx,x] = ode45(@prostateON, tEstimate, [HSCinit;CRCinit;androgenInit]);
%                     before_e2 = abs((y(i,1)-y(i-2,1))/(ty(i)-ty(i-2)));
%                     after_e2 = abs((x(3,1)-x(1,1))/(tx(3)-tx(1)));
%                     delta_f1_e2 = before_e2 - after_e2;

                    delta_f1_e2 = (alpha1 + mu1_OFF) - (beta1 + m1_OFF);
                    % delta_f1_e2 = (alpha1 + mu1) - (beta1 + m1);
                    dx11_OFF_ON = dx11_final_OFF + (delta_f1_e2*dtau_dtheta1);
                    dx12_OFF_ON = dx12_final_OFF + (delta_f1_e2*dtau_dtheta2);

                    
                    % State variable x2
%                     firstTerm = ((alpha2*d)/x30)*(x3_event - hOFF);
%                     secondTerm = (m1/x30)*(hOFF - x3_event);
%                     delta_f2_e2 = (firstTerm*x2_event) - (secondTerm*x1_event);

%                     before_e2 = abs((y(i,2)-y(i-2,2))/(ty(i)-ty(i-2)));
%                     after_e2 = abs((x(3,2)-x(1,2))/(tx(3)-tx(1)));
%                     delta_f2_e2 = before_e2 - after_e2;

                    delta_f2_e2 = beta2 - (alpha2 + m1_OFF);
                    % delta_f2_e2 = beta2 - (alpha2 + m1);
                    dx21_OFF_ON = dx21_final_OFF + (delta_f2_e2*dtau_dtheta1);
                    dx22_OFF_ON = dx22_final_OFF + (delta_f2_e2*dtau_dtheta2);

                    % Store state derivatives
                    eventDynamics(5,eventCounter) = dx11_OFF_ON;
                    eventDynamics(6,eventCounter) = dx12_OFF_ON;
                    eventDynamics(7,eventCounter) = dx21_OFF_ON;
                    eventDynamics(8,eventCounter) = dx22_OFF_ON;

                    % Perform state transition
                    mode = 1;

                    % Store end time of this OFF period
                    OFF_period(2,eventCounter-1) = clock(end)+ty(i);

                    % Updates global clock
                    clock = [clock; clock(end)+ty(1:i)];

                    % Flag to indicate that upper threshold was reached
                    flagHigh = 1;

                    % Exits the for loop
                    break;
                end
            end

            if flagHigh == 0 && clock(end) < T
                % Updates initial conditions
                HSCinit = y(end,1);
                CRCinit = y(end,2);
                androgenInit = y(end,3);

                % Stores ODE results
                HSCpop = [HSCpop; y(:,1)];
                CRCpop = [CRCpop; y(:,2)];
                androgen = [androgen; y(:,3)];
                timeODE = [timeODE; ty(1:end)];

                % No state transition (patient remains off therapy)
                % Repeat OFF therapy cycle
                mode = 0;
                %mode = 1;

                % Store end time of this OFF period
                OFF_period(2,eventCounter) = clock(end)+ty(end);
                eventCounter = eventCounter+1;

                % Updates global clock
                clock = [clock; clock(end)+ty(1:end)];
            end
        end
    end

    % Eliminate first row of duplicate zeros
    % Store only data until time T
    timeODE = timeODE(2:end);
    clock = clock(2:find(clock<T,1,'last'));
    HSCpop = HSCpop(2:length(clock)+1);
    CRCpop = CRCpop(2:length(clock)+1);
    androgen = androgen(2:length(clock)+1);
    ON_period = ON_period(:,1:eventCounter-1);
    OFF_period = OFF_period(:,1:eventCounter-1);
    eventDynamics = eventDynamics(:,1:eventCounter);

    % Parse population dynamics to substitute negative values by zero
    HSCpop(find(HSCpop < 0)) = 0.0;
    CRCpop(find(CRCpop < 0)) = 0.0;
    androgen(find(androgen < 0)) = 0.0;
    eventDynamics(9,find(eventDynamics(9,:) < 0)) = 0.0;
    eventDynamics(10,find(eventDynamics(10,:) < 0)) = 0.0;
    eventDynamics(11,find(eventDynamics(11,:) < 0)) = 0.0;

    %%%%%%%%%%%%%%%%%%
    % Cost over [0,T]
    %%%%%%%%%%%%%%%%%%

    % Cost metric defined as weighted sum of:
    % 1 - Average PSA level over [0,T]
    % 2 - Average duration of ON period over [0,T]

    % Penalty factor for each term of the cost function
    W1 = 0.9;
    W2 = 0.1;

    % PSA level
    delta_t = 0.0;
    mysum_upper = 0.0;
    mysum_lower = 0.0;
    for k = 1:length(clock)-1
        % Integral of sum of HSC and CRC populations over [0,T]
        % Approximate by Riemann sum of discretized population trajectories
        delta_t = clock(k+1) - clock(k);
        mysum_upper = mysum_upper + (HSCpop(k+1)+CRCpop(k+1))*delta_t;
        mysum_lower = mysum_lower + (HSCpop(k)+CRCpop(k))*delta_t;
    end
    mysum_PSA = (mysum_upper+mysum_lower)/2;

    % Duration of ON period
    mysum_ON = 0.0;
    for j = 1:size(ON_period,2)
        if (ON_period(1,j) > T)
            break;
        else
            if (ON_period(2,j) > T)
                mysum_ON = mysum_ON + (T - ON_period(1,j));
            else
                mysum_ON = mysum_ON + (ON_period(2,j) - ON_period(1,j));
            end
        end
    end

    % Total cost over [0,T]
    mysum_J = W1*(1/PSAinit)*mysum_PSA + W2*(1/T)*mysum_ON;
    cost = (1/T)*mysum_J;


    %%%%%%%%%%%%%%%%%%
    % Cost derivative
    %%%%%%%%%%%%%%%%%%

    % First term
    % Evaluate integral of x'1(t)+x'2(t) over [0,T]
    % Approximate integral by Riemann Sum
    
    % ON periods
    ONloop1_lower = 0;
    ONloop1_upper = 0;
    ONloop2_lower = 0;
    ONloop2_upper = 0;
    for i = 1:length(dx1dtheta1_ON)-1
        ONloop1_upper = ONloop1_upper + (dx1dtheta1_ON(i+1) + dx2dtheta1_ON(i+1))*tStep;
        ONloop1_lower = ONloop1_lower + (dx1dtheta1_ON(i) + dx2dtheta1_ON(i))*tStep;
    end

    ONloop1 = (ONloop1_lower + ONloop1_upper)/2; % dxdtheta1

    for i = 1:length(dx1dtheta2_ON)-1
        ONloop2_upper = ONloop2_upper + (dx1dtheta2_ON(i+1) + dx2dtheta2_ON(i+1))*tStep;
        ONloop2_lower = ONloop2_lower + (dx1dtheta2_ON(i) + dx2dtheta2_ON(i))*tStep;
    end

    ONloop2 = (ONloop2_lower + ONloop2_upper)/2; % dxdtheta2

    % OFF periods
    OFFloop1_lower = 0;
    OFFloop1_upper = 0;
    OFFloop2_lower = 0;
    OFFloop2_upper = 0;
    for i = 1:length(dx1dtheta1_OFF)-1
        OFFloop1_upper = OFFloop1_upper + (dx1dtheta1_OFF(i+1) + dx2dtheta1_OFF(i+1))*tStep;
        OFFloop1_lower = OFFloop1_lower + (dx1dtheta1_OFF(i) + dx2dtheta1_OFF(i))*tStep;
    end

    OFFloop1 = (OFFloop1_lower + OFFloop1_upper)/2; % dxdtheta1

    for i = 1:length(dx1dtheta2_OFF)-1
        OFFloop2_upper = OFFloop2_upper + (dx1dtheta2_OFF(i+1) + dx2dtheta2_OFF(i+1))*tStep;
        OFFloop2_lower = OFFloop2_lower + (dx1dtheta2_OFF(i) + dx2dtheta2_OFF(i))*tStep;
    end

    OFFloop2 = (OFFloop2_lower + OFFloop2_upper)/2; % dxdtheta2

    % First term in cost derivative
    dL1dtheta1 = (1/PSAinit)*(ONloop1 + OFFloop1);
    dL1dtheta2 = (1/PSAinit)*(ONloop2 + OFFloop2);

    % Second term
    % Evaluate duration of ON periods * difference of event time derivatives
    dL2dtheta1 = 0;
    dL2dtheta2 = 0;
    
    % First ON period is not coupled with previous OFF period, so it can be
    % computed outside of the loop
    dL2dtheta1 = dL2dtheta1 + ((ON_period(2,1) - ON_period(1,1))*ON_period(3,1));
    dL2dtheta2 = dL2dtheta2 + ((ON_period(2,1) - ON_period(1,1))*ON_period(4,1));
    
    for j = 2:size(ON_period,2)
        
        % Last term only exists if sample path ends with incomplete ON period
        if (ON_period(1,j) > T || ON_period(2,j) > T)
            dL2dtheta1 = dL2dtheta1 - ((T - ON_period(1,j))*OFF_period(3,j-1));
            dL2dtheta2 = dL2dtheta2 - ((T - ON_period(1,j))*OFF_period(4,j-1));    
        
        else
            dL2dtheta1 = dL2dtheta1 + ((ON_period(2,j) - ON_period(1,j))*(ON_period(3,j) - OFF_period(3,j-1)));
            dL2dtheta2 = dL2dtheta2 + ((ON_period(2,j) - ON_period(1,j))*(ON_period(4,j) - OFF_period(4,j-1)));
        end
    end

    % Cost derivatives
    dLdtheta1 = (1/T)*((W1*dL1dtheta1) + (1/T)*(W2*dL2dtheta1));
    dLdtheta2 = (1/T)*((W1*dL1dtheta2) + (1/T)*(W2*dL2dtheta2));

    % Display results
    display(iter);
    display(theta1);
    display(theta2);
    display(cost);
    display((1/T)*W1*dL1dtheta1);
    display((1/T)*W1*dL1dtheta2);
    display((1/T)*W2*dL2dtheta1);
    display((1/T)*W2*dL2dtheta2);
    display(dLdtheta1);
    display(dLdtheta2);
    
    % Stores results before gradient-based update
    history(1,iter) = iter;
    history(2,iter) = cost;
    history(3,iter) = dLdtheta1;
    history(4,iter) = dLdtheta2;
    history(5,iter) = theta1;
    history(6,iter) = theta2;
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Update controllable parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step sequences
    step1 = 300/(iter^(3/2));
    step2 = 100/(iter^(3/2));
    
    % Update controllable parameters
    theta1 = theta1 - (step1*dLdtheta1);
    theta2 = theta2 - (step2*dLdtheta2);
    
    theta1 = max(theta1, 1.0);
    theta1 = min(theta1, 7.5);
    theta2 = max(theta2, 8.0);
    theta2 = min(theta2, 15.0);
end


% Plots results
figure(1)
subplot(1,2,1);
plot(history(1,1:end), history(2,1:end), '-r*');
xlabel('Iterations');
ylabel('Average Cost');
subplot(1,2,2);
plot(history(1,1:end), history(5,1:end), '-g*');
hold on;
plot(history(1,1:end), history(6,1:end), '-b+');
xlabel('Iterations');
ylabel('PSA threshold values');
legend('\theta_{1}','\theta_{2}');
set(1, 'Color', 'w');


% % Plot simulation profiles
% figure(index)
% plot(clock,HSCpop+CRCpop,'-r')
% title('Patient #15');
% xlabel('Time (days)');
% ylabel('PSA level (nmol/l)');
% %legend('Road 1','Road 2');
% %axis([0 clock(end) 4 30]);
% set(index, 'Color', 'w');
% index = index+1;
% 
% figure(index)
% plot(clock,androgen,'-b')
% title('Patient #15');
% xlabel('Time (days)');
% ylabel('Androgen level (nmol/l)');
% %legend('Road 1','Road 2');
% %axis([0 clock(end) 0 10]);
% set(index, 'Color', 'w');
% index = index+1;
% 
% figure(index)
% plot(clock,HSCpop,'-m')
% title('Patient #15');
% xlabel('Time (days)');
% ylabel('HSC level (nmol/l)');
% %legend('Road 1','Road 2');
% %axis([0 clock(end) 4 30]);
% set(index, 'Color', 'w');
% index = index+1;
% 
% figure(index)
% plot(clock,CRCpop,'-g')
% title('Patient #15');
% xlabel('Time (days)');
% ylabel('CRC level (nmol/l)');
% %legend('Road 1','Road 2');
% %axis([0 clock(end) 0 5]);
% set(index, 'Color', 'w');
% index = index+1;
% 
% 
% % Save simulation results to file
% fid = fopen('simulation.txt','w');
% for j = 1:length(HSCpop)
%    fprintf(fid,'%f\t%f\t%f\t%f\t%f\r\n', clock(j),HSCpop(j),CRCpop(j),androgen(j),HSCpop(j)+CRCpop(j));
% end
% fclose(fid);