sbioloadproject('insulindemo', 'm1')

warnSettings = warning('off', 'SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless');

mealDose = sbioselect(m1, 'Name', 'Single Meal');
get(mealDose)

configset = getconfigset(m1,'active');
configset.StopTime = 7;

configset.TimeUnits

normalMealSim = sbiosimulate(m1, configset, [], mealDose);

diabeticVar = sbioselect(m1, 'Name', 'Type 2 diabetic')

diabeticMealSim = sbiosimulate(m1, configset, diabeticVar, mealDose);

outputNames = {'Plasma Glu Conc', 'Plasma Ins Conc', 'Glu Prod', ...
    'Glu Appear Rate', 'Glu Util', 'Ins Secr'};
figure;
for i = 1:numel(outputNames)
    subplot(2, 3, i);

    [tNormal, yNormal  ]  = normalMealSim.selectbyname(outputNames{i});
    [tDiabetic, yDiabetic]  = diabeticMealSim.selectbyname(outputNames{i});
    
    plot( tNormal    , yNormal   , '-'       , ... 
          tDiabetic  , yDiabetic , '--'      );
  
    % Annotate figures 
    outputParam = sbioselect(m1, 'Name', outputNames{i});  
    title(outputNames{i});
    xlabel('time (hour)');
    if strcmp(outputParam.Type, 'parameter')
        ylabel(outputParam.ValueUnits);
    else
        ylabel(outputParam.InitialAmountUnits);
    end
    xlim([0 7]);
    
    % Add legend
    if i == 3
        legend({'Normal', 'Diabetic'}, 'Location', 'Best');
    end
    
end

configset.StopTime = 24;

dayDose  = sbioselect(m1, 'Name', 'Daily Life');

normalDaySim = sbiosimulate(m1, configset, [], dayDose);

impairVars{1} =  sbioselect(m1, 'Name', 'Low insulin sensitivity'    ) ;
impairVars{2} = [impairVars{1}, ...
                    sbioselect(m1, 'Name', 'High beta cell responsivity')];
impairVars{3} =  sbioselect(m1, 'Name', 'Low beta cell responsivity' ) ;
impairVars{4} = [impairVars{3}, ...
                    sbioselect(m1, 'Name', 'High insulin sensitivity'   )];

for i = 1:4
    impairSims(i) = sbiosimulate(m1, configset, impairVars{i}, dayDose);
end

figure;
outputNames = {'Plasma Glu Conc', 'Plasma Ins Conc'};

legendLabels = {{'Normal'}, ...
    {'-Ins =\beta', '-Ins +\beta'}, ...
    {'=Ins -\beta', '+Ins -\beta'}};
yLimits = [80 240; 0 500];

for i = 1:numel(outputNames)
     
    [tNormal, yNormal] = selectbyname(normalDaySim , outputNames{i} );
    [tImpair, yImpair] = selectbyname(impairSims   , outputNames{i} );
    
    % Plot Normal
    subplot(2, 3, 3*i-2 );
    plot(tNormal, yNormal, 'b-');
    xlim([0 24]);
    ylim(yLimits(i,:));
    xlabel('time (hour)');
    legend(legendLabels{1}, 'Location', 'NorthWest');
    
    % Plot Low Insulin
    subplot(2, 3, 3*i-1 );
    plot(tImpair{1}, yImpair{1}, 'g--', tImpair{2}, yImpair{2}, 'r:');
    xlim([0 24]); 
    ylim(yLimits(i,:));
    xlabel('time (hour)');
    legend(legendLabels{2}, 'Location', 'NorthWest');
    title(outputNames{i});
    
    % Plot Low Beta
    subplot(2, 3, 3*i   );
    plot(tImpair{3}, yImpair{3}, 'c-.', tImpair{4}, yImpair{4}, 'm-');
    xlim([0 24]); 
    ylim(yLimits(i,:));
    xlabel('time (hour)');
    legend(legendLabels{3}, 'Location', 'NorthWest');
    
end

% Load the experimental data
fitData = groupedData(readtable('GlucoseData.csv', 'Delimiter', ','));

% Set the units on the data
fitData.Properties.VariableUnits = {...
    'hour', ...                % Time units
    'milligram/minute', ...    % GluRate units
    'milligram/deciliter', ... % PlasmaGluConc units
    'milligram/minute', ...    % GluUtil units
    };

% Identify which model components corresponds to observed data variables.
gastroFitObs = '[Glu Appear Rate] = GluRate';

% Estimate the value of the following parameters:
gastroFitEst = estimatedInfo({'kmax', 'kmin', 'kabs', 'Dose'});

% Ensure the parameter estimates are always positive during estimation by
% using a log transform on all parameters.
[gastroFitEst.Transform] = deal('log');

% Set the initial estimate for Dose to the reported meal dose amount. The
% remaining initial estimates will be taken from the parameter values in
% the model.
gastroFitEst(4).InitialValue = mealDose.Amount;

% Generate simulation data with the initial parameter estimates
configset.StopTime = 7;
gastroInitSim = sbiosimulate(m1, mealDose);
 
% Fit the data using |nlinfit|, displaying output at each iteration
fitOptions = statset('Display', 'iter');
[gastroFitResults, gastroFitSims] = sbiofit(m1, fitData, ...
    gastroFitObs, gastroFitEst, [], 'nlinfit', fitOptions);

% Fit the data, plotting the objective function at each iteration
fitOptions2 = optimoptions('fminunc', 'PlotFcns', @optimplotfval);
[gastroFitResults(2), gastroFitSims(2)] = sbiofit(m1, fitData, ...
    gastroFitObs, gastroFitEst, [], 'fminunc', fitOptions2);

