شبیه سازی پاسخ گلوکز و انسولین در انسان عادی و دیابتی با متلب

جمعه ۵ اسفند ۱۴۰۱
بازدید ۶۴۹ نفر
AmRo
(10 امتیاز از 2 رای)
Loading...

شبیه سازی پاسخ گلوکز و انسولین در انسان عادی و دیابتی با MATLAB

این مثال نحوه شبیه سازی و تجزیه و تحلیل یک مدل را در SimBiology با استفاده از یک مدل مبتنی بر فیزیولوژیک سیستم گلوکز و انسولین در انسان های عادی و دیابتی نشان می دهد. این مثال به تولباکس Statistics and Machine Learning Toolbox و تولباکس Optimization Toolbox نیاز دارد. اهداف کلی در این مثل به این شکل است که یک مدل SimBiology پاسخ گلوکز و انسولین را اجرا می کند. سپس پاسخ گلوکز و انسولین را به یک یا چند وعده غذایی برای افراد عادی و آسیب دیده (دیابتی) شبیه سازی می کند و تخمین پارامتر را با استفاده از sbiofit با استراتژی تابع اجباری انجام می دهد.

پیش زمینه

دالا من (Dalla Man) و همکاران در مقاله خود در سال 2007 الگویی برای پاسخ گلوکز و انسولین انسان بعد از غذا معرفی کردند. این مدل دینامیک سیستم را با استفاده از معادلات دیفرانسیل معمولی توصیف می کند. نویسندگان از مدل خود برای شبیه سازی پاسخ گلوکز و انسولین پس از یک یا چند وعده غذایی، برای افراد عادی و برای انسان های با انواع مختلف اختلالات انسولین استفاده کردند. اختلالات به عنوان مجموعه های متناوب مقادیر پارامترها و شرایط اولیه نشان داده شده است.

معرفی مدل شبیه سازی شده

در این مثال آموزشی متلب مدل SimBiology و m1 برای ترجمه معادلات مدل در دالا من و همکاران (2007) به واکنش ها، قوانین و رویدادها استفاده شده است. سازماندهی مدل به دو بخش صورت گرفته که یکی برای گونه ها و واکنش های مرتبط با گلوکز (به نام ظاهر گلوکز) و دیگری برای گونه ها و واکنش های مرتبط با انسولین (به نام ترشح انسولین) می باشد. همچنین با استفاده از مقادیر پارامتر و شرایط اولیه از معادلات مدل و از جدول 1 و شکل 1 استفاده شده و شامل یک معادله برای نرخ تخلیه معده همانطور که در مقاله Dala Man و همکاران (۲۰۰۶) ارایه شده، می باشد. تنظیم واحدها برای همه گونه ها، محفظه ها و پارامترها همانطور که توسط دالا من و همکاران در مقاله 2007 مشخص شده، انجام گرفته که به مدل SimBiology اجازه می دهد تا با استفاده از تبدیل واحد شبیه سازی شود (توجه داشته باشید که SimBiology از پارامترهای بدون بعد با تنظیم ویژگی ValueUnits آنها بر روی dimensionless پشتیبانی می کند).

پیکربندی TimeUnits نیز باید بر روی ساعت تنظیم گردد، زیرا شبیه سازی ها طی 7 یا 24 ساعت انجام می شود. برای تبدیل گونه ها و پارامترهایی که با وزن بدن در مدل اصلی نرمال شده اند نیز از مبنای 1 کیلوگرم وزن بدن استفاده شده است. با انجام این کار، واحدهای گونه در مقدار یا غلظت همانطور که توسط SimBiology نیاز است، انجام می شود. در این مثال آموزشی اختلالات انسولین در مدل سیمبیولوژی به صورت اشیای متغیر با نام های زیر نشان داده شده است:

دیابت نوع دوم
انسولین کم
واکنش پذیری بالای سلول های بتا
واکنش پذیری کم سلول های بتا
حساسیت به انسولین زیاد

وعده های غذایی در مدل SimBiology به عنوان دوز اشیاء نشان داده شده است:

دوزی به نام Single Meal نشان دهنده یک وعده ۷۸ گرمی گلوکز در آغاز یک شبیه سازی است.
دوزی به نام Daily Life نشان دهنده ارزش غذایی یک روز نسبت به شبیه سازی که از نیمه شب شروع می شود، است. صبحانه ۴۵ گرم گلوکز در ۸ ساعت، زمان شبیه سازی (۸ صبح)، ناهار ۷۰ گرم گلوکز در ۱۲ ساعت (۹ شب)، و شام ۷۰ گرم گلوکز در ۲۰ ساعت (۸ شب) است.

