Delay differential equation of glucose insulin model
6 views (last 30 days)
Show older comments
%% Minimal Glucose dde model
% Name: Aniruddha
% Date: 16.10.2019
% delay differential equation:
% DDE1
% dG(t)/dt= Gin-f2(G(t))-f3(G(t))f4(I(t))+f5(I(t-T2))
% d(I)/dt=f1(G(t))-I(t)/t1
% DDE2
%dG(t)/dt= Gin-f2(G(t))-f3(G(t))f4(I(t))+f5(I(t-T2))
% d(I(t))/dt= Qf1(G(t))-I(t)/t1 +(1-Q)f1(G(t-T1)
%f1(G)=209/(1+e^(-G/300V3 +6.6)
%f2(G)=72(1-e^(-G/144V3))
%f3(G)=0.01G/V3
%f4(I)=[90/(1+e^(-1.772log(I/V1) + 7.76))] +4
%f5(I)=180/(1+e^(0.29I/V1 -7.5))
%Mapping: G(t)=y(1), I(t)=y(2) I1(t)=y(3)
%% clean up:
clc;clear all;
%% parameter values:
p.V1=3;
p.Gin=216
p.Q=0.5;
p.t1=6;
p.V3=10;
p.T2=50;
p.Eg=180;
p.T1=5;
lags=[p.T1 p.T2];
tf=5;
sol=dde23(@ddemodel,lags,[0,tf]);
t=linspace(0,tf,100);
y=ddeval(sol,t);
figure;
plot(t,y)
function dX=ddemodel(t,y,p)
G=y(1);
I=y(2);
I1=y(3);
f1=@(G) 209/(1+e^(-G/300*p.V3 +6.6));
f2=@(G) 72*(1-e^(-G/144*p.V3));
f3=@(G) 0.01*G/p.V3;
f4=@(I) [90/(1+e^(-1.772*log(I/p.V1) + 7.76))] +4;
f5=@(I) 180/(1+e^(0.29*I/p.V1 -7.5));
dG(t)/dt= p.Gin-f2(G)-f3(G)*f4(I)+f5(I(t-p.T2));
d(I)/dt=f1(G)-I(t)/p.t1;
d(I1(t))/dt= p.Q*f1(G)-I(t)/p.t1 +(1-p.Q)*f1(G(t-T1))
df=[dG dI dI1]';
end
Answers (1)
Stephan
on 27 Oct 2019
Another try, after reading a little bit about dde's:
[t,y] = delayed;
% plot results
subplot(3,1,1)
plot(t,y(1,:),'LineWidth',2)
title('G(t)')
subplot(3,1,2)
plot(t,y(2,:),'LineWidth',2)
title('I(t)')
subplot(3,1,3)
plot(t,y(3,:),'LineWidth',2)
title('I1(t)')
%% solve system
function [t,y] = delayed
% parameter values:
p.V1=3;
p.Gin=216;
p.Q=0.5;
p.t1=6;
p.V3=10;
p.T2=50;
p.Eg=180;
p.T1=5;
lags=[p.T1 p.T2];
tf=5;
sol=dde23(@ddemodel,lags,zeros(3,1),[0,tf]);
t=linspace(0,tf,100);
y=deval(sol,t);
function df=ddemodel(~,y,lags)
df = zeros(3,1);
% functions
G=y(1);
I=y(2);
I1=y(3);
f1=@(G) 209/(1+exp(-G/300*p.V3 +6.6));
f2=@(G) 72*(1-exp(-G/144*p.V3));
f3=@(G) 0.01*G/p.V3;
f4=@(I) 90/(1+exp(-1.772*log(I/p.V1) + 7.76)) +4;
f5=@(I) 180/(1+exp(0.29*I/p.V1 -7.5));
df(1) = p.Gin-f2(G)-f3(G)*f4(I)+f5(lags(2));
df(2) =f1(G)-I/p.t1;
df(3) = p.Q*f1(G)-I/p.t1 +(1-p.Q)*f1(lags(1));
end
end
gives the following result:
I think this may be a bit more helpful. If not let me know - but also let me know if it helps...
0 Comments
See Also
Categories
Find more on Model Predictive Control Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!