Bueno-Orovio–Cherry–Fenton model explained

The Bueno-Orovio–Cherry–Fenton model, also simply called Bueno-Orovio model, is a minimal ionic model for human ventricular cells.^[1] It belongs to the category of phenomenological models, because of its characteristic of describing the electrophysiological behaviour of cardiac muscle cells without taking into account in a detailed way the underlying physiology and the specific mechanisms occurring inside the cells.^{[2] [3]}

This mathematical model reproduces both single cell and important tissue-level properties, accounting for physiological action potential development and conduction velocity estimations.^[1] It also provides specific parameters choices, derived from parameter-fitting algorithms of the MATLAB Optimization Toolbox, for the modeling of epicardial, endocardial and myd-myocardial tissues. In this way it is possible to match the action potential morphologies, observed from experimental data, in the three different regions of the human ventricles.The Bueno-Orovio–Cherry–Fenton model is also able to describe reentrant and spiral wave dynamics, which occurs for instance during tachycardia or other types of arrhythmias.

From the mathematical perspective, it consists of a system of four differential equations. One PDE, similar to the monodomain model, for an adimensional version of the transmembrane potential, and three ODEs that define the evolution of the so-called gating variables, i.e. probability density functions whose aim is to model the fraction of open ion channels across a cell membrane.^[4]

Mathematical modeling

The system of four differential equations reads as follows:

\begin{cases} \dfrac{\partial u}{\partial t} = \nabla ⋅ (D \nabla u) - ( J _{fi} + J _{so} + J _{si} ) & in \Omega x ( 0 , T) \\[ 5 pt] \dfrac{\partial v}{\partial t} = \dfrac{( 1 - H(u - \theta _{v))} (v _{infty} - v)}{\tau ^{- }} _{v} - \dfrac{H(u - \theta _{v)} v}{\tau ^{+ }} _{v} & in \Omega x ( 0 , T) \\[ 5 pt] \dfrac{\partial w}{\partial t} = \dfrac{( 1 - H(u - \theta _{w))} (w _{infty} - w)}{\tau ^{- }} _{w} - \dfrac{H(u - \theta _{w)} w}{\tau ^{+ }} _{w} & in \Omega x ( 0 , T) \\[ 5 pt] \dfrac{\partial s}{\partial t} = \dfrac{ 1 }{\tau _{s}} \left( \dfrac{ 1 + \tanh(k _{s} (u - u _{s))}{2}-s} \right) & in \Omega x ( 0 , T) \\[ 5 pt] ( D \nabla u ) ⋅ \boldsymbol{N} = 0 & on \partial \Omega x ( 0 , T) \\[ 5 pt] u = u _{0,} v = v _{0,} w = w _{0,} s = s _{0} & in \Omega x \{ 0 \} \end{cases}

where

\Omega is the spatial domain and T is the final time. The initial conditions are u _{0} = 0 , v _{0} = 1 , w _{0} = 1 , s _{0} = 0 .H(x - x _{0)} refers to the Heaviside function centered in x _{0} . The adimensional transmembrane potential u can be rescaled in mV by means of the affine transformation V _{mV} = 85 . 7 u - 84 . v , w and s refer to gating variables, where in particular s can be also used as an indication of intracellular calcium {{Ca} ^{2+} _{i}} concentration (given in the adimensional range [0, 1] instead of molar concentration).^[5] J _{fi} , J _{so} and J _{si} are the fast inward, slow outward and slow inward currents respectively, given by the following expressions:J _{fi} =- (v H(u - \theta _{v)} (u - \theta _{v)} (u _{u} - u))/(\tau _{fi}) , J _{so} = ((u - u _{o)} ( 1- H(u - \theta _{w))})/(\tau _{o}) + (H(u - \theta _{w)})/(\tau _{so}) , J _{si} =- (H(u - \theta _{w)} w s)/(\tau _{si}) ,

All the above-mentioned ionic density currents are partially adimensional and are expressed in

\dfrac{ 1 }{ seconds }.

Different parameters sets, as shown in Table 1, can be used to reproduce the action potential development of epicardial, endocardial and mid-myocardial human ventricular cells. There are some constants of the model, which are not located in Table 1, that can be deduced with the following formulas:

\tau ^{-} _{v} = ( 1 - H(u - \theta ^{- ))} _{v} \tau ^{-} _{v _{1}} + H(u - \theta ^{- )} _{v} \tau ^{- ,} _{v _{2}} \tau ^{-} _{w} = \tau ^{-} _{w _{1}} + (\tau ^{-} _{w _{2}} - \tau ^{- )( 1} _{w _{1}} + \tanh(k ^{- (u} _{w} - u ^{- )))/ 2 ,} _{w} \tau _{so} = \tau _{{so} _{1}} + (\tau _{{so} _{2}} - \tau _{{so} _{1})(1} + \tanh(k _{so} (u - u _{so} )))/ 2 , \tau _{s} = ( 1 - H(u - \theta _{w))} \tau _{s _{1}} + H(u - \theta _{w)} \tau _{s _{2}} , \tau _{o} = ( 1 - H(u - \theta _{o))} \tau _{o _{1}} + H(u - \theta _{o)} \tau _{o _{2}} , v _{infty} = 1 - H(u - \theta ^{- ),} _{v} w _{infty} = ( 1 - H(u - \theta _{o))(1} - u/\tau _{w _{infty}} ) + H(u - \theta _{o)} w ^{*.} _{inf ty}

where the temporal constants, i.e.

\tau ^{- ,} _{v} \tau ^{- ,} _{w} \tau _{so} , \tau _{s,} \tau _{o} are expressed in seconds, whereas v _{infty} and w _{infty} are adimensional.

The diffusion coefficient

D results in a value of 1 . 171 \pm 0 . 221 \dfrac{ cm ^{2}{seconds} }, which comes from experimental tests on human ventricular tissues.

In order to trigger the action potential development in a certain position of the domain

\Omega , a forcing term J _{app(\boldsymbol} x, t) , which represents an externally applied density current, is usually added at the right hand side of the PDE and acts for a short time interval only.

Table 1: values of the parameters for different positions of the human heart

Parameter	u _{o}	u _{u}	\theta _{v}	\theta _{w}	\theta ^{-} _{v}	\theta _{o}	\tau ^{-} _{v _{1}}	\tau ^{-} _{v _{2}}	\tau ^{+} _{v}	\tau ^{-} _{w _{1}}	\tau ^{-} _{w _{2}}	k ^{-} _{w}	u ^{-} _{w}	\tau ^{+} _{w}	\tau _{fi}	\tau _{o _{1}}	\tau _{o _{2}}	\tau _{{so} _{1}}	\tau _{{so} _{2}}	k _{so}	u _{so}	\tau _{s _{1}}	\tau _{s _{2}}	k _{s}	u _{s}	\tau _{si}	\tau _{w _{infty}}	w ^{*} _{inf ty}
Unity of measure	-	-	-	-	-	-	seconds	seconds	seconds	seconds	seconds	-	-	seconds	seconds	seconds	seconds	seconds	seconds	-	-	seconds	seconds	-	-	seconds	-	-
Epicardium	0	1.55	0.3	0.13	0.006	0.006	60e-3	1150e-3	1.4506e-3	60e-3	15e-3	65	0.03	200e-3	0.11e-3	400e-3	6e-3	30.0181e-3	0.9957e-3	2.0458	0.65	2.7342e-3	16e-3	2.0994	0.9087	1.8875e-3	0.07	0.94
Endocardium	0	1.56	0.3	0.13	0.2	0.006	75e-3	10e-3	1.4506e-3	6e-3	140e-3	200	0.016	280e-3	0.1e-3	470e-3	6e-3	40e-3	1.2e-3	2	0.65	2.7342e-3	2e-3	2.0994	0.9087	2.9013e-3	0.0273	0.78
Myocardium	0	1.61	0.3	0.13	0.1	0.005	80e-3	1.4506e-3	1.4506e-3	70e-3	8e-3	200	0.016	280e-3	0.078e-3	410e-3	7e-3	91e-3	0.8e-3	2.1	0.6	2.7342e-3	4e-3	2.0994	0.9087	3.3849e-3	0.01	0.5

Weak formulation

Assume that

\boldsymbol{z} refers to the vector containing all the gating variables, i.e. \boldsymbol{z} = [v, w, s] ^{T} , and \boldsymbol{F}: R ^{4} → R ^{3} contains the corresponding three right hand sides of the ionic model. The Bueno-Orovio–Cherry–Fenton model can be rewritten in the compact form:^[6] \begin{cases} \dfrac{\partial u}{\partial t} - \nabla ⋅ (D \nabla u) + ( J _{fi} + J _{so} + J _{si} ) = 0 & in \Omega x ( 0 , T) \\[ 5 pt] \dfrac{\partial \boldsymbol{z}}{\partial t} = \boldsymbol{F}(u, \boldsymbol{z}) & in \Omega x ( 0 , T) \\ \end{cases}

Let

