www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/insulindemo.m

    %% Simulating the Glucose-Insulin Response
%
% This example shows how to simulate and analyze a model in SimBiology(R)
% using a physiologically based model of the glucose-insulin system in
% normal and diabetic humans.
%
% This example requires Statistics and Machine Learning Toolbox(TM) and
% Optimization Toolbox(TM).

% Copyright 2011-2014 The MathWorks, Inc.

%% References
% # Meal Simulation Model of the Glucose-Insulin System. C. Dalla Man, R.A.
% Rizza, and C. Cobelli. _IEEE Transactions on Biomedical Engineering_
% (2007) 54(10), 1740-1749.
% # A System Model of Oral Glucose Absorption: Validation on Gold Standard
% Data. C. Dalla Man, M. Camilleri, and C. Cobelli. _IEEE Transactions on
% Biomedical Engineering_ (2006) 53(12), 2472-2478.

%% Aims
%
% * Implement a SimBiology model of the glucose-insulin response.
% * Simulate the glucose-insulin response to one or more meals for normal
% and impaired (diabetic) subjects.
% * Perform parameter estimation using |sbiofit| with a forcing function
% strategy.

%% Background
%
% In their 2007 publication, Dalla Man et al. develop a model for the human
% glucose-insulin response after a meal. This model describes the dynamics
% of the system using ordinary differential equations. The authors used
% their model to simulate the glucose-insulin response after one or more
% meals, for normal human subjects and for human subjects with various
% kinds of insulin impairments. The impairments were represented as
% alternate sets of parameter values and initial conditions.
%
% We implemented the SimBiology model, |m1|, by:
%
% * Translating the model equations in Dalla Man et al. (2007) into
% reactions, rules, and events.
% * Organizing the model into two compartments, one for glucose-related
% species and reactions (named |Glucose appearance|) and one for
% insulin-related species and reactions (named |Insulin secretion|).
% * Using the parameter values and initial conditions from the model
% equations and from Table 1 and Figure 1.
% * Including an equation for the gastric emptying rate as presented in
% Dalla Man et al. (2006).
% * Setting the units for all species, compartments, and parameters as
% specified by Dalla Man et al. (2007), which allows the SimBiology model
% to be simulated using unit conversion. (Note that SimBiology also
% supports the use of dimensionless parameters by setting their |ValueUnits|
% property to |dimensionless|.)
% * Setting the configuration set |TimeUnits| to |hour|, since simulations
% were conducted over 7 or 24 hours.
% * Using a basis of 1 kilogram of body weight to transform species and
% parameters that were normalized by body weight in the original model.
% Doing so made species units in amount or concentration, as required by
% SimBiology.
%
% We represented the insulin impairments in the SimBiology model as variant
% objects with the following names:
%
% * Type 2 diabetic
% * Low insulin sensitivity
% * High beta cell responsivity
% * Low beta cell responsivity
% * High insulin sensitivity
%
% We represented the meals in the SimBiology model as dose objects:
%
% * A dose named |Single Meal| represents a single meal of 78 grams of
% glucose at the start of a simulation.
% * A dose named |Daily Life| represents one day's worth of meals, relative
% to a simulation starting at midnight: breakfast is 45 grams of glucose at
% 8 hours of simulation time (8 a.m.), lunch is 70 grams of glucose at 12
% hours (noon), and dinner is 70 grams of glucose at 20 hours (8 p.m.).
%
% A diagram of the SimBiology model is shown below:
%
% <<../insulindemo_diagram.png>>

%% Setup
% Load the model.
sbioloadproject('insulindemo', 'm1')

%%
% Suppress an informational warning that is issued during simulations.
warnSettings = warning('off', 'SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless');

%% Simulating the Glucose-Insulin Response for a Normal Subject

%%
% Select the |Single Meal| dose object and display its properties.
mealDose = sbioselect(m1, 'Name', 'Single Meal');
get(mealDose)

%%
% Simulate for 7 hours.
configset = getconfigset(m1,'active');
configset.StopTime = 7;

%%
% Display the simulation time units (and |StopTime| units).
configset.TimeUnits

%%
% Simulate a single meal for a normal subject.
normalMealSim = sbiosimulate(m1, configset, [], mealDose);

%% Simulating the Glucose-Insulin Response for a Type 2 Diabetic

%%
% Select the Type 2 diabetic variant and display its properties.
diabeticVar = sbioselect(m1, 'Name', 'Type 2 diabetic')

%%
% Simulate a single meal for a Type 2 diabetic.
diabeticMealSim = sbiosimulate(m1, configset, diabeticVar, mealDose);

%%
% Compare the results for the most important outputs of the simulation.
%
% * Plasma Glucose (species |Plasma Glu Conc|)
% * Plasma Insulin (species |Plasma Ins Conc|)
% * Endogenous Glucose Production (parameter |Glu Prod|)
% * Glucose Rate of Appearance (parameter |Glu Appear Rate|)
% * Glucose Utilization (parameter |Glu Util|)
% * Insulin Secretion (parameter |Ins Secr|)

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

%%
% Note the much higher concentrations of glucose and insulin in the plasma,
% as well as the prolonged duration of glucose utilization and insulin
% secretion.

%% Simulating One Day with Three Meals for a Normal Subject

%%
% Set the simulation |StopTime| to 24 hours.
configset.StopTime = 24;

