MS03.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/MS03.hpp"
 
MS03::MS03(
  Unknown*                      u, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  AbstractIonicModel(u,attrs)
{
 
  m_name          = "Mitchell-Schaeffer 2003";
  m_nStateVars    = 1;
  m_stateVarNames = {"h [adim]"};
 
  // ================================================================================
  // Default parameters
 
  m_vMin        = -90; 
  m_vMax        =  50;
  m_cellSurface =  1.;
  
  m_constants = new CepsReal[5];  
  m_constants[_tauIn   ] = 0.3;
  m_constants[_tauOut  ] = 6;
  m_constants[_tauOpen ] = 130;
  m_constants[_tauClose] = 150;
  m_constants[_vGate   ] = 0.13;
 
  m_vInit = 0.00000820413566106744;
  m_hInit = 0.8789655121804799; 
 
  m_paperCm   = 1.;
  m_paperStim = 0.2;
 
  // ================================================================================
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params,nullptr);
  
}
 
void 
MS03::getInitialCondition(CepsReal* v, CepsReal* y) const
{
  // *v = 0.0;
  // y[_h] = 1.0;
  *v    = m_vInit;
  y[_h] = m_hInit;
}
 
void
MS03::computeRates(
  CepsReal         t,
  CepsReal*        y, 
  CepsReal*        vm, 
  CepsReal*        dtyLin, 
  CepsReal*        dtyNLin, 
  CepsReal*        dtv, 
  DegreeOfFreedom* dof
) const
{ 
 
  CepsReal* c = m_constants;
  CepsReal  v = *vm;
 
  CepsReal tau  = v<c[_vGate] ? c[_tauOpen] : c[_tauClose];
  CepsReal wInf = v<c[_vGate] ? 1:0;
  dtyLin [_h] = -1/tau;
  dtyNLin[_h] = wInf/tau;
 
  CepsReal cm    = getCmInternal(t,dof);
  CepsReal iStim = getStimulation(dof,t);
  CepsReal iion  = y[_h]*(1-v)*v*v/c[_tauIn] - v/c[_tauOut];
  #ifdef DEBUG
    ceps::checkNanOrInf(iion,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
 
  *dtv  = iion + iStim;
  *dtv /= cm;
 
}
 
CepsReal
MS03::convertVoltageToCepsUnit(const CepsReal& v) const
{
  return m_vMin + v*(m_vMax-m_vMin);
}
 
CepsReal
MS03::convertVoltageFromCepsUnit(const CepsReal& v) const
{
  return (v-m_vMin)/(m_vMax-m_vMin);
}
 
CepsReal
MS03::convertCurrentFromCepsUnit(const CepsReal& i) const
{
  return i/(m_vMax-m_vMin);
}
 
CepsReal
MS03::convertDtvToCepsUnit(const CepsReal& dtv) const
{
  return dtv*(m_vMax-m_vMin);
}
 
CepsReal
MS03::convertCmFromCepsUnit(const CepsReal& cm) const
{
  return cm;
}
 
CepsReal
MS03::convertCmToCepsUnit(const CepsReal& cm) const
{
  return cm;
}
 
void
MS03::setupWithParameters(InputParameters* params, FunctionDictionary* dico)
{
 
  m_constants[_tauIn   ] = params->getReal("MS03 tauIn   ",m_constants[_tauIn   ]); 
  m_constants[_tauOut  ] = params->getReal("MS03 tauOut  ",m_constants[_tauOut  ]); 
  m_constants[_tauOpen ] = params->getReal("MS03 tauOpen ",m_constants[_tauOpen ]); 
  m_constants[_tauClose] = params->getReal("MS03 tauClose",m_constants[_tauClose]); 
  m_constants[_vGate   ] = params->getReal("MS03 vGate   ",m_constants[_vGate   ]); 
 
 
  m_vMin        = params->getReal("MS03 vMin"        ,m_vMin       );
  m_vMax        = params->getReal("MS03 vMax"        ,m_vMax       );
  m_vInit       = params->getReal("MS03 v init"      ,m_vInit      );
  m_hInit       = params->getReal("MS03 h init"      ,m_hInit      );
  m_cellSurface = params->getReal("MS03 cell surface",m_cellSurface);
 
  return;
}
 
 
// =================================================================================================
// Modded MS03 (Regularized+bounded currents)
 
MS03Modified::MS03Modified(
  Unknown*                      u, 
  const CepsSet<CepsAttribute>& attrs,
  InputParameters*              params
):
  MS03(u,attrs,params)
{
 
  m_name       = "Mitchell-Schaeffer 2003";
  m_nStateVars = 1;
  m_stateVarNames = {"h [adim]"};
 
  // ================================================================================
  // Default parameters
  delete[] m_constants; 
  m_constants = new CepsReal[8];  
  m_constants[_tauIn   ] = 0.3;
  m_constants[_tauOut  ] = 6;
  m_constants[_tauOpen ] = 130;
  m_constants[_tauClose] = 150;
  m_constants[_vGate   ] = 0.13;
  m_constants[_etaGate ] = 0.02;
  m_constants[_vB      ] = 5;
  m_constants[_vA      ] = -5;
 
  // ================================================================================
  // Parameters from text inputs, if any
  if (ceps::isValidPtr(params))
    setupWithParameters(params,nullptr);
  
}
 
void
MS03Modified::computeRates(
  CepsReal         t,
  CepsReal*        y, 
  CepsReal*        vm, 
  CepsReal*        dtyLin, 
  CepsReal*        dtyNLin, 
  CepsReal*        dtv, 
  DegreeOfFreedom* dof
) const
{ 
 
  CepsReal* c = m_constants;
  CepsReal  v = *vm;
  CepsReal  a = m_constants[_vA];
  CepsReal  b = m_constants[_vB];
 
  CepsReal hinfty = 0.5*(1-std::tanh((v-c[_vGate])/c[_etaGate]));
  CepsReal coeff  = 1/c[_tauClose] + hinfty*(c[_tauClose]-c[_tauOpen])/(c[_tauClose]*c[_tauOpen]);
 
  dtyLin [_h] = -coeff;
  dtyNLin[_h] = coeff*hinfty;
 
  CepsReal cm    = getCmInternal(t,dof);
  CepsReal iStim = getStimulation(dof,t);
 
  CepsReal iIonB = y[_h]*(1-b)*b*b/c[_tauIn] - b/c[_tauOut];
  CepsReal iIonA = y[_h]*(1-a)*a*a/c[_tauIn] - a/c[_tauOut];
 
  // Current limiter
  CepsReal iIon;
  if (ceps::greaterThan(v,b) and not ceps::equals(iIonB,0.))
  {
    CepsReal kB = -1/iIonB*(y[_h]/c[_tauIn]*b*(3*b-2)+1/c[_tauOut]);
    iIon = iIonB*(1+std::tanh(kB*(v-b)));
  }
  else if (ceps::lesserThan(v,a) and not ceps::equals(iIonA,0.))
  {
    CepsReal kA = 1/iIonA*(y[_h]/c[_tauIn]*a*(3*a-2)+1/c[_tauOut]);
    iIon = iIonA*(1-std::tanh(kA*(v-a)));
  }
  else
    iIon = y[_h]*(1-v)*v*v/c[_tauIn] - v/c[_tauOut];
 
  #ifdef DEBUG
    ceps::checkNanOrInf(iIon,
      "Ionic current on DoF " + CepsString(dof->getGlobalIndex() + " at time" + CepsString(t))
    );
  #endif
 
  *dtv  = iIon + iStim;
  *dtv /= cm;
 
}
 
 
void
MS03Modified::setupWithParameters(InputParameters* params, FunctionDictionary* dico)
{
 
  MS03::setupWithParameters(params,dico);
  m_constants[_etaGate] = params->getReal("MS03 etaGate",m_constants[_etaGate]); 
  m_constants[_vA     ] = params->getReal("MS03 a"      ,m_constants[_vA     ]); 
  m_constants[_vB     ] = params->getReal("MS03 b"      ,m_constants[_vB     ]); 
 
  return;
}