IonicModel Class Reference

Abstract ionic model class. More...

#include <ionic_model.h>

Inheritance diagram for IonicModel:

Public Types
using	InitialStates = std::vector< std::pair< std::string, double > >

Public Member Functions
	IonicModel (const InitialStates &initial_X_, const InitialStates &initial_Xg_, const double Vrest_)
	Constructor.
	IonicModel (const InitialStates &initial_X_, const InitialStates &initial_Xg_, const double Vrest_, const double Vscale_, const double Tscale_, const double Voffset_)
	Constructor with scaling factors.
virtual	~IonicModel ()=default
	Virtual destructor.
virtual std::unique_ptr< IonicModelParameters >	get_parameters () const
	Construct an instance of model parameters for this model.
virtual void	read_parameters (const IonicModelParameters &params)
	Read model parameters from a parameter object.
virtual void	distribute_parameters (const CmMod &cm_mod, const cmType &cm)
	Distribute model parameters to all parallel processes.
void	init (Vector< double > &X, Vector< double > &Xg) const
	Setup model initial conditions.
void	integ (const odeType &ode_solver_params, const int zone_id, const double t, const double dt, const double Istim, const double Ksac, Vector< double > &X, Vector< double > &Xg) const
	Integrate over one time step.
unsigned int	nX () const
	Get the number of state variables.
unsigned int	nG () const
	Get the number of gating variables.
virtual unsigned int	get_calcium_index () const =0
	Get the index of the intracellular calcium concentration in the state vector.
virtual std::vector< std::pair< std::string, int > >	get_output_variables () const
	Get a list of state variables to export to VTU.
std::vector< outputType >	get_registered_outputs () const
	Get output variable information for output registration.

Protected Member Functions
virtual void	update_g (const unsigned int zone_id, const double dt, const Vector< double > &X, Vector< double > &Xg) const =0
	Update variables with analytical solution.
virtual Vector< double >	getf (const unsigned int zone_id, const Vector< double > &X, const Vector< double > &Xg, const double I_stim, const double I_sac) const =0
	Model right hand side.
virtual Array< double >	getj (const unsigned int zone_id, const Vector< double > &X, const Vector< double > &Xg, const double Ksac) const
	Model jacobian.
Integration methods.
void	integ_fe (const unsigned int zone_id, Vector< double > &X, Vector< double > &Xg, const double Ts, const double Ti, const double Istim, const double Ksac) const
	Integrate the model with the forward Euler method.
void	integ_rk (const unsigned int zone_id, Vector< double > &X, Vector< double > &Xg, const double Ts, const double Ti, const double Istim, const double Ksac) const
	Integrate the model with the fourth-order explicit Runge-Kutta method.
void	integ_cn2 (const unsigned int zone_id, Vector< double > &X, Vector< double > &Xg, const double Ts, const double Ti, const double Istim, const double Ksac, const unsigned int max_iter, const double rtol, const double atol) const
	Integrate the model with the Crank-Nicolson method.

Protected Attributes
InitialStates	initial_X
	Initial states.
InitialStates	initial_Xg
	Initial gating variables.
const double	Vrest
Scaling factors.
Individual ionic models may need to rescale the time or voltage variable, e.g. to bring them into dimensionless form. These are the factors used for that purpose. They are assigned in the constructor of this class.
const double	Vscale
	Voltage scaling [mV].
const double	Tscale
	Time scaling [ms].
const double	Voffset
	Voltage offset parameter [mV].

Detailed Description

Abstract ionic model class.

Mathematical model

This class represents an abstract ionic model, i.e. an ODE system in the form

\[ \begin{aligned} \frac{\partial v}{\partial t} + I_\text{ion}(v, \mathbf{w}) &= I_\text{ext}(t) \\ \frac{\partial \mathbf{w}}{\partial t} &= \mathbf{f}(v, \mathbf{w}) \end{aligned} \]

where \(v\) is the transmembrane potential, \(\mathbf{w}\) is a vector of ionic model variables (which may be gating variables or ionic concentrations), \(I_\text{ion}\) is the ionic current, and \(I_\text{ext}\) is an externally applied current.

Individual models differ in the number of state variables \(\mathbf{w}\) and in the expressions of \(I_\text{ion}\) and \(\mathbf{F}\). These are specified by classes derived from this.

References:

• Colli Franzone, Pavarino, Scacchi. Mathematical Cardiac Electrophysiology. Springer, 2014.
• Goktepe, Kuhl. Computational Modeling of Cardiac Electrophysiology: a Novel Finite Element Approach. International Journal for Numerical Methods in Biomedical Engineering, 2009.

Numerical methods

The vector of state variables \(\mathbf{w}\) is partitioned into two groups, the state variables \(\mathbf{x}\) (including the transmembrane potential \(v\)) and the gating variables \(\mathbf{x}_g\). The ODE system is advanced through a Rush-Larsen scheme, in which the gating variables are updated through an analytical expression and the state variables through a time-stepping scheme.

This class currently supports forward Euler (FE), fourth-order explicit Runge-Kutta (RK) and Crank-Nicolson (CN) for advancing the state variables. Refer to the documentation of the functions integ_fe, integ_rk and integ_cn2 for details on the individual methods. Beware that the CN method is only supported for ionic models that implement the getj method, and an exception is raised otherwise.

Implementing concrete ionic models

To implement a new ionic model, the following steps need to be taken:

1. Create a class derived from IonicModel.
2. Override the methods init, update_g, getf for that ionic model. If implicit solvers are desired, also override getj.
3. Create a class derived from IonicModelParameters to manage the parameters specific to the new ionic model.
4. Override the methods get_parameters, read_parameters and distribute_parameters to manage the parameters of the new ionic model.
5. Override the get_calcium_index method to return the index of the calcium proxy variable.
6. Register the new class into the ionic model factory by using the macro REGISTER_IONIC_MODEL. The macro should be called in a .cpp file, not in a header file.
7. Edit the files CepMod.h and CepMod.cpp to add the label for the new ionic model type.

Member Typedef Documentation

InitialStates

using IonicModel::InitialStates = std::vector<std::pair<std::string, double> >

Alias for initial states vector. Each initial state is a pair of a label for that state variable and its initial value.

Todo:

[michelebucelli] This would work better with a struct, due to the fields having meaningful names instead of first and second.

Constructor & Destructor Documentation

IonicModel() [1/2]

IonicModel::IonicModel ( const InitialStates &  initial_X_, const InitialStates &  initial_Xg_, const double  Vrest_  )

inline

Constructor.

IonicModel() [2/2]

IonicModel::IonicModel ( const InitialStates &  initial_X_, const InitialStates &  initial_Xg_, const double  Vrest_, const double  Vscale_, const double  Tscale_, const double  Voffset_  )

inline

Constructor with scaling factors.

~IonicModel()

virtual IonicModel::~IonicModel ( )

virtualdefault

Virtual destructor.

Member Function Documentation

distribute_parameters()

void IonicModel::distribute_parameters ( const CmMod &  cm_mod, const cmType &  cm  )

virtual

Distribute model parameters to all parallel processes.

By default, this method only takes care of the parameters related to initial states. If a derived ionic model has other parameters, then this method should be overridden to distribute those as well.

Reimplemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

get_calcium_index()

virtual unsigned int IonicModel::get_calcium_index ( ) const

pure virtual

Get the index of the intracellular calcium concentration in the state vector.

This is the index of the variable to be used for electromechanics coupling.

Implemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

get_output_variables()

virtual std::vector< std::pair< std::string, int > > IonicModel::get_output_variables ( ) const

inlinevirtual

Get a list of state variables to export to VTU.

By default, this function returns the calcium index as the only exported state variable. Derived classes can override this to export additional states if needed. Beware that, for phenomenological models, the calcium variable might actually be a proxy for calcium, rather than the actual concentration (and in particular it might be non-dimensional or have different units than a molar concentration).

Returns

A vector of pairs {variable_name, state_vector_index}. The vector need not include the transmembrane potential V, which is exported by other means.

get_parameters()

virtual std::unique_ptr< IonicModelParameters > IonicModel::get_parameters ( ) const

inlinevirtual

Construct an instance of model parameters for this model.

Reimplemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

get_registered_outputs()

std::vector< outputType > IonicModel::get_registered_outputs

(

)

const

Get output variable information for output registration.

getf()

virtual Vector< double > IonicModel::getf ( const unsigned int  zone_id, const Vector< double > &  X, const Vector< double > &  Xg, const double  I_stim, const double  I_sac  ) const

protectedpure virtual