p \in U = H ^{1(\Omega)} and \boldsymbol{q} \in \boldsymbol{W} = [L ^{2(\Omega)]} ^{3} be two generic test functions.

To obtain the weak formulation:

• multiply by

p \in U the first equation of the model and by \boldsymbol{q} \in \boldsymbol{W} the equations for the evolution of the gating variables. Integrating both members of all the equations in the domain \Omega :\begin{cases} \displaystyle \int _{\Omega} \dfrac{du(t)}{dt} p d\Omega - \int _{\Omega} \nabla ⋅ (D \nabla u(t)) p d\Omega + \int _{\Omega} (J _{fi} + J _{so} + J _{si} ) p d\Omega = 0 & \forall p \in U \\[ 5 pt] \displaystyle \int _{\Omega} \dfrac{d\boldsymbol{z}(t)}{dt} \boldsymbol{q} d\Omega = \int _{\Omega} \boldsymbol{F}(u(t), \boldsymbol{z}(t)) \boldsymbol{q} d\Omega & \forall \boldsymbol{q} \in \boldsymbol{W} \end{cases}

• perform an integration by parts of the diffusive term of the first monodomain-like equation:

- \int _{\Omega} \nabla ⋅ (D \nabla u(t)) p d\Omega = \int _{\Omega} D \nabla u(t) \nabla p d\Omega - \cancelto{ NeumannB.C. }

Finally the weak formulation reads:

Find

u \in L ^{2(0,} T; H ^{1(\Omega))} and \boldsymbol{z} \in L ^{2(0,} T; [L ^{2(\Omega)]} ^{3)} , \forall t \in ( 0 , T) , such that^[6] \begin{cases} \displaystyle \int _{\Omega} (du(t))/(dt) p d\Omega + \int _{\Omega} D \nabla u(t) ⋅ \nabla p d\Omega + \int _{\Omega} (J _{fi} + J _{so} + J _{si} ) p d\Omega = 0 & \forall p \in U \\[ 5 pt] \displaystyle \int _{\Omega} (d\boldsymbol{z)/((t)}{dt}) \boldsymbol{q} d\Omega = \int _{\Omega} \boldsymbol{F}(u(t), \boldsymbol{z}(t)) \boldsymbol{q} d\Omega & \forall \boldsymbol{q} \in \boldsymbol{W} \\[ 5 pt] u( 0 ) = u _{0} \\[ 5 pt] \boldsymbol{z}( 0 ) = \boldsymbol{z} _{0 \end{cases}}

Numerical discretization

There are several methods to discretize in space this system of equations, such as the finite element method (FEM) or isogeometric analysis (IGA).^{[7] [8]}

Time discretization can be performed in several ways as well, such as using a backward differentiation formula (BDF) of order

\sigma or a Runge–Kutta method (RK).

Space discretization with FEM

Let

l{T} _{h} be a tessellation of the computational domain \Omega by means of a certain type of elements (such as tetrahedrons or hexahedra), with h representing a chosen measure of the size of a single element K \in l{T} _{h} . Consider the set P ^{r(K)} of polynomials with degree smaller or equal than r over an element K . Define l{X} ^{r} _{h} = \{f \in C ^{0(\bar} \Omega): f| _{K} \in P ^{r(K)} \forall K \in l{T} _{h} \} as the finite dimensional space, whose dimension is N _{h=\dim(l{X}} ^{r)} ^{h} . The set of basis functions of l{X} ^{r} _{h} is referred to as \{ \phi _{i} \} ^{N _{h}} _{i =1} .

The semidiscretized formulation of the first equation of the model reads: find

u _{h} = u _{h(t)} = \sum ^{N _{h}} _{j =1} \bar{u} _{j(t)} \phi _{j} projection of the solution u(t) on l{X} ^{r} _{h} , \forall t \in ( 0 , T) , such that^[5] \int _{\Omega} { • u }_ \phi_i \, d \Omega + \int_\Omega (D \nabla u_h) \cdot \nabla \phi_i \, d\Omega + \int_\Omega J_\text(u_h, \boldsymbol_h) \phi_i \, d \Omega = 0 \quad \text i = 1, \ldots, N_h

with

u _{h(0)} = \sum ^{N _{h}} _{j =1} \left( \int _{\Omega} u _{0} \phi _{j} d\Omega \right) \phi _{j} , \boldsymbol{z} _{h} = \boldsymbol{z} _{h(t)} = [v _{h(t),} w _{h(t),} s ^{T} _{h(t)]} semidiscretized version of the three gating variables, and J _{ion(u} _{h,} \boldsymbol{z} _{h} ) = J _{fi} (u _{h,} \boldsymbol{z} _{h)} + J _{so} (u _{h,} \boldsymbol{z} _{h)} + J _{si} (u _{h,} \boldsymbol{z} _{h)} is the total ionic density current.^[5]

