Bueno-Orovio–Cherry–Fenton model

Short description: Phenomenological ionic model for human ventricular cells

The Bueno-Orovio–Cherry–Fenton model, also simply called Bueno-Orovio model, is a minimal ionic model for human ventricular cells. 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.

This mathematical model reproduces both single cell and important tissue-level properties, accounting for physiological action potential development and conduction velocity estimations. 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.

Mathematical modeling

Evolution in time only (i.e. case of a single cardiac cell) of the Bueno-Orovio ionic model variables [math]\displaystyle{ u, v, w, s }[/math]

The system of four differential equations reads as follows:

[math]\displaystyle{ \begin{cases} \dfrac{\partial u}{\partial t} = \nabla \cdot (D \nabla u) - ( J_{fi} + J_{so} + J_{si}) & \text{in } \Omega \times (0, T) \\[5pt] \dfrac{\partial v}{\partial t} = \dfrac{(1 - H(u - \theta_v)) (v_\infty - v)}{\tau_v^-} - \dfrac{H(u - \theta_v) v}{\tau_v^+} & \text{in } \Omega \times (0, T) \\[5pt] \dfrac{\partial w}{\partial t} = \dfrac{(1 - H(u - \theta_w)) (w_\infty - w)}{\tau_w^-} - \dfrac{H(u - \theta_w) w}{\tau_w^+} & \text{in } \Omega \times (0, T) \\[5pt] \dfrac{\partial s}{\partial t} = \dfrac{1}{\tau_s} \left( \dfrac{1 + \tanh(k_s (u - u_s))}{2}-s \right) & \text{in } \Omega \times (0, T) \\[5pt] \big( D \nabla u \big) \cdot \boldsymbol{N} = 0 & \text{on } \partial \Omega \times (0, T) \\[5pt] u = u_0, \, v = v_0, \, w = w_0, \, s = s_0 & \text{in } \Omega \times \{0\} \end{cases} }[/math]

where [math]\displaystyle{ \Omega }[/math] is the spatial domain and [math]\displaystyle{ T }[/math] is the final time. The initial conditions are [math]\displaystyle{ u_0 = 0 }[/math], [math]\displaystyle{ v_0 = 1 }[/math], [math]\displaystyle{ w_0 = 1 }[/math], [math]\displaystyle{ s_0 = 0 }[/math].

[math]\displaystyle{ H(x - x_0) }[/math] refers to the Heaviside function centered in [math]\displaystyle{ x_0 }[/math]. The adimensional transmembrane potential [math]\displaystyle{ u }[/math] can be rescaled in mV by means of the affine transformation [math]\displaystyle{ V_{mV} = 85.7 u - 84 }[/math]. [math]\displaystyle{ v }[/math], [math]\displaystyle{ w }[/math] and [math]\displaystyle{ s }[/math] refer to gating variables, where in particular [math]\displaystyle{ s }[/math] can be also used as an indication of intracellular calcium [math]\displaystyle{ {{Ca}_i}^{2+} }[/math] concentration (given in the adimensional range [0, 1] instead of molar concentration).

[math]\displaystyle{ J_{fi}, J_{so} }[/math] and [math]\displaystyle{ J_{si} }[/math] are the fast inward, slow outward and slow inward currents respectively, given by the following expressions:

[math]\displaystyle{ J_{fi} = -\frac{v H(u - \theta_v) (u - \theta_v) (u_u - u)}{\tau_{fi}}, }[/math]

[math]\displaystyle{ J_{so} = \frac{(u - u_o) (1 - H(u - \theta_w))}{\tau_o} + \frac{H(u - \theta_w)}{\tau_{so}}, }[/math]

[math]\displaystyle{ J_{si} = -\frac{H(u - \theta_w) w s}{\tau_{si}}, }[/math]

All the above-mentioned ionic density currents are partially adimensional and are expressed in [math]\displaystyle{ \dfrac{1}{\text{seconds}} }[/math].

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:

[math]\displaystyle{ \tau_v^- = (1 - H(u - \theta_v^-)) \tau_{v_1}^- + H(u - \theta_v^-) \tau_{v_2}^-, }[/math]

[math]\displaystyle{ \tau_w^- = \tau_{w_1}^- + (\tau_{w_2}^- - \tau_{w_1}^-)(1 + \tanh(k_w^-(u - u_w^-)))/2, }[/math]

[math]\displaystyle{ \tau_{so} = \tau_{{so}_1} + (\tau_{{so}_2} - \tau_{{so}_1})(1 + \tanh(k_{so} (u - u_{so})))/2, }[/math]

[math]\displaystyle{ \tau_s = (1 - H(u - \theta_w)) \tau_{s_1} + H(u - \theta_w) \tau_{s_2}, }[/math]