Model right hand side.

Defines the right-hand side function for the potential and ionic equations. Must be ovverridden in derived classes.

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in]	X	Vector of state variables.
[in]	Xg	Vector of gating variables.
[in]	I_stim	Applied current.
[in]	I_sac	Stretch-activated current.

Returns

A vector containing the right-hand side of the model equations.

Implemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

getj()

virtual Array< double > IonicModel::getj ( const unsigned int  zone_id, const Vector< double > &  X, const Vector< double > &  Xg, const double  Ksac  ) const

inlineprotectedvirtual

Model jacobian.

Defines the jacobian matrix of the model equations, that is the matirx of derivatives of the function evaulated by getf.

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in]	X	Vector of state variables.
[in]	Xg	Vector of gating variables.
[in]	Ksac	Stretch-activated current coefficient.

Returns

A matrix containing the jacobian of the model equations.

Reimplemented in AlievPanfilov, BuenoOrovio, and FitzHughNagumo.

init()

void IonicModel::init	(	Vector< double > &	X,
		Vector< double > &	Xg
	)		const

Setup model initial conditions.

Parameters

[out]	X	Vector of state variables to be initialized.
[out]	Xg	Vector of gating variables to be initialized.

integ()

void IonicModel::integ	(	const odeType &	ode_solver_params,
		const int	zone_id,
		const double	t,
		const double	dt,
		const double	Istim,
		const double	Ksac,
		Vector< double > &	X,
		Vector< double > &	Xg
	)		const

Integrate over one time step.

Parameters

[in]	ode_solver_params	ODE solver parameters structure, including the solver method and the stopping criterion information.
[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in]	t	Current time.
[in]	dt	Integration time step.
[in]	Istim	Amplitude of the stimulus current at the current time.
[in]	Ksac	Amplitude of the stretch-activated current.
[in,out]	X	Vector of state variables to be updated.
[in,out]	Xg	Vector of gating variables to be updated.

integ_cn2()

void IonicModel::integ_cn2 ( const unsigned int  zone_id, Vector< double > &  X, Vector< double > &  Xg, const double  Ts, const double  Ti, const double  Istim, const double  Ksac, const unsigned int  max_iter, const double  rtol, const double  atol  ) const

protected

Integrate the model with the Crank-Nicolson method.

The state variables \(\mathbf{x}\) are updated by solving the following non-linear problem:

\[ \frac{\mathbf{x}^{n+1} - \mathbf{x}^n}{\Delta t} = \frac{1}{2} \left( \mathbf{f}(\mathbf{x}^n, \mathbf{x}_g^n) + \mathbf{f}(\mathbf{x}^{n+1}, \mathbf{x}_g^n) \right) \]

The problem is solved through Newton's method. Subsequently, the gating variables are updated:

\[ \mathbf{x}_g^{n+1} = \texttt{update_g}(\Delta t, \mathbf{x}^{n+1}, \mathbf{x}_g^n) \]

Since this method involves solving a nonlinear system of equations, it is only available for those derived classes that implement the Jacobian matrix of the system through getj. An exception will be raised otherwise.

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in,out]	X	Vector of state variables to be updated.
[in,out]	Xg	Vector of gating variables to be updated.
[in]	Ts	Current time.
[in]	Ti	Time step.
[in]	Istim	Applied current.
[in]	Ksac	Stretch-activated current coefficient.
[in]	max_iter	Maximum number of Newton iterations. Beware that no error is raised if the maximum number of iterations is reached, so that the solution might be affected by the truncation of the iterations.
[in]	rtol	Relative tolerance for the Newton method.
[in]	atol	Absolute tolerance for the Newton method.

integ_fe()

void IonicModel::integ_fe ( const unsigned int  zone_id, Vector< double > &  X, Vector< double > &  Xg, const double  Ts, const double  Ti, const double  Istim, const double  Ksac  ) const

protected

Integrate the model with the forward Euler method.

\[ \begin{aligned} \mathbf{x}^{n+1} &= \mathbf{x}^n + \Delta t \mathbf{f}(\mathbf{x}^n, \mathbf{x}_g^n) \\ \mathbf{x}_g^{n+1} &= \texttt{update_g}(\Delta t, \mathbf{x}^n, \mathbf{x}_g^n) \end{aligned} \]

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in,out]	X	Vector of state variables to be updated.
[in,out]	Xg	Vector of gating variables to be updated.
[in]	Ts	Current time.
[in]	Ti	Time step.
[in]	Istim	Applied current.
[in]	Ksac	Stretch-activated current coefficient.

