% DOM_10.m

% This MATLAB m.file contains the program for the Dual Oscillator Model,
% with equations included for oxidative phosphorylation.
% This was used to make Figs. 1 and 2 in Merrins, Fendler, Zhang, Sherman, 
% Bertram. Satin, "Metabolic Oscillations in Pancreatic
% Islets Depend on the Intracellular Calcium Level but Not Calcium
% Oscillations", Biophysical Journal, vol. 99, pp. 76-84, 2010.

% Variables:
%    v -- voltage
%    n -- activation variable for a delayed rectifier
%    c -- free cytosolic calcium concentration
%    cer -- concentration of free calcium in the endoplasmic reticulum
%    Ins -- insulin concentration
%    g6p -- glucose 6-phosphate concentration
%    fbp -- fructose bis-phosphate concentration
%    adpm -- mitochondrial ADP concentration
%    nadhm -- mitochondrial NADH concentration
%    cam -- mitochondrial free calcium concentration
%    psim -- mitochondrial inner membrane potential
%    adp -- cytosolic ADP concentration


function [t,x] = call_DOM();
   tmax = 600000;  % stop before Dz is added
%   tmax = 1800000;  % continue past Dz application
   tspan= [0 tmax];
   v0=-38; n0=0; c0=0.34; cer0=578; cam0=1.7; g6p0=189; fbp0=2.9; adpm0=11.5;
   nadhm0=0.6; psim0=161; adp0=1427; Ins0=0.8;
   x0 = [v0 n0 c0 cer0 cam0 g6p0 fbp0 adpm0 nadhm0 psim0 adp0 Ins0 ];
%   [t,x] = ode113(@DOM,tspan,x0);
   [t,x] = ode15s(@DOM,tspan,x0);

%  change defaults for plotting
   set(0,'DefaultAxesFontSize',18)

   tsec=t/1000;
   tmin=tsec/60;
   tmax_min=tmax/60000;

%  set up variable vectors
   v=x(:,1); n=x(:,2); c=x(:,3); cer=x(:,4); cam=x(:,5); g6p=x(:,6); fbp=x(:,7);
   adpm=x(:,8); nadhm=x(:,9); psim=x(:,10); adp=x(:,11); Ins=x(:,12);

%  Time course figure
   figure;

%  Top left panel
   subplot(2,2,1);
   h=plot(tmin,v,'black');
   set(h,'LineWidth',2);
   axis([0 tmax_min -80 -10]);
   box off;
   ylabel('V (mV)');

%  Top right panel
   subplot(2,2,2);
   h=plot(tmin,c,'black');
   set(h,'LineWidth',2);
   axis([0 tmax_min 0.1 0.4]);
   box off;
   ylabel('Ca (\muM)');

%  Bottom left panel
   subplot(2,2,3);
   h=plot(tmin,nadhm,'black');
   set(h,'LineWidth',2);
   axis([0 tmax_min 0 4]);
   box off;
   xlabel('t (m)');
   ylabel('NADH_m (mM)');

%  Bottom right panel
   subplot(2,2,4);
   h=plot(tmin,fbp,'black');
   set(h,'LineWidth',2);
   axis([0 tmax_min 0 10]);
   box off;
   xlabel('t (m)');
   ylabel('FBP (\muM)');

%  Print data set into text file
%   set up a matrix of output. Must transpose the column vectors.
    out = [tmin.'; v.'; n.'; c.'; cer.'; cam.'; g6p.'; fbp.'; adpm.';
          nadhm.'; psim.'; adp.'; Ins.'];
    fid = fopen('data.dat','w');
    fprintf(fid,'%9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n',out);
    fclose(fid);

   function dxdt = DOM(t,x);

%       set up variable vectors
        v=x(1); n=x(2); c=x(3); cer=x(4); cam=x(5); g6p=x(6); fbp=x(7);
        adpm=x(8); nadhm=x(9); psim=x(10); adp=x(11); Ins=x(12);

%       Parameters most likely to vary:
        Jgk=0.47;   % MOs persist in Dz
%        Jgk=0.53;   % MOs stop in Dz
        gkatpbar=11000;
        gkca=250;

%       Diazoxide and KCl 
        Dzamp=1; Dzstart=600000; KClstart=3600000;
        vk0=-75; vk1=-62.5; 
        Dz = Dzamp*heaviside(t - Dzstart);
        vk = vk0 + (vk1 - vk0)*heaviside(t - KClstart);

%       Membrane capacitance
        cm=5300;

%       Ikca
        kd=0.5;
        Ikca = gkca/(1+(kd/c)^2)*(v-vk);       