%%
% Select daily meal dose.
dayDose  = sbioselect(m1, 'Name', 'Daily Life');

%%
% Simulate three meals for a normal subject.
normalDaySim = sbiosimulate(m1, configset, [], dayDose);

%% Simulating One Day with Three Meals for Impaired Subjects
% Simulate the following combinations of impairments:
%
% * Impairment 1: Low insulin sensitivity
% * Impairment 2: Impairment 1 with high beta cell responsivity
% * Impairment 3: Low beta cell responsivity
% * Impairment 4: Impairment 3 with high insulin sensitivity

%%
% Store the impairments in a cell array.
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'   )];

%%
% Simulate each impairment.
for i = 1:4
    impairSims(i) = sbiosimulate(m1, configset, impairVars{i}, dayDose);
end

%%
% Compare the plasma glucose and plasma insulin results.
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

%%
% Note that either low insulin sensitivity (dashed green line, 
% $-Ins =\beta$) or low beta-cell sensitivity (dashed-dotted cyan line,
% $=Ins -\beta$) lead to increased and prolonged plasma glucose
% concentrations (top row of plots). Low sensitivity in one system can be
% partially compensated by high sensitivity in another system. For example,
% low insulin sensitivity and high beta-cell sensitivity (dotted red line,
% $-Ins +\beta$) results in relatively normal plasma glucose concentrations
% (top row of plots). However, in this case, the resulting plasma insulin
% concentration is extremely high (bottom row of plots).

%%  Parameter Estimation Methodology
% Rather than simultaneously estimating parameters for the entire model,
% the authors perform parameter estimation for different subsystems of the
% model using a forcing function strategy. This approach requires
% additional experimental data for the "inputs" of the submodel. During
% fitting, the input data determine the dynamics of the inputs species. (In
% the full model, the dynamics of the inputs are determined from the
% differential equations.) In SimBiology terms, you can implement a forcing
% function as a repeated assignment rule that controls the value of a
% species or parameter that serves as an input for a subsystem of the
% model. In the following sections, we use the parameter fitting
% capabilities of SimBiology to refine the authors' reported parameter
% values.

%% Fitting the Gastrointestinal Model of Glucose Appearance using |nlinfit|
% The gastrointestinal model represents how glucose in a meal is
% transported through the stomach, gut, and intestine, and then absorbed
% into the plasma. The input to this subsystem is the amount of glucose in
% a meal, and the output is the rate of appearance of glucose in the
% plasma. However, we also estimate the meal size since the value reported
% by the authors is inconsistent with the parameters and simulation
% results. Because this input only occurs at the start of the simulation,
% no forcing function is required.
%
% The function |sbiofit| supports the estimation of parameters in
% SimBiology models using several different algorithms from MATLAB(TM),
% Statistics and Machine Learning Toolbox, Optimization Toolbox, and Global
% Optimization Toolbox. First, estimate the parameters using Statistics and
% Machine Learning Toolbox function |nlinfit|.

% 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);

%% Fitting the Data Using |fminunc|
% Now, estimate the parameters using the Optimization Toolbox function
% |fminunc|.

% 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);

%%
% Compare the simulation before and after fitting.
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');

%%
% Plot the change in parameter values, relative to reported values.
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')

%%
% Note that the model fits the experimental data significantly better if
% the meal size (|Dose|) is significantly larger than reported, the
% parameter |kmax| is significantly larger than reported, and |kabs| is
% smaller than reported.

%% Fitting the Muscle and Adipose Tissue Model of Glucose Utilization
% The muscle and adipose tissue model represents how glucose is utilized in
% the body. The "inputs" to this subsystem are the concentration of insulin
% in the plasma (|Plasma Ins Conc|), the endogenous glucose production
% (|Glu Prod|), and the rate of appearance of glucose (|Glu Appear Rate|). The
% "outputs" are the concentration of glucose in the plasma (|Plasma Glu Conc|) 
% and the rate of glucose utilization (|Glu Util|).
%
% Because the inputs are a function of time, they need to be implemented as
% forcing functions. Specifically, the values of |Plasma Ins Conc|, |Glu Prod|,
% and |Glu Appear Rate| are controlled by repeated assignments that call
% functions to do linear interpolation of the reported experimental values.
% When using these functions to control a species or parameter, you must
% make inactive any other rule that is used to set its value. To facilitate
% the selection of these rules, the rule Name properties contain meaningful
% names.

% 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);

%%
% Plot the change in parameter values, relative to reported values.
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');

%%
% Clean up the changes to the model.
PlasmaInsRule.Active = true;
GluProdRule.Active   = true;
GluRateRule.Active   = true;

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

%%
% Compare the simulation before and after fitting
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');

%%
% Note that significantly increasing some parameters, such as |Vmx|, allows
% a much better fit of late-time plasma glucose concentrations.

%% Cleanup
% Restore warning settings.
warning(warnSettings);

%% Conclusions
% SimBiology contains several features that facilitate the implementation
% and simulation of a complex model of the glucose-insulin system.
% Reactions, events, and rules provide a natural way to describe the
% dynamics of the system. Unit conversion allows species and parameters to
% be specified in convenient units and ensures the dimensional consistency
% of the model. Dose objects are a simple way to describe recurring inputs
% to a model, such as the daily meal schedule in this example. SimBiology
% also provides built-in support for analysis tasks like simulation and
% parameter estimation.