Delay differential equation of glucose insulin model

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

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Stephan]

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