BR77.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/BR77.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
 
BR77::BR77(
  Unknown*                      u, 
  AbstractPdeProblem*           pb, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  AbstractIonicModel(u,attrs)
{
  
  CEPS_ABORT_IF(ceps::isNullPtr(pb),
    "Passed pointer to cardiac problem is null" 
  );
 
  m_name          = "Beeler-Reuter 1977";
  m_nStateVars    = 7;
  m_stateVarNames = {"m [adim]",
                     "h [adim]",
                     "j [adim]",
                     "CaI [concentration units]",
                     "d [adim]",
                     "f [adim]",
                     "x1 [adim]"};
 
  // ================================================================================
  // Default parameters
  m_constants = ceps::newArray<CepsReal>(2);
  m_paperCm         = 1.e-2;
  m_paperStim       = 0.5;
  m_constants[_gS ] = 9.e-4;
  m_constants[_ENa] = 5.e+1;
 
  // Create entries in the function dictionary
  FunctionDictionary* dico = pb->getFunctionDictionary();
  dico->add("br77gNa"  ,"CONSTANT 4.e-2");  
  dico->add("br77gNaCa","CONSTANT 3.e-5");  
 
  // ================================================================================
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params,dico,pb);
 
  // ================================================================================
  // Initialize scalar fields for gNa and gNaCa
  AbstractDiscretization*             sd   = pb->getSpatialDiscretization();
  const CepsVector<DegreeOfFreedom*>& dofs = sd->getDegreesOfFreedomForUnknown(m_unknown);
 
  m_gNa   = ceps::getNew<ScalarField<DegreeOfFreedom>>(dico->getScalar("br77gNa"  ),&dofs,true);
  m_gNaCa = ceps::getNew<ScalarField<DegreeOfFreedom>>(dico->getScalar("br77gNaCa"),&dofs,true);
 
  // Restrict support if the model is defined on selected attributes
  if (not m_attributes.empty())
  {
    AttributesSelector<CepsVector<DegreeOfFreedom*>::const_iterator> selec;
    selec.setAttributes(m_attributes.begin(),m_attributes.end());
    m_gNa  ->computeSupport(&selec);
    m_gNaCa->computeSupport(&selec);
  }
  
}
 
void 
BR77::getInitialCondition(CepsReal* v, CepsReal* y) const
{
  *v = -84.624; 
  y[_m  ] = 0.011;  
  y[_h  ] = 0.988;  
  y[_j  ] = 0.975;  
  y[_CaI] = 1.e-4; 
  y[_d  ] = 3.e-3; 
  y[_f  ] = 0.994;
  y[_x1 ] = 1.e-4; 
}
 
