function insulin_ODE
param = [
0.0477985900779305
7.78160761103111	
1.61147523779118	
0.00792717614041237	
3.88248960751442e-05	
0.284425974469310	
9.93311225447353	
1.00000000282413e-05	
1.00277786609339e-05 
0.363030286526969	  
0.408127912886804	
0.0948960328385619  
0.00919578911309774 
7.70618517548016
46.2222520609226
4.33812187331668
0.419682384304397
0.124330492920416
0.00752464611370891
1.95497593092361	
0.00105342379833469
0.00145811601430322
0.0959163701057347
14.9913345914433
2.77699437841580
2.97538195830595
0.924599254651007	
10.5641549514879
0.151729454062921	
4.59698375609346	
4.86146143292690
0.0549645715836672
38.5402920478707
83.8672192461257
0.111097316860158	
0.0363622452066626
]';
y0 = [param(15);0;0;0;0;0;0;0;0;param(32)*param(33)/param(31);param(33)];
time_steady = [0,10000];

Ins = 0.01;  % basal Ins conc.
[t,time_course0]= ode15s(@(t,y) ODEs(t,y,param,Ins),time_steady,y0);
Ins = 1;  % Ins conc.
[t,time_course]= ode15s(@(t,y) ODEs(t,y,param,Ins),0:480,time_course0(end,:));

figure(1)
plot(t,time_course(:,5)) %pAKT
hold on
figure(2)
plot(t,param(35)*time_course(:,8)) %pGSK
hold on
figure(3)
plot(t,param(34)*time_course(:,6)) %pS6K
hold on
figure(4)
plot(t,param(36)*time_course(:,10)) %pG6Pase
hold on
end

function dydt = ODEs(t,y,param,Ins)
reaction1 = (param(15)-y(1))*param(1);          % sIRS <-> IRS
reaction2 = y(1).*Ins*param(2) - y(2)*param(3); % IRS + Ins <-> p1IRS
reaction3 = y(2)*param(4);                      % p1IRS -> (Degradation)
reaction4 = y(2).*y(11).*param(5);              % p1IRS + pmTOR -> p1p2IRS + pmTOR 
reaction5 = y(3).*param(6);                     % p1p2IRS -> p1IRS
reaction6 = y(1).*y(11).*param(7);              % IRS + pmTOR  ->  p2IRS + pmTOR 
reaction7 = y(4).*param(8);                     % p2IRS -> IRS
reaction8 = y(4).*param(9);                     % p2IRS -> (Degradation)
reaction9 = y(4).*Ins*param(10)-y(3)*param(11); % p2IRS + Ins <-> p1p2IRS
reaction10 = y(3).*param(12);                   %p1p2IRS -> (Degradation)
reaction11 = y(2).*(param(16)-y(5))*param(13);  % Akt  -> pAkt
reaction12 = y(5).*param(14);                   % pAkt  -> Akt

reaction13 = y(5).*(param(23)-y(11))*param(17); % pAkt + mTOR  -> pAkt + pmTOR
reaction14 = y(11).*param(18);                  % pmTOR  -> mTOR
reaction15 = y(11).*(param(25)-y(6))*param(19); % pmTOR + S6K -> pmTOR +pS6K
reaction16 = y(6).*y(7)*param(20);              % pS6K + pX -> S6K+pX
reaction17 = y(11).*(param(24)-y(7))*param(21); % pmTOR + X -> pmTOR +pX
reaction18 = y(7).*param(22);                   % pX -> X

reaction19 = y(5).*(param(28)-y(8))*param(26);  % pAkt + GSK3b -> pAkt + pGSK3b
reaction20 = y(8)*param(27);                    % pGSK3b -> GSK3b

reaction21 = ((param(32)*param(33)/param(31))-y(9))*param(29); % pFox01 -> Fox01
reaction22 = y(5).*y(9)*param(30);              % pAkt + Fox01 -> pAkt + pFox01
reaction23 = y(9)*param(31);                    % Fox01 -> G6Pase
reaction24 = y(10)*param(32);                   % G6Pase -> degradation
%%
dydt(1,:) = reaction1 - reaction2 - reaction6 + reaction7;  % IRS
dydt(2,:) = reaction2 - reaction3 - reaction4 + reaction5;  % p1IRS
dydt(3,:) = reaction4 - reaction5 + reaction9 - reaction10; % p1p2IRS
dydt(4,:) = reaction6 - reaction7 - reaction8 - reaction9;  % p2IRS
dydt(5,:) = reaction11 - reaction12;  % pAkt

dydt(6,:) = reaction15 - reaction16;  % pS6K
dydt(7,:) = reaction17 - reaction18;  % pX
dydt(11,:) = reaction13 - reaction14; % pmTOR

dydt(8,:) = reaction19 - reaction20;  % pGSK3b

dydt(9,:) = reaction21 - reaction22;  % Fox01
dydt(10,:) = reaction23 - reaction24; % G6Pase
end