TTP06.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/TTP06.hpp"
 
TTP06::TTP06(
  const TTP06::Type&            type,
  Unknown*                      u, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  AbstractIonicModel(u,attrs),
  m_type(type)
{
  if (m_type == Type::Endo) 
    m_name = "TenTusscher Panfilov 2006 (endocardium)";
  else if (m_type == Type::MidMyo)
    m_name = "TenTusscher Panvilov 2006 (mid myocardium)";
  else 
    m_name = "TenTusscher Panvilov 2006 (epicardium)";
 
  m_nStateVars    = 18;
  m_stateVarNames = {
    "K_i [millimolar]",
    "Na_i [millimolar]",
    "Ca_i [millimolar]",
    "xr1 [adim]",
    "xr2 [adim]",
    "xs [adim]",
    "m [adim]",
    "h [adim]",
    "j [adim]",
    "Ca_ss [millimolar]",
    "d [adim]",
    "f [adim]",
    "f2 [adim]",
    "f_Ca_ss [adim]",
    "s [adim]",
    "r [adim]",
    "Ca_sr [millimolar]",
    "R_prime [adim]"
  };
 
  // ================================================================================
  // Default parameters
  m_constants = new CepsReal[53];
  m_constants[_R      ] = 8314.472;
  m_constants[_T      ] = 310;
  m_constants[_F      ] = 96485.3415;
  m_constants[_vc     ] = 0.016404;
  m_constants[_pkna   ] = 0.03;
  m_constants[_ko     ] = 5.4;
  m_constants[_nao    ] = 140;
  m_constants[_cao    ] = 2;
  m_constants[_gk1    ] = 5.405;
  m_constants[_gkr    ] = 0.153;
  m_constants[_gna    ] = 14.838;
  m_constants[_gbna   ] = 0.00029;
  m_constants[_gcal   ] = 0.0000398;
  m_constants[_gbca   ] = 0.000592;
  m_constants[_pnak   ] = 2.724;
  m_constants[_kmk    ] = 1;
  m_constants[_kmna   ] = 40;
  m_constants[_knaca  ] = 1000;
  m_constants[_ksat   ] = 0.1;
  m_constants[_alpha  ] = 2.5;
  m_constants[_gamma  ] = 0.35;
  m_constants[_kmca   ] = 1.38;
  m_constants[_kmnai  ] = 87.5;
  m_constants[_gpca   ] = 0.1238;
  m_constants[_kpca   ] = 0.0005;
  m_constants[_gpk    ] = 0.0146;
  m_constants[_k1prime] = 0.15;
  m_constants[_k2prime] = 0.045;
  m_constants[_k3     ] = 0.06;
  m_constants[_k4     ] = 0.005;
  m_constants[_ec     ] = 1.5;
  m_constants[_maxsr  ] = 2.5;
  m_constants[_minsr  ] = 1;
  m_constants[_vrel   ] = 0.102;
  m_constants[_vxfer  ] = 0.0038;
  m_constants[_kup    ] = 0.00025;
  m_constants[_vleak  ] = 0.00036;
  m_constants[_vmaxup ] = 0.006375;
  m_constants[_bufc   ] = 0.2;
  m_constants[_kbufc  ] = 0.001;
  m_constants[_bufsr  ] = 10;
  m_constants[_kbuf   ] = 0.3;
  m_constants[_bufss  ] = 0.4;
  m_constants[_kbufss ] = 0.00025;
  m_constants[_vsr    ] = 0.001094;
  m_constants[_vss    ] = 0.00005468;
  m_constants[_gks    ] = m_type == Type::MidMyo ? 0.098 : 0.392;
  m_constants[_gto    ] = m_type == Type::Endo   ? 0.073 : 0.294;
 
  // ================================================================================
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params,nullptr);
 
  // The cited paper mentions a surface to volume ratio of 0.2 um-1
  m_cellSurface  = 0.2*m_constants[_vc]*1.e-8; // from um2 to cm2
  // However the model mentions "per surface unit" for all constants. Surface unit is unclear.
  // so m_cellSurface is not used...
 
  m_paperCm   = 0.185;
  m_paperStim = 52*m_paperCm; // 
 
 
}
 