%       ICa
        vca=25; gca=1000;
        vm=-20; sm=12;
        minf = 1/(1+exp((vm-v)/sm));
        Ica = gca*minf*(v-vca);

%       Ik
        gk=2000;
        vn=-16; sn=5; taun=20;
        ninf = 1/(1+exp((vn-v)/sn));
        Ik = gk*n*(v-vk);

%       IKatp
        kdd=17; ktd=26; ktt=1; 

%       Cytosolic nucleotide concentrations
        khyd=0.00005; JhydSS=0.00005;
        amp=500; atot=2500;
        atp = atot-adp;
        Jhyd=(khyd*c+JhydSS)*atp;
        mgadp = 0.165*adp;
        adp3m = 0.135*adp;
        atp4m = 0.05*atp;
        RATc = atp/adp;
      
%       KATP channel open probability (from Magnus and Keizer, 1998)
        topo = 0.08*(1+2*mgadp/kdd) + 0.89*(mgadp/kdd)^2;
        bottomo = (1+mgadp/kdd)^2 * (1+adp3m/ktd+atp4m/ktt);
        katpo = (topo/bottomo);
        Ikatp = (1-Dz)*gkatpbar*katpo*(v-vk)+Dz*gkatpbar*(v-vk);

%       calcium handling, sigmav=cyt volume/ER volume
        alpha=4.50e-6; fcyt=0.01; fer=0.01; kpmca=0.1;
        per=0.0002; kserca=0.4; lambdaer=1; epser=1;
        sigmav=31; cbas=0.05;

%       calcium fluxes
        Jmem = -(alpha*Ica + kpmca*(c-cbas));
        Jserca = kserca*c;
        Jleak = per*(cer - c);
        Jer = epser*(Jleak - Jserca)/lambdaer;

%      glycolytic parameters
%        Rgk--glucokinase rate
%        atot--total adenine nucleotide concentration (micromolar)
%        k1--Kd for AMP binding
%        k2--Kd for FBP binding
%        k3--Kd for F6P binding
%        k4--Kd for ATP binding
%        famp,etc--Kd amplification factors for heterotropic binding
%        Rgpdh--glyceraldehyde phosphate dehydrogenase rate

        phigk=0.3; kgpdh=0.0005;
        k1=30; k2=1; k3=50000; k4=330;
        f13=0.02; f43=20; f23=0.2; f42=20; f41=20; lambda=0.06;
        VmaxPFK=5; kappa=0.001; Jgpdh_bas=0.0005;

%     Glycolytic expressions
        f6p = phigk*g6p;
        Jgpdh = kgpdh*sqrt(fbp);
        Jgk_ms=kappa*Jgk;
 
%        Iterative calculation of PFK
%        alpha=1 -- AMP bound
%        beta=1 -- FBP bound
%        gamma=1 -- F6P bound
%        delta=1 -- ATP bound

%      (alpha,beta,gamma,delta)
%      (0,0,0,0)
       weight1=1;
       topa1=0;
       bottom1=1;
          
%      (0,0,0,1)
       weight2=atp^2/k4;
       topa2=topa1;
       bottom2=bottom1+weight2;

%      (0,0,1,0)
       weight3=f6p^2/k3;
       topa3=topa2+weight3;
       bottom3=bottom2+weight3;

%      (0,0,1,1)
       weight4=(f6p*atp)^2/(f43*k3*k4);
       topa4=topa3+weight4;
       bottom4=bottom3+weight4;

%      (0,1,0,0)
       weight5=fbp/k2;
       topa5=topa4;
       bottom5=bottom4+weight5;

%      (0,1,0,1)
       weight6=(fbp*atp^2)/(k2*k4*f42);
       topa6=topa5;
       bottom6=bottom5+weight6;

%      (0,1,1,0)
       weight7=(fbp*f6p^2)/(k2*k3*f23);
       topa7=topa6+weight7;
       bottom7=bottom6+weight7;

%      (0,1,1,1)
       weight8=(fbp*f6p^2*atp^2)/(k2*k3*k4*f23*f42*f43);
       topa8=topa7+weight8;
       bottom8=bottom7+weight8;

%      (1,0,0,0)
       weight9=amp/k1;
       topa9=topa8;
       bottom9=bottom8+weight9;

%      (1,0,0,1)
       weight10=(amp*atp^2)/(k1*k4*f41);
       topa10=topa9;
       bottom10=bottom9+weight10;

%     (1,0,1,0)
      weight11=(amp*f6p^2)/(k1*k3*f13);
      topa11=topa10+weight11;
      bottom11=bottom10+weight11;