[math]\displaystyle{ \tau_o = (1 - H(u - \theta_o)) \tau_{o_1} + H(u - \theta_o) \tau_{o_2}, }[/math]

[math]\displaystyle{ v_{\infty} = 1 - H(u - \theta_v^-), }[/math]

[math]\displaystyle{ w_\infty = (1 - H(u - \theta_o))(1 - u/\tau_{w_\infty}) + H(u - \theta_o) w_\infty^*. }[/math]

where the temporal constants, i.e. [math]\displaystyle{ \tau_v^-, \tau_w^-, \tau_{so}, \tau_s, \tau_o }[/math] are expressed in seconds, whereas [math]\displaystyle{ v_{\infty} }[/math] and [math]\displaystyle{ w_{\infty} }[/math] are adimensional.

The diffusion coefficient [math]\displaystyle{ D }[/math] results in a value of [math]\displaystyle{ 1.171 \pm 0.221 \dfrac{\text{cm}^2}{\text{seconds}} }[/math], which comes from experimental tests on human ventricular tissues.

In order to trigger the action potential development in a certain position of the domain [math]\displaystyle{ \Omega }[/math], a forcing term [math]\displaystyle{ J_\text{app}(\boldsymbol x, t) }[/math], 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	[math]\displaystyle{ u_o }[/math]	[math]\displaystyle{ u_u }[/math]	[math]\displaystyle{ \theta_v }[/math]	[math]\displaystyle{ \theta_w }[/math]	[math]\displaystyle{ \theta_v^- }[/math]	[math]\displaystyle{ \theta_o }[/math]	[math]\displaystyle{ \tau_{v_1}^- }[/math]	[math]\displaystyle{ \tau_{v_2}^- }[/math]	[math]\displaystyle{ \tau_{v}^+ }[/math]	[math]\displaystyle{ \tau_{w_1}^- }[/math]	[math]\displaystyle{ \tau_{w_2}^- }[/math]	[math]\displaystyle{ k_w^- }[/math]	[math]\displaystyle{ u_w^- }[/math]	[math]\displaystyle{ \tau_w^+ }[/math]	[math]\displaystyle{ \tau_{fi} }[/math]	[math]\displaystyle{ \tau_{o_1} }[/math]	[math]\displaystyle{ \tau_{o_2} }[/math]	[math]\displaystyle{ \tau_{{so}_1} }[/math]	[math]\displaystyle{ \tau_{{so}_2} }[/math]	[math]\displaystyle{ k_{so} }[/math]	[math]\displaystyle{ u_{so} }[/math]	[math]\displaystyle{ \tau_{s_1} }[/math]	[math]\displaystyle{ \tau_{s_2} }[/math]	[math]\displaystyle{ k_s }[/math]	[math]\displaystyle{ u_s }[/math]	[math]\displaystyle{ \tau_{si} }[/math]	[math]\displaystyle{ \tau_{w_{\infty}} }[/math]	[math]\displaystyle{ w_{\infty}^* }[/math]
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 [math]\displaystyle{ \boldsymbol{z} }[/math] refers to the vector containing all the gating variables, i.e. [math]\displaystyle{ \boldsymbol{z} = [v, w, s]^T }[/math], and [math]\displaystyle{ \boldsymbol{F}:\mathbb{R}^4 \rightarrow \mathbb{R}^3 }[/math] contains the corresponding three right hand sides of the ionic model. The Bueno-Orovio–Cherry–Fenton model can be rewritten in the compact form:

[math]\displaystyle{ \begin{cases} \dfrac{\partial u}{\partial t} - \nabla \cdot (D \nabla u) + ( J_{fi} + J_{so} + J_{si} ) = 0 & \text{in } \Omega \times (0, T) \\[5pt] \dfrac{\partial \boldsymbol{z}}{\partial t} = \boldsymbol{F}(u, \boldsymbol{z}) & \text{in } \Omega \times (0, T) \\ \end{cases} }[/math]

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

To obtain the weak formulation:

• multiply by [math]\displaystyle{ p \in U }[/math] the first equation of the model and by [math]\displaystyle{ \boldsymbol{q} \in \boldsymbol{W} }[/math] the equations for the evolution of the gating variables. Integrating both members of all the equations in the domain [math]\displaystyle{ \Omega }[/math]:

[math]\displaystyle{ \begin{cases} \displaystyle \int_\Omega \dfrac{du(t)}{dt} p \,d\Omega - \int_\Omega \nabla \cdot (D \nabla u(t)) p \,d\Omega + \int_\Omega (J_{fi} + J_{so} + J_{si}) p \,d\Omega = 0 & \forall p \in U \\[5pt] \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} }[/math]

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