void 
TTP06::getInitialCondition(CepsReal* v, CepsReal* y) const
{
 
  if (m_type == Type::Endo)
  {  
    *v = -86.709;
    y[_ki    ] = 138.4;
    y[_nai   ] = 10.355;
    y[_cai   ] = 0.00013;
    y[_xr1   ] = 0.00448;
    y[_xr2   ] = 0.476;
    y[_xs    ] = 0.0087;
    y[_m     ] = 0.00155;
    y[_h     ] = 0.7573;
    y[_j     ] = 0.7225;
    y[_cass  ] = 0.00036;
    y[_d     ] = 3.164e-5;
    y[_f     ] = 0.8009;
    y[_f2    ] = 0.9778;
    y[_fcass ] = 0.9953;
    y[_s     ] = 0.3212;
    y[_r     ] = 2.235e-8;
    y[_casr  ] = 3.715;
    y[_rprime] = 0.9068;
  }
  else if (m_type == Type::MidMyo)
  {
    *v = -85.423;  
    y[_ki    ] = 138.52;
    y[_nai   ] = 10.132;
    y[_cai   ] = 0.000153;
    y[_xr1   ] = 0.0165;
    y[_xr2   ] = 0.473;
    y[_xs    ] = 0.0174;
    y[_m     ] = 0.00165;
    y[_h     ] = 0.749;
    y[_j     ] = 0.6788;
    y[_cass  ] = 0.00042;
    y[_d     ] = 3.288e-5;
    y[_f     ] = 0.7026;
    y[_f2    ] = 0.9526;
    y[_fcass ] = 0.9942;
    y[_s     ] = 0.999998;
    y[_r     ] = 2.347e-8;
    y[_casr  ] = 4.272;
    y[_rprime] = 0.8978;
  }
  else 
  {
    *v = -85.23;
    y[_ki    ] = 136.89;
    y[_nai   ] = 8.604;
    y[_cai   ] = 0.000126;
    y[_xr1   ] = 0.00621;
    y[_xr2   ] = 0.4712;
    y[_xs    ] = 0.0095;
    y[_m     ] = 0.00172;
    y[_h     ] = 0.7444;
    y[_j     ] = 0.7045;
    y[_cass  ] = 0.00036;
    y[_d     ] = 3.373e-5;
    y[_f     ] = 0.7888;
    y[_f2    ] = 0.9755;
    y[_fcass ] = 0.9953;
    y[_s     ] = 0.9999998;
    y[_r     ] = 2.42e-8;
    y[_casr  ] = 3.64;
    y[_rprime] = 0.9073;
  }
 
}
 