void
BR77::computeRates(
  CepsReal         t,
  CepsReal*        y, 
  CepsReal*        vm, 
  CepsReal*        dtyLin, 
  CepsReal*        dtyNLin, 
  CepsReal*        dtv,
  DegreeOfFreedom* dof
) const
{ 
 
  CepsReal* cst = m_constants;
  CepsReal  v   = *vm;
 
  // Gate variables 
  CepsReal alpha;
  CepsReal beta;
 
  alpha = -(v+47.)/std::expm1(-0.1*(v+47.));
  beta  = 40*std::exp(-0.056*(v+72.));
  dtyLin [_m] = -alpha-beta;
  dtyNLin[_m] = alpha;
 
  alpha = 0.126*std::exp(-0.25*(v+77.));
  beta  = 1.7*sigmoid(v,-22.5,-1./0.082); 
  dtyLin [_h] = -alpha-beta;
  dtyNLin[_h] = alpha;
 
  // The sigmoids version refer to the original paper and can also be found in Myokit.
  // In comments, the cellML version, with divisions
 
  alpha = 0.055*std::exp(-0.25*(v+78))*sigmoid(v,-78,-1./0.2); // /(1+std::exp(-0.2*(v+78))); 
  beta  = 0.3                         *sigmoid(v,-32,-1./0.1); // /(1+std::exp(-0.1*(v+32)));
  dtyLin [_j] = -alpha-beta;
  dtyNLin[_j] = alpha;
 
  alpha = 0.095*std::exp(-(v-5. )/100.)*sigmoid(v,+ 5,-1./0.072); // /(1+std::exp(-(v-5. )/13.89));
  beta  = 0.07 *std::exp(-(v+44.)/ 59.)*sigmoid(v,-44,+1./0.05 ); // /(1+std::exp( (v+44.)/20.  ));
  dtyLin [_d] = -alpha-beta;
  dtyNLin[_d] = alpha;
 
  alpha = 0.012 *std::exp(-(v+28.)/125.)*sigmoid(v,-28,+1./0.15); // /(1+std::exp( (v+28.)/6.67));
  beta  = 0.0065*std::exp(-(v+30.)/50. )*sigmoid(v,-30,-1./0.2 ); // /(1+std::exp(-(v+30.)/5.  ));
  dtyLin [_f] = -alpha-beta;
  dtyNLin[_f] = alpha;
 
  alpha = 0.0005*std::exp( (v+50.)/12.10)*sigmoid(v,-50, 1./0.057); // /(1+std::exp( (v+50.)/17.5));
  beta  = 0.0013*std::exp(-(v+20.)/16.67)*sigmoid(v,-20,-1./0.04 ); // /(1+std::exp(-(v+20.)/25.0));
  dtyLin [_x1] = -alpha-beta;
  dtyNLin[_x1] = alpha;
 
  // Last state variable caI
  CepsReal Es  = -82.3-13.0287*std::log(y[_CaI]*0.001);
  CepsReal iS  = cst[_gS]*y[_d]*y[_f]*(v-Es);
  dtyNLin[_CaI] = 7.e-6 -0.01*iS;
  dtyLin [_CaI] = -0.07;
  
  // Fetch Ionic current
  CepsReal gNa   = m_gNa  ->getValue(dof,t);
  CepsReal gNaCa = m_gNaCa->getValue(dof,t);
 
  CepsReal iNa = this->getINa(y,gNaCa,gNa,v);
  CepsReal iX1 = y[_x1] *0.008*std::expm1(0.04*(v+77)) / std::exp(0.04*(v+35.));
  CepsReal tmp = std::exp(0.04*(v+53.));
  CepsReal iK1 = 0.0035*(4.*std::expm1(0.04*(v+85.))/(tmp*tmp+tmp) 
                          + 0.2/0.04*uOverExpm1u(-0.04*(v+23.)));
 
  CepsReal cm    = getCmInternal(t,dof);
  CepsReal iStim = getStimulation(dof,t);
  CepsReal iion  = -(iNa+iS+iK1+iX1);
  #ifdef DEBUG
    ceps::checkNanOrInf(iion,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
 
  *dtv= iion + iStim;
  *dtv /= cm;
 
}
 
CepsReal
BR77::convertCurrentFromCepsUnit(const CepsReal& i) const
{
  return i/100.;
}
 
CepsReal
BR77::convertCmToCepsUnit(const CepsReal& cm) const
{
  return 100*cm;
}
 
CepsReal
BR77::convertCmFromCepsUnit(const CepsReal& cm) const
{
  return 0.01*cm;
}
 
void
BR77::setupWithParameters(InputParameters* params, FunctionDictionary* dico, AbstractPdeProblem* pb)
{
  m_constants[_gS ] = params->getReal("BR77 gS" ,m_constants[_gS ]);
  m_constants[_ENa] = params->getReal("BR77 ENa",m_constants[_ENa]);
 
  for (CepsString key: {"BR77gNa","BR77gNaCa"})
    if (params->hasKey(key))
    {
      dico->deleteScalar(key);
      dico->add(key,params->getString(key),pb->getGeometry());
    }
 
}
 
CepsReal
BR77::getINa(CepsReal* y, CepsReal gNaCa, CepsReal gNa, CepsReal v) const
{
  return (gNa*std::pow(y[_m],3.)*y[_h]*y[_j]+gNaCa)*(v-m_constants[_ENa]);
}
 
// ==========================================================================
 
BR77Modified::BR77Modified(
  Unknown*                      u, 
  AbstractPdeProblem*           pb, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  BR77(u,pb,attrs,params)
{
 
  m_name       = "Beeler-Reuter 1977 Modified";
 
  // ================================================================================
  // Default parameters
  delete[] m_constants;
  m_constants = ceps::newArray<CepsReal>(4);
  m_constants[_gS ] = 9.e-4;
  m_constants[_ENa] = 5.e+1;
  m_constants[_vA ] = -85*5;
  m_constants[_vB ] = 30+5*115;
 
  // ================================================================================
  // Parameters from text inputs, if any
  FunctionDictionary* dico = pb->getFunctionDictionary();
  if (ceps::isValidPtr(params))
    setupWithParameters(params,dico,pb);
 
}
 
void
BR77Modified::setupWithParameters(InputParameters* params, FunctionDictionary* dico, AbstractPdeProblem* pb)
{
  BR77::setupWithParameters(params,dico,pb);
  m_constants[_vA] = params->getReal("BR77 a" ,m_constants[_vA]);
  m_constants[_vB] = params->getReal("BR77 b" ,m_constants[_vB]);
}
 
CepsReal
BR77Modified::getINa(CepsReal* y, CepsReal gNaCa, CepsReal gNa, CepsReal v) const
{
 
  CepsReal iNa;
  CepsReal a = m_constants[_vA];
  CepsReal b = m_constants[_vB];
  
  CepsReal iNaA = BR77::getINa(y,gNaCa,gNa,a);
  CepsReal iNaB = BR77::getINa(y,gNaCa,gNa,b);
 
  if (ceps::greaterThan(v,b) and not ceps::equals(iNaB,0.))
  {
    CepsReal kB = -1/iNaB*(gNa*std::pow(y[_m],3.)*y[_h]*y[_j]+gNaCa);
    iNa = iNaB*(1+std::tanh(kB*(v-b)));
  }
  else if (ceps::lesserThan(v,a) and not ceps::equals(iNaA,0.))
  {
    CepsReal kA = 1/iNaA*(gNa*std::pow(y[_m],3.)*y[_h]*y[_j]+gNaCa);
    iNa = iNaA*(1-std::tanh(kA*(v-a)));
  }
  else
    iNa = BR77::getINa(y,gNaCa,gNa,v);
 
  return iNa; 
}