ORdmD16.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "ode/ionicModels/ORdmD16.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
 
using std::exp;
using std::expm1;
using std::pow;
using std::log;
using std::sqrt;
using ceps::approxEquals;
 
ORdmD16::ORdmD16(
  const ORdmD16::Type&          type,
  Unknown*                      u, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  AbstractIonicModel(u,attrs),
  m_type(type)
{
  
  m_name          = "O'Hara-Rudy modified 2016";
  m_nStateVars    = 40;
  m_stateVarNames = {
    "Na_i []",
    "Na_ss []",
    "K_i []",
    "K_ss []",
    "Ca_i []",
    "Ca_ss []",
    "Ca_nsr []",
    "Ca_jsr []",
    "m [adim]",
    "h_fast [adim]",
    "h_slow [adim]",
    "j [adim]",
    "h_slow_p [adim]",
    "jp []",
    "m_late [adim]",
    "h_late [adim]",
    "h_late_p [adim]",
    "a [adim]",
    "i_fast [adim]",
    "i_slow [adim]",
    "a_p [adim]",
    "i_fast_p [adim]",
    "i_slow_p [adim]",
    "d [adim]",
    "f_fast [adim]",
    "f_slow [adim]",
    "fca_fast [adim]",
    "fca_slow [adim]",
    "jca [adim]",
    "nca [adim]",
    "f_fast_p [adim]",
    "fca_fast_p [adim]",
    "xr_fast [adim]",
    "xr_slow [adim]",
    "xs1 [adim]",
    "xs2 [adim]",
    "xk1 [adim]",
    "Jrelnp []",
    "Jrelp []",
    "CaMKt []",
  };
 
  // ================================================================================
  // Default parameters
  m_constants  = ceps::newArray<CepsReal>(32); 
  m_dconstants = ceps::newArray<CepsReal>(15);
  CepsReal* c  = m_constants ;
  CepsReal* dc = m_dconstants;
  c[_IKr_factor   ] = 1.119;
  c[_IKs_factor   ] = 1.648;
  c[_IK1_factor   ] = 1.414;
  c[_ICaL_factor  ] = 1.018;
  c[_INaL_factor  ] = 2.274;
  c[_nao          ] = 140.0;
  c[_cao          ] = 1.8;
  c[_ko           ] = 5.4;
  c[_R            ] = 8314.0;
  c[_T            ] = 310.0;
  c[_F            ] = 96485.0;
  c[_L            ] = 0.01;
  c[_rad          ] = 0.0011;
  c[_vmyo_fraction] = 0.68;
  c[_vnsr_fraction] = 0.0552;
  c[_vjsr_fraction] = 0.0048;
  c[_vss_fraction ] = 0.02;
  c[_drug         ] = 1448;
  c[_ICNa         ] = 5800;
  c[_hNa          ] = 1;
  c[_ICNaL        ] = 18870;
  c[_hNaL         ] = 0.6;
  c[_ICto         ] = 9266;
  c[_hto          ] = 0.7;
  c[_ICCal        ] = 26349.5;
  c[_hCal         ] = 1.185;
  c[_ICKr         ] = 1500;
  c[_hKr          ] = 0.88;
  c[_ICKs         ] = 20000;
  c[_hKs          ] = 1;
  c[_ICK1         ] = std::nan("");
  c[_hK1          ] = std::nan("");
  c[_PKNa         ] = 0.018333;      
 
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params);
 
  dc[_vcell   ] = 1000*3.14*c[_rad]*c[_rad]*c[_L];   // 3.14... seriously...
  dc[_ageo    ] = (2.*3.14*c[_rad])*(c[_rad]+c[_L]); // But we stick to the original
  dc[_acap    ] = 2*dc[_ageo];
  dc[_vmyo    ] = c[_vmyo_fraction]*dc[_vcell];
  dc[_vnsr    ] = c[_vnsr_fraction]*dc[_vcell];
  dc[_vjsr    ] = c[_vjsr_fraction]*dc[_vcell];
  dc[_vss     ] = c[ _vss_fraction]*dc[_vcell];
  dc[_FRT     ] = c[_F]/(c[_R]*c[_T]);
  dc[_gNaFrac ] = (std::isnan(c[_ICNa ]) or std::isnan(c[_hNa ])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICNa ],c[_hNa ]));
  dc[_gNaLFrac] = (std::isnan(c[_ICNaL]) or std::isnan(c[_hNaL])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICNaL],c[_hNaL]));
  dc[_gtoFrac ] = (std::isnan(c[_ICto ]) or std::isnan(c[_hto ])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICto ],c[_hto ]));
  dc[_gCalFrac] = (std::isnan(c[_ICCal]) or std::isnan(c[_hCal])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICCal],c[_hCal]));
  dc[_gKrFrac ] = (std::isnan(c[_ICKr ]) or std::isnan(c[_hKr ])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICKr ],c[_hKr ]));
  dc[_gKsFrac ] = (std::isnan(c[_ICKs ]) or std::isnan(c[_hKs ])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICKs ],c[_hKs ]));
  dc[_gK1Frac ] = (std::isnan(c[_ICK1 ]) or std::isnan(c[_hK1 ])) ? 1. : 1/(1.+pow(c[_drug]/c[_ICK1 ],c[_hK1 ]));
 
 
  // Other constants are deduced from these ones, almost
  m_paperCm   = 1.;
  m_paperStim = -80;
 
}
 
