В этом примере показано, как симулировать и анализировать модель в 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 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 также оказывает встроенную поддержку для аналитических задач как оценка параметра и симуляция.