On boundary stimulation and optimal boundary control of the bidomain equations

Highlights

• Investigation of boundary stimulation and boundary control of bidomain equations.

• Development and implementation of optimal control approach to bidomain equations.

• Optimal control approach to successful defibrillation.

• Parallel finite element based algorithm is devised and its efficiency is demonstrated.

Abstract

The bidomain equations with Neumann boundary stimulation and optimal control of these stimuli are investigated. First an analytical framework for boundary control is provided. Then a parallel finite element based algorithm is devised and its efficiency is demonstrated not only for the direct problem but also for the optimal control problem. The computations realize a model configuration corresponding to optimal boundary defibrillation of a reentry phenomenon by applying current density stimuli.

Introduction

The heart is an electrically-controlled mechanical pump which drives blood through the circulatory system with remarkable efficiency. Under healthy conditions its electrical activation is highly organized, in disease, however, disturbances in the formation and/or propagation of electrical impulses may induce reentrant activation patterns which precipitate its rhythm significantly. Ultimately, such fast rhythms may evolve to highly disorganized almost chaotic activation patterns, an electrical state referred to as fibrillation. Under such conditions the heart looses its capacity to pump a sufficiently large mass of blood through the circulatory system. Without therapy, death would ensue within minutes. The only reliable therapy to terminate this otherwise lethal condition is the delivery of a strong electrical shock. This therapy, referred to as electrical defibrillation, is nowadays reliably achieved in a large patient population via the implantation of devices, so-called implantable cardioverters defibrillators (ICD), which monitor the heart rate and, if needed, deliver a discharge to restore a normal rhythm.

Although ICD therapy has proved to be efficient and reliable in preventing sudden cardiac death [2], it is far from ideal. There are several known adverse effects secondary to the administration of strong electrical shocks which are caused by the high field strengths required to terminate fibrillation with a sufficiently high probability. Moreover, psychological effects play an important role as well since shock delivery is perceived as extremely painful by conscious patients, leading to traumatization and a reduced quality of life. The link between the high shock strengths required and adverse effects provides the motivation for posing the defibrillation process as an optimization problem, where one aims to achieve defibrillation with minimal energy and, consequently, with minimal detrimental side effects.

The optimal control approach to defibrillation is to determine an applied electrical field in such a way that it optimizes a given design objective, which is, in our case, the restoration of a tissue state in which fibrillatory propagation cannot be maintained. This can be achieved by driving the whole tissue to a resting state, or equivalently, to an excited state. In both cases the main ingredients for maintaining fibrillation, namely the presence of both propagating wavefronts and a sufficient mass of excitable tissue at rest, referred to as “excitable gap”, in which these wavefronts can travel, are missing. Achieving these objectives is challenging since, on biophysical grounds, shock-induced changes in polarization of both polarities are always present during shock delivery [26], [22].

In a previous work [17], [16], [15] we addressed these points by modeling the controller action representing the current delivered by the electrodes as distributed forces. One of the main objectives of the current work consists in analyzing the case when the action of the electrodes is modeled as Neumann boundary conditions.

From a methodological point of view, in most, if not all, recent finite element modeling studies the effect of extracellularly applied electric fields has been accounted for either by imposing inhomogeneous Dirichlet boundary conditions to model extracellular potential stimuli, or, by using current volume sources to model current stimuli.

The use of inhomogeneous Neumann boundary conditions for modeling current stimuli has been the method of choice in early pioneering monodomain modeling work where the finite difference method was employed to model impulse propagation in 1D strands or 2D sheets [27]. Surprisingly, to the best of our knowledge, in the bidomain literature Neumann boundary conditions for modeling current injection via electrodes have not been rigorously stated yet, neither in bidomain forward models nor in the context of optimal control. While equivalence between both formulations, i.e. Neumann boundary conditions and volume sources, can be achieved for any given setup, in the latter case where currents are injected via shape functions into 2D or 3D elements, the total injected current depends on spatial discretization and choice of weighting function. Thus, in the present work we aim to investigate the suitability of using inhomogeneous Neumann boundary conditions in a bidomain model, specifically, the feasibility of optimal boundary control for the bidomain equations.

A second important issue of the current work consists in comparing the nonlinear conjugate gradient and the Newton method as iterative solution processes to solve the resulting optimization problems. Here we point out that due to the complicated dynamical systems behavior of the bidomain equations including wave phenomena which change significantly during the iterations of the optimal control scheme, it is out front not evident that the strength of improved rate of Newton’s method over gradient based methods can be achieved. Unless it is possible to obtain sufficiently accurate approximations to the Hessian the Newton method will provide no improvement or may fail. This behavior has been addressed in the passed for control of fluids, see e.g. [8], [10], but there is much less experience with optimal control of reaction diffusion systems.

Finally any numerical optimization approach requires repeated solution of the bidomain equations and the associated adjoint equations. While an efficient solution strategy is already of paramount importance for the direct numerical simulation of bidomain equation it is indispensable in the context of optimal control. For this reason our numerical realization relies on parallelization. We report on the parallel efficiency both for the direct simulation and for the optimization algorithms.

The optimal pacing of the cardiac tissue is expressed by optimal control with partial differential equations as constraints. Let $\mathrm{\Omega }\subset {\mathbb{R}}^{\mathbb{d}}\text{,}\mathbb{d}\in \{\mathbb{2}\text{,}\mathbb{3}\}$, denote a bounded connected domain with Lipschitz continuous boundary $\partial \mathrm{\Omega }$. The space–time domain and its lateral boundary are denoted by $Q=\mathrm{\Omega }\times (0\text{,}T]$ and $\mathrm{\Sigma }=\partial \mathrm{\Omega }\times (0\text{,}T]$, respectively. Also we denote the observation domain by ${\mathrm{\Omega }}_{\mathit{obs}}\subset \mathrm{\Omega }$. The standard form control problem is expressed as:

$$(\text{P})\left\{\begin{array}{c}\min J(v\text{,}{I}_e)\text{,}\hfill \\ e(u\text{,}v\text{,}w\text{,}{I}_e)=0\text{,}\hfill \end{array}\right.$$

where $u\text{,}v$ and w are the state variables, and $I_e$ is the extracellular current which is utilized as a control variable in the optimal control problem and $e=0$ stands formally for the dynamical system constraint. The dynamical behavior of the intra- and extracellular potentials is described by the coupled system of reaction–diffusion equations which can be expressed as follows

$$0=\nabla ·({\sigma }_i+{\sigma }_e)\nabla u+\nabla ·{\sigma }_i\nabla v\text{in}Q$$

$$\frac{\partial v}{\partial t}=\nabla ·{\sigma }_i\nabla v+\nabla ·{\sigma }_i\nabla u-I_{\mathit{ion}}(v\text{,}w)+I_{\mathit{tr}}(x\text{,}t)\text{in}Q$$

$$\frac{\partial w}{\partial t}=G(v\text{,}w)\text{in}Q\text{,}$$

where $u:Q\to \mathbb{R}$ is the extracellular potential, $v:Q\to \mathbb{R}$ is the transmembrane voltage, $w:Q\to {\mathbb{R}}^{\mathbb{n}}$ represents the ionic current variables, ${\sigma }_i:\mathrm{\Omega }\to {\mathbb{R}}^{\mathbb{d}\times \mathbb{d}}$ and ${\sigma }_e:\mathrm{\Omega }\to {\mathbb{R}}^{\mathbb{d}\times \mathbb{d}}$ are respectively the intracellular and extracellular conductivity tensors. The term $I_{\mathit{tr}}(x\text{,}t)$ is the transmembrane current density stimulus as delivered by an intracellular electrode. The $I_{\mathit{ion}}(v\text{,}w)$ is the current density flowing through the ionic channels and the function $G(v\text{,}w)$ determines the evolution of the gating variables, which is determined by an electrophysiological cell model, see e.g. [1] for more description on these models. Eq. (1.2) above is an elliptic type equation, Eq. (1.3) is a parabolic type equation and Eq. (1.4) is a set of ordinary differential equations which can be solved independently for each node. Typically, the conductivity tensors, which were considered in our computations, are expressed in the following form,

$${\sigma }_c=\left(\begin{array}{cc}{\sigma }_{\mathit{cl}}\hfill & 0\hfill \\ 0\hfill & {\sigma }_{\mathit{ct}}\hfill \end{array}\right)\text{,}\text{where}c=i\text{,}e\text{,}$$

where ${\sigma }_{\mathit{cl}}$ and ${\sigma }_{\mathit{ct}}$ are longitudinal and transverse fiber conductivities, respectively.

The membrane model for the ionic activity is described by a set of ordinary differential equations. The dimension of the ODE system is a consequence of the ionic model. In our numerical computations, we used a modified FitzHugh–Nagumo (FHN) model, called Rogers–McCulloch model [24], which consists of only two state variables and has a cubic non-linearity in the transmembrane potential

$$I_{\mathit{ion}}(v\text{,}w)=\mathit{gv}\left(1-\frac{v}{v_{\mathit{th}}}\right)\left(1-\frac{v}{v_p}\right)+{\eta }_1\mathit{vw}\text{.}$$

$$G(v\text{,}w)={\eta }_2\left(\frac{v}{v_p}-{\eta }_{3}w\right)$$

where $g\text{,}{\eta }_1\text{,}{\eta }_2\text{,}{\eta }_3$ are prescribed positive real coefficients, $v_{\mathit{th}}>0$ is the threshold potential and $v_p>v_{\mathit{th}}$ is the peak potential.

In the absence of a conductive bath both intracellular and extracellular domains are electrically isolated along the tissue boundaries and homogeneous Neumann boundary conditions are appropriate to reflect this fact, except for those parts of the boundary where extracellular stimuli are applied. The initial values of the transmembrane voltage and state variables are prescribed by given values $v_0\in L^2(\mathrm{\Omega })$ and $w_0\in L^2(\mathrm{\Omega })$. The initial and boundary conditions are therefore prescribed as

$$\eta ·({\sigma }_i\nabla v+{\sigma }_i\nabla u)=0\text{on}\mathrm{\Sigma }$$

$$\eta ·{\sigma }_e\nabla u=I_e\text{on}{\mathrm{\Gamma }}_{12}\times (0\text{,}T]$$

$$\eta ·{\sigma }_e\nabla u=0\text{on}{\mathrm{\Gamma }}_3\times (0\text{,}T]$$

$$v(x\text{,}0)=v_0\text{and}w(x\text{,}0)=w_0\text{on}\mathrm{\Omega }\text{,}$$

where $\eta$ denotes the outwards normal to the boundary of $\mathrm{\Omega }$. Here $I_e$ is the extracellular current density stimulus which acts as control along the boundary ${\mathrm{\Gamma }}_{12}={\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2$, where ${\mathrm{\Gamma }}_i\text{,}I=1\text{,}2\text{,}3$ are mutually disjoint and satisfy ${\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\cup {\mathrm{\Gamma }}_3=\partial \mathrm{\Omega }$. For compatibility reasons it is assumed throughout that

$${\int }_{\partial \mathrm{\Omega }}{I}_e(t\text{,}·)ds=0$$

for almost every $t\in (0\text{,}T)$. In the numerical experiments $I_e$ will be only temporally dependent and will be of the form

$$I_e={\widehat{I}}_e(t)({\chi }_{{\mathrm{\Gamma }}_1}-{\chi }_{{\mathrm{\Gamma }}_2})\text{,}$$

where ${\chi }_{{\mathrm{\Gamma }}_i}$ is the characteristic function of the set ${\mathrm{\Gamma }}_i\text{,}i=1\text{,}2$. Then condition (1.12) is satisfied if $|{\mathrm{\Gamma }}_1|=|{\mathrm{\Gamma }}_2|$. The support regions ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ can be considered to represent a cathode and an anode, respectively.

In terms of optimal control we recall that the optimality system involves the primal as well as the adjoint equations. Each of these two systems has similar complexity and must be solved within an iterative solution process. Also, the linearized primal and dual equations need to be solved frequently in each iteration of Newton–Krylov methods. The computational costs involved in this solution process are significant, rendering the use of efficient numerical approaches a key ingredient to obtain results within reasonable time frames. Parallelization techniques, as they are routinely used in virtually any bidomain modeling study, suggest themselves quite naturally also as a mean to speed up the solution process of the optimal control problem. In this regard, we parallelized our optimization codes based on the public domain package DUNE [4]. In the numerical section we compare the improvements in terms of computational cost that can be offered by Newton–Krylov optimization algorithm over gradient type algorithms.

The remainder of the paper is organized as follows. In Section 2, the boundary control formulation of the bidomain model is presented. There we discuss the necessary first order optimality conditions for the boundary control problem. Also, we lay out the necessary steps to implement Newton’s method in our numerical realization. The numerical discretization of the optimality system is described in Section 3. At the end of that section a brief overview of parallelization results is given. In Section 4 numerical results are presented for the different test cases. A summary is given then in the last section.