فلوچارت مدل SimBiology در زیر نشان داده شده است:

شروع کد و شبیه سازی

بارگذاری مدل

sbioloadproject('insulindemo', 'm1')

1	sbioloadproject('insulindemo', 'm1')

هشدارهای اطلاعاتی که در طول شبیه سازی ها صادر می شوند را غیرفعال کنید.

warnSettings = warning('off', 'SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless');

1	warnSettings = warning('off', 'SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless');

شبیه سازی پاسخ گلوکز و انسولین برای یک فرد عادی

شی Single Meal dose را انتخاب کنید و خواص آن را نمایش دهید.

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: []

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;

1 2	configset = getconfigset(m1,'active'); configset.StopTime = 7;

نمایش واحدهای زمان شبیه سازی (و واحدهای StopTime).

configset.TimeUnits

1	configset.TimeUnits

ans = 
'hour'

1 2	ans = 'hour'

شبیه سازی یک وعده غذایی برای یک موضوع عادی.

normalMealSim = sbiosimulate(m1, configset, [], mealDose);

1	normalMealSim = sbiosimulate(m1, configset, [], mealDose);

شبیه سازی پاسخ گلوکز و انسولین برای دیابت نوع دوم

دیابت نوع 2 را انتخاب کنید و خواص آن را نمایش دهید.

diabeticVar = sbioselect(m1, 'Name', 'Type 2 diabetic')

1	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

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

شبیه سازی یک وعده غذایی برای یک دیابتی نوع دوم.

diabeticMealSim = sbiosimulate(m1, configset, diabeticVar, mealDose);

1	diabeticMealSim = sbiosimulate(m1, configset, diabeticVar, mealDose);

نتایج را برای مهم ترین خروجی های شبیه سازی مقایسه کنید.

گلوکز پلاسما (species Plasma Glu Conc)
انسولین پلاسما (species Plasma Ins Conc)
تولید درون زا گلوکز (پارامتر Glu Prod)
نرخ ظاهری گلوکز (پارامتر نرخ ظاهری گلوکز)
استفاده از گلوکز (پارامتر 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

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

به غلظت بسیار بالاتر گلوکز و انسولین در پلاسما و همچنین استفاده طولانی مدت از گلوکز و ترشح انسولین توجه کنید.

شبیه سازی یک روز با سه وعده غذایی برای یک سوژه عادی

در این مرحله StopTime شبیه سازی را روی 24 ساعت تنظیم کنید.

configset.StopTime = 24;

1	configset.StopTime = 24;

دوز وعده غذایی روزانه را انتخاب کنید.

dayDose  = sbioselect(m1, 'Name', 'Daily Life');

1	dayDose = sbioselect(m1, 'Name', 'Daily Life');

سه وعده غذایی را برای یک موضوع عادی شبیه سازی کنید.

normalDaySim = sbiosimulate(m1, configset, [], dayDose);

1	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'   )];

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

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

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=β) یا حساسیت کم به سلول‌های بتا (خط فیروزه‌ای نقطه چین، −Ins=β) منجر به افزایش و طولانی‌تر شدن غلظت گلوکز پلاسما می‌شود (ردیف بالای نمودارها). حساسیت کم در یک سیستم می تواند تا حدی با حساسیت بالا در سیستم دیگر جبران شود. برای مثال، حساسیت کم به انسولین و حساسیت بالای سلول های بتا (خط قرمز نقطه چین، −Ins+β) منجر به غلظت نسبتاً طبیعی گلوکز پلاسما (ردیف بالای نمودارها) می شود. با این حال، در این مورد، غلظت انسولین پلاسما به دست آمده بسیار زیاد است (ردیف پایین نمودارها).

روش تخمین پارامتر

