%-------------------------------------------------------------------------%
% Author          :   Athil Althaf A Z, Urvi Mishra and Seeraja Konnur    %
% Date            :   15-June-2018                                        %
% Version         :   2.0                                                 %
% Description     :   Pancreatic beta cell Second Messenger interactions  %
% Copyright       :   IISER Kolkata,NMIMS Mumbai,IBAB Bangalore           %
%-------------------------------------------------------------------------%
Ts = 6000;              % Total simulation time (s)
dt = .01;                 % Timestep (s)
t = 0:dt:Ts;            % Defining the array for time
fNc = 1;          % Conversion factor for converting microm^-2 to microM (microM micro m^2)   Doubt
fNc1 = 21.1;
a=zeros(1,length(t));
t1 = 3000;
t2 = 5000;
t3 = 5000;
t4 = 3000;
tpde=4000; % same as above
t6 = 5000;              % AM3 input time 
t7 = 1000;              % AR7 input time
t9 = 3000;              % kATP Blocker time
t10 =5000;              % VmPLP input time
%-------------------------------------------------------------------------%
%                   Glucose Metabolism                                    %
%-------------------------------------------------------------------------%
Cac1=0.35;
ATD_in = 1.5;           % Initial value of ATD eq.2
tATD = 100;             % Time constant (s) eq.2
ATDm = 32;              % Saturated ADT value eq.4
KGF = 5200;             % Coefficient of half saturated glucose (microM) eq.4
kGFr = 7.4;             % Coefficient of coversion FFA to glucose eq.4
hgla = 5;               % Hill coefficient eq.4
FFA = 0;                % Free Fatty acid concentration (microM) eq.4
AT = 4000;              % Total free ATP and free ADP concentration (microM) eq.5
Glui = 3000;            % Initial low glucose concentration (microM) eq.7
Glu1 = 8000;            % Increased Glucose concentration (microM) eq.7
Glu = Glui*ones(1,length(t)); % Initial Glucose level array
ATD = a;   % Defining the Array for ATP -ADP ratio eq.2
ATD(1) = ATD_in;            % Initializing the ATD array
ADPf = a;  % Initializing the free ADP array
ATP = a;  % Initializing the free ATP array
ADPf(1) = AT/(1+ATD(1));   % Initializing first element of ADPf
ATP(1) = AT - ADPf(1);     % Initializing first element of ATP
%-------------------------------------------------------------------------%
%                      Membrane Potential,channels and Pumps              %
%-------------------------------------------------------------------------%
Cm = 6158;              % Membrane capacitance (fF) eq.8
Vp_in = -62;            % Initial membrane potential (mV) eq.8
gmKATPi = 24000;        % Maximum conductance for IKATP (pS) eq.9
gmKATP1 = 9000;       % Action of KATP blocker at t9 (pS)  eq.14
gmKATP = gmKATPi*ones(1,length(t));       % Normal level of KATP conductance   eq.14
EK = -75;               % Reversal potential of pottassium (mV) eq.9
KKPI2 = 1125*fNc1;       % Activation constant (microM  converted by fNc) eq.10     Doubt
kdd = 17;               % Dissosiation coefficient (microM) eq.11
ktd = 26;               % Dissosiation coefficient (microM) eq.11
ktt = 50;               % Dissosiation coefficient (microM) eq.11
kKPKa = 1;              % Inhibition constant (microM^-1) eq.13                  doubt
gmKr = 45000;           % Maximum conductance for IKr (pS) eq.15
VdKr = -9;              % Half activaton potential (mV) eq.16
kdKr = 5;               % Slope of half activation potential (mV) eq.16
fKri = 1;               % Inactivating gating variable eq.17
gmCa = 900;             % Maximum conductance for ICa (pS) eq.18
ECa = 100;              % Reversal potential for calcium (mV) eq.18
VdCa = -19;             % Half activation potential (mV) eq.19
kdCa = 9.5;             % Slope of half activation potential (mV) eq.19
kdCap = .01;            % Coefficient of constitutive channel activity eq.19
VfCa = -9;              % Half inactivation potential (mV) eq.20
kfCa = 8;               % Slope of half activation potential (mV) eq.20
gmSOC = 10;             % Maximum whole cell conductance (pS) eq.21
KNS = 200;              % Calcium ER inhibition constant (microM) eq.22
PmCa = 6000;            % Maximum IPCa current (fA) eq.23
KpCa = .2;              % Cyosolic Calcium activation constant eq.23
ENa = 70;               % Reversal potential of sodium (mV) eq.24
gmNb = 10;              % Permanent conductance for INab (pS) eq.25
gmNM3 = 0;              % Maximum conductance of NALCN channels (pS) eq.25
Vp = a;      % Defining the Array for membrane potential eq.8
Vp(1) = Vp_in;          % Adding initial value of the Membrane potential


MgADPf = a; % Initializing array for MgADP eq.12
OKATP = a;  % Initializing array for OKATP eq.11
IKATP = a;  % Initializing array for IKATP eq.9
fKPI = a;   % Initialising array for fKPI eq.9
dKri = a;   % Initializing array for dKri eq.16
IKr = a;    % Initializing array for IKr eq.15
dCai = a;   % Initializing array for dCai eq.19
fCai = a;   % Initializing array for fCai eq.19
ICa = a;    % Initializing array for ICa eq.19
fSOC = a;   % Initializing array for fSOC eq.22
ISOC = a;   % Initializing array for ISOC eq.21
IPCa = a;   % Initializing array for IPCa eq.23
INab = a;   % Initializing array for INab eq.24
gNab =a;    % Initializing array for gNab eq.24