void 
ORdmD16::getInitialCondition(CepsReal* v, CepsReal* y) const
{
  *v         = -87.84      ; // mV
  y[_nai   ] = 7.23        ; // mMol
  y[_nass  ] = 7.23        ; // mMol
  y[_ki    ] = 143.79      ; // mMol
  y[_kss   ] = 143.79      ; // mMol
  y[_cai   ] = 8.54e-05    ; // mMol
  y[_cass  ] = 8.43e-05    ; // mMol
  y[_cansr ] = 1.61        ; // mMol
  y[_cajsr ] = 1.56        ; // mMol
  y[_m     ] = 0.0074621   ; // 
  y[_hf    ] = 0.692591    ; //
  y[_hs    ] = 0.692574    ; //
  y[_j     ] = 0.692477    ; //
  y[_hsp   ] = 0.448501    ; //
  y[_jp    ] = 0.692413    ; //
  y[_mL    ] = 0.000194015 ; //
  y[_hL    ] = 0.496116    ; //
  y[_hLp   ] = 0.265885    ; //
  y[_a     ] = 0.00101185  ; //
  y[_iF    ] = 0.999542    ; //
  y[_iS    ] = 0.589579    ; //
  y[_ap    ] = 0.000515567 ; //
  y[_iFp   ] = 0.999542    ; //
  y[_iSp   ] = 0.641861    ; //
  y[_d     ] = 2.43015e-09 ; //
  y[_ff    ] = 1           ; //
  y[_fs    ] = 0.916071    ; //
  y[_fcaf  ] = 1           ; //
  y[_fcas  ] = 0.99982     ; //
  y[_jca   ] = 0.999977    ; //
  y[_nca   ] = 0.00267171  ; //
  y[_ffp   ] = 1           ; //
  y[_fcafp ] = 1           ; //
  y[_xrf   ] = 8.26608e-06 ; //
  y[_xrs   ] = 0.453268    ; //
  y[_xs1   ] = 0.270492    ; //
  y[_xs2   ] = 0.0001963   ; //
  y[_xk1   ] = 0.996801    ; //
  y[_Jrelnp] = 2.53943e-05 ; //
  y[_Jrelp ] = 3.17262e-07 ; //
  y[_CaMKt ] = 0.0124065   ; //
}
 