[math]\displaystyle{ - \int_\Omega \nabla \cdot (D \nabla u(t)) p \,d\Omega = \int_\Omega D \nabla u(t) \nabla p \, d\Omega - \cancelto{\text{Neumann B.C.}}{\int_{\partial \Omega} (D \nabla u(t)) \cdot \boldsymbol{N} p \,d\Omega} }[/math]

Finally the weak formulation reads:

Find [math]\displaystyle{ u \in L^2(0, T; H^1(\Omega)) }[/math] and [math]\displaystyle{ \boldsymbol{z} \in L^2(0, T; [L^2(\Omega)]^3) }[/math], [math]\displaystyle{ \forall t \in (0, T) }[/math], such that

[math]\displaystyle{ \begin{cases} \displaystyle \int_\Omega \frac{du(t)}{dt} p \,d\Omega + \int_{\Omega} D \nabla u(t) \cdot \nabla p \,d\Omega + \int_\Omega (J_{fi} + J_{so} + J_{si}) p \,d\Omega = 0 & \forall p \in U \\[5pt] \displaystyle \int_\Omega \frac{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} \\[5pt] u(0) = u_0 \\[5pt] \boldsymbol{z}(0) = \boldsymbol{z}_0 \end{cases} }[/math]

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).

Time discretization can be performed in several ways as well, such as using a backward differentiation formula (BDF) of order [math]\displaystyle{ \sigma }[/math] or a Runge–Kutta method (RK).

Space discretization with FEM

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

The semidiscretized formulation of the first equation of the model reads: find [math]\displaystyle{ u_h = u_h(t) = \sum_{j=1}^{N_h} \bar{u}_j(t) \phi_j }[/math] projection of the solution [math]\displaystyle{ u(t) }[/math] on [math]\displaystyle{ \mathcal{X}_{h}^r }[/math], [math]\displaystyle{ \forall t \in (0, T) }[/math], such that

[math]\displaystyle{ \int_\Omega {\dot{u}}_{h} \phi_i \, d \Omega + \int_\Omega (D \nabla u_h ) \cdot \nabla \phi_i \, d\Omega + \int_\Omega J_\text{ion}(u_h, \boldsymbol{z}_h) \phi_i \, d \Omega = 0 \quad \text{for } i = 1, \ldots, N_h }[/math]

with [math]\displaystyle{ u_h(0) = \sum_{j=1}^{N_h} \left( \int_\Omega u_0 \phi_j \,d\Omega \right) \phi_j }[/math], [math]\displaystyle{ \boldsymbol{z}_{h} = \boldsymbol{z}_h(t) = [v_h(t), w_h(t), s_h(t)]^T }[/math] semidiscretized version of the three gating variables, and [math]\displaystyle{ J_\text{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) }[/math] is the total ionic density current.

The space discretized version of the first equation can be rewritten as a system of non-linear ODEs by setting [math]\displaystyle{ \boldsymbol{U}_h(t) = \{ \bar{u}_j(t) \}_{j=1}^{N_h} }[/math] and [math]\displaystyle{ \boldsymbol{Z}_h(t) = \{ \bar{\boldsymbol{z}}_j(t) \}_{j=1}^{N_h} }[/math]:

[math]\displaystyle{ \begin{cases} \mathbb{M} \dot{\boldsymbol{U}}_h(t) + \mathbb{K} \boldsymbol{U}_h(t) + \boldsymbol{J}_\text{ion}(\boldsymbol{U}_h(t), \boldsymbol{Z}_h(t)) = 0 & \forall t \in (0, T) \\ \boldsymbol{U}_h(0) = \boldsymbol{U}_{0, h} \end{cases} }[/math]

where [math]\displaystyle{ \mathbb{M}_{ij} = \int_\Omega \phi_j \phi_i \,d \Omega }[/math], [math]\displaystyle{ \mathbb{K}_{ij} = \int_\Omega D \nabla \phi_j \cdot \nabla \phi_i \,d \Omega }[/math] and [math]\displaystyle{ \left( \boldsymbol{J}_\text{ion}(\boldsymbol{U}_{h}(t), \boldsymbol{z}_{h}(t)) \right)_i = \int_\Omega J_\text{ion}(u_h, \boldsymbol{z}_{h}) \phi_i \,d \Omega }[/math].

The non-linear term [math]\displaystyle{ \boldsymbol{J}_\text{ion}(\boldsymbol{U}_{h}(t), \boldsymbol{Z}_{h}(t)) }[/math] can be treated in different ways, such as using state variable interpolation (SVI) or ionic currents interpolation (ICI). In the framework of SVI, by denoting with [math]\displaystyle{ \{ \boldsymbol{x}_s^K \}_{s=1}^{N_q} }[/math] and [math]\displaystyle{ \{ \omega_s^K \}_{s=1}^{N_q} }[/math] the quadrature nodes and weights of a generic element of the mesh [math]\displaystyle{ K \in \mathcal{T}_{h} }[/math], both [math]\displaystyle{ u_h }[/math] and [math]\displaystyle{ \boldsymbol{z}_h }[/math] are evaluated at the quadrature nodes:

[math]\displaystyle{ \int_\Omega J_\text{ion}(u_h, \boldsymbol{z}_h) \phi_i \, d \Omega \approx \sum_{K \in \mathcal{T}_h} \left( \sum_{s=1}^{N_q} J_\text{ion} \left( \sum_{j=1}^{N_h} \bar{u}_j(t) \phi_j(\boldsymbol{x}_s^K), \sum_{j=1}^{N_h} \bar{\boldsymbol{z}}_j(t) \phi_j(\boldsymbol{x}_s^K) \right) \phi_i(\boldsymbol{x}_s^K) \omega_s^K \right) }[/math]

The equations for the three gating variables, which are ODEs, are directly solved in all the degrees of freedom (DOF) of the tessellation [math]\displaystyle{ \mathcal{T}_h }[/math] separately, leading to the following semidiscrete form:

[math]\displaystyle{ \begin{cases} \dot{\boldsymbol{Z}}_h(t) = \boldsymbol{F}(\boldsymbol{U}_h(t), \boldsymbol{Z}_h(t)) & \forall t \in (0, T) \\ \boldsymbol{Z}_h(0) = \boldsymbol{Z}_{0, h} \end{cases} }[/math]

Time discretization with BDF (implicit scheme)

With reference to the time interval [math]\displaystyle{ (0, T] }[/math], let [math]\displaystyle{ \Delta t = \dfrac{T}{N} }[/math] be the chosen time step, with [math]\displaystyle{ N }[/math] number of subintervals. A uniform partition in time [math]\displaystyle{ [t_0 = 0, t_1 = \Delta t, \ldots, t_k, \ldots, t_{N-1}, t_N = T] }[/math] is finally obtained.

At this stage, the full discretization of the Bueno-Orovio ionic model can be performed both in a monolithic and segregated fashion. With respect to the first methodology (the monolithic one), at time [math]\displaystyle{ t = t^k }[/math], the full problem is entirely solved in one step in order to get [math]\displaystyle{ (\boldsymbol{U}_h^{k + 1}, \boldsymbol{Z}_h^{k + 1}) }[/math] by means of either Newton method or Fixed-point iterations:

[math]\displaystyle{ \begin{cases} \mathbb{M} \alpha \dfrac{ \boldsymbol{U}_{h}^{k + 1} - \boldsymbol{U}_{\text{ext}, h}^{k}}{\Delta t} + \mathbb{K} \boldsymbol{U}_h^{k + 1} + \boldsymbol{J}_\text{ion}(\boldsymbol{U}_{h}^{k + 1}, \boldsymbol{Z}_{h}^{k + 1}) = 0 \\ \alpha \dfrac{\boldsymbol{Z}_h^{k + 1} - \boldsymbol{Z}_{\text{ext}, h}^{k}}{\Delta t} = \boldsymbol{F}(\boldsymbol{U}_h^{k + 1}, \boldsymbol{Z}_h^{k + 1} ) \end{cases} }[/math]

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

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, [math]\displaystyle{ \boldsymbol{Z}_h^{k + 1} }[/math] is calculated, using an extrapolation at previous timesteps [math]\displaystyle{ \boldsymbol{U}_{\text{ext}, h}^k }[/math] for the transmembrane potential at the right hand side:

[math]\displaystyle{ \alpha \dfrac{\boldsymbol{Z}_h^{k + 1} - \boldsymbol{Z}_{\text{ext}, h}^{k}}{\Delta t} = \boldsymbol{F}(\boldsymbol{U}_{\text{ext}, h}^k, \boldsymbol{Z}_h^{k + 1} ) }[/math]

• Secondly, [math]\displaystyle{ \boldsymbol{U}_h^{k + 1} }[/math] is computed, exploiting the value of [math]\displaystyle{ \boldsymbol{Z}_h^{k + 1} }[/math] that has just been calculated:

[math]\displaystyle{ \mathbb{M} \alpha \dfrac{ \boldsymbol{U}_h^{k + 1} - \boldsymbol{U}_{\text{ext}, h}^k}{\Delta t} + \mathbb{K} \boldsymbol{U}_h^{k + 1} + \boldsymbol{J}_\text{ion}(\boldsymbol{U}_h^{k + 1}, \boldsymbol{Z}_h^{k + 1}) = 0 }[/math]

Another possible segregated scheme would be the one in which [math]\displaystyle{ \boldsymbol{U}_{h}^{k + 1} }[/math] is calculated first, and then it is used in the equations for [math]\displaystyle{ \boldsymbol{Z}_h^{k + 1} }[/math].

See also

• Monodomain model
• Bidomain model