В этом примере показано, как симулировать и анализировать модель в SimBiology® с помощью физиологически основанной модели системы инсулина глюкозы в нормальных и страдающих от диабета людях.
Этот пример требует Statistics and Machine Learning Toolbox™ и Optimization Toolbox™.
Имитационная модель еды системы инсулина глюкозы. Человек К. Даллы, Р.А. Ризза, и К. Кобелли. Транзакции IEEE на биоинженерии (2007) 54 (10), 1740-1749.
Системная модель устного поглощения глюкозы: валидация на данных о золотом стандарте. Человек К. Даллы, М. Камиллери, и К. Кобелли. Транзакции IEEE на биоинженерии (2006) 53 (12), 2472-2478.
Реализуйте модель SimBiology ответа инсулина глюкозы.
Симулируйте ответ инсулина глюкозы на одну или несколько еды для нормального, и повредил (диабетические) предметы.
Выполните оценку параметра с помощью sbiofit со стратегией функции принуждения.
В их 2 007 публикациях Человек Даллы и др. разрабатывает модель для человеческого ответа инсулина глюкозы после еды. Эта модель описывает динамику системы с помощью обыкновенных дифференциальных уравнений. Авторы использовали свою модель, чтобы симулировать ответ инсулина глюкозы после одной или нескольких еды для нормальных испытуемых людей и для испытуемых людей с различными видами ухудшений инсулина. Ухудшения были представлены как альтернативные наборы значений параметров и начальных условий.
Мы реализовали модель SimBiology, m1:
Перевод уравнений модели в Человеке Даллы и др. (2007) в реакции, правила и события.
Организовывая модель в два отсека, один для связанных с глюкозой разновидностей и реакций (названный Glucose appearance) и один для связанных с инсулином разновидностей и реакций (названный Insulin secretion).
Используя значения параметров и начальные условия от уравнений модели и из Таблицы 1 и рисунка 1.
Включая уравнение для уровня освобождения желудка, как представлено в Человеке Даллы и др. (2006).
Устанавливая модули для всех разновидностей, отсеков и параметров, как задано Человеком Даллы и др. (2007), который позволяет модели SimBiology быть симулированной с помощью модульного преобразования. (Обратите внимание на то, что SimBiology также поддерживает использование безразмерных параметров путем установки их ValueUnits свойство к dimensionless.)
Установка конфигурации модели TimeUnits к hour, поскольку симуляции проводились более чем 7 или 24 часа.
Используя базис 1 килограмма массы тела, чтобы преобразовать разновидности и параметры, которые были нормированы на массу тела в исходной модели. Выполнение так сделанного модулями разновидностей в сумме или концентрации, как требуется SimBiology.
Мы представляли ухудшения инсулина в модели SimBiology, когда вариант возражает со следующими именами:
Диабетик типа 2
Низкая чувствительность инсулина
Высокая бета чувствительность ячейки
Низкая бета чувствительность ячейки
Высокая чувствительность инсулина
Мы представляли еду в модели SimBiology, когда доза возражает:
Доза под названием Single Meal представляет одну еду 78 граммов глюкозы в начале симуляции.
Доза под названием Daily Life представляет ценность одного дня еды, относительно симуляции, запускающейся в полночь: завтрак составляет 45 граммов глюкозы в 8 часов времени симуляции (8 a.m.), ланч составляет 70 граммов глюкозы в 12 часов (полдень), и ужин составляет 70 граммов глюкозы в 20 часов (8 p.m.).
Схему модели 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 kp1 Value 3.09 21 parameter kp2 Value 0.0007 22 parameter kp3 Value 0.005 23 parameter kp4 Value 0.0786 24 parameter ki Value 0.0066 25 parameter [Ins Ind Glu U... Value 1 26 parameter Vm0 Value 4.65 27 parameter Vmx Value 0.034 28 parameter Km Value 466.21 29 parameter p2U Value 0.084 30 parameter K Value 0.99 31 parameter alpha Value 0.013 32 parameter beta Value 0.05 33 parameter gamma Value 0.5 34 parameter ke1 Value 0.0007 35 parameter ke2 Value 269 36 parameter Basal Plasma G... Value 164.18 37 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 2.23537 22.9971 0.596828
2 1.65253 1.36198 0.104842
3 1.64911 0.0938221 0.0264691
4 1.64908 0.000752512 0.00253512
5 1.64908 0.0837661 1.2591e-16
Iterations terminated: relative norm of the current step is less than OPTIONS.TolX
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.51
1 682.537 15418.9 0.477459
2 603.109 25644.9 0.879582
3 503.599 16447.4 0.802291
4 376.519 12343.9 0.324297
5 368.782 8923.16 0.0700439
6 332.709 5653.32 0.0479788
7 306.419 851.074 0.0813683
8 271.443 6011.07 0.156832
9 251.228 3163.88 0.0428459
10 234.782 1700.19 0.0854005
11 233.126 1159.08 0.00724257
12 229.584 2401.43 0.0256606
13 224.26 1594.28 0.038696
14 224.013 1272.59 0.00210287
15 222.788 890.28 0.00323053
16 222.319 1034.28 0.00140676
17 221.957 767.18 0.00268741
18 221.953 445.383 1.35982e-07
19 221.953 945.749 5.04347e-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 также оказывает встроенную поддержку для аналитических задач как оценка параметра и симуляция.