The space discretized version of the first equation can be rewritten as a system of non-linear ODEs by setting

\boldsymbol{U} _{h(t)} = \{ \bar{u} _{j(t)} \} ^{N _{h}} _{j =1} and \boldsymbol{Z} _{h(t)} = \{ \bar{\boldsymbol{z}} _{j(t)} \} ^{N _{h}} _{j =1} :\begin{cases} M • \boldsymbol{U }_h(t) + \mathbb \boldsymbol_h(t) + \boldsymbol_\text(\boldsymbol_h(t), \boldsymbol_h(t)) = 0 & \forall t \in (0, T) \\\boldsymbol_h(0) = \boldsymbol_\end

where

M _{ij} = \int _{\Omega} \phi _{j} \phi _{i} d \Omega , K _{ij} = \int _{\Omega} D \nabla \phi _{j} ⋅ \nabla \phi _{i} d \Omega and\left( \boldsymbol{J} _{ion(\boldsymbol{U}} _{h} (t), \boldsymbol{z} _{h} (t)) \right) _{i} = \int _{\Omega} J _{ion(u} _{h,} \boldsymbol{z} _{h} ) \phi _{i} d \Omega .^[5]

The non-linear term

\boldsymbol{J} _{ion(\boldsymbol{U}} _{h} (t), \boldsymbol{Z} _{h} (t)) can be treated in different ways, such as using state variable interpolation (SVI) or ionic currents interpolation (ICI).^{[9] [10]} In the framework of SVI, by denoting with \{ \boldsymbol{x} ^{K} _{s} \} ^{N _{q}} _{s =1} and \{ \omega ^{K} _{s} \} ^{N _{q}} _{s =1} the quadrature nodes and weights of a generic element of the mesh K \in l{T} _{h} , both u _{h} and \boldsymbol{z} _{h} are evaluated at the quadrature nodes:\int _{\Omega} J _{ion(u} _{h,} \boldsymbol{z} _{h)} \phi _{i} d \Omega ≈ \sum _{K} _{h}} \left( \sum ^{N _{q}} _{s =1} J _{ion} \left( \sum ^{N _{h}} _{j =1} \bar{u} _{j(t)} \phi _{j(\boldsymbol{x}} ^{K),} ^{s} \sum ^{N _{h}} _{j =1} \bar{\boldsymbol{z}} _{j(t)} \phi _{j(\boldsymbol{x}} ^{K)} ^{s} \right) \phi _{i(\boldsymbol{x}} ^{K)} ^{s} \omega ^{K} _{s} \right)

The equations for the three gating variables, which are ODEs, are directly solved in all the degrees of freedom (DOF) of the tessellation

l{T} _{h} separately, leading to the following semidiscrete form:• \begin{cases} \boldsymbol{Z }_h(t) = \boldsymbol(\boldsymbol_h(t), \boldsymbol_h(t)) & \forall t \in (0, T) \\\boldsymbol_h(0) = \boldsymbol_\end

Time discretization with BDF (implicit scheme)

With reference to the time interval

( 0 , T] , let \Delta t = \dfrac{T}{N} be the chosen time step, with N number of subintervals. A uniform partition in time [t _{0} = 0 , t _{1} = \Delta t, \ldots, t _{k,} \ldots, t _{N-1} , t _{N} = T] is finally obtained.^[7]

At this stage, the full discretization of the Bueno-Orovio ionic model can be performed both in a monolithic and segregated fashion.^[11] With respect to the first methodology (the monolithic one), at time

t = t ^{k} , the full problem is entirely solved in one step in order to get (\boldsymbol{U} ^{k +1} _{h} , \boldsymbol{Z} ^{k +1} _{h} ) by means of either Newton method or Fixed-point iterations:\begin{cases} M \alpha \dfrac{ \boldsymbol{U} ^{k +1} _{h} - \boldsymbol{U} ^{k} _{ext , h} } + \mathbb \boldsymbol_h^ + \boldsymbol_\text(\boldsymbol_^, \boldsymbol_^) = 0 \\\alpha \dfrac = \boldsymbol(\boldsymbol_h^, \boldsymbol_h^)\end

where

\boldsymbol{U} ^{k} _{ext , h} and \boldsymbol{Z} ^{k} _{ext , h} are extrapolations of transmembrane potential and gating variables at previous timesteps with respect to t ^{k+1} , considering as many time instants as the order \sigma of the BDF scheme. \alpha is a coefficient which depends on the BDF order \sigma .