%     (1,0,1,1)
      weight12=(amp*f6p^2*atp^2)/(k1*k3*k4*f13*f41*f43);
      topa12=topa11+weight12;
      bottom12=bottom11+weight12;

%     (1,1,0,0)
      weight13=(amp*fbp)/(k1*k2);
      topa13=topa12;
      bottom13=bottom12+weight13;

%     (1,1,0,1)
      weight14=(amp*fbp*atp^2)/(k1*k2*k4*f42*f41);
      topa14=topa13;
      bottom14=bottom13+weight14;

%     (1,1,1,0) -- the most active state of the enzyme
      weight15=(amp*fbp*f6p^2)/(k1*k2*k3*f23*f13);
      topa15=topa14;
      topb=weight15;
      bottom15=bottom14+weight15;

%     (1,1,1,1)
      weight16=(amp*fbp*f6p^2*atp^2)/(k1*k2*k3*k4*f23*f13*f42*f41*f43);
      topa16=topa15+weight16;
      bottom16=bottom15+weight16;

%     Phosphofructokinase rate
      pfk=(lambda*VmaxPFK*topa16 + VmaxPFK*topb)/bottom16;
      pfk_ms=kappa*pfk;

%     Mitochondrial Ca2+ handling
%        delta transforms  (mito volume) -> (cyto volume)/ms
%        gamma transforms uM to mM        

      fmito=0.01;
      delta=3.9/53.2; gamma=0.001;

%     Uniporter (uM/ms)
      p21=0.04; p22=1.1;
      Juni=(p21*psim-p22)*c^2;

%     Na/Ca exchanger (uM/ms)
      p23=0.01; p24=0.016;
      JNaCa=(p23/c)*cam*exp(p24*psim);

%     mito -> cyto (uM/ms)
      Jmito=JNaCa-Juni;

%     Mito fluxes (MK98 and Cortassa et al. 2003)

%     Mitochondrial adenine and pyridine nucleotide conservation. All in mM.
      Amtot=15; NADmtot=10;
      NADm=NADmtot-nadhm; 
      ATPm=Amtot-adpm;
      RATm=ATPm/adpm;

%     Pyruvate dehydrogenase (PDH) (uM/ms)
      p1=400; p2=1; p3=0.01;
      Jpdh=(p1/(p2+nadhm/NADm))*(cam/(p3+cam))*(Jgpdh+Jgpdh_bas);

%     Mitochondrial inner membrane capacitance 
      Cmito=1.8;

%     H+ leakage through mitochondrial inner membrane (uM/ms)
      p17=0.002; p18=-0.03;
      JHleak=p17*psim+p18;

%     Respiration (uM/ms)
      p4=0.6; p5=0.1; p6=177; p7=5;
      MM1=p4*nadhm/(p5+nadhm);
      JO=MM1/(1+exp((psim-p6)/p7));

%     Proton pumping due to respiration  (uM/ms)
      p8=7; p9=0.1; p10=177; p11=5;
      MM2=p8*nadhm/(p9+nadhm);
      JHres=MM2/(1+exp((psim-p10)/p11));

%     Phosphorylation (uM/ms)
      p13=10; p14=190; p15=8.5; p16=35;
      b2=(p16*p13)/(p13+ATPm);
      JF1F0=b2/(1.0+exp((p14-psim)/p15));

%     Proton flux due to ATPase (uM/ms)
      JHatp=3*JF1F0;

%     ADP/ATP translocator (uM/ms)
      p19=0.35; p20=2;
      FRT=96480/(310.16*8315);
      Jant=p19*(RATm/(RATm+p20))/exp(-0.5*FRT*psim);
       
%     Insulin secretion
      taui=1000; ki=0.15;
      Iinf = (c^2)/(ki^2 + c^2);

%     Differential equations
      vdot = -(Ik + Ica + Ikca + Ikatp)/cm;
      ndot = (ninf-n)/taun;
      cdot = fcyt*(Jmem + Jer + delta*Jmito);
      cerdot = -fer*sigmav*Jer;
      camdot = -fmito*Jmito;
      g6pdot = Jgk_ms - pfk_ms;
      fbpdot = pfk_ms - 0.5*Jgpdh;
      adpmdot = gamma*(Jant-JF1F0);
      nadhmdot = gamma*(Jpdh-JO);
      psimdot = (JHres-JHatp-Jant-JHleak-JNaCa-2*Juni)/Cmito;
      adpdot = -delta*Jant + Jhyd;
      Insdot = (Iinf-Ins)/taui;
      dxdt=[vdot;ndot;cdot;cerdot;camdot;g6pdot;fbpdot;adpmdot;nadhmdot;psimdot;adpdot;Insdot];
    end
end
     