void
ORdmD16::computeRates(
  CepsReal         t,
  CepsReal*        y, 
  CepsReal*        vm, 
  CepsReal*        dtyLin, 
  CepsReal*        dtyNLin, 
  CepsReal*        dtv,
  DegreeOfFreedom* dof
) const
{ 
 
  CepsReal* c   = m_constants;
  CepsReal* d   = m_dconstants;
  CepsReal  v   = *vm;
  CepsReal  cm  = getCmInternal(t,dof);
 
  // CaMK 
 
  CepsReal KmCaMK   = 0.15;
  CepsReal aCaMK    = 0.05;
  CepsReal bCaMK    = 6.8e-4;
  CepsReal CaMKo    = 0.05;
  CepsReal KmCaM    = 1.5e-3;
  CepsReal CaMKb    = CaMKo*(1.-y[_CaMKt])/(1.+KmCaM/y[_cass]);
  CepsReal CaMKa    = CaMKb + y[_CaMKt];
  CepsReal fraction = 1./(1.+KmCaMK/CaMKa);
 
  dtyLin [_CaMKt] = -bCaMK;
  dtyNLin[_CaMKt] = aCaMK*CaMKb*(CaMKb+y[_CaMKt]);
 
  // Fast Na+ current, late Na+ current
  CepsReal INa,INaL;
  computeINa(&INa,&INaL,v,y,dtyLin,dtyNLin,fraction);
 
  // Transient outward K+ current
  CepsReal Ito;
  computeIto(&Ito,v,y,dtyLin,dtyNLin,fraction);
  
  // L-type Ca2+ channel 
  CepsReal ICaL, ICaNa, ICaK;
  computeICa(&ICaL,&ICaNa,&ICaK,v,y,dtyLin,dtyNLin,fraction);
 
  // Rapid delayed rectifier K+ current
  CepsReal IKr;
  computeIKr(&IKr,v,y,dtyLin,dtyNLin);
 
  // Slow delayed rectifier K+ current
  CepsReal IKs;
  computeIKs(&IKs,v,y,dtyLin,dtyNLin);
 
  // Inward Rectifier K+ Current 
  CepsReal IK1;
  computeIK1(&IK1,v,y,dtyLin,dtyNLin);
 
  // Na+/Ca2+ exchanger
  CepsReal INaCa_i, INaCa_ss;
  computeINaCa(&INaCa_i,&INaCa_ss,v,y);
 
  // Sodium/Potassium ATPase Current
  CepsReal INaK;
  computeINaK(&INaK,v,y);
 
  // Background currents
  CepsReal IKb,INab,ICab;
  computeBgCurrents(&IKb,&INab,&ICab,v,y);
 
  CepsReal IpCa;
  computeIPCa(&IpCa,y);
 
  // Potential
  CepsReal IStim = getStimulation(dof,t);
  CepsReal Iion  = INa+INaL+Ito+ICaL+ICaNa+ICaK+IKr+IKs+IK1+INaCa_i+INaCa_ss+INaK+INab+IKb+IpCa+ICab;
  #ifdef DEBUG
    ceps::checkNanOrInf(Iion,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
  *dtv = -(Iion+IStim)/cm;
 
  // Diffusion fluxes
  CepsReal JdiffNa = (y[_nass]-y[_nai])/2.;
  CepsReal JdiffK  = (y[_kss ]-y[_ki ])/2.;
  CepsReal JdiffCa = (y[_cass]-y[_cai])/0.2;
  //ryanodione receptor calcium induced calcium release from the jsr
  CepsReal Jrel;
  computeJrel(&Jrel,ICaL,y,dtyLin,dtyNLin,fraction);
 
  //  Calcium Uptake via SERCA Pump
  CepsReal Jup;
  computeJup(&Jup,y,fraction);
  
  CepsReal Jtr = (y[_cansr]-y[_cajsr])/100.;
 
  // update concentration rates
  {  
    CepsReal cmdnmax = 0.05;  
    if (m_type==ORdmD16::Type::Epi)
      cmdnmax *= 1.3;
 
    dtyLin [_nai] = 0.;
    dtyNLin[_nai] = -(INa+INaL+3.0*INaCa_i+3.0*INaK+INab)*d[_acap]/(c[_F]*d[_vmyo])+JdiffNa*d[_vss]/d[_vmyo];
 
    dtyLin [_nass] = 0.;
    dtyNLin[_nass] = -(ICaNa+3.0*INaCa_ss)*d[_acap]/(c[_F]*d[_vss])-JdiffNa;
 
    dtyLin [_ki] = 0.;
    dtyNLin[_ki] = -(Ito+IKr+IKs+IK1+IKb+IStim-2.0*INaK)*d[_acap]/(c[_F]*d[_vmyo])+JdiffK*d[_vss]/d[_vmyo];
 
    dtyLin [_kss] = 0.;
    dtyNLin[_kss] = -ICaK*d[_acap]/(c[_F]*d[_vss])-JdiffK;
 
    CepsReal kmcmdn  = 0.00238;
    CepsReal trpnmax = 0.07;
    CepsReal kmtrpn  = 0.0005;
    CepsReal csqnmax = 10.0;
    CepsReal kmcsqn  = 0.8;
    CepsReal Bcai    = 1.0/(1.0+cmdnmax*kmcmdn/std::pow(kmcmdn+y[_cai],2.0)
                               +trpnmax*kmtrpn/std::pow(kmtrpn+y[_cai],2.0));
    dtyLin [_cai]    = 0.;
    dtyNLin[_cai]    = Bcai*(-(IpCa+ICab-2.0*INaCa_i)*d[_acap]/(2.0*c[_F]*d[_vmyo])
                             -Jup*d[_vnsr]/d[_vmyo]+JdiffCa*d[_vss]/d[_vmyo]);
 
    CepsReal BSRmax  = 0.047;
    CepsReal KmBSR   = 0.00087;
    CepsReal BSLmax  = 1.124;
    CepsReal KmBSL   = 0.0087;
    CepsReal Bcass   = 1.0/(1.0+BSRmax*KmBSR/std::pow(KmBSR+y[_cass],2.0)
                               +BSLmax*KmBSL/std::pow(KmBSL+y[_cass],2.0));
    dtyLin [_cass]   = 0.;
    dtyNLin[_cass]   = Bcass*(-(ICaL-2.0*INaCa_ss)*d[_acap]/(2.0*c[_F]*d[_vss])+Jrel*d[_vjsr]/d[_vss]-JdiffCa);
 
    dtyLin [_cansr] = 0.;
    dtyNLin[_cansr] = Jup - Jtr*d[_vjsr]/d[_vnsr];
 
    CepsReal Bcajsr = 1/(1+csqnmax*kmcsqn/std::pow(kmcsqn+y[_cajsr],2.0));
    dtyLin [_cajsr] = 0.;
    dtyNLin[_cajsr] = Bcajsr*(Jtr-Jrel);
 
  }
 
}
 
void
ORdmD16::setupWithParameters(InputParameters* params)
{
 
  CepsString key = "ORDMD16";
 
  m_constants[_IKr_factor   ] = params->getReal(key+"IKr_factor   ",m_constants[_IKr_factor   ]);
  m_constants[_IKs_factor   ] = params->getReal(key+"IKs_factor   ",m_constants[_IKs_factor   ]);
  m_constants[_IK1_factor   ] = params->getReal(key+"IK1_factor   ",m_constants[_IK1_factor   ]);
  m_constants[_ICaL_factor  ] = params->getReal(key+"ICaL_factor  ",m_constants[_ICaL_factor  ]);
  m_constants[_INaL_factor  ] = params->getReal(key+"INaL_factor  ",m_constants[_INaL_factor  ]);
  m_constants[_nao          ] = params->getReal(key+"nao          ",m_constants[_nao          ]);
  m_constants[_cao          ] = params->getReal(key+"cao          ",m_constants[_cao          ]);
  m_constants[_ko           ] = params->getReal(key+"ko           ",m_constants[_ko           ]);
  m_constants[_R            ] = params->getReal(key+"R            ",m_constants[_R            ]);
  m_constants[_T            ] = params->getReal(key+"T            ",m_constants[_T            ]);
  m_constants[_F            ] = params->getReal(key+"F            ",m_constants[_F            ]);
  m_constants[_L            ] = params->getReal(key+"L            ",m_constants[_L            ]);
  m_constants[_rad          ] = params->getReal(key+"rad          ",m_constants[_rad          ]);
  m_constants[_vmyo_fraction] = params->getReal(key+"vmyo_fraction",m_constants[_vmyo_fraction]);
  m_constants[_vnsr_fraction] = params->getReal(key+"vnsr_fraction",m_constants[_vnsr_fraction]);
  m_constants[_vjsr_fraction] = params->getReal(key+"vjsr_fraction",m_constants[_vjsr_fraction]);
  m_constants[_vss_fraction ] = params->getReal(key+"vss_fraction ",m_constants[_vss_fraction ]);
  m_constants[_drug         ] = params->getReal(key+"drug         ",m_constants[_drug         ]);
  m_constants[_ICNa         ] = params->getReal(key+"ICNa         ",m_constants[_ICNa         ]);
  m_constants[_hNa          ] = params->getReal(key+"hNa          ",m_constants[_hNa          ]);
  m_constants[_ICNaL        ] = params->getReal(key+"ICNaL        ",m_constants[_ICNaL        ]);
  m_constants[_hNaL         ] = params->getReal(key+"hNaL         ",m_constants[_hNaL         ]);
  m_constants[_ICto         ] = params->getReal(key+"ICto         ",m_constants[_ICto         ]);
  m_constants[_hto          ] = params->getReal(key+"hto          ",m_constants[_hto          ]);
  m_constants[_ICCal        ] = params->getReal(key+"ICCal        ",m_constants[_ICCal        ]);
  m_constants[_hCal         ] = params->getReal(key+"hCal         ",m_constants[_hCal         ]);
  m_constants[_ICKr         ] = params->getReal(key+"ICKr         ",m_constants[_ICKr         ]);
  m_constants[_hKr          ] = params->getReal(key+"hKr          ",m_constants[_hKr          ]);
  m_constants[_ICKs         ] = params->getReal(key+"ICKs         ",m_constants[_ICKs         ]);
  m_constants[_hKs          ] = params->getReal(key+"hKs          ",m_constants[_hKs          ]);
  m_constants[_ICK1         ] = params->getReal(key+"ICK1         ",m_constants[_ICK1         ]);
  m_constants[_hK1          ] = params->getReal(key+"hK1          ",m_constants[_hK1          ]);
  // Other constants are deduced from these ones
 
}
 
void
ORdmD16::computeINa(
  CepsReal* INa, 
  CepsReal* INaL, 
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin,
  CepsReal  fraction
) const
{
  CepsReal wInfm = sigmoid(v,-(39.57+9.4),-7.5);
  CepsReal taum  = 1./(6.765*exp((v+11.64)/34.77) +8.552*exp(-(v+77.42)/5.955));
  dtyLin [_m] = -1./taum;
  dtyNLin[_m] = wInfm/taum;
 
  CepsReal wInfh = sigmoid(v,-78.5,6.22);
  CepsReal tauhf = 1./(1.432e-5*exp(-(v+1.196)/6.285)+6.149*exp((v+0.5096)/20.27));
  CepsReal tauhs = 1./(0.009794*exp(-(v+17.95)/28.05)+0.3343*exp((v+5.730)/56.66));
  dtyLin [_hf] = -1./tauhf;
  dtyNLin[_hf] = wInfh/tauhf;
  dtyLin [_hs] = -1./tauhs;
  dtyNLin[_hs] = wInfh/tauhs;
 
  CepsReal Ahf   = 0.99;
  CepsReal h     = Ahf*y[_hf] +(1.-Ahf)*y[_hs];
  CepsReal wInfj = wInfh;
  CepsReal tauj = 2.038+1.0/(0.02136*exp(-(v+100.6)/8.281)+0.3052*exp((v+0.9941)/38.45));
  dtyLin [_j] = -1./tauj;
  dtyNLin[_j] = wInfj/tauj;
 
  CepsReal wInfhsp = sigmoid(v,-(78.5+6.2),6.22);
  CepsReal tauhsp  = 3.0*tauhs;
  dtyLin [_hsp] = -1./tauhsp;
  dtyNLin[_hsp] = wInfhsp/tauhsp;
 
  CepsReal hp    = Ahf*y[_hf] + (1.-Ahf)*y[_hsp];
  CepsReal taujp = 1.46*tauj;
  dtyLin [_jp] = -1./taujp;
  dtyNLin[_jp] = wInfj/taujp;
 
  CepsReal gNa   = 31*m_dconstants[_gNaFrac];
  CepsReal ENa   = log(m_constants[_nao]/y[_nai])/m_dconstants[_FRT];
  *INa = gNa*(v-ENa)*pow(y[_m],3.)*((1-fraction)*h*y[_j]+fraction*hp*y[_jp]);
 
  CepsReal wInfmL = sigmoid(v,-42.85,-5.264);
  CepsReal taumL  = taum;
  dtyLin [_mL]    = -1./taumL;
  dtyNLin[_mL]    = wInfmL/taumL;
  
  CepsReal wInfhL = sigmoid(v,-87.61,7.488);
  CepsReal tauhL  = 200.;
  dtyLin [_hL]    = -1./tauhL;
  dtyNLin[_hL]    = wInfhL/tauhL;
  
  CepsReal wInfhLp = sigmoid(v,-93.81,7.488);
  CepsReal tauhLp  = 3.*tauhL;
  dtyLin [_hLp]    = -1./tauhLp;
  dtyNLin[_hLp]    = wInfhLp/tauhLp;
 
  CepsReal gNaL = 0.0075*m_constants[_INaL_factor]*m_dconstants[_gNaLFrac];  
  if (m_type==ORdmD16::Type::Epi)
    gNaL *= 0.6;
  *INaL = gNaL*(v-ENa)*y[_mL]*((1.-fraction)*y[_hL]+fraction*y[_hLp]);
 
  return;
 
}
 
void
ORdmD16::computeIto(
  CepsReal* Ito, 
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin,
  CepsReal  fraction
) const
{
 
  CepsReal wInfa = sigmoid(v,14.34,-14.82);
  CepsReal taua  = 1.0515/(sigmoid(v,18.4099,-29.3814)/1.2089 + sigmoid(v,-100,29.3814)*3.5); 
  dtyLin [_a]    = -1./taua;
  dtyNLin[_a]    = wInfa/taua;
  
  CepsReal wInfi    = sigmoid(v,-43.95,5.711);
  CepsReal deltaEpi = 1.;
  if (m_type==ORdmD16::Type::Epi)
    deltaEpi -= 0.95*sigmoid(v,-70,5.);
  CepsReal tauiF = deltaEpi*(4.562+1/(0.3933*exp((-(v+100.0))/100.0)+0.08004*exp((v+50.0)/16.59)));
  CepsReal tauiS = deltaEpi*(23.62+1/(0.001416*exp((-(v+96.52))/59.05)+1.780e-8*exp((v+114.1)/8.079)));
  dtyLin [_iF]   = -1./tauiF;
  dtyNLin[_iF]   = wInfi/tauiF;
  dtyLin [_iS]   = -1./tauiS;
  dtyNLin[_iS]   = wInfi/tauiS;
  
  CepsReal AiF    = sigmoid(v,213.6,151.2);
  CepsReal i      = AiF*y[_iF]+(1.-AiF)*y[_iS];
  CepsReal wInfap = sigmoid(v,24.34,-14.82);
  dtyLin [_ap]    = -1./taua;
  dtyNLin[_ap]    = wInfap/taua;
  
  CepsReal dtiDevelop = 1.354+1.0e-4/(exp((v-167.4)/15.89)+exp(-(v-12.23)/0.2154));
  CepsReal dtiRecover = 1.0-0.5*sigmoid(v,-70,20);
  CepsReal tauiFp     = dtiDevelop*dtiRecover*tauiF;
  CepsReal tauiSp     = dtiDevelop*dtiRecover*tauiS;
  dtyLin [_iFp]       = -1./tauiFp;
  dtyNLin[_iFp]       = wInfi/tauiFp;
  dtyLin [_iSp]       = -1./tauiSp;
  dtyNLin[_iSp]       = wInfi/tauiSp;
  
  CepsReal ip  = AiF*y[_iFp]+(1.-AiF)*y[_iSp];
  CepsReal gto = 0.02*m_dconstants[_gtoFrac];
  CepsReal EK  = log(m_constants[_ko ]/y[_ki ])/m_dconstants[_FRT];
  if (m_type!=ORdmD16::Type::Endo)
    gto *= 4.0;
  *Ito   = gto*(v-EK)*((1-fraction)*i*y[_a]+fraction*ip*y[_ap]);
 
  return;
 
}
 
 
void
ORdmD16::computeICa(
  CepsReal* ICaL,
  CepsReal* ICaNa,
  CepsReal* ICaK,
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin,
  CepsReal  fraction
) const
{
 
  CepsReal wInfd = sigmoid(v,-3.940,-4.23);
  CepsReal taud  = 0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0)));
  dtyLin [_d]    = -1./taud;
  dtyNLin[_d]    = wInfd/taud;
 
  CepsReal wInff = sigmoid(v,-19.58,3.696);
  CepsReal tauff = 7.0+1.0/(0.0045*exp(-(v+20.0)/10.0)+0.0045*exp((v+20.0)/10.0));
  CepsReal taufs = 1000.0+1.0/(0.000035*exp(-(v+5.0)/4.0)+0.000035*exp((v+5.0)/6.0));
  dtyLin [_ff]   = -1./tauff;
  dtyNLin[_ff]   = wInff/tauff;
  dtyLin [_fs]   = -1./taufs;
  dtyNLin[_fs]   = wInff/taufs;
 
  CepsReal wInffca  = wInff;
  CepsReal taufcaf  = 7.0+1.0/(0.04*exp(-(v-4.0)/7.0)+0.04*exp((v-4.0)/7.0));
  CepsReal taufcas  = 100.0+1.0/(0.00012*exp(-v/3.0)+0.00012*exp(v/7.0));
  dtyLin [_fcaf]    = -1./taufcaf;
  dtyNLin[_fcaf]    = wInffca/taufcaf;
  dtyLin [_fcas]    = -1./taufcas;
  dtyNLin[_fcas]    = wInffca/taufcas;
 
  CepsReal taujca   = 75.;
  dtyLin [_jca]     = -1./taujca;
  dtyNLin[_jca]     = wInffca/taujca;
 
  CepsReal tauffp   = 2.5*tauff;
  dtyLin [_ffp]     = -1./tauffp;
  dtyNLin[_ffp]     = wInff/tauffp;
 
  CepsReal taufcafp = 2.5*taufcaf;
  dtyLin [_fcafp]   = -1./taufcafp;
  dtyNLin[_fcafp]   = wInffca/taufcafp;
 
  CepsReal Aff   = 0.6;
  CepsReal f     = Aff*y[_ff]+(1.-Aff)*y[_fs];
  CepsReal Afcaf = 0.3+0.6*sigmoid(v,10.,10.);
  CepsReal fca   = Afcaf*y[_fcaf ]+(1.-Afcaf)*y[_fcas];
  CepsReal fp    = Aff  *y[_ffp  ]+(1.-Aff  )*y[_fs  ];
  CepsReal fcap  = Afcaf*y[_fcafp]+(1.-Afcaf)*y[_fcas];
  CepsReal Kmn   = 0.002;
  CepsReal k2n   = 1000.0;
  CepsReal km2n  = y[_jca];
  CepsReal anca  = 1./(k2n/km2n+pow(1+Kmn/y[_cass],4.));
  dtyLin [_nca]  = -km2n;
  dtyNLin[_nca]  = anca*k2n;
 
  CepsReal vfrt    = v*m_dconstants[_FRT];
  CepsReal vffrt   = vfrt*m_constants[_F];
  CepsReal e1vfrt  = exp(vfrt);
  CepsReal e2vfrt  = e1vfrt*e1vfrt;
  CepsReal phiCaL  = 4.0*vffrt*(     y[_cass]*e2vfrt-0.341*m_constants[_cao])/(e2vfrt-1);
  CepsReal phiCaNa = 1.0*vffrt*(0.75*y[_nass]*e1vfrt-0.750*m_constants[_nao])/(e1vfrt-1);
  CepsReal phiCaK  = 1.0*vffrt*(0.75*y[_kss ]*e1vfrt-0.750*m_constants[_ko ])/(e1vfrt-1);
  CepsReal PCa     = 1.e-4;
  if (m_type==ORdmD16::Type::Epi)
    PCa *= 1.2;
  else if (m_type==ORdmD16::Type::MidMyo)
    PCa *= 2.5;
  CepsReal PCap   = 1.100e+0*PCa;
  CepsReal PCaNa  = 1.250e-3*PCa;
  CepsReal PCaK   = 3.574e-4*PCa;
  CepsReal PCaNap = 1.250e-3*PCap;
  CepsReal PCaKp  = 3.574e-4*PCap;
                   
  *ICaNa = (1.0-fraction)*PCaNa *phiCaNa*y[_d]*(f *(1.0-y[_nca])+y[_jca]*fca *y[_nca])
         +      fraction *PCaNap*phiCaNa*y[_d]*(fp*(1.0-y[_nca])+y[_jca]*fcap*y[_nca]);
  *ICaK  = (1.0-fraction)*PCaK  *phiCaK *y[_d]*(f *(1.0-y[_nca])+y[_jca]*fca *y[_nca])
         +      fraction *PCaKp *phiCaK *y[_d]*(fp*(1.0-y[_nca])+y[_jca]*fcap*y[_nca]);
  *ICaL  = (1.0-fraction)*PCa   *phiCaL *y[_d]*(f *(1.0-y[_nca])+y[_jca]*fca *y[_nca])
         +      fraction *PCap  *phiCaL *y[_d]*(fp*(1.0-y[_nca])+y[_jca]*fcap*y[_nca]);
  *ICaL *= m_dconstants[_gCalFrac]*m_constants[_ICaL_factor];
 
  return;
 
}
 