If a segregated method is employed, the equation for the evolution in time of the transmembrane potential and the ones of the gating variables are numerically solved separately:

• Firstly,

\boldsymbol{Z} ^{k +1} _{h} is calculated, using an extrapolation at previous timesteps \boldsymbol{U} ^{k} _{ext , h} for the transmembrane potential at the right hand side:\alpha \dfrac{\boldsymbol{Z} ^{k +1} _{h} - \boldsymbol{Z} ^{k} _{ext , h} } = \boldsymbol(\boldsymbol_^k, \boldsymbol_h^)

• Secondly,

\boldsymbol{U} ^{k +1} _{h} is computed, exploiting the value of \boldsymbol{Z} ^{k +1} _{h} that has just been calculated:M \alpha \dfrac{ \boldsymbol{U} ^{k +1} _{h} - \boldsymbol{U} ^{k}{\Delta} _{ext , h} t} + K \boldsymbol{U} ^{k +1} _{h} + \boldsymbol{J} _{ion(\boldsymbol{U}} ^{k +1} ^{h} , \boldsymbol{Z} ^{k +1} _{h} ) = 0

Another possible segregated scheme would be the one in which

\boldsymbol{U} ^{k +1} _{h} is calculated first, and then it is used in the equations for \boldsymbol{Z} ^{k +1} _{h} .

See also

• Monodomain model
• Bidomain model

Notes and References

1. Bueno-Orovio . A. . Cherry . E.M. . Fenton . F. H. . Minimal model for human ventricular action potentials in tissue . Journal of Theoretical Biology . August 2008 . 253 . 3 . 544–560 . 10.1016/j.jtbi.2008.03.029. 18495166 .
2. Book: Colli Franzone . P. . Pavarino . L.F. . Scacchi . S. . Mathematical cardiac electrophysiology . 30 October 2014 . 978-3-319-04801-7.
3. Book: Sundnes . J. . Lines . G.T. . Cai . X. . Nielsen . B.F. . Mardal . K.-A. . Tveito . A. . Computing the electrical activity in the heart . 26 June 2007 . Springer . 978-3-540-33437-8.
4. Book: Keener . J. . Sneyd . J. . Mathematical physiology . 27 October 2008 . Springer . 978-0-387-79387-0 . 2nd.
5. Gerbi . A. . Dede’ . L. . Quarteroni . A. . A monolithic algorithm for the simulation of cardiac electromechanics in the human left ventricle . Mathematics in Engineering. 2018 . 1 . 1 . 1–37 . 10.3934/Mine.2018.1.1. free . 11311/1066015 . free .
6. Pegolotti . L. . Dedè . L. . Quarteroni . A. . Isogeometric Analysis of the electrophysiology in the human heart: Numerical simulation of the bidomain equations on the atria . Computer Methods in Applied Mechanics and Engineering . January 2019 . 343 . 52–73 . 10.1016/j.cma.2018.08.032. 2019CMAME.343...52P .
7. Book: Quarteroni . A. . Numerical models for differential problems . 25 April 2014 . 978-88-470-5522-3 . Second.
8. Book: Cottrell . J. . Hughes . T.J.R. . Bazilevs . Y. . Isogeometric analysis: toward integration of CAD and FEA . 15 September 2009 . Wiley . 978-0-470-74873-2.
9. Pathmanathan . P. . Mirams . G.R. . Southern . J. . Whiteley . J.P. . The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations . International Journal for Numerical Methods in Biomedical Engineering . November 2011 . 27 . 11 . 1751–1770 . 10.1002/cnm.1438.
10. Pathmanathan . P. . Bernabeu . M.O. . Niederer . S.A. . Gavaghan . D.J. . Kay . D. . Computational modelling of cardiac electrophysiology: explanation of the variability of results from different numerical solvers . International Journal for Numerical Methods in Biomedical Engineering . August 2012 . 28 . 8 . 890–903 . 10.1002/cnm.2467. 25099569 .
11. Gerbi . A. . Dede' . L. . Quarteroni . A. . Numerical approximation of cardiac electro-fluid-mechanical models: coupling strategies for large-scale simulation . PhD Thesis . École Polytechnique Fédérale de Lausanne.

This article is licensed under the GNU Free Documentation License. It uses material from the Wikipedia article "Bueno-Orovio–Cherry–Fenton model".

Except where otherwise indicated, Everything.Explained.Today is © Copyright 2009-2024, A B Cryer, All Rights Reserved. Cookie policy.