%-------------------------------------------------------------------------%
%                      Calcium Dynamics                                   %
%-------------------------------------------------------------------------%
Cac_in = .09;           % Initial cytosolic calcium concentration (microM) eq.86
ksg = .00001;           % Coefficient of calcium sequestration rate (s^-1) eq.86
CaER_in = 160;          % Initial ER calcium concentration (microM) eq.87
PCaER = .026;              % Maximum SERCA pump rate (microMs^-1) eq.88
Kser = .5;              % Calcium activation constant (microM) eq.88
kmIP = 7;               % Maximum permeability of IP3 activated channel (microMs^-1) eq.89
kleak = .002;           % Coefficient of passive leak from ER (s^-1) eq.89
KRPCa = .35;            % Cytosolic calcium activation constant (microM) eq.90
dinact = .4;            % Cytosolic inactivation constant (microM) eq.90 %% obtained from ref.203
KIP3 = 3.2;             % IP3 activation constant (microM) eq.90
KIP3R = .3;             % PKA activation coefficient 90
Cac = a;     % Defining array for cytosolic calcium eq.86
fcf = .01;              % Fraction of free calcium in cytoplasm eq.87
fER = .03;              % Fraction of free calcium at ER eq.87
F = 96.487;             % Faraday constant CmicroM^-1 eq.86 
Cac(1) = Cac_in;        % Adding initial value for cytosolic calcium
CaER  = a;   % Defining array for ER calcium eq.87
VER = 280;                    % Effective volume of ER microm^3 eq.87
Vc = 764;                     % Effective volume of cytoplasm eq.86
CaER(1) = CaER_in;      % Adding initial value for ER calcium
Jser = a;    % Initializing Jser array eq.88
Jrel = a;    % Initializing Jrel array eq.89
fIPKA = a;   % Initializing fIPKA array eq.90
PRIP3 = a;   % Initializing Jrel array eq.90



%-----------------------------------------------49--------------------------------------------------------%

