Этот пример показывает, как моделировать и анализировать модель в SimBiology ® с использованием физиологически основанной модели глюкозно-инсулиновой системы у нормальных и диабетических людей.
В этом примере требуются Toolbox™ статистики и машинного обучения, а также Toolbox™ оптимизации.
Имитационная модель еды системы глюкоза-инсулин. К. Далла Ман, Р. А. Ризза и К. Кобелли. IEEE Transactions on Biomedical Engineering (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);
Выберите диабетический вариант типа 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
Обратите внимание, что либо низкая чувствительность к инсулину (пунктирная зеленая линия, -β), либо низкая чувствительность к бета-клеткам (пунктирная голубая линия, Ins-β) приводят к увеличению и увеличению концентраций глюкозы в плазме (верхний ряд графиков). Низкая чувствительность в одной системе может быть частично компенсирована высокой чувствительностью в другой системе. Например, низкая чувствительность к инсулину и высокая чувствительность к бета-клеткам (пунктирная красная линия+ β) приводит к относительно нормальным концентрациям глюкозы в плазме (верхний ряд графиков). Однако в этом случае полученная концентрация инсулина в плазме чрезвычайно высока (нижний ряд графиков).
Вместо одновременной оценки параметров для всей модели авторы выполняют оценку параметров для различных подсистем модели с использованием стратегии форсирующей функции. Этот подход требует дополнительных экспериментальных данных для «входов» подмодели. Во время подгонки входные данные определяют динамику входных видов. (В полной модели динамика входов определяется из дифференциальных уравнений.) В терминах 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Теперь оцените параметры с помощью функции «Панель инструментов оптимизации» 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 управляются повторными назначениями, которые вызывают функции для выполнения линейной интерполяции сообщенных экспериментальных значений. При использовании этих функций для управления видом или параметром необходимо сделать неактивным любое другое правило, которое используется для установки его значения. Чтобы облегчить выбор этих правил, свойства правила Имя содержат значимые имена.
% 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 также предоставляет встроенную поддержку для таких задач анализа, как моделирование и оценка параметров.