ODE 6 Model equation

%%  Minimal Glucose insulin 6 dimension ode model
% Name: Aniruddha
% Date: 16.09.2019
%  6 model ode differential equation:
%dG(t)/dt=Gin-f2(G(t))-f3(G(t))f4(Ii(t))+f5(x(3)) %G(t)= mass of glucose

%dIp(t)/dt=f1(G(t))-E[(Ip(t)/Vp)-(Ii(t)/Vi)]-Ip(t)/tp  
%Ip(t)=mass of insulin in plasma

%dIi(t)/dt=E[(Ip(t)/Vp)-(Ii(t)/Vi)]-Ip(t)/ti 
%Ii(t)= mass of insulin in intercellulr space

%dx1(t)/dt=3/td(Ip(t)-x1(t)) 
%dx2(t)/dt=3/td(x1(t)-x2(t))
%dx3(t)/dt=3/td(x2(t)-x3(t))
% Mapping G=y(1), Ip=y(2),Ii=y(3),x1=y(4),x2=y(5),x3=y(6)
%% clean up:
clc;
%% parameter values
p.Vp=3; % plasma insulin distribution volume
p.Vi=11;% Effective volume of intercellular space
p.Vg=10;
p.E=0.2;%Diffusion transfer rate
p.tp=6; % insulin degradation  time constant in plasma
p.ti=100;% insulin degradation  time constant in intercellular spaces
p.td=36;
p.Rm=210;
p.a1=300;
p.C1=2000;
p.Ub=72;
p.C2=144;
p.C3=1000;
p.U0=40;
p.Um=940;
p.b8=1.77;
p.C4=80;
p.Rg=180;
p.a2=0.29;
p.C5=26;
p.Gin=3;
%%  Solve the ode:
tspan=[0 200];
y0=[0;0;0;0;0;0]';
[T,Y]=ode45(@model,tspan,y0,[],p);
% Rewrite the state variables:
G=Y(:,1);
Ip=Y(:,2);
Ii=Y(:,3);
x1=Y(:,4);
x2=Y(:,5);
x3=Y(:,6);
%% Illustration:
 figure(1)
  subplot(1,6,1)
  plot(T,G)
  xlabel('time [min]')
  ylabel('mass of glucose ')
  grid on;
  
  subplot(1,6,2)
  plot(T,Ip)
  xlabel('time [min]')
  ylabel('mass of insulin in plasma ')
  grid on;
  
  subplot(1,6,3)
  plot(T,Ii)
  xlabel('time [min]')
  ylabel('mass of insulin in intercellular spaces ')
  grid on;
  
 subplot(1,6,4)
  plot(T,x1)
  xlabel('time [min]')
  ylabel('x1 ')
  grid on;
  
  subplot(1,6,5)
  plot(T,x2)
  xlabel('time [min]')
  ylabel('x2 ')
  grid on;
  
  subplot(1,6,6)
  plot(T,x3)
  xlabel('time [min]')
  ylabel('x3 ')
  grid on;
  
  
  
  
  
  function ds=model(t,y,p)
G=y(1);
Ip=y(2);
Ii=y(3);
x1=y(4);
x2=y(5);
x3=y(6);
%% 6 model equation:
f1(G)=p.Rm/(1+exp((p.C1-G/p.Vg)/p.a1));
f2(G)=(p.Ub)*(1-exp(-G/p.C2*p.Vg));
f3(G)=G/p.C3*p.Vg;
f4(Ii)=p.U0 + (p.Um-p.U0)/1+exp(-p.b8(Ii/p.C4(1/Vi + 1/Eti)));
f5(x3)= p.Rg/1+exp(p.a2(x3/p.Vp -p.C5));
dG=  p.Gin-f2(G)-f3(G)*f4(Ii)+f5(x3);  %%ODE for mass of glucose
dIp=  f1(G)-p.E*((Ip/p.Vp)-(Ii/p.Vi))-Ip/p.tp; %%ODE for mass of insulin in plasma
dIi=p.E*((Ip/p.Vp)-(Ii/p.Vi))-Ip/p.ti;   %% ODE for mass of insulin in intercellular spaces
dx1=(3/p.td)*(Ip-x1);
dx2=(3/p.td)*(x1-x2);
dx3=(3/p.td)*(x2-x3);
ds=[dG,dIp,dIi,dx1,dx2,dx3]';
end
  
  
  
  
  

I am getting error in ode45

Answer 1

[Star Strider]

See the documentation section on Anonymous Functions to understand how to write them correctly.

Your code is also missing several multiplication operators (if I am guessing correctly). MATLAB does not recognise implicit multiplication. (I supplied them.)

There are also missing constants in your ‘p’ structure. I supplied them to test to code, however you need to supply the correct values for them.

Change to these (and supply all the necessary constants):

f1 = @(G) p.Rm/(1+exp((p.C1-G/p.Vg)/p.a1));
f2 = @(G) (p.Ub)*(1-exp(-G/p.C2*p.Vg));
f3 = @(G) G/p.C3*p.Vg;
f4 = @(Ii) p.U0 + (p.Um-p.U0)/1+exp(-p.b8*(Ii/p.C4*(1/p.Vi + 1/p.Eti)));
f5 = @(x3) p.Rg/1+exp(p.a2*(x3/p.Vp -p.C5)); 

and your code runs with out error.

You must determine if it is giving the correct results.

[Star Strider]

Aniruddh Deshmukh’s Answer moved here:

It is still showing error on ode45 line

[Star Strider]

When I changed your anonymous functions to the ones that I posted, and suppliued values and assignments for the missing constants, the code ran without error.

What error is the ode45 call throwing when you run it?

[Star Strider]

You have not provided enough information for me to respond.

Change ‘f1’ though ‘f5’ in the code you posted to the code I provided, supply the missing constants in your ‘p’ structure, and the code should run without error, and do all the plots.

[Star Strider]

The function declaration for ‘Modeleqn’ needs to change to:

function df = Modeleqn(t,y,p)

so the call to it in the ode45 call needs to change to:

[T,Y]=ode45(@(t,y)Modeleqn(t,y,p),tspan,y0);

With those changes, your code ran without error when I tested it.