void
TTP06::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 a7  = sigmoid(v,-20,7);
  CepsReal a20 = 1102.5*gaussian(v,-27,225) +200*sigmoid(v,13,-10) +180*sigmoid(v,-30,10) + 20;
  dtyLin [_f] = -1/a20;
  dtyNLin[_f] = a7/a20;
 
  CepsReal a8  = 0.67*sigmoid(v,-35,7) + 0.33;
  CepsReal a21 = 562.*gaussian(v,-27,240) +31*sigmoid(v,25,-10) +80*sigmoid(v,-30,10);
  dtyLin [_f2] = -1/a21;
  dtyNLin[_f2] = a8/a21;
 
  CepsReal a9  = 0.6/(1+pow(y[_cass]/0.05,2)) +0.4;
  CepsReal a22 = 80./(1+pow(y[_cass]/0.05,2)) +2.0;
  dtyLin [_fcass] = -1/a22;
  dtyNLin[_fcass] = a9/a22;
 
  CepsReal a10 = m_type == Type::Endo ? sigmoid(v,-28,5) : sigmoid(v,-20.,5.);  
  CepsReal a23 = m_type == Type::Endo ?  8. + 1000.*gaussian(v,-67.,1000.) 
                                      : 85.*gaussian(v,-45.,320.) +5.*sigmoid(v,20.,5.) +3;
  dtyLin [_s] = -1/a23;
  dtyNLin[_s] = a10/a23;
 
  CepsReal a11 = sigmoid(v,20,-6);
  CepsReal a24 = 9.5*gaussian(v,-40,1800) + 0.8;
  dtyLin [_r] = -1/a24;
  dtyNLin[_r] = a11/a24;
 
  CepsReal a0  = sigmoid(v,-26,-7);
  CepsReal a13 = 450*sigmoid(v,-45,-10);
  CepsReal a26 = 6*sigmoid(v,-30,11.5);
  CepsReal a34 = a13*a26;
  dtyLin [_xr1] = -1/a34;
  dtyNLin[_xr1] = a0/a34;
 
  CepsReal a1  = sigmoid(v,-88,24);
  CepsReal a14 = 3*sigmoid(v,-60,-20);
  CepsReal a27 = 1.12*sigmoid(v,60,20);
  CepsReal a35 = a14*a27;
  dtyLin [_xr2] = -1/a35;
  dtyNLin[_xr2] = a1/a35;
 
  CepsReal a2  = sigmoid(v,-5,-14);
  CepsReal a15 = 1400*std::sqrt(sigmoid(v,5,-6));
  CepsReal a28 = sigmoid(v,35,15);
  CepsReal a36 = a15*a28+80;
  dtyLin [_xs] = -1/a36;
  dtyNLin[_xs] = a2/a36;
 
  CepsReal a3  = pow(sigmoid(v,-56.86,-9.03),2);
  CepsReal a16 = sigmoid(v,-60,-5);
  CepsReal a29 = 0.1*sigmoid(v,-35,5) + 0.1*sigmoid(v,50,200);
  CepsReal a37 = a16*a29;
  dtyLin [_m] = -1/a37;
  dtyNLin[_m] = a3/a37;
 
  CepsReal a4  = pow(sigmoid(v,-71.55,7.43),2);
  CepsReal a17 = v<-40 ? 0.057*exp(-(v+80)/6.8) : 0;
  CepsReal a30 = v<-40 ? 2.7*exp(0.079*v) +3.1e5*exp(0.3485*v) : 0.77/0.13*sigmoid(v,-10.66,-11.1);
  CepsReal a38 = 1/(a17+a30);
  dtyLin [_h] = -1/a38;
  dtyNLin[_h] = a4/a38;
 
  CepsReal a5  = a4;
  CepsReal a18 = v<-40 ? (-25428*exp(0.2444*v)-6.948e-6*exp(-0.04391*v))*(v+37.78)*sigmoid(v,-79.23,1/0.311) : 0;
  CepsReal a31 = v<-40 ? 0.02424*exp(-0.01052*v)*sigmoid(v,-40.14,-1/0.1378) : 0.6*exp(0.057*v)*sigmoid(v,-32,-10);
  CepsReal a39 = 1/(a18+a31);
  dtyLin [_j] = -1/a39;
  dtyNLin[_j] = a5/a39;
 
  CepsReal a6  = sigmoid(v,-8,-7.5);
  CepsReal a19 = 1.4*sigmoid(v,-35,-13) + 0.25;
  CepsReal a32 = 1.4*sigmoid(v,-5,5);
  CepsReal a40 = sigmoid(v,50,-20);
  CepsReal a42 = a19*a32+a40;
  dtyLin [_d] = -1/a42;
  dtyNLin[_d] = a6/a42;
 
  CepsReal texp = exp(-v/RTF);
  CepsReal a55  = c[_pnak]*c[_ko]/(c[_ko]+c[_kmk]) *y[_nai]/(y[_nai]+c[_kmna])
                / (1 +0.1245*pow(texp,0.1) +0.0353*texp);
  CepsReal a25  = RTF*std::log(c[_nao]/y[_nai]);
  CepsReal a50  = c[_gna]*y[_m]*y[_m]*y[_m]*y[_h]*y[_j]*(v-a25);
  CepsReal a51  = c[_gbna]*(v-a25);
  CepsReal a56a = c[_cao]*pow(texp,-c[_gamma])*pow(y[_nai],3); 
  CepsReal a56b = c[_alpha]*y[_cai]*pow(texp,1-c[_gamma])*pow(c[_nao],3);
  CepsReal a56c = (pow(c[_kmnai],3)+pow(c[_nao],3))*(c[_kmca]+c[_cao]);
  CepsReal a56d = 1+c[_ksat]*pow(texp,1-c[_gamma]);
  CepsReal a56  = c[_knaca]*(a56a-a56b)/a56c/a56d;
  dtyLin [_nai] = 0.;
  dtyNLin[_nai] = -cm*(a50+a51+3*a55+3*a56)/c[_vc]/c[_F];
 
  CepsReal a33  = RTF*std::log(c[_ko]/y[_ki]);
  CepsReal a44  = 0.1*sigmoid(v,a33+200,1/0.06);
  CepsReal a45  = (3*exp(2e-4*(v-a33+100)) +exp(0.1*(v-a33-10)))*sigmoid(v,a33,-2);
  CepsReal a46  = a44/(a44+a45);
  CepsReal a47  = c[_gk1]*a46*std::sqrt(c[_ko]/5.4)*(v-a33);
  CepsReal a54  = c[_gto]*y[_r]*y[_s]*(v-a33);
  CepsReal a48  = c[_gkr]*std::sqrt(c[_ko]/5.4)*y[_xr1]*y[_xr2]*(v-a33);
 
  CepsReal a41  = RTF*std::log( (c[_ko]+c[_pkna]*c[_nao])/(y[_ki]+c[_pkna]*y[_nai]) );
  CepsReal a49  = c[_gks]*y[_xs]*y[_xs]*(v-a41);
  //CepsReal a52a = c[_gcal]*y[_d]*y[_f]*y[_f2]*y[_fcass]*4*(v-15)*c[_F]/RTF;
  //CepsReal a52b = 0.25*y[_cass]*exp(2*(v-15)/RTF)-c[_cao];
  //CepsReal a52c = expm1(2*(v-15)/RTF);
  //CepsReal a52  = a52a*a52b/a52c;
  CepsReal a52a = c[_gcal]*y[_d]*y[_f]*y[_f2]*y[_fcass]*2*c[_F]*uOverExpm1u(2*(v-15)/RTF);
  CepsReal a52b = 0.25*y[_cass]*exp(2*(v-15)/RTF)-c[_cao];
  CepsReal a52  = a52a*a52b;
 
  CepsReal a43  = 0.5*RTF*std::log(c[_cao]/y[_cai]);
  CepsReal a53  = c[_gbca]*(v-a43);
  CepsReal a58  = c[_gpk]*(v-a33)*sigmoid(v,25,-5.98);
  CepsReal a57  = c[_gpca]*y[_cai]/(y[_cai]+c[_kpca]);
 
  CepsReal iStim = getStimulation(dof,t)/cm;
  CepsReal iIon  = -(a47+a54+a48+a49+a52+a55+a50+a51+a56+a53+a58+a57);
  #ifdef DEBUG
    ceps::checkNanOrInf(iIon,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
 
  *dtv = iIon + iStim; 
 
  dtyLin [_ki] = 0.;
  dtyNLin[_ki] = -cm*(a47+a54+a48+a49+a58+iStim-2*a55)/(c[_vc]*c[_F]);
 
  CepsReal a59 = c[_vmaxup]/(1+c[_kup]*c[_kup]/y[_cai]/y[_cai]);
  CepsReal a60 = c[_vleak]*(y[_casr]-y[_cai]);
  CepsReal a61 = c[_vxfer]*(y[_cass]-y[_cai]);
  CepsReal a63 = 1/(1+c[_bufc]*c[_kbufc]/pow(y[_cai]+c[_kbufc],2));
  dtyLin [_cai] = 0.;
  dtyNLin[_cai] = a63*((a61+(a60-a59)*c[_vsr]/c[_vc]) - cm*(a53+a57-2*a56)/(2*c[_vc]*c[_F]));
 
  CepsReal a62 = c[_maxsr]-(c[_maxsr]-c[_minsr])/(1+pow(c[_ec]/y[_casr],2));
  CepsReal a65 = c[_k2prime]*a62;
  dtyLin [_rprime] = -c[_k4];
  dtyNLin[_rprime] = c[_k4]-a65*y[_cass]*y[_rprime];
 
  CepsReal a64 = c[_k1prime]/a62;
  CepsReal a66 = a64*y[_cass]*y[_cass]*y[_rprime]/(c[_k3]+a64*y[_cass]*y[_cass]);
  CepsReal a67 = c[_vrel]*a66*(y[_casr]-y[_cass]);
  CepsReal a68 = 1/(1+c[_bufsr]*c[_kbuf]/pow(y[_casr]+c[_kbuf],2));
  dtyLin [_casr] = 0.;
  dtyNLin[_casr] = a68*(a59-a67-a60);
 
  CepsReal a69 = 1/(1+c[_bufss]*c[_kbufss]/pow(y[_cass]+c[_kbufss],2));
  dtyLin [_cass] = 0.;
  dtyNLin[_cass] = a69*( (-a52*cm/(2*c[_vss]*c[_F]) +a67*c[_vsr]/c[_vss]) -a61*c[_vc]/c[_vss]);
 
}
 