نویسندگان به جای تخمین همزمان پارامترها برای کل مدل، تخمین پارامتر را برای زیرسیستم های مختلف مدل با استفاده از یک استراتژی تابع اجباری انجام می دهند. این رویکرد به داده‌های تجربی اضافی برای «ورودی‌های» مدل فرعی نیاز دارد. در طول برازش، داده های ورودی، دینامیک گونه های ورودی را تعیین می کند. (در مدل کامل، دینامیک ورودی ها از معادلات دیفرانسیل تعیین می شود.) در قوانین SimBiology، شما می توانید یک تابع اجباری را به عنوان یک قانون تخصیص مکرر اجرا کنید که مقدار یک گونه یا پارامتری را که به عنوان ورودی برای یک ورودی عمل می کند، کنترل می کند. زیر سیستم مدل در بخش‌های بعدی، از قابلیت‌های برازش پارامتر SimBiology برای اصلاح مقادیر پارامتر گزارش شده توسط نویسندگان استفاده می‌کنیم.

تناسب مدل گوارشی ظاهر گلوکز با استفاده از nlinfit

مدل گوارشی نشان می دهد که چگونه گلوکز موجود در یک وعده غذایی از طریق معده و روده منتقل می شود و سپس به پلاسما جذب می شود. ورودی این زیرسیستم مقدار گلوکز در یک وعده غذایی است و خروجی سرعت ظاهر شدن گلوکز در پلاسما می باشد. با این حال، ما اندازه وعده غذایی را نیز تخمین می زنیم زیرا مقدار گزارش شده توسط نویسندگان با پارامترها و نتایج شبیه سازی ناسازگار است. از آنجا که این ورودی فقط در شروع شبیه سازی رخ می دهد، هیچ تابع اجباری مورد نیاز نیست. تابع sbiofit از تخمین پارامترها در مدل های SimBiology با استفاده از چندین الگوریتم مختلف از MATLAB، Statistics and Machine Learning Toolbox، Optimization Toolbox و Global Optimization Toolbox پشتیبانی می کند. ابتدا پارامترها را با استفاده از تابع Statistics و تابع nlinfit به کمک Machine Learning Toolbox تخمین بزنید.

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

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

پایان تکرارها: تغییر نسبی در SSE کم تر از Options.tolfun می باشد.

برازش داده ها با استفاده از fminunc

اکنون پارامترها را با استفاده از تابع fminunc به کمک Optimization Toolbox تخمین بزنید.

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

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

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')

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')

توجه داشته باشید که اگر اندازه وعده غذایی (دوز) به طور قابل توجهی بزرگتر از آنچه گزارش شده باشد، پارامتر 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')

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

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')

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])

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')

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])

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

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

پایان تکرارها: هنجار نسبی مرحله فعلی کمتر از 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');

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;

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

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

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

تصویر glucose-insulin-response-simbiology-matlab_5633_12 شبیه سازی پاسخ گلوکز و انسولین در انسان عادی و دیابتی با متلب

توجه داشته باشید که افزایش قابل توجه برخی از پارامترها، مانند Vmx، به تناسب بسیار بهتری از غلظت دیررس گلوکز پلاسما اجازه می دهد.

پاک سازی

تنظیمات هشدار را بازیابی کنید.

warning(warnSettings);

1	warning(warnSettings);

نتیجه گیری

SimBiology شامل چندین ویژگی است که پیاده سازی و شبیه سازی یک مدل پیچیده از سیستم گلوکز و انسولین را تسهیل می کند. واکنش ها، رویدادها و قوانین روشی طبیعی را برای توصیف پویایی سیستم فراهم می کنند. تبدیل واحد اجازه می دهد تا گونه ها و پارامترها در واحدهای مناسب مشخص شوند و سازگاری ابعادی مدل را تضمین می کند. اشیاء دوز یک راه ساده برای توصیف ورودی های تکرار شونده به یک مدل هستند، مانند برنامه غذایی روزانه در این مثال. همچنین SimBiology پشتیبانی داخلی برای کارهای تجزیه و تحلیل و شبیه سازی و تخمین پارامتر ارائه می دهد.

نکته پایانی: این مقاله آموزشی متلب از مثال های آماده خود نرم افزار MATLAB است که اگر تولباکس SimBiology ، تولباکس بهینه سازی (Optimization Toolbox) و تولباکس آمار و یادگیری ماشین (Statistics and Machine Learning Toolbox) بر روی متلب شما نصب باشد، می توانید با تایپ دستور زیر در پنجره کامند متلب به این کدها و دستورات دسترسی داشته باشید.

openExample('simbio/insulindemo')

1	openExample('simbio/insulindemo')

منابع

Mathworks
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.