void
ORdmD16::computeIKr(
  CepsReal* IKr, 
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin
) const
{
 
  CepsReal wInfxr = sigmoid(v,-8.337,-6.789);
  CepsReal tauxrf = 12.98+1.0/(0.3652*exp((v-31.66)/3.869)+4.123e-5*exp((-(v-47.78))/20.38));
  CepsReal tauxrs = 1.865+1.0/(0.06629*exp((v-34.70)/7.355)+1.128e-5*exp((-(v-29.74))/25.94));
  dtyLin [_xrf]   = -1./tauxrf;
  dtyNLin[_xrf]   = wInfxr/tauxrf;
  dtyLin [_xrs]   = -1./tauxrs;
  dtyNLin[_xrs]   = wInfxr/tauxrs;
 
  CepsReal Axrf = sigmoid(v,-54.81,38.21);
  CepsReal xr   = Axrf*y[_xrf] + (1-Axrf)*y[_xrs];
  CepsReal rkr  = sigmoid(v,-55,75)*sigmoid(v,10,30);
  CepsReal EK   = log(m_constants[_ko ]/y[_ki ])/m_dconstants[_FRT];
  CepsReal GKr  = 0.046*m_constants[_IKr_factor]*m_dconstants[_gKrFrac];
  if (m_type==ORdmD16::Type::Epi)
    GKr *= 1.3;
  else if (m_type==ORdmD16::Type::MidMyo)
    GKr *= 0.8;
  *IKr = GKr*sqrt(m_constants[_ko]/5.4)*xr*rkr*(v-EK);
 
  return;
 
}
 