integ_rk()

void IonicModel::integ_rk ( const unsigned int  zone_id, Vector< double > &  X, Vector< double > &  Xg, const double  Ts, const double  Ti, const double  Istim, const double  Ksac  ) const

protected

Integrate the model with the fourth-order explicit Runge-Kutta method.

\[ \begin{aligned} \mathbf{x}^{(1)} &= \mathbf{x}^n \\ \mathbf{f}^{(1)} &= \mathbf{f}(\mathbf{x}^{(1)}, \mathbf{x}_g^n) \\ \mathbf{x}_g^{(1)} &= \texttt{update_g}(\Delta t/2, \mathbf{x}^{(1)}, \mathbf{x}_g^n) \\ \\ \mathbf{x}^{(2)} &= \mathbf{x}^n + \frac{\Delta t}{2}\mathbf{f}^{(1)} \\ \mathbf{f}^{(2)} &= \mathbf{f}(\mathbf{x}^{(2)}, \mathbf{x}_g^{(1)}) \\ \\ \mathbf{x}^{(3)} &= \mathbf{x}^n + \frac{\Delta t}{2}\mathbf{f}^{(2)} \\ \mathbf{f}^{(3)} &= \mathbf{f}(\mathbf{x}^{(3)}, \mathbf{x}_g^{(1)}) \\ \\ \mathbf{x}_g^{(4)} &= \texttt{update_g}(\Delta t, \mathbf{x}^{(1)}, \mathbf{x}_g^{(1)}) \\ \mathbf{x}^{(4)} &= \mathbf{x} + \Delta t \mathbf{f}^{(3)} \\ \mathbf{f}^{(4)} &= \mathbf{f}(\mathbf{x}^{(4)}, \mathbf{x}_g^{(4)}) \\ \\ \mathbf{x}^{n+1} &= \mathbf{x} + \frac{\Delta t}{6} \left( \mathbf{f}^{(1)} + 2\mathbf{f}^{(2)} + 2\mathbf{f}^{(3)} + \mathbf{f}^{(4)} \right) \\ \mathbf{x}_g^{n+1} &= \mathbf{x}_g^{(4)} \end{aligned} \]

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in,out]	X	Vector of state variables to be updated.
[in,out]	Xg	Vector of gating variables to be updated.
[in]	Ts	Current time.
[in]	Ti	Time step.
[in]	Istim	Applied current.
[in]	Ksac	Stretch-activated current coefficient.

nG()

unsigned int IonicModel::nG ( ) const

inline

Get the number of gating variables.

nX()

unsigned int IonicModel::nX ( ) const

inline

Get the number of state variables.

read_parameters()

void IonicModel::read_parameters ( const IonicModelParameters &  params )

virtual

Read model parameters from a parameter object.

By default, this method only takes care of the parameters related to initial states. If a derived ionic model has other parameters, then this method should be overridden to read those as well.

Reimplemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

update_g()

virtual void IonicModel::update_g ( const unsigned int  zone_id, const double  dt, const Vector< double > &  X, Vector< double > &  Xg  ) const

protectedpure virtual

Update variables with analytical solution.

Parameters

[in]	zone_id	Identifier for the transmural zone (epicardium, endocardium, myocardium).
[in]	dt	Time step.
[in]	X	Vector of state variables.
[out]	Xg	Vector of gating variables.

Implemented in AlievPanfilov, BuenoOrovio, FitzHughNagumo, and TTP.

Member Data Documentation

initial_X

InitialStates IonicModel::initial_X

protected

Initial states.

initial_Xg

InitialStates IonicModel::initial_Xg

protected

Initial gating variables.

Tscale

const double IonicModel::Tscale

protected

Time scaling [ms].

Voffset

const double IonicModel::Voffset

protected

Voltage offset parameter [mV].

Vrest

const double IonicModel::Vrest

protected

Resting transmembrane potential. It is used to define the stretch-activated current.

Vscale

const double IonicModel::Vscale

protected

Voltage scaling [mV].

---

The documentation for this class was generated from the following files:

• solver/ionic_model.h
• solver/ionic_model.cpp