AM3i = 0.0022; %Ln (microM)
RG6_in = 0.01*fNc; %RnGN(#per micro m^2)
ReM3d_in = 0.1*fNc; %Rnd(#per micro m^2)
LRG6_in = 0.01*fNc; %LnRnGn(#per micro m^2)
GaT6_in = 0.01*fNc; %Gan(GTP)(#per micro m^2)
GaD6_in = 0.01*fNc; %Gan(GDT)(#per micro m^2)
PLCM3_in = 0.1*fNc; %Gan(GTP)En(#per micro m^2)
ReM3t = 2*fNc; %Rnt(#per micro m^2)
G6t = 20*fNc; %Gnt(#per micro m^2)
KAM3 = 0.22; %K(M)(microM)
R61 = 10; %k1n(microM)
R62 = 0.68/fNc; %k2n(per# microm^2 per s)
R62r = 6.8; %k2nr(per s)
R63 = 0.65; %k3n(per s)
R64 = 1/fNc; %k4n(per# microm^2 per s)
R65 = 0.026; %k5n(per s)
R66 = 0.4; %k6n(per s)
R67 = 1/fNc; %k7n(per# microm^2 per s)
R68 = 0.007; %k8n(per s)
R69 = 0.00105; %k9n(per s)
AM31 = 2.2;

%----------------------------------------------------------50-----------------------------------------------%

AR7i = 0.0506; %Ln (microM)
ReR40d_in = 0.1*fNc; %Rnd(#per micro m^2)
LR7_in = 0.01*fNc; %LnRn(#per micro m^2)
LRG7_in = 0.01*fNc; %LnRnGn(#per micro m^2)
GaT7_in = 0.01*fNc; %Gan(GTP)(#per micro m^2)
GaD7_in = 0.01*fNc; %Gan(GDT)(#per micro m^2)
PLCR40_in = 0.1*fNc; %Gan(GTP)En(#per micro m^2)
ReR40t = 2*fNc; %Rnt(#per micro m^2)
G7t = 20*fNc; %Gnt(#per micro m^2)
PLCpt = 3*fNc; %Ent(#per micro m^2)
KAR7 = 4.6; %K(M)(microM)
R71 = 7; %k1n(microM)
R72 = 1/fNc; %k2n(per# microm^2 per s)
R72r = 0.68; %k2nr(per s)
R73 = 0.7; %k3n(per s)
R74 = 1/fNc; %k4n(per# microm^2 per s)
R75 = 0.026; %k5n(per s)
R76 = 0.4; %k6n(per s)
R77 = 1/fNc; %k7n(per# microm^2 per s)
R78 = 5*10^-5; %k8n(per s)
R79 = 2.83*10^-4; %k9n(per s)
PLCR40 = a;
ReR40 = a;
GaD7 = a;
GaT7 = a;
LRG7 = a;
G7 = a;
Gbg7 = a;
AR7 = AR7i*ones(1,length(t));
ReR40d = a;
LR7 = a;
PLCR40(1) = PLCR40_in;
GaD7(1) = GaD7_in ;
GaT7(1) = GaT7_in ;
LRG7(1) = LRG7_in ;
%AR7(1) = AR7i;
ReR40d(1) = ReR40d_in ;
LR7(1) = LR7_in ;
AR71 = 82;

%--------------------------------------------------73-86-----------------------------------------------------%

P4P_in = 4000*fNc; % #per micro m^2
kPI = 0.0015; %per s
kPIr = 0.006; %per s
kP4P = 0.02; %per s
kP4Pr = 0.014; %per s
PI = 140000*fNc; % #per micro m^2
KCaPI = 0.3;
KP4PK = 0.5;
kcpi = 0.2;
PIP2_in = 4200*fNc; %#per micro m^2
KPIP2 = 2370*fNc; % #per micro m^2 converted to micro M
kpPL = 15;% micro mol per s
PKCa_in = 0.1;
KCaPL = 0.4; % micro M
KCCaPL = 0.2; % micro M
IP3_in = 1; % micro M
kdIP3  = 0.04; % per s
DAG_in = 23*fNc; % #per micro m^2
kdDAG = 0.05; % per s
kPKC = 3*10^-6; % per s
kPKCr = 0.0034;% per s
PKCt = 1;

%-----------------------------------------------------M3R6-------------------------------------%

GaD6 = a;
PLCM3 = a;
GaT6 = a;
ReM3d = a;
LRG6 = a;
RG6 = a;
Gbg6 = a;
G6 = a;
AM3 = AM3i*ones(1,length(t));
GaD6(1) = GaD6_in;
PLCM3(1) = PLCM3_in;
GaT6(1) = GaT6_in;
ReM3d(1) = ReM3d_in;
LRG6(1) = LRG6_in;
RG6(1) = RG6_in;
%AM3(1) = AM3i;

P4P = a;
PIP2 = a;
IP3 = a;
DAG = a;
PKCa = a;
PLCpa = a;
fPIP2 = a;
PKCi = a;
fPI = a;
PLCpf = a;

VPLP = a;
VPLC = a;
VPL = a;
P4P(1)=P4P_in;
PIP2(1)=PIP2_in;
Cac_in = 0.09;% micro M
PKCa(1)= PKCa_in;
VmPLPi = 700; % micro mol per s
VmPLP1 = 100000;  % micro mol per s
VmPLP =VmPLPi*ones(1,length(t));
VmPLC = 50; % micro mol per s
IP3(1)=IP3_in;
DAG(1) = DAG_in; 
%----------------------------------incretin-------------------------------%
% for GLP1
GLP1i=(3.2)*(10^-7);
GLP11=(6.2)*(10^-4);
ReGLd_in=0.1*fNc;
LR1_in=0.01*fNc;
LRG1_in=0.01*fNc;
GaT1_in=0.01*fNc;
GaD1_in=0.01*fNc;
ACGLP_in=0.3*fNc;
ReGLt=2*fNc;
G1t=20*fNc;
ACpt=3*fNc;
KGLP1=(3.1)*(10^-5);
R11=7;
R12=1/fNc;
R12r=0.68;
R13=0.7;
R14=1/fNc;
R15=0.026;
R16=0.4;
R17=1/fNc;
R18=0.00005;
R19=0.000283;

LR1=a;
LR1(1)=LR1_in;

LRG1=a;
LRG1(1)=LRG1_in;

ReGLd=a;
ReGLd(1)=ReGLd_in;

GaT1=a;
GaT1(1)=GaT1_in;

ACGLP=a;

ACGLP(1)=ACGLP_in;

GaD1=a;
GaD1(1)=GaD1_in;

GLP1=GLP1i*ones(1,length(t));
ACpf=a;


ReGL=a;
G1=a;
Gbg1=a;

%.....for GIP...%

GIPi=(1.37)*(10^-6);
GIP1=(3.42)*(10^-3);
ReGld_in=0.1*fNc;
LR2_in=0.01*fNc;
LRG2_in=0.01*fNc;
GaT2_in=0.01*fNc;
GaD2_in=0.01*fNc;
ACGIP_in=0.3*fNc;
ReGlt=2.5*fNc;
G2t=25*fNc;
KGIP=(1.71)*(10^-4);
R21=7;
R22=1/fNc;
R22r=0.68;
R23=0.7;
R24=1/fNc;
R25=0.026;
R26=0.4;
R27=1/fNc;
R28=0.00005;
R29=0.000283;
%...for catecholamines....%
AR3i=0.002;
AR31=3.0;
ReARd_in=0.1*fNc;
LR3_in=0.01*fNc;
LRG3_in=0.01*fNc;
GaT3_in=0.01*fNc;
GaD3_in=0.01*fNc;
ACAR_in=0.3*fNc;
ReARt=3*fNc;
G3t=20*fNc;
KAR3=0.3;
R31=7;
R32=1/fNc;
R32r=0.68;
R33=0.7;
R34=1/fNc;
R35=0.026;
R36=0.4;
R37=1/fNc;
R38=0.00005;
R39=0.000283;

LR2=a;
LR2(1)=LR2_in;

LRG2=a;
LRG2(1)=LRG2_in;

ReGld=a;
ReGld(1)=ReGld_in;

GaT2=a;
GaT2(1)=GaT2_in;

ACGIP=a;

ACGIP(1)=ACGIP_in;

GaD2=a;
GaD2(1)=GaD2_in;

GIP=GIPi*ones(1,length(t));



ReGl=a;
G2=a;
Gbg2=a;

%...for catev...%

LR3=a;
LR3(1)=LR3_in;

LRG3=a;
LRG3(1)=LRG3_in;

ReARd=a;
ReARd(1)=ReARd_in;

GaT3=a;
GaT3(1)=GaT3_in;

ACAR=a;

ACAR(1)=ACAR_in;

GaD3=a;
GaD3(1)=GaD3_in;

AR3=AR3i*ones(1,length(t));



ReAR=a;
G3=a;
Gbg3=a;
%...calmodulin...%

CaCaM_in=0.42; %uM;
% 10^3 factor is not considered here!!!
k1f=(2.3); % uM-1 s-1
k1b=(2.4); % ms-1; 
k2f=(2.3); % (units:uM-1 s-1) (forward rate constant eq 52)
k2b=(2.4); %(units: ms-1)  (backward rate constant eq 52)
k3f=(160); %(units:uM-1 s-1)  (forward rate constant eq 53)
k3b=(405); %(units:ms-1)    (backward rate constant eq 53)
k4f=(160) ; %(units:uM-1 s-1)  (forward rate constant eq 54)
k4b=(405); %(units:ms-1)    (backward rate constant eq 54)
CaM0=11.25; % (units:uM)          (total amount of calmodulin eq 55)

%...for cAMP pka Epac.....%
kda=0.01; %(units:  (kda is the coefficient of PDE-independent cAMP degradation)
Vmcam=2; % (units: umol/s) (maximum ACp acitivity)
Kpcam=0.348; %  (Units:uM)(CaM(active) activation constant)
KNCa=75; % (Units:uM)([Ca2+]c inhibition constant.)
VmACc=0.2 ;%  (Units:umol/s)(maximum ACc activity)
KmAACS=1030; %  (Units:uM) ATP activation constant
KmCACS=0.5 ;%  (Units:uM) [Ca2+]c activation constant
kACS=0.01 ; %  (Units:umol/s)coefficient of the constitutive ACc activity.
kipdei=1; % (activation or inhibition coefficient)
kipde1=-0.4; % don't know value dummy input
Kdpe=0.348; %  (Units:uM)CaMa activation constant
Kpde=3; %  (Units:uM)cAMP activation constant
PKAa_in=0.24; % conc of PKA;
kak=1; % scaling factor,
tpka=900; %  (Units:s)time constant;
Kpcm=2.9; %  (Units:uM)cAMP activation coefficient
hpca=1.4; % Hill coefficient
EPa_in=0.01;
tep=900; %  (Units:s)time constant
Kmep=20.2; %  (Units:uM)cAMP activation constant for Epac,
hce=2; %  Hill coefficient


CaCaM=a;
CaCaM(1)=CaCaM_in;
Ca2CaM=a;
Ca3CaM=a;
Ca4CaM=a;
CaM=a;
CaMa=a;

%...cAMP..%
ACpa=a;
ACpar=a;
facca=a;
VACp=a;
VACc=a;
Vpde=a;
VAC=a;
cAMP=a;
PKAi=a;
PKAa=a;
PKAa(1)=PKAa_in;
fKPKa = a;
EPi=a;
EPa=a;
EPa(1)=EPa_in;
fpca=a;
fEcA=a;
kipde=kipdei*ones(1,length(t));
Vgpde=0.04*ones(1,length(t));
Vcpde=1.4*ones(1,length(t));
%-------------------------------------------------------------------------%
%                      Modeling Part                                      %
%-------------------------------------------------------------------------%
for i = 1:length(t)-1                     
    if t(i)>t1
        % Place for inputs and dummy variables
        % Dummy inputs mimicks the behavior of those variables
        Glu(i) =  Glu1;          % Increasing the Glucose concentration eq.7 
        Cac(i) = Cac1; 
        
    end
    if t(i)>t2
      GLP1(i)=GLP11;
    end
%     if t(i)>t3
%           GIP(i)=GIP1;
%     end
%     if t(i)>t4
%         AR3(i)=AR31;
%     end
%     if t(i)> t6
%          AM3(i) = AM31;
%     end
%     if t(i)> t7;
%          AR7(i) = AR71;
%    end
%    if t(i)>t9                       % Blocking KATP eq.14
%         gmKATP(i) = gmKATP1;
%     end
%    Membrane potential part
    if t(i)> t10;
        VmPLP(i+1) = VmPLP(i+1) + VmPLPi;
    end
    if t(i)>tpde
    kipde(i+1)=kipde(i+1)+kipde1;
    end
           
    % Place for solving the differential equations and the allied parameter

    % Glucose metabolism part
    ATD0 = ATDm*(Glu(i)+kGFr*FFA)^hgla/((Glu(i)+kGFr*FFA)^hgla+KGF^hgla);     % eq.4   
    %dATD(i) = dt*(ATD0 - ATD(i))/tATD;                                  % eq.2
    %ATD(i) = ATD(i) + dATD(i);      % Appending the ATD by euler method eq.2
    ATD(i+1) = ATD(i) + dt*(ATD0 - ATD(i))/tATD;
    ADPf(i+1) = AT/(1+ATD(i));        % eq.6
    ATP(i+1) = AT - ADPf(i);          % eq.5
 
    fKPI(i) = PIP2(i)/(PIP2(i)+KKPI2);       % PIP2 is dummy eq.10  %% Membrane Mechanism
    fKPKa(i) = 1/(1+kKPKa*PKAa(i));       % PKAa is dummy eq.13
    MgADPf(i) = 0.055*fKPKa(i)*ADPf(i);                 % eq.12 
    OKATP(i) = (.088*(1+2*MgADPf(i)/kdd)+.89*(MgADPf(i)/kdd)^2)/((1+MgADPf(i)/kdd)^2*(1+.45*MgADPf(i)/ktd +ATP(i)/ktt));    % eq.11
    IKATP(i) = gmKATP(i)*fKPI(i)*OKATP(i)*(Vp(i)-EK);      % eq.9
    dKri(i) = 1/(1+exp((VdKr - Vp(i))/kdKr));        % eq.16
    IKr(i)  = gmKr*dKri(i)^2*fKri*(Vp(i)-EK);        % eq.15
    dCai(i) = 1/(1+exp((VdCa-Vp(i))/kdCa))+kdCap;    % eq.19
    fCai(i) = 1/(1+exp(-(VfCa-Vp(i))/kfCa));         % eq.20
    ICa(i) = gmCa*dCai(i)*fCai(i)*(Vp(i)-ECa);       % eq.18
    fSOC(i) = KNS^4/(KNS^4+CaER(i)^4);               % eq.22
    ISOC(i) = gmSOC*fSOC(i)*(Vp(i)-ECa);             % eq.21
    IPCa(i) = PmCa*Cac(i)^2/(KpCa^2+Cac(i)^2);       % eq.23
    gNab(i) = gmNb + gmNM3*LRG6(i);                     % eq.24
    INab(i) = gNab(i)*(Vp(i)-ENa);                   % eq.24  
%     dVp(i) = -dt*(IKATP(i)+IKr(i)+ICa(i)+ISOC(i)+IPCa(i)+INab(i))/Cm;   % eq.8
    
    Vp(i+1) = Vp(i)-dt*(IKATP(i)+IKr(i)+ICa(i)+ISOC(i)+IPCa(i)+INab(i))/Cm;  % eq.8                       % Appending the Vp array
    % Part for Calcium dynamics
    Jser(i) = PCaER*Cac(i)^2/(Kser^2 + Cac(i)^2);    % eq.88
    fIPKA(i) = 1/(PKAa(i) + KIP3R);                      % eq.90
    PRIP3(i) = (Cac(i)/(KRPCa+Cac(i)))^3*(IP3(i)/(fIPKA(i)*KIP3+IP3(i)))^3*(dinact/(dinact+Cac(i)))^3; %  eq.90
    Jrel(i) = (kmIP*PRIP3(i) + kleak)*(CaER(i) - Cac(i));      % eq.89
%    Cac(i+1) = Cac(i) + dt*(fcf*((-ICa(i)-ISOC(i)-2*IPCa(i))/(2*F*Vc)-Jser(i)+Jrel(i)/Vc) - ksg*Cac(i)); % eq.86 
    CaER(i+1)= CaER(i) + dt*(fER*(Jser(i)*Vc-Jrel(i))/VER);      % eq.87
    
    ACpf(i) = ACpt- ACGLP(i)-ACGIP(i)-ACAR(i);
    ReGL(i)=ReGLt-ReGLd(i)-LR1(i)-LRG1(i);
    G1(i)=G1t-LRG1(i)-GaT1(i)-GaD1(i)-ACGLP(i);
    
    Gbg1(i)=GaT1(i)+GaD1(i)+ACGLP(i);
    
    LR1(i+1)=LR1(i) + dt*(R11*ReGL(i)*(GLP1(i)/(KGLP1+GLP1(i))) + R12r*LRG1(i)-R12*G1(i)*LR1(i)-R19*LR1(i));
    LRG1(i+1)=LRG1(i)+ dt*(R12*G1(i)*LR1(i)-R12r*LRG1(i)-R13*LRG1(i));
    ReGLd(i+1)=ReGLd(i) + dt*(R19*LR1(i)-R18*ReGLd(i));
    GaT1(i+1)=GaT1(i) + dt*(R13*LRG1(i)-R14*GaT1(i)*ACpf(i)-R15*GaT1(i));
    ACGLP(i+1)= ACGLP(i)+dt*(R14*GaT1(i)*ACpf(i)-R16*ACGLP(i));
    GaD1(i+1)=GaD1(i) + dt*(R15*GaT1(i)+R16*ACGLP(i)-R17*GaD1(i)*Gbg1(i));
    
    %......for GIP....%
    
      
    ReGl(i)=ReGlt-ReGld(i)-LR2(i)-LRG2(i);
    G2(i)=G2t-LRG2(i)-GaT2(i)-GaD2(i)-ACGIP(i);
    
    Gbg2(i)=GaT2(i)+GaD2(i)+ACGIP(i);
    
    LR2(i+1)=LR2(i) + dt*(R21*ReGl(i)*(GIP(i)/(KGIP+GIP(i))) + R22r*LRG2(i)-R22*G2(i)*LR2(i)-R29*LR2(i));
    LRG2(i+1)=LRG2(i)+ dt*(R22*G2(i)*LR2(i)-R22r*LRG2(i)-R23*LRG2(i));
    ReGld(i+1)=ReGld(i) + dt*(R29*LR2(i)-R28*ReGld(i));
    GaT2(i+1)=GaT2(i) + dt*(R23*LRG2(i)-R24*GaT2(i)*ACpf(i)-R25*GaT2(i));
    ACGIP(i+1)= ACGIP(i)+dt*(R24*GaT2(i)*ACpf(i)-R26*ACGIP(i));
    GaD2(i+1)=GaD2(i) + dt*(R25*GaT2(i)+R26*ACGIP(i)-R27*GaD2(i)*Gbg2(i));
    
    
    %.....for catecholamines...%
    
    ReAR(i)=ReARt-ReARd(i)-LR3(i)-LRG3(i);
    G3(i)=G3t-LRG3(i)-GaT3(i)-GaD3(i)-ACAR(i);
    
    Gbg3(i)=GaT3(i)+GaD3(i)+ACAR(i);
    
    LR3(i+1)=LR3(i) + dt*(R31*ReAR(i)*(AR3(i)/(KAR3+AR3(i))) + R32r*LRG3(i)-R32*G3(i)*LR3(i)-R39*LR3(i));
    LRG3(i+1)=LRG3(i)+ dt*(R32*G3(i)*LR3(i)-R32r*LRG3(i)-R33*LRG3(i));
    ReARd(i+1)=ReARd(i) + dt*(R39*LR3(i)-R38*ReARd(i));
    GaT3(i+1)=GaT3(i) + dt*(R33*LRG3(i)-R34*GaT3(i)*ACpf(i)-R35*GaT3(i));
    ACAR(i+1)= ACAR(i)+dt*(R34*GaT3(i)*ACpf(i)-R36*ACAR(i));
    GaD3(i+1)=GaD3(i) + dt*(R35*GaT3(i)+R36*ACAR(i)-R37*GaD3(i)*Gbg3(i));
    %....calmodulin...%
  
  Ca2CaM(i)=(k2f/k2b)*Cac(i)*CaCaM(i);
  Ca3CaM(i)=(k3f/k3b)*Cac(i)*Ca2CaM(i);
  Ca3CaM(i)=(k4f/k4b)*Cac(i)*Ca3CaM(i);
  
  CaM(i)=CaM0-CaCaM(i)-Ca2CaM(i)-Ca3CaM(i)-Ca4CaM(i);
  CaMa(i)=Ca3CaM(i)+Ca4CaM(i);
  
  CaCaM(i+1)=CaCaM(i) + dt*(k1f*Cac(i)*CaM(i)-k1b*CaCaM(i));

%...for cAMP...%

ACpf(i)=ACpt-ACGLP(i)-ACGIP(i)-ACAR(i);
ACpa(i)=ACGLP(i)+ACGIP(i);
ACpar(i)=ACpa(i)/ACpt;

facca(i)=((CaMa(i)/CaMa(i)+Kpcam))*(KNCa/(KNCa+Cac(i)));
VACp(i)=(Vmcam)*(facca(i))*ACpar(i);
VACc(i)=(((VmACc)*ATP(i))/(ATP(i)+KmAACS))*((Cac(i))/(Cac(i)+KmCACS))+kACS;
Vpde(i)=kipde(i)*(Vgpde(i)+Vcpde(i)*(CaMa(i)/(CaMa(i)+Kdpe)))*(cAMP(i)/(cAMP(i)+Kpde));

VAC(i)=VACp(i) + VACc(i);
cAMP(i+1)=cAMP(i) + dt*(VAC(i)-Vpde(i)-kda*cAMP(i));

    %...for PKA...%
    PKAi(i)=1-PKAa(i);
    fpca(i)=(cAMP(i)^1.4)/((Kpcm^1.4)+(cAMP(i)^1.4));
    PKAa(i+1)=PKAa(i) + (dt*(kak*(fpca(i)*PKAi(i)-PKAa(i))/tpka));

%...for Epac...%

EPi(i)=1-EPa(i);
fEcA(i)=(cAMP(i)^1.4)/((Kmep^1.4)+(cAMP(i)^1.4));
EPa(i+1)=EPa(i) + (dt*(kak*(fEcA(i)*EPi(i)-EPa(i))/tep));

    
    % IP3,PIP2 and PLC pathways
    PKCi(i) = PKCt - PKCa(i);                                                                          %% eq 85
    PKCa(i+1) =PKCa(i)+ dt*(kPKC*DAG(i)*PKCi(i) - kPKCr*PKCa(i));                                      %% eq 84
    fPI(i) = (Cac(i)^2/((Cac(i)^2)+(KCaPI^2)))*(PKCa(i)^2/((PKCa(i)^2)+(KP4PK^2)))+kcpi;               %% eq 74
    P4P(i+1) = P4P(i)+ dt*(kPI*fPI(i)*PI+kP4Pr*PIP2(i) - kPIr*P4P(i) - kP4P*P4P(i));                   %% eq 73
    fPIP2(i) = PIP2(i)/(KPIP2+PIP2(i));                                                                %% eq 76
    PLCpa(i) = PLCM3(i) + PLCR40(i);                                                                   %% eq 79
    PLCpf(i) = PLCpt - PLCpa(i);                                                                       %% eq 80
    VPLP(i) = (VmPLP(i))*(PLCpa(i)/PLCpt)*(Cac(i)^2/((KCaPL^2)+Cac(i)^2));                             %% eq 78
    VPLC(i) = VmPLC*(Cac(i)^2)/(KCCaPL^2 + Cac(i)^2);                                                  %% eq 81
    VPL(i) = VPLP(i) + VPLC(i) + kpPL;                                                                 %% eq 77
    PIP2(i+1) = PIP2(i)+ dt*(kP4P*P4P(i) - fPIP2(i)*VPL(i) - kP4Pr*PIP2(i));                           %% eq 75
    IP3(i+1) =IP3(i)+ dt*(VPL(i)*fNc*fPIP2(i) - kdIP3*IP3(i));                                         %% eq 82
    DAG(i+1) =DAG(i)+ dt*(VPL(i)*fPIP2(i) - kdDAG*DAG(i));   
%------------------------------------------------------------ M3R ----------------------------------------------------------%
    G6(i) = G6t - RG6(i) - LRG6(i) - GaT6(i) - GaD6(i) + PLCM3(i);                                  %% eq 49
    Gbg6(i) = GaT6(i) + GaD6(i) + PLCM3(i);                                                         %% '' 
    RG6(i+1) = RG6(i)+ dt*(((R62*ReM3t*G6(i)) -(R61*RG6(i)*AM3(i)))/(KAM3+AM3(i)) - R62r*RG6(i));   %% ''
    LRG6(i+1) = LRG6(i)+ dt*(R61*RG6(i)*AM3(i)/(KAM3+AM3(i)) - R63*LRG6(i) - R69*LRG6(i));          %% ''
    ReM3d(i+1) = ReM3d(i) + dt*(R69*LRG6(i) - R68*ReM3d(i));                                        %% ''
    GaT6(i+1) = GaT6(i)+ dt*(R63*LRG6(i) - R64*GaT6(i)*(PLCpf(i)) - R65*GaT6(i));                      %% ''
    PLCM3(i+1) = PLCM3(i)+ dt*(R64*GaT6(i)*PLCpf(i) - R66*PLCM3(i));                                   %% ''
    GaD6(i+1) = GaD6(i)+ dt*(R65*GaT6(i) +R66*PLCM3(i) - R67*GaD6(i)*Gbg6(i));                      %% ''
    ReM3 = ReM3t - ReM3d(i) - RG6(i) - LRG6(i); 

%----------------------------------------------------------FFAR1/GPR40-----------------------------------------------------%
                                                                                 
    ReR40(i) = ReR40t - ReR40d(i) - LRG7(i) - LR7(i);                                              %% eq 50
    G7(i) = G7t - LRG7(i) - GaT7(i) - GaD7(i) - PLCR40(i);                                         %% ''  
    Gbg7(i+1) = GaT7(i) + GaD7(i) + PLCR40(i);                                                     %% ''
    LR7(i+1) = LR7(i) + dt*(R71*ReR40(i)*AR7(i)/(KAR7 + AR7(i)) + R72r*LRG7(i) - R72*G7(i)*LR7(i) - R79*LR7(i));  %% ''
    LRG7(i+1) = LRG7(i)+ dt*(R72*G7(i)*LR7(i)*R72r*LRG7(i) - R73*LRG7(i));                         %% ''
    ReR40d(i+1) = ReR40d(i)+ dt*(R79*LR7(i) - R78*ReR40d(i));                                      %% ''
    GaT7(i+1) = GaT7(i)+ dt*(R73*LRG7(i) - R74*GaT7(i)*PLCpf(i) - R75*GaT7(i));                       %% ''
    PLCR40(i+1) = PLCR40(i)+ dt*(R74*GaT7(i)*PLCpf(i) - R76*PLCR40(i));                               %% ''
    GaD7(i+1) = GaD7(i)+ dt*(R75*GaT7(i) + R76*PLCR40(i) - R77*GaD7(i)*Gbg7(i));                   %% ''
end

%~PLOTS~%
% Section of plots
% figure(1);
% plot(t,ATD0)
% xlabel('time s')
% ylabel('Units')
% title('ATD_0')
% 
% figure(2);
% plot(t,ATD)
% xlabel('time s')
% ylabel('Units')
% title('ATD')
% 
% figure(3);
% plot(t,ADPf)
% xlabel('time s')
% ylabel('Concentration \muM')
% title('free ADP (ADPf)')
% 
% figure(4);
% plot(t,ATP)
% xlabel('time s')
% ylabel('Concentration \muM')
% title('ATP')
% 
% figure(5);
% plot(t,fKPI)
% xlabel('time s')
% ylabel('Units')
% title('f_{KPI}')
% 
% figure(6);
% plot(t,fKPKa)
% xlabel('time s')
% ylabel('Concentration^{-1} \muM^{-1}')        %Check this part 
% title('f_{KPKa}')
% 
% figure(7);
% plot(t,MgADPf)
% xlabel('time s')
% ylabel('Concentration \muM')
% title('Mg Bound ADP MgADP_f')
% 
% figure(8);
% plot(t,OKATP)
% xlabel('time s')
% ylabel('Probability')
% title('O_{KATP} fraction of open channels')
% 
% figure(9);
% plot(t,IKATP)
% xlabel('time s')
% ylabel('Current fA')
% title('ATP Dependant  pottassium channel current I_{KATP}')
% 
% figure(10);
% plot(t,dKri)
% xlabel('time s')
% ylabel('Probability')
% title('d_{Kri} channel activation probability')
% 
% figure(11);
% plot(t,IKr)
% xlabel('time s')
% ylabel('Current fA')
% title('Voltage Dependant  pottassium channel current  I_{Kr}')
% 
% figure(12);
% plot(t,dCai)
% xlabel('time s')
% ylabel('Probability')
% title('d_{Cai} channel activation probability')
% 
% figure(13);
% plot(t,fCai)
% xlabel('time s')
% ylabel('Probability')
% title('f_{Cai} channel inactivation probability')
% 
% figure(14);
% plot(t,ICa)
% xlabel('time s')
% ylabel('Current fA')
% title('Voltage Dependant  pottassium channel current I_{Ca}')
% 
% figure(15);
% plot(t,fSOC)
% xlabel('time s')
% ylabel('Probability')
% title('f_{SOC} channel inactivation probability')
% 
% figure(16);
% plot(t,ISOC)
% xlabel('time s')
% ylabel('Current fA')
% title('Store operated channel current I_{SOC}')
% 
% figure(17);
% plot(t,IPCa)
% xlabel('time s')
% ylabel('Current fA')
% title('Calcium pump current I_{PCa}')
% 
% figure(18);
% plot(t,gNab)
% xlabel('time s')
% ylabel(' Conductance pS')
% title('Sodium Conductance g_{Nab}')
% 
% figure(19);
% plot(t,INab)
% xlabel('time s')
% ylabel('Current fA')
% title('Sodium background current I_{Nab}')
% 
% figure(20);
% plot(t,Vp)
% xlabel('time s')
% ylabel('Voltage mV')
% title('Membrane potential V_p')
% 
% figure(21);
% plot(t,Jser)
% xlabel('time s')
% ylabel('flux \mumols^-1 ')
% title('Flux of Calcium from cytosol to ER J_{SER}')
% 
% figure(22);
% plot(t,fIPKA)
% xlabel('time s')
% ylabel('fraction')
% title(' f_{IPKA}')
% 
% figure(23);
% plot(t,PRIP3)
% xlabel('time s')
% ylabel('Probability')
% title('Channel open probability P_{RIP3}')
% 
% figure(24);
% plot(t,Jrel)
% xlabel('time s')
% ylabel('flux \mumols^-1 ')
% title('Flux of Calcium from ER to cytosol J_{rel}')
% 
% figure(25);
% plot(t,Cac)
% xlabel('time s')
% ylabel('Concentration \muM')
% title('Cytosolic calcium concentration')
% 
% 
% figure(26);
% plot(t,CaER)
% xlabel('time s')
% ylabel('Concentration \muM')
% title('ER calcium concentration')

%....Incretin catecholamines graph...%

% % figure(27);
% % plot(t,LR1)
% % xlabel('time s')
% % ylabel('LR1')
% % title('LR1')
% % 
% % figure(28);
% % plot(t,LR2)
% % xlabel('time s')
% % ylabel('LR2')
% % title('LR2')
% % 
% % figure(29);
% % plot(t,LR3)
% % xlabel('time s')
% % ylabel('LR3')
% % title('LR3')
% % 
% % figure(30);
% % plot(t,LRG1)
% % xlabel('time s')
% % ylabel('LRG1')
% % title('LRG1')
% % 
% % figure(31);
% % plot(t,LRG2)
% % xlabel('time s')
% % ylabel('LRG2')
% % title('LRG2')
% % 
% % figure(32);
% % plot(t,LRG3)
% % xlabel('time s')
% % ylabel('LRG3')
% % title('LRG3')
% % 
% % figure(33);
% % plot(t,G1)
% % xlabel('time s')
% % ylabel('G1')
% % title('G1')
% % 
% % figure(34);
% % plot(t,G2)
% % xlabel('time s')
% % ylabel('G2')
% % title('G2')
% % 
% % figure(35);
% % plot(t,G3)
% % xlabel('time s')
% % ylabel('G3')
% % title('G3')
% % 
% % figure(36);
% % plot(t,GaD1)
% % xlabel('time s')
% % ylabel('GaD1')
% % title('GaD1')
% % 
% % figure(37);
% % plot(t,GaD2)
% % xlabel('time s')
% % ylabel('GaD2')
% % title('GaD2')
% % 
% % figure(38);
% % plot(t,GaD3)
% % xlabel('time s')
% % ylabel('GaD3')
% % title('GaD3')
% % 
% % figure(39);
% % plot(t,ACGLP)
% % xlabel('time s')
% % ylabel('ACGLP')
% % title('ACGLP')
% % 
% % figure(40);
% % plot(t,ACGIP)
% % xlabel('time s')
% % ylabel('ACGIP')
% % title('ACGIP')
% % 
% % figure(41);
% % plot(t,ACAR)
% % xlabel('time s')
% % ylabel('ACAR')
% % title('ACAR')
% %
% %figure(42);
% %plot(t,PKAa)
% %xlabel('time s')
% %ylabel('PKAa')
% %title('PKAa')
% %
% %figure(43);
% %plot(t,CaMa)
% %xlabel('time s')
% %ylabel('CaMa')
% %title('CaMa')