void
ORdmD16::computeIKs(
  CepsReal* IKs, 
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin
) const
{
 
  CepsReal wInfxs1 = sigmoid(v,-11.60,-8.932);
  CepsReal tauxs1  = 817.3+1.0/(2.326e-4*exp((v+48.28)/17.80)+0.001292*exp((-(v+210.0))/230.0));
  dtyLin [_xs1]    = -1./tauxs1;
  dtyNLin[_xs1]    = wInfxs1/tauxs1;
 
  CepsReal wInfxs2 = wInfxs1;
  CepsReal tauxs2  = 1.0/(0.01*exp((v-50.0)/20.0)+0.0193*exp((-(v+66.54))/31.0));
  dtyLin [_xs2]    = -1./tauxs2;
  dtyNLin[_xs2]    = wInfxs2/tauxs2;
 
  CepsReal KsCa = 1.0+0.6/(1.0+std::pow(3.8e-5/y[_cai],1.4));
  CepsReal EKs  = log((m_constants[_ko]+m_constants[_PKNa]*m_constants[_nao])
                       /(y[_ki ]+m_constants[_PKNa]*y[_nai]))/m_dconstants[_FRT];
  CepsReal GKs  = 0.0034*m_constants[_IKs_factor]*m_dconstants[_gKsFrac];
  if (m_type==ORdmD16::Type::Epi)
    GKs *= 1.4;
  *IKs = GKs*KsCa*y[_xs1]*y[_xs2]*(v-EKs);
 
  return;
 
}
 