gastroSims = selectbyname([gastroInitSim gastroFitSims], 'Glu Appear Rate');

figure;
plot(gastroSims(1).Time , gastroSims(1).Data , '-'  , ...
     gastroSims(2).Time , gastroSims(2).Data , '--' , ...
     gastroSims(3).Time , gastroSims(3).Data , '-.' , ...
     fitData.Time       , fitData.GluRate, 'x'  );
 
xlabel('Time (hour)');
ylabel('mg/min');
legend('Reported', 'Estimated (nlinfit)', ...
    'Estimated (fminunc)', 'Experimental');
title('Glucose Appearance Fit');

figure;
fitResults = [gastroFitResults(1).ParameterEstimates.Estimate ...
    gastroFitResults(2).ParameterEstimates.Estimate];
% The initial values for kmax, kmin, and kabs come from the model.
gastroFitInitValues = [
    get(sbioselect(m1, 'Name', 'kmax'), 'Value')
    get(sbioselect(m1, 'Name', 'kmin'), 'Value')
    get(sbioselect(m1, 'Name', 'kabs'), 'Value')
    gastroFitEst(4).InitialValue
    ];
relFitChange = fitResults ./ [gastroFitInitValues gastroFitInitValues] - 1;
bar(relFitChange);
ax = gca;
ax.XTickLabel = {gastroFitEst.Name};
ylabel('Relative change in estimated values');
title('Comparing Reported and Estimated Gastrointestinal Parameter Values');
legend({'nlinfit', 'fminunc'}, 'Location', 'North')

% Create forcing functions for the "inputs":
% Plasma Insulin
PlasmaInsRule       = sbioselect(m1, 'Name', 'Plasma Ins Conc definition');
PlasmaInsForcingFcn = sbioselect(m1, 'Name', 'Plasma Ins Conc forcing function')

PlasmaInsRule.Active       = false;
PlasmaInsForcingFcn.Active = true;

% Endogenous Glucose Production (Glu Prod)
GluProdRule       = sbioselect(m1, 'Name', 'Glu Prod definition');
GluProdForcingFcn = sbioselect(m1, 'Name', 'Glu Prod forcing function')

GluProdRule.Active       = false;
GluProdForcingFcn.Active = true;

% Glucose Rate of Appearance (Glu Appear Rate)
GluRateRule       = sbioselect(m1, 'Name', 'Glu Appear Rate definition');
GluRateForcingFcn = sbioselect(m1, 'Name', 'Glu Appear Rate forcing function')

GluRateRule.Active       = false;
GluRateForcingFcn.Active = true;

% Simulate with the initial parameter values
muscleInitSim = sbiosimulate(m1);

% Identify which model components corresponds to observed data variables.
muscleFitObs = {'[Plasma Glu Conc] = PlasmaGluConc', ...
    '[Glu Util] = GluUtil'};

% Estimate the value of the following parameters:
muscleFitEst = estimatedInfo({'[Plasma Volume (Glu)]', 'k1', 'k2', ...
    'Vm0', 'Vmx', 'Km', 'p2U'});

% Ensure the parameter estimates are always positive during estimation by
% using a log transform on all parameters.
[muscleFitEst.Transform] = deal('log');

% Fit the data, displaying output at each iteration 
[muscleFitResults, muscleFitSim] = sbiofit(m1, fitData, ...
    muscleFitObs, muscleFitEst, [], 'nlinfit', fitOptions);

figure;
muscleFitInitValues = [
    get(sbioselect(m1, 'Name', 'Plasma Volume (Glu)'), 'Value')
    get(sbioselect(m1, 'Name', 'k1'), 'Value')
    get(sbioselect(m1, 'Name', 'k2'), 'Value')
    get(sbioselect(m1, 'Name', 'Vm0'), 'Value')
    get(sbioselect(m1, 'Name', 'Vmx'), 'Value')
    get(sbioselect(m1, 'Name', 'Km'), 'Value')
    get(sbioselect(m1, 'Name', 'p2U'), 'Value')
    ];

bar(muscleFitResults.ParameterEstimates.Estimate ./ muscleFitInitValues - 1);
ax = gca;
ax.XTickLabel = {muscleFitEst.Name};
ylabel('Relative change in estimated values');
title('Comparing Reported and Estimated Glucose Parameter Values');

PlasmaInsRule.Active = true;
GluProdRule.Active   = true;
GluRateRule.Active   = true;

PlasmaInsForcingFcn.Active = false;
GluProdForcingFcn.Active   = false;
GluRateForcingFcn.Active   = false;

muscleSims = selectbyname([muscleInitSim muscleFitSim], ...
    {'Plasma Glu Conc', 'Glu Util'});
figure;
plot(muscleSims(1).Time, muscleSims(1).Data(:,1), '-', ...
    muscleSims(2).Time, muscleSims(2).Data(:,1), '--', ...
    fitData.Time, fitData.PlasmaGluConc, 'x');
xlabel('Time (hour)');
ylabel('mg/dl');
legend('Initial (Simulation)', 'Estimated (Simulation)', 'Experimental');
title('Plasma Glucose Fit');

figure;
plot(muscleSims(1).Time, muscleSims(1).Data(:,2), '-', ...
    muscleSims(2).Time, muscleSims(2).Data(:,2), '--', ...
    fitData.Time, fitData.GluUtil, 'x');
xlabel('Time (hour)');
ylabel('mg/min');
legend('Initial (Simulation)', 'Estimated (Simulation)', 'Experimental');
title('Glucose Utilization Fit');

warning(warnSettings);