Этот пример показывает, как моделировать и анализировать модель в SimBiology ® с использованием физиологической модели системы глюкозо-инсулина у нормальных и диабетических людей.
Этот пример требует Statistics and Machine Learning Toolbox™ и Optimization Toolbox™.
Симуляция питания системы глюкозы-инсулина. К. Далла Ман, Р. А. Рицца и К. Кобелли. Транзакции IEEE по биомедицинской инженерии (2007) 54 (10), 1740-1749.
Системная модель поглощения глюкозы в полости рта: валидация на данных Золотого стандарта. К. Далла Ман, М. Камильери и К. Кобелли. Транзакции IEEE по биомедицинской инженерии (2006) 53 (12), 2472-2478.
Реализуйте SimBiology модель глюкозо-инсулинового ответа.
Симулируйте ответ глюкозо-инсулина на один или несколько приемов пищи для нормальных и нарушенных (диабетических) субъектов.
Выполните оценку параметра, используя sbiofit со стратегией принудительной функции.
В своей публикации 2007 года Dalla Man et al. разработать модель ответа глюкозы-инсулина человека после еды. Эта модель описывает динамику системы, используя обыкновенные дифференциальные уравнения. Авторы использовали свою модель для моделирования ответа глюкозо-инсулина после одного или нескольких приемов пищи, для нормальных людей и для людей с различного рода нарушениями инсулина. Искажения были представлены как альтернативные наборы значений параметров и начальных условий.
Мы реализовали модель SimBiology, m1около:
Перемещение уравнений модели в Dalla Man et al. (2007) в реакции, правила и события.
Организация модели в два отделения, одно для связанных с глюкозой видов и реакций (названо Glucose appearance) и один для связанных с инсулином видов и реакций (названный Insulin secretion).
Использование значений параметров и начальных условий из уравнений модели и из таблицы 1 и фигуры 1.
Включая уравнение для скорости опорожнения желудка, представленное в Dalla Man et al. (2006).
Установка модулей для всех видов, отсеков и параметров, как задано Dalla Man et al. (2007), которая позволяет моделировать модель SimBiology с помощью единичного преобразования. (Обратите внимание, что SimBiology также поддерживает использование безразмерных параметров, устанавливая их ValueUnits свойство к dimensionless.)
Установка конфигурации модели TimeUnits на hour, поскольку симуляции проводились в течение 7 или 24 часов.
Использование базиса из 1 килограмма массы тела для преобразования видов и параметров, которые были нормированы по массе тела в исходной модели. Это делает видовые модули в количестве или концентрации, как того требует SimBiology.
Мы представляли нарушения инсулина в модели SimBiology как варианты объектов со следующими именами:
Диабетик 2 типа
Низкая чувствительность к инсулину
Высокая бета- камера быстродействие
Низкая бета- камера быстродействие
Высокая чувствительность к инсулину
Мы представляли еду в модели SimBiology в качестве объектов дозы:
Доза по наименованию Single Meal представляет собой один прием пищи из 78 граммов глюкозы в начале симуляции.
Доза по наименованию Daily Life представляет собой ценность одного дня еды, относительно симуляции начиная с полуночи: завтрак - 45 граммов глюкозы в 8 часа времени симуляции (8 утра), обед - 70 граммов глюкозы в 12 часа (полдень), а ужин - 70 граммов глюкозы в 20 часа (8 вечера).
Схема модели SimBiology показана ниже:
Загрузите модель.
sbioloadproject('insulindemo', 'm1')
Подавить информационное предупреждение, которое выдается во время симуляций.
warnSettings = warning('off', 'SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless');
Выберите Single Meal дозировать объект и отобразить его свойства.
mealDose = sbioselect(m1, 'Name', 'Single Meal'); get(mealDose)
ans = struct with fields:
Amount: 78
Interval: 0
Rate: 0
RepeatCount: 0
StartTime: 0
Active: 0
AmountUnits: 'gram'
DurationParameterName: ''
EventMode: 'stop'
LagParameterName: ''
RateUnits: ''
TargetName: 'Dose'
TimeUnits: 'hour'
Name: 'Single Meal'
Parent: [1x1 SimBiology.Model]
Notes: ''
Tag: ''
Type: 'repeatdose'
UserData: []
Симулируйте в течение 7 часов.
configset = getconfigset(m1,'active');
configset.StopTime = 7;Отображение временных модулей симуляции (и StopTime модулей).
configset.TimeUnits
ans = 'hour'
Симулируйте один прием пищи для нормального субъекта.
normalMealSim = sbiosimulate(m1, configset, [], mealDose);
Выберите диабетический вариант Type 2 и отобразите его свойства.
diabeticVar = sbioselect(m1, 'Name', 'Type 2 diabetic')
diabeticVar = SimBiology Variant - Type 2 diabetic (inactive) ContentIndex: Type: Name: Property: Value: 1 parameter Plasma Volume ... Value 1.49 2 parameter k1 Value 0.042 3 parameter k2 Value 0.071 4 parameter Plasma Volume ... Value 0.04 5 parameter m1 Value 0.379 6 parameter m2 Value 0.673 7 parameter m4 Value 0.269 8 parameter m5 Value 0.0526 9 parameter m6 Value 0.8118 10 parameter Hepatic Extrac... Value 0.6 11 parameter kmax Value 0.0465 12 parameter kmin Value 0.0076 13 parameter kabs Value 0.023 14 parameter kgri Value 0.0465 15 parameter f Value 0.9 16 parameter a Value 6e-05 17 parameter b Value 0.68 18 parameter c Value 0.00023 19 parameter d Value 0.09 20 parameter Stomach Glu Af... Value 125 21 parameter kp1 Value 3.09 22 parameter kp2 Value 0.0007 23 parameter kp3 Value 0.005 24 parameter kp4 Value 0.0786 25 parameter ki Value 0.0066 26 parameter [Ins Ind Glu U... Value 1 27 parameter Vm0 Value 4.65 28 parameter Vmx Value 0.034 29 parameter Km Value 466.21 30 parameter p2U Value 0.084 31 parameter K Value 0.99 32 parameter alpha Value 0.013 33 parameter beta Value 0.05 34 parameter gamma Value 0.5 35 parameter ke1 Value 0.0007 36 parameter ke2 Value 269 37 parameter Basal Plasma G... Value 164.18 38 parameter Basal Plasma I... Value 54.81
Симулируйте один прием пищи для диабетика 2 типа.
diabeticMealSim = sbiosimulate(m1, configset, diabeticVar, mealDose);
Сравните результаты для наиболее важных выходов симуляции.
Плазма Глюкоза (вид Plasma Glu Conc)
Плазменный инсулин (вид Plasma Ins Conc)
Эндогенное производство глюкозы (параметр Glu Prod)
Скорость внешнего вида глюкозы (параметр Glu Appear Rate)
Использование глюкозы (параметр Glu Util)
Секреция инсулина (параметр 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
Обратите внимание на намного более высокие концентрации глюкозы и инсулина в плазме, а также на длительную длительность использования глюкозы и секреции инсулина.
Установите StopTime симуляции до 24 часов.
configset.StopTime = 24;
Выберите ежедневную дозу еды.
dayDose = sbioselect(m1, 'Name', 'Daily Life');
Симулируйте три приема пищи для нормального субъекта.
normalDaySim = sbiosimulate(m1, configset, [], dayDose);
Симулируйте следующие комбинации нарушений:
Нарушение 1: Низкая чувствительность к инсулину
Ухудшение 2: Ухудшение 1 с высокой бета камеры чувствительностью
Нарушение 3: Низкая бета камеры чувствительность
Нарушение 4: Нарушение 3 с высокой чувствительностью к инсулину
Сохраните нарушения в массиве ячеек.
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
Обратите внимание, что либо низкая чувствительность к инсулину (штриховая зеленая линия, ) или низкая чувствительность бета-клеток (пунктирно-пунктирная голубая линия, ) приводят к повышенным и длительным концентрациям глюкозы в плазме (верхняя строка графиков). Низкая чувствительность в одной системе может быть частично компенсирована высокой чувствительностью в другой системе. Например, низкая чувствительность к инсулину и высокая чувствительность к бета-клеткам (пунктирная красная линия, ) приводит к относительно нормальным концентрациям глюкозы в плазме (верхняя строка графиков). Однако в этом случае полученная концентрация инсулина в плазме чрезвычайно высока (нижняя строка графиков).
Вместо того, чтобы одновременно оценивать параметры для всей модели, авторы выполняют оценку параметра для различных подсистем модели, используя стратегию принудительной функции. Этот подход требует дополнительных экспериментальных данных для «входов» подмодели. Во время подбора кривой входные данные определяют динамику видов входов. (В полной модели динамика входов определяется из дифференциальных уравнений.) В терминах SimBiology можно реализовать принудительную функцию как повторное правило назначения, которое управляет значением вида или параметра, который служит входом для подсистемы модели. В следующих разделах мы используем возможности подгонки параметров SimBiology, чтобы уточнить сообщенные значения параметров авторов.
nlinfitЖелудочно-кишечная модель представляет, как глюкоза в еде переносится через желудок, кишечник и кишечник, а затем всасывается в плазму. Входом в эту подсистему является количество глюкозы в еде, а выходом - скорость внешнего вида глюкозы в плазме. Тем не менее, мы также оцениваем размер еды, поскольку значение, сообщенное авторами, не соответствует параметрам и результатам симуляции. Поскольку этот вход происходит только в начале симуляции, принудительная функция не требуется.
Функция sbiofit поддерживает оценку параметров в моделях SimBiology, используя несколько различных алгоритмов из MATLAB™, Statistics and Machine Learning Toolbox, Optimization Toolbox и Global Optimization Toolbox. Сначала оцените параметры с помощью функции Statistics and Machine Learning Toolbox 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);
Norm of Norm of
Iteration SSE Gradient Step
-----------------------------------------------------------
0 43.798
1 3.1179 32.3923 0.462126
2 2.13847 1.04093 0.0510667
3 2.13751 0.000736706 0.0045159
4 2.13749 5.76665e-05 0.000866769
5 2.13749 5.87539e-06 0.000194306
6 2.13749 4.05387e-07 4.77789e-05
7 2.13749 2.0556e-08 1.31015e-05
Iterations terminated: relative change in SSE less than OPTIONS.TolFun
fminuncТеперь оцените параметры с помощью функции Optimization Toolbox 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);
Сравните симуляцию до и после подбора кривой.
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')
Обратите внимание, что модель подходит для экспериментальных данных значительно лучше, если размер еды (Dose) значительно больше, чем сообщалось, параметр kmax значительно больше, чем сообщалось, и kabs меньше, чем сообщалось.
Модель мышцы и жировой ткани представляет, как глюкоза используется в организме. «входами» в эту подсистему являются концентрации инсулина в плазме (Plasma Ins Conc), эндогенное производство глюкозы (Glu Prod), и скорость внешнего вида глюкозы (Glu Appear Rate). «выходами» являются концентрации глюкозы в плазме (Plasma Glu Conc) и скорость использования глюкозы (Glu Util).
Поскольку входы являются функцией времени, они должны быть реализованы как принудительные функции. В частности, значения Plasma Ins Conc, Glu Prod, и Glu Appear Rate управляются повторными назначениями, которые вызывают функции для выполнения линейной интерполяции сообщенных экспериментальных значений. При использовании этих функций для управления видом или параметром необходимо сделать неактивным любое другое правило, которое используется для установки его значения. Чтобы облегчить выбор этих правил, свойства Name правила содержат значимые имена.
% 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')
PlasmaInsForcingFcn = SimBiology Rule Array Index: RuleType: Rule: 1 repeatedAssignment [Plasma Ins Conc] = [picomole per liter]*PlasmaInsulin(time/[one hour])
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')
GluProdForcingFcn = SimBiology Rule Array Index: RuleType: Rule: 1 repeatedAssignment [Glu Prod] = [milligram per minute]*EndogenousGlucoseProduction(time/[one hour])
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')
GluRateForcingFcn = SimBiology Rule Array Index: RuleType: Rule: 1 repeatedAssignment [Glu Appear Rate] = [milligram per minute]*GlucoseAppearanceRate(time/[one hour])
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);
Norm of Norm of
Iteration SSE Gradient Step
-----------------------------------------------------------
0 2181.52
1 1085.53 39020.2 0.610029
2 1059.87 34535.8 0.741592
3 743.846 18217.2 0.752823
4 544.725 12620.2 0.738337
5 328.744 2879.49 0.383848
6 321.805 3860.01 0.159715
7 269.919 4675.63 0.200186
8 265.575 2548.53 0.0536833
9 263.684 2060.19 0.0293006
10 262.848 652.619 0.0294674
11 262.346 1087.79 0.00407125
12 261.897 674.32 0.00347999
13 261.886 950.085 0.00199629
14 261.523 1327.79 0.00115098
15 261.475 723.834 0.0025953
16 261.28 312.505 0.00210966
17 257.985 625.055 0.0336004
18 256.732 633.64 0.0251723
19 256.535 526.063 0.0837353
20 252.174 2095.9 0.167036
21 251.642 398.83 0.0677341
22 250.661 575.847 0.0124814
23 250.601 1371.22 0.00134259
24 250.509 670.542 0.00182034
25 250.341 881.41 0.00336496
26 250.205 601.85 0.0028711
27 250.136 588.918 0.00241238
28 250.096 405.799 0.0012213
29 250.026 253.856 0.0196512
30 249.718 973.329 0.0262556
31 249.718 2398.16 3.39399e-17
Iterations terminated: relative norm of the current step is less than OPTIONS.TolX
Постройте график изменения значений параметров относительно сообщаемых значений.
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');
Обратите внимание, что значительное увеличение некоторых параметров, таких как Vmx, позволяет намного лучше подгонка поздние концентрации глюкозы в плазме.
Восстановите параметры предупреждения.
warning(warnSettings);
SimBiology содержит несколько функций, которые облегчают реализацию и симуляцию сложной модели системы глюкоза-инсулин. Реакции, события и правила обеспечивают естественный способ описания динамики системы. Преобразование модулей измерения позволяет задавать виды и параметры в удобных модулях и обеспечивает размерную согласованность модели. Объекты дозы являются простым способом описания повторяющихся входов в модель, таких как дневное расписание приема пищи в этом примере. SimBiology также обеспечивает встроенную поддержку задач анализа, таких как симуляция и оценка параметра.