void
ORdmD16::computeIK1(
  CepsReal* IK1, 
  CepsReal  v,
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin
) const
{
 
  CepsReal wInfxk1 = sigmoid(v,-2.5538*m_constants[_ko]-144.59,-1.5692*m_constants[_ko]-3.8115);
  CepsReal tauxk1  = 122.2/(exp((-(v+127.2))/20.36)+exp((v+236.8)/69.33));
  dtyLin [_xk1] = -1./tauxk1;
  dtyNLin[_xk1] = wInfxk1/tauxk1;
 
  CepsReal rk1 = sigmoid(v,-105.8+2.6*m_constants[_ko],9.493);
  CepsReal EK  = log(m_constants[_ko ]/y[_ki ])/m_dconstants[_FRT];
  CepsReal GK1 = 0.1908*m_constants[_IK1_factor]*m_dconstants[_gK1Frac];
  if (m_type==ORdmD16::Type::Epi)
    GK1 *= 1.2;
  else if (m_type==ORdmD16::Type::MidMyo)
    GK1 *= 1.3;
  *IK1 = GK1*sqrt(m_constants[_ko])*rk1*y[_xk1]*(v-EK);
 
  return;
 
}
 
void
ORdmD16::computeINaCa(CepsReal* INaCa_i, CepsReal* INaCa_ss, CepsReal  v, CepsReal* y) const
{
  *INaCa_i  = 0.8*computeINaCaInternal(v,y[_nai ],y[_cai ]);
  *INaCa_ss = 0.2*computeINaCaInternal(v,y[_nass],y[_cass]);
  return;
}
 
CepsReal
ORdmD16::computeINaCaInternal(CepsReal v, CepsReal na, CepsReal ca) const
{
 
  CepsReal kna1    = 15.0;
  CepsReal kna2    = 5.0;
  CepsReal kna3    = 88.12;
  CepsReal kasymm  = 12.5;
  CepsReal wna     = 6.e4;
  CepsReal wca     = 6.e4;
  CepsReal wnaca   = 5.e3;
  CepsReal kcaon   = 1.5e6;
  CepsReal kcaoff  = 5.e3;
  CepsReal qca     = 0.1670;
  CepsReal qna     = 0.5224;
  CepsReal hca     = exp(qca*v*m_dconstants[_FRT]);
  CepsReal hna     = exp(qna*v*m_dconstants[_FRT]);
 
  CepsReal h1      = 1+na/kna3*(1+hna);
  CepsReal h2      = na*hna/(kna3*h1);
  CepsReal h3      = 1.0/h1;
  CepsReal h4      = 1.0+na/kna1*(1+na/kna2);
  CepsReal h5      = na*na/(h4*kna1*kna2);
  CepsReal h6      = 1.0/h4;
  CepsReal h7      = 1.0+m_constants[_nao]/kna3*(1.0+1.0/hna);
  CepsReal h8      = m_constants[_nao]/(kna3*hna*h7);
  CepsReal h9      = 1.0/h7;
  CepsReal h10     = kasymm+1.0+m_constants[_nao]/kna1*(1.0+m_constants[_nao]/kna2);
  CepsReal h11     = m_constants[_nao]*m_constants[_nao]/(h10*kna1*kna2);
  CepsReal h12     = 1.0/h10;
  CepsReal k1      = h12*m_constants[_cao]*kcaon;
  CepsReal k2      = kcaoff;
  CepsReal k3p     = h9*wca;
  CepsReal k3pp    = h8*wnaca;
  CepsReal k3      = k3p+k3pp;
  CepsReal k4p     = h3*wca/hca;
  CepsReal k4pp    = h2*wnaca;
  CepsReal k4      = k4p+k4pp;
  CepsReal k5      = kcaoff;
  CepsReal k6      = h6*ca*kcaon;
  CepsReal k7      = h5*h2*wna;
  CepsReal k8      = h8*h11*wna;
  CepsReal x1      = k2*k4*(k7+k6)+k5*k7*(k2+k3);
  CepsReal x2      = k1*k7*(k4+k5)+k4*k6*(k1+k8);
  CepsReal x3      = k1*k3*(k7+k6)+k8*k6*(k2+k3);
  CepsReal x4      = k2*k8*(k4+k5)+k3*k5*(k1+k8);
  CepsReal E1      = x1/(x1+x2+x3+x4);
  CepsReal E2      = x2/(x1+x2+x3+x4);
  CepsReal E3      = x3/(x1+x2+x3+x4);
  CepsReal E4      = x4/(x1+x2+x3+x4);
  CepsReal KmCaAct = 150.0e-6;
  CepsReal allo    = 1.0/(1.0+pow(KmCaAct/ca,2.0));
  CepsReal JncxNa  = 3.0*(E4*k7-E1*k8)+E3*k4pp-E2*k3pp;
  CepsReal JncxCa  = E2*k2-E1*k1;
  CepsReal zna     = 1.0;
  CepsReal zca     = 2.0;
  CepsReal Gncx    = 0.0008;
  if (m_type==ORdmD16::Type::Epi)
    Gncx *= 1.1;
  else if (m_type==ORdmD16::Type::MidMyo)
    Gncx *= 1.4;
  return Gncx*allo*(zna*JncxNa+zca*JncxCa);
 
}
 
void
ORdmD16::computeINaK(CepsReal* INaK, CepsReal  v, CepsReal* y) const
{
 
  CepsReal k1p     = 949.5;
  CepsReal k1m     = 182.4;
  CepsReal k2p     = 687.2;
  CepsReal k2m     = 39.4;
  CepsReal k3p     = 1899.0;
  CepsReal k3m     = 79300.0;
  CepsReal k4p     = 639.0;
  CepsReal k4m     = 40.0;
  CepsReal Knai0   = 9.073;
  CepsReal Knao0   = 27.78;
  CepsReal delta   = -0.1550;
  CepsReal evfrt   = exp(v*m_dconstants[_FRT]/3.0);
  CepsReal Knai    = Knai0*pow(evfrt,  delta);
  CepsReal Knao    = Knao0*pow(evfrt,1-delta);
  CepsReal Kki     = 0.5;
  CepsReal Kko     = 0.3582;
  CepsReal MgADP   = 0.05;
  CepsReal MgATP   = 9.8;
  CepsReal Kmgatp  = 1.698e-7;
  CepsReal H       = 1.0e-7;
  CepsReal eP      = 4.2;
  CepsReal Khp     = 1.698e-7;
  CepsReal Knap    = 224.0;
  CepsReal Kxkur   = 292.0;
  CepsReal P       = eP/(1.0+H/Khp+y[_nai]/Knap+y[_ki]/Kxkur);
  CepsReal a1      = k1p*pow(y[_nai]/Knai,3.0)/(pow(1.0+y[_nai]/Knai,3.0)+pow(1.0+y[_ki]/Kki,2.0)-1.0);
  CepsReal b1      = k1m*MgADP;
  CepsReal a2      = k2p;
  CepsReal b2      = k2m*pow(m_constants[_nao]/Knao,3.0)/(pow(1.0+m_constants[_nao]/Knao,3.0)+pow(1.0+y[_ko]/Kko,2.0)-1.0);
  CepsReal a3      = k3p*pow(m_constants[_ko ]/Kko ,2.0)/(pow(1.0+m_constants[_nao]/Knao,3.0)+pow(1.0+y[_ko]/Kko,2.0)-1.0);
  CepsReal b3      = k3m*P*H/(1.0+MgATP/Kmgatp);
  CepsReal a4      = k4p*MgATP/Kmgatp/(1.0+MgATP/Kmgatp);
  CepsReal b4      = k4m*pow(y[_ki]/Kki,2.0)/(pow(1.0+y[_nai]/Knai,3.0)+pow(1.0+y[_ki]/Kki,2.0)-1.0);
  CepsReal x1      = a4*a1*a2+b2*b4*b3+a2*b4*b3+b3*a1*a2;
  CepsReal x2      = b2*b1*b4+a1*a2*a3+a3*b1*b4+a2*a3*b4;
  CepsReal x3      = a2*a3*a4+b3*b2*b1+b2*b1*a4+a3*a4*b1;
  CepsReal x4      = b4*b3*b2+a3*a4*a1+b2*a4*a1+b3*b2*a1;
  CepsReal E1      = x1/(x1+x2+x3+x4);
  CepsReal E2      = x2/(x1+x2+x3+x4);
  CepsReal E3      = x3/(x1+x2+x3+x4);
  CepsReal E4      = x4/(x1+x2+x3+x4);
  CepsReal zk      = 1.0;
  CepsReal JnakNa  = 3.0*(E1*a3-E2*b3);
  CepsReal JnakK   = 2.0*(E4*b1-E3*a1);
  CepsReal Pnak    = 30;
  CepsReal zna     = 1.0;
  if (m_type==ORdmD16::Type::Epi)
    Pnak *= 0.9;
  else if (m_type==ORdmD16::Type::MidMyo)
    Pnak *= 0.7;
 
  *INaK = Pnak*(zna*JnakNa+zk*JnakK);
 
  return;
}
 
void
ORdmD16::computeBgCurrents(CepsReal* IKb, CepsReal* INab, CepsReal* ICab, CepsReal  v, CepsReal* y) const
{
  CepsReal xkb = sigmoid(v,14.48,-18.34);
  CepsReal EK  = log(m_constants[_ko ]/y[_ki ])/m_dconstants[_FRT];
  CepsReal GKb = 0.003;
  if (m_type==ORdmD16::Type::Epi)
    GKb *= 0.6;
  *IKb = GKb*xkb*(v-EK);
 
  CepsReal vfrt   = v*m_dconstants[_FRT];
  CepsReal vffrt  = vfrt*m_constants[_F];
  CepsReal e1vfrt = exp(vfrt);
  CepsReal e2vfrt = e1vfrt*e1vfrt;
 
  CepsReal PNab = 3.75e-10;
  CepsReal PCab = 2.5e-8;
  *INab = PNab*1.0*vffrt*(y[_nai]*e1vfrt-      m_constants[_nao])/(e1vfrt-1);
  *ICab = PCab*4.0*vffrt*(y[_cai]*e2vfrt-0.341*m_constants[_cao])/(e2vfrt-1);
  return;
}
 
void
ORdmD16::computeIPCa(CepsReal* IPCa, CepsReal* y) const
{
  CepsReal GpCa = 5.e-4;
  *IPCa = GpCa*y[_cai]/(5.e-4+y[_cai]);
}
 
void
ORdmD16::computeJrel(
  CepsReal* Jrel,
  CepsReal  ICaL, 
  CepsReal* y, 
  CepsReal* dtyLin, 
  CepsReal* dtyNLin,
  CepsReal  fraction
) const
{
  CepsReal bt       = 4.75;
  CepsReal a_rel    = 0.5*bt;
  CepsReal wInfJrel = -a_rel*ICaL/(1.0+pow(1.5/y[_cajsr],8.0));
  if (m_type==ORdmD16::Type::MidMyo)
    wInfJrel *= 1.7;
  CepsReal tauJrel = std::max(0.001,bt/(1.0+0.0123/y[_cajsr]));
  dtyLin [_Jrelnp] = -1/tauJrel;
  dtyNLin[_Jrelnp] = wInfJrel/tauJrel;
 
  CepsReal btp       = 1.25*bt;
  CepsReal a_relp    = 0.5*btp;
  CepsReal wInfJrelp = -a_relp*ICaL/(1.0+pow(1.5/y[_cajsr],8.0));
  if (m_type==ORdmD16::Type::MidMyo)
    wInfJrelp *= 1.7;
  CepsReal tauJrelp = std::max(0.001,btp/(1.0+0.0123/y[_cajsr]));
  dtyLin [_Jrelp]   = -1/tauJrelp;
  dtyNLin[_Jrelp]   = wInfJrelp/tauJrelp;
  *Jrel   = (1-fraction)*y[_Jrelnp] + fraction*y[_Jrelp];
 
  return;
}
 
void
ORdmD16::computeJup(CepsReal* Jup, CepsReal* y, CepsReal fraction) const
{
  CepsReal Jupnp =      0.004375*y[_cai]/(y[_cai]+0.00092);
  CepsReal Jupp  = 2.75*0.004375*y[_cai]/(y[_cai]+0.00092-0.00017);
  if (m_type==ORdmD16::Type::Epi)
  {
    Jupnp *= 1.3;
    Jupp  *= 1.3;
  }
  CepsReal Jleak = 0.0039375*y[_cansr]/15.;
  *Jup = (1-fraction)*Jupnp + fraction*Jupp - Jleak;
 
  return;
}