void
TTP06::setupWithParameters(InputParameters* p, FunctionDictionary* dico)
{
 
  CepsString key = "TTP06_";
  if (m_type==Type::Endo)
    key += "endo ";
  if (m_type==Type::Epi)
    key += "epi ";
  if (m_type==Type::MidMyo)
    key += "midmyo ";
 
  m_constants[_R      ] = p->getReal(key+"R      ",m_constants[_R      ]);
  m_constants[_T      ] = p->getReal(key+"T      ",m_constants[_T      ]);
  m_constants[_F      ] = p->getReal(key+"F      ",m_constants[_F      ]);
  m_constants[_vc     ] = p->getReal(key+"vc     ",m_constants[_vc     ]);
  m_constants[_pkna   ] = p->getReal(key+"pkna   ",m_constants[_pkna   ]);
  m_constants[_ko     ] = p->getReal(key+"ko     ",m_constants[_ko     ]);
  m_constants[_nao    ] = p->getReal(key+"nao    ",m_constants[_nao    ]);
  m_constants[_cao    ] = p->getReal(key+"cao    ",m_constants[_cao    ]);
  m_constants[_gk1    ] = p->getReal(key+"gk1    ",m_constants[_gk1    ]);
  m_constants[_gkr    ] = p->getReal(key+"gkr    ",m_constants[_gkr    ]);
  m_constants[_gks    ] = p->getReal(key+"gks    ",m_constants[_gks    ]);
  m_constants[_gna    ] = p->getReal(key+"gna    ",m_constants[_gna    ]);
  m_constants[_gbna   ] = p->getReal(key+"gbna   ",m_constants[_gbna   ]);
  m_constants[_gcal   ] = p->getReal(key+"gcal   ",m_constants[_gcal   ]);
  m_constants[_gbca   ] = p->getReal(key+"gbca   ",m_constants[_gbca   ]);
  m_constants[_gto    ] = p->getReal(key+"gto    ",m_constants[_gto    ]);
  m_constants[_pnak   ] = p->getReal(key+"pnak   ",m_constants[_pnak   ]);
  m_constants[_kmk    ] = p->getReal(key+"kmk    ",m_constants[_kmk    ]);
  m_constants[_kmna   ] = p->getReal(key+"kmna   ",m_constants[_kmna   ]);
  m_constants[_knaca  ] = p->getReal(key+"knaca  ",m_constants[_knaca  ]);
  m_constants[_ksat   ] = p->getReal(key+"ksat   ",m_constants[_ksat   ]);
  m_constants[_alpha  ] = p->getReal(key+"alpha  ",m_constants[_alpha  ]);
  m_constants[_gamma  ] = p->getReal(key+"gamma  ",m_constants[_gamma  ]);
  m_constants[_kmca   ] = p->getReal(key+"kmca   ",m_constants[_kmca   ]);
  m_constants[_kmna   ] = p->getReal(key+"kmna   ",m_constants[_kmna   ]);
  m_constants[_gpca   ] = p->getReal(key+"gpca   ",m_constants[_gpca   ]);
  m_constants[_kpca   ] = p->getReal(key+"kpca   ",m_constants[_kpca   ]);
  m_constants[_gpk    ] = p->getReal(key+"gpk    ",m_constants[_gpk    ]);
  m_constants[_k1prime] = p->getReal(key+"k1prime",m_constants[_k1prime]);
  m_constants[_k2prime] = p->getReal(key+"k2prime",m_constants[_k2prime]);
  m_constants[_k3     ] = p->getReal(key+"k3     ",m_constants[_k3     ]);
  m_constants[_k4     ] = p->getReal(key+"k4     ",m_constants[_k4     ]);
  m_constants[_ec     ] = p->getReal(key+"ec     ",m_constants[_ec     ]);
  m_constants[_maxsr  ] = p->getReal(key+"maxsr  ",m_constants[_maxsr  ]);
  m_constants[_minsr  ] = p->getReal(key+"minsr  ",m_constants[_minsr  ]);
  m_constants[_vrel   ] = p->getReal(key+"vrel   ",m_constants[_vrel   ]);
  m_constants[_vxfer  ] = p->getReal(key+"vxfer  ",m_constants[_vxfer  ]);
  m_constants[_kup    ] = p->getReal(key+"kup    ",m_constants[_kup    ]);
  m_constants[_vleak  ] = p->getReal(key+"vleak  ",m_constants[_vleak  ]);
  m_constants[_vmaxup ] = p->getReal(key+"vmaxup ",m_constants[_vmaxup ]);
  m_constants[_bufc   ] = p->getReal(key+"bufc   ",m_constants[_bufc   ]);
  m_constants[_kbufc  ] = p->getReal(key+"kbufc  ",m_constants[_kbufc  ]);
  m_constants[_bufsr  ] = p->getReal(key+"bufsr  ",m_constants[_bufsr  ]);
  m_constants[_kbuf   ] = p->getReal(key+"kbuf   ",m_constants[_kbuf   ]);
  m_constants[_bufss  ] = p->getReal(key+"bufss  ",m_constants[_bufss  ]);
  m_constants[_kbufss ] = p->getReal(key+"kbufss ",m_constants[_kbufss ]);
  m_constants[_vsr    ] = p->getReal(key+"vsr    ",m_constants[_vsr    ]);
  m_constants[_vss    ] = p->getReal(key+"vss    ",m_constants[_vss    ]);
}