CRN98.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/CRN98.hpp"
 
 
CRN98::CRN98(
  Unknown*                      u, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  AbstractIonicModel(u,attrs)
{
 
  m_name          = "Courtemanche-Ramirez-Nattel 1998";
  m_nStateVars    = 20;
  m_stateVarNames = {
    "Na_i [millimolar]",
    "m [adim]",
    "h [adim]",
    "j [adim]",
    "K_i [millimolar]",
    "oa [adim]",
    "oi [adim]",
    "ua [adim]",
    "ui [adim]",
    "xr [adim]",
    "xs [adim]",
    "Ca_i [millimolar]",
    "d [adim]",
    "f [adim]",
    "f_Ca [adim]",
    "Ca_rel [millimolar]",
    "u [adim]",
    "v [adim]",
    "w [adim]",
    "Ca_up [millimolar]"
  };
 
  // ================================================================================
  // Default parameters
  m_constants = ceps::newArray<CepsReal>(43);
  m_constants[_R       ] = 8.3143;
  m_constants[_T       ] = 310;
  m_constants[_F       ] = 96.4867;
  m_constants[_vcell   ] = 20100;
  m_constants[_vi      ] = 0.68  *m_constants[_vcell];
  m_constants[_vup     ] = 0.0552*m_constants[_vcell];
  m_constants[_vrel    ] = 0.0048*m_constants[_vcell];
  m_constants[_ko      ] = 5.4;
  m_constants[_nao     ] = 140;
  m_constants[_cao     ] = 1.8;
  m_constants[_gna     ] = 7.8;
  m_constants[_gk1     ] = 0.09;
  m_constants[_gto     ] = 0.1652;
  m_constants[_gkr     ] = 0.029411765;
  m_constants[_gks     ] = 0.12941176;
  m_constants[_gcal    ] = 0.12375;
  m_constants[_gbca    ] = 0.001131;
  m_constants[_gbna    ] = 0.0006744375;
  m_constants[_inakmax ] = 0.59933874;
  m_constants[_inacamax] = 1600;
  m_constants[_icapmax ] = 0.275;
  m_constants[_iupmax  ] = 0.005;
  m_constants[_kq10    ] = 3;
  m_constants[_sigma   ] = 1./7.*std::expm1(m_constants[_nao]/67.3);
  m_constants[_kmnai   ] = 10;
  m_constants[_kmko    ] = 1.5;
  m_constants[_kmna    ] = 87.5;
  m_constants[_kmca    ] = 1.38;
  m_constants[_ksat    ] = 0.1;
  m_constants[_krel    ] = 30;
  m_constants[_kup     ] = 0.00092;
  m_constants[_caupmax ] = 15;
  m_constants[_cmdnmax ] = 0.05;
  m_constants[_trpnmax ] = 0.07;
  m_constants[_csqnmax ] = 10;
  m_constants[_kmcmdn  ] = 0.00238;
  m_constants[_kmtrpn  ] = 0.0005;
  m_constants[_kmcsqn  ] = 0.8;
  m_constants[_taufca  ] = 2;
  m_constants[_tautr   ] = 180;
  m_constants[_gbk     ] = 0;
  m_constants[_gamma   ] = 0.35;
  m_constants[_tauu    ] = 8;
 
  // ================================================================================
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params,nullptr);
  
  // Assuming a radius of 10um. Vcell is in um3
  m_cellSurface  = 2*m_constants[_vcell]/10 + 2*M_PI*100; // in um2
  m_cellSurface *= 1.e-8; // in cm2
 
  m_paperCm   = 100; //pF
  m_paperStim = 2000; //pA
 
}
 
void 
CRN98::getInitialCondition(CepsReal* v, CepsReal* y) const
{
  *v = -81.18;
  y[_nai  ] = 1.117e+1;
  y[_m    ] = 2.908e-3;
  y[_h    ] = 9.649e-1;
  y[_j    ] = 9.775e-1;
  y[_ki   ] = 1.390e+2;
  y[_oa   ] = 3.043e-2;
  y[_oi   ] = 9.992e-1;
  y[_ua   ] = 4.966e-3;
  y[_ui   ] = 9.986e-1;
  y[_xr   ] = 3.296e-5;
  y[_xs   ] = 1.869e-2;
  y[_cai  ] = 1.013e-4;
  y[_d    ] = 1.367e-4;
  y[_f    ] = 9.996e-1;
  y[_fca  ] = 7.755e-1;
  y[_carel] = 1.488e+0;
  y[_caup ] = 1.488e+0;
  y[_u    ] = 2.35e-112;
  y[_v    ] = 1;
  y[_w    ] = 0.9992;
}
 
void
CRN98::computeRates(
  CepsReal         t,
  CepsReal*        y, 
  CepsReal*        vm, 
  CepsReal*        dtyLin, 
  CepsReal*        dtyNLin, 
  CepsReal*        dtv, 
  DegreeOfFreedom* dof
) const
{ 
 
  using std::exp;
  using std::expm1;
  using std::pow;
  using ceps::approxEquals;
 
  CepsReal* c = m_constants;
  CepsReal  v = *vm;
 
  CepsReal cm  = getCmInternal(t,dof);
  CepsReal RTF = c[_R]*c[_T]/c[_F];
 
  CepsReal fCaInf = 1/(1+y[_cai]/0.00035);
  dtyLin [_fca] = -1/c[_taufca];
  dtyNLin[_fca] = fCaInf/c[_taufca];
 
  CepsReal dInf = 1/(1+exp(-0.125*(v+10)));
  CepsReal texp = expm1(-(v+10.)/6.24);
  CepsReal taud = approxEquals(v,-10,1.e-10) ? 4.579/(2+texp) : -texp/(0.035*(v+10.)*(2+texp));
  dtyLin [_d] = -1/taud;
  dtyNLin[_d] = dInf/taud;
 
  CepsReal fexp = exp(-(v+28.)/6.9);
  CepsReal fInf = fexp/(1+fexp);
  CepsReal tauf = 9./(0.02+0.0197*exp(-pow(0.0337*(v+10.),2)));
  dtyLin [_f] = -1/tauf;
  dtyNLin[_f] = fInf/tauf;
 
  CepsReal wexp = expm1(-(v-7.9)/5.);
  CepsReal tauw = approxEquals(v,7.9,1.e-10) ? 1.2/1.3 : -6.*wexp/(1.3+0.3*wexp)/(v-7.9);
  CepsReal wInf = 1-1/(1+exp(-(v-40)/17.));
  dtyLin [_w] = -1/tauw;
  dtyNLin[_w] = wInf/tauw;
 
  CepsReal alpham = approxEquals(v,-47.13,1.e-10) ? 3.2 : -0.32*(v+47.13)/expm1(-0.1*(v+47.13));
  CepsReal betam  = 0.08*exp(-v/11.);
  dtyLin [_m] = -alpham-betam;
  dtyNLin[_m] = alpham;
 
  CepsReal alphah = v<-40 ? 0.135*exp(-(v+80)/6.8) : 0.;
  CepsReal betah  = v<-40 ? 3.56*exp(0.079*v) + 3.1e5*exp(0.35*v) : 1/(0.13*(1+exp(-(v+10.66)/11.1)));
  dtyLin [_h] = -alphah-betah;
  dtyNLin[_h] = alphah;
 
  CepsReal alphaj = v<-40 ? (-127140.*exp(0.2444*v)-3.474e-5*exp(-0.04391*v))*(v+37.78)
                            /(1+exp(0.311*(v+79.23))) : 0.;
  CepsReal betaj  = v<-40 ? 0.1212*exp(-0.01052*v)/(1+exp(-0.1378*(v+40.14)))
                          : 0.3*exp(-2.35e-7*v)/(1+exp(-0.1*(v+32)));
  dtyLin [_j] = -alphaj-betaj;
  dtyNLin[_j] = alphaj;
 
  CepsReal alphaoa = 0.65/(exp(-(v+10)/8.5)+exp(-(v-30)/59.));
  CepsReal betaoa  = 0.65/(2.5+exp((v+82)/17.));
  CepsReal tauoa   = 1/(alphaoa+betaoa) /c[_kq10];
  CepsReal oaInf   = 1/(1+exp(-(v+20.47)/17.54));
  dtyLin [_oa] = -1/tauoa;
  dtyNLin[_oa] = oaInf/tauoa;
 
  CepsReal alphaoi = 1/(18.53+exp((v+113.7)/10.95));
  CepsReal betaoi  = 1/(35.56+exp(-(v+1.26)/7.44));
  CepsReal tauoi   = 1/(alphaoi+betaoi) /c[_kq10];
  CepsReal oiInf   = 1/(1+exp((v+43.1)/5.3));
  dtyLin [_oi] = -1/tauoi;
  dtyNLin[_oi] = oiInf/tauoi;
 
  CepsReal tauua   = tauoa;
  CepsReal uaInf   = 1/(1+exp(-(v+30.3)/9.6));
  dtyLin [_ua] = -1/tauua;
  dtyNLin[_ua] = uaInf/tauua;
 
  CepsReal alphaui = 1/(21+exp(-(v-185.)/28.));
  CepsReal betaui  = exp((v-158)/16.);
  CepsReal tauui   = 1/(alphaui+betaui) /c[_kq10];
  CepsReal uiInf   = 1/(1+exp((v-99.45)/27.48));
  dtyLin [_ui] = -1/tauui;
  dtyNLin[_ui] = uiInf/tauui;
 
  CepsReal alphaxr = approxEquals(v,-14.10,1.e-10) ? 0.0015
                                                   : -3.e-4*(v+14.10)/expm1(-0.2*(v+14.10));
  CepsReal betaxr  = approxEquals(v,3.3328,1.e-10) ? 0.000378361
                                                   : 7.3898e-5*(v-3.3328)/expm1((v-3.3328)/5.1237);
  CepsReal tauxr = 1/(alphaxr+betaxr);
  CepsReal xrInf = 1/(1+exp(-(v+14.1)/6.5));
  dtyLin [_xr] = -1/tauxr;
  dtyNLin[_xr] = xrInf/tauxr;
  
  CepsReal alphaxs = approxEquals(v,19.9,1.e-10) ? 6.80e-4 : -4e-5*(v-19.9)/expm1(-(v-19.9)/17);
  CepsReal betaxs  = approxEquals(v,19.9,1.e-10) ? 3.15e-4 : 3.5e-5*(v-19.9)/expm1((v-19.9)/9);
  CepsReal tauxs   = 0.5/(alphaxs+betaxs);
  CepsReal xsInf   = 1/std::sqrt(1+exp(-(v-19.9)/12.7));
  dtyLin [_xs] = -1/tauxs;
  dtyNLin[_xs] = xsInf/tauxs;
 
  CepsReal Ek   = RTF*std::log(c[_ko]/y[_ki]);
  CepsReal ik1  = cm*c[_gk1]*(v-Ek)/(1+exp(0.07*(v+80)));
  CepsReal ito  = cm*c[_gto]*y[_oa]*y[_oa]*y[_oa]*y[_oi]*(v-Ek);
  CepsReal gkur = 0.005+0.05/(1+exp(-(v-15)/13));
  CepsReal ikur = cm*gkur*y[_ua]*y[_ua]*y[_ua]*y[_ui]*(v-Ek);
  CepsReal ikr  = cm*c[_gkr]*y[_xr]*(v-Ek)/(1+exp((v+15)/22.4));
  CepsReal iks  = cm*c[_gks]*y[_xs]*y[_xs]*(v-Ek);
  CepsReal fnak = 1/(1+0.1245*exp(-0.1/RTF*v)+0.0365*c[_sigma]*exp(-v/RTF));
  CepsReal inak = cm*fnak*c[_inakmax]/(1+pow(c[_kmnai]/y[_nai],1.5)) * c[_ko]/(c[_ko]+c[_kmko]);
  CepsReal ibk  = cm*c[_gbk]*(v-Ek);
  dtyLin [_ki] = 0.; 
  dtyNLin[_ki] = (2*inak-ik1-ito-ikur-ikr-iks-ibk)/(c[_vi]*c[_F]);
 
  CepsReal ena   = RTF*std::log(c[_nao]/y[_nai]);
  CepsReal ina   = cm*c[_gna]*y[_m]*y[_m]*y[_m]*y[_h]*y[_j]*(v-ena);
  CepsReal inaca = cm*c[_inacamax]*
                  (
                    exp(c[_gamma]*v/RTF)*y[_nai]*y[_nai]*y[_nai]*c[_cao]
                   -exp((c[_gamma]-1)*v/RTF)*c[_nao]*c[_nao]*c[_nao]*y[_cai]
                ) / (
                    (c[_kmna]*c[_kmna]*c[_kmna] + c[_nao]*c[_nao]*c[_nao])*(c[_kmca]+c[_cao])
                   *(1+c[_ksat]*exp((c[_gamma]-1)*v/RTF))
                  );
  CepsReal ibna  = cm*c[_gbna]*(v-ena);
  dtyLin [_nai] = 0.;
  dtyNLin[_nai] = (-3*inak-3*inaca-ibna-ina)/(c[_vi]*c[_F]);
 
  CepsReal ical = cm*c[_gcal]*y[_d]*y[_f]*y[_fca]*(v-65);
  CepsReal icap = cm*c[_icapmax]*y[_cai]/(5.e-4+y[_cai]);
  CepsReal eca  = 0.5*RTF*std::log(c[_cao]/y[_cai]);
  CepsReal ibca = cm*c[_gbca]*(v-eca);
 
  CepsReal iStim = getStimulation(dof,t);
  CepsReal iIon  = -(ina+ik1+ito+ikur+ikr+iks+ibna+ibca+inak+icap+inaca+ical);
  #ifdef DEBUG
    ceps::checkNanOrInf(iIon,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
 
  *dtv = iIon + iStim;
  *dtv /= cm;
 
  CepsReal irel = c[_krel]*y[_u]*y[_u]*y[_v]*y[_w]*(y[_carel]-y[_cai]);
  CepsReal itr  = (y[_caup]-y[_carel])/c[_tautr];
  dtyLin [_carel] = 0.;
  dtyNLin[_carel] = (itr-irel)/(1+(c[_csqnmax]*c[_kmcsqn])/pow(y[_carel]+c[_kmcsqn],2));
 
  CepsReal fn   = 1.e-12*(c[_vrel]*irel-0.5/c[_F]*(0.5*ical-0.2*inaca));
  CepsReal uInf = 1/(1+exp(-(fn-3.4175e-13)/1.367e-15));
  dtyLin [_u] = -1/c[_tauu];
  dtyNLin[_u] = uInf/c[_tauu];
 
  CepsReal tauv = 1.91+2.09/(1+exp(-(fn-3.4175e-13)/1.367e-15));
  CepsReal vInf = 1-1/(1+exp(-(fn-6.835e-14)/1.367e-15));
  dtyLin [_v] = -1/tauv;    
  dtyNLin[_v] = vInf/tauv;
 
  CepsReal iup     = c[_iupmax]/(1+c[_kup]/y[_cai]);
  CepsReal iupleak = (c[_iupmax]*y[_caup])/c[_caupmax];
  dtyLin [_caup] = 0.;
  dtyNLin[_caup] = iup-(iupleak+(itr*c[_vrel])/c[_vup]);
 
  CepsReal b1 = (2*inaca-(icap+ical+ibca)) / (2*c[_vi]*c[_F])
              + (c[_vup]*(iupleak-iup)+irel*c[_vrel])/c[_vi];
  CepsReal b2 = 1 + c[_trpnmax]*c[_kmtrpn]/pow(y[_cai]+c[_kmtrpn],2)
                  + c[_cmdnmax]*c[_kmcmdn]/pow(y[_cai]+c[_kmcmdn],2);
  dtyLin [_cai] = 0.;
  dtyNLin[_cai] = b1/b2;
 
}
 
CepsReal
CRN98::convertCurrentFromCepsUnit(const CepsReal& i) const
{
  return i*m_cellSurface*1.e6;
}
 
CepsReal
CRN98::convertCmFromCepsUnit(const CepsReal& cm) const
{
  return 1.e6*m_cellSurface*cm;
}
 
CepsReal
CRN98::convertCmToCepsUnit(const CepsReal& cm) const
{
  return cm/m_cellSurface*1.e-6;
}
 
void
CRN98::setupWithParameters(InputParameters* params, FunctionDictionary* dico)
{
 
  m_constants[_R       ] = params->getReal("CRN98 R       ",m_constants[_R       ]); 
  m_constants[_T       ] = params->getReal("CRN98 T       ",m_constants[_T       ]); 
  m_constants[_F       ] = params->getReal("CRN98 F       ",m_constants[_F       ]); 
  m_constants[_vcell   ] = params->getReal("CRN98 vcell   ",m_constants[_vcell   ]); 
  m_constants[_vi      ] = params->getReal("CRN98 vi      ",m_constants[_vi      ]); 
  m_constants[_vup     ] = params->getReal("CRN98 vup     ",m_constants[_vup     ]); 
  m_constants[_vrel    ] = params->getReal("CRN98 vrel    ",m_constants[_vrel    ]); 
  m_constants[_ko      ] = params->getReal("CRN98 ko      ",m_constants[_ko      ]); 
  m_constants[_nao     ] = params->getReal("CRN98 nao     ",m_constants[_nao     ]); 
  m_constants[_cao     ] = params->getReal("CRN98 cao     ",m_constants[_cao     ]); 
  m_constants[_gna     ] = params->getReal("CRN98 gna     ",m_constants[_gna     ]); 
  m_constants[_gk1     ] = params->getReal("CRN98 gk1     ",m_constants[_gk1     ]); 
  m_constants[_gto     ] = params->getReal("CRN98 gto     ",m_constants[_gto     ]); 
  m_constants[_gkr     ] = params->getReal("CRN98 gkr     ",m_constants[_gkr     ]); 
  m_constants[_gks     ] = params->getReal("CRN98 gks     ",m_constants[_gks     ]); 
  m_constants[_gcal    ] = params->getReal("CRN98 gcal    ",m_constants[_gcal    ]); 
  m_constants[_gbca    ] = params->getReal("CRN98 gbca    ",m_constants[_gbca    ]); 
  m_constants[_gbna    ] = params->getReal("CRN98 gbna    ",m_constants[_gbna    ]); 
  m_constants[_inakmax ] = params->getReal("CRN98 inakmax ",m_constants[_inakmax ]); 
  m_constants[_inacamax] = params->getReal("CRN98 inacamax",m_constants[_inacamax]); 
  m_constants[_icapmax ] = params->getReal("CRN98 icapmax ",m_constants[_icapmax ]); 
  m_constants[_iupmax  ] = params->getReal("CRN98 iupmax  ",m_constants[_iupmax  ]); 
  m_constants[_kq10    ] = params->getReal("CRN98 kq10    ",m_constants[_kq10    ]); 
  m_constants[_sigma   ] = params->getReal("CRN98 sigma   ",m_constants[_sigma   ]); 
  m_constants[_kmnai   ] = params->getReal("CRN98 kmnai   ",m_constants[_kmnai   ]); 
  m_constants[_kmko    ] = params->getReal("CRN98 kmko    ",m_constants[_kmko    ]); 
  m_constants[_kmna    ] = params->getReal("CRN98 kmna    ",m_constants[_kmna    ]); 
  m_constants[_kmca    ] = params->getReal("CRN98 kmca    ",m_constants[_kmca    ]); 
  m_constants[_ksat    ] = params->getReal("CRN98 ksat    ",m_constants[_ksat    ]); 
  m_constants[_krel    ] = params->getReal("CRN98 krel    ",m_constants[_krel    ]); 
  m_constants[_kup     ] = params->getReal("CRN98 kup     ",m_constants[_kup     ]); 
  m_constants[_caupmax ] = params->getReal("CRN98 caupmax ",m_constants[_caupmax ]); 
  m_constants[_cmdnmax ] = params->getReal("CRN98 cmdnmax ",m_constants[_cmdnmax ]); 
  m_constants[_trpnmax ] = params->getReal("CRN98 trpnmax ",m_constants[_trpnmax ]); 
  m_constants[_csqnmax ] = params->getReal("CRN98 csqnmax ",m_constants[_csqnmax ]); 
  m_constants[_kmcmdn  ] = params->getReal("CRN98 kmcmdn  ",m_constants[_kmcmdn  ]); 
  m_constants[_kmtrpn  ] = params->getReal("CRN98 kmtrpn  ",m_constants[_kmtrpn  ]); 
  m_constants[_kmcsqn  ] = params->getReal("CRN98 kmcsqn  ",m_constants[_kmcsqn  ]); 
  m_constants[_taufca  ] = params->getReal("CRN98 taufca  ",m_constants[_taufca  ]); 
  m_constants[_tautr   ] = params->getReal("CRN98 tautr   ",m_constants[_tautr   ]); 
  m_constants[_gbk     ] = params->getReal("CRN98 gbk     ",m_constants[_gbk     ]); 
  m_constants[_gamma   ] = params->getReal("CRN98 gamma   ",m_constants[_gamma   ]); 
  m_constants[_tauu    ] = params->getReal("CRN98 tauu    ",m_constants[_tauu    ]); 
  m_cellSurface          = params->getReal("CRN98 cell surface", m_cellSurface    ); 
 
  return;
}