The Effects of Control Domain Position on Optimal Control of Cardiac Arrhythmia

Abstract

We study the effects of control domain position on optimal control problem of cardiac arrhythmia. The well-established mathematical model, namely, the monodomain model, is used to represent the electrical behavior of the cardiac tissue. Two test cases are considered in this paper: the control domain for Test Case 1 is close to the excitation domain, while the control domain for Test Case 2 is a little bit far away from the excitation domain. The modified Dai-Yuan nonlinear conjugate gradient method is employed for computing the optimal applied extracellular current. Our experiment results show that the excitation wavefronts are successfully dampened out for both test cases. Also, experiment results show that the position of the control domain for Test Case 1 is preferable than that of Test Case 2 during the defibrillation process.

1. Introduction

Cardiac arrhythmia refers to the condition in which normal heart rhythm is disrupted. During arrhythmia, the heart may beat too rapidly, too slowly, or even irregularly. This abnormal and irregular heart rhythm causes the heart muscle to fibrillate and stops pumping blood to the body and the brain. Eventually, death can occur within minutes unless the normal heart rhythm is restored. For prevention of this sudden cardiac death, the only effective therapy is through electrical defibrillation [1–3].

The optimal control approach to cardiac arrhythmia was first proposed by Nagaiah et al. [4] with the control objective being dampen the excitation wavefront of the transmembrane potential using optimal applied extracellular current. Since the monodomain model is chosen by Nagaiah et al. [4] for representing the electrical behavior of the cardiac tissue, we, therefore, call it optimal control problem of monodomain model in the rest of the paper.

The optimal control problem of monodomain model is a PDE-constrained optimization problem that is constrained by the monodomain model. The monodomain model is a well-established mathematical model extensively used for simulation of cardiac electrical activity [5–7]. It consists of a parabolic partial differential equation (PDE) coupled to a system of nonlinear ordinary differential equations (ODEs) representing cell ionic activity.

Three classical nonlinear conjugate gradient methods have been employed by Nagaiah et al. [4] for solving the optimal control problem of monodomain model, namely, the Dai-Yuan (DY) method [8], the Hager-Zhang (HZ) method [9], and the Polak-Ribière-Polyak (PRP) method [10, 11]. Later, the optimal control problem of monodomain model has been subsequently studied by researchers with attempt to develop efficient numerical solving techniques [12–14].

Recently, the effects of the size of the control domain on the optimal control problem of monodomain model have been studied by Ng and Rohanin [15]. Two test cases with different sizes of the control domain are considered, and their numerical results show that higher extracellular current is required in the dampening process when the size of the control domain changed to a smaller one. Also, longer time is needed to dampen out the excitation wavefront when small size control domain is used.

In this paper, we are concerned with the effects of the control domain position on the optimal control problem of monodomain model. Again, two test cases are considered with each one of them located at different positions. The indirect method, also known as the optimize-then-discretize approach is used to discretize the optimal control problem. The discretized optimal control problem is then solved by the modified nonlinear conjugate gradient method, namely, the modified Dai-Yuan (MDY) method. As shown in [13], the MDY method is well performed under Armijo line search.

The outline of this paper is as follows. In Section 2, we introduce the optimal control problem of monodomain model with Rogers-modified FitzHugh-Nagumo ion kinetic. In Section 3, we present the numerical approach used to discretize the optimal control problem. In Section 4, we discuss the optimization procedure used to solve the discretized optimal control problem. In Section 5, we present and compare the experiment results for both test cases. In Section 6, we conclude our paper with a short discussion.

2. Optimal Control Problem of Monodomain Model

Let 𝑇 be the final simulation time, Ω⊂ℜ² the computational domain with Lipschitz boundary 𝜕Ω, Ω_𝑐⊂Ω the control domain, and Ω_𝑜⊂Ω the observation domain. We further set 𝐻=Ω×[0,𝑇] and 𝜕𝐻=𝜕Ω×[0,𝑇]. The optimal control problem of monodomain model with Rogers-modified FitzHugh-Nagumo ion kinetic can be summarized as follows: min𝐽𝑉,𝑤,𝐼_𝑒=12^𝑇_0Ω𝑜|𝑉|²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐𝜆𝑑𝑡,s.t.𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒=0,in𝐻,𝜕𝑤()𝐷𝜕𝑡−𝑓𝑉,𝑤=0,in𝐻,_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻,𝑥,0=𝑉⁰,𝑤(𝑥,0)=𝑤⁰,onΩ, (2.1)(2.2)(2.3) where 𝐼_ion(𝑉,𝑤)=𝑐₁𝑉𝑉1−𝑉_th𝑉1−𝑉_𝑝+𝑐₂𝑓()𝑤𝑉,𝑉,𝑤=𝑐₃𝑉𝑉_𝑝−𝑐₄𝑤. (2.4) Here, 𝛼 is the regularization parameter, 𝜂 is the outer normal to Ω, 𝐷_𝑖 is the intracellular conductivity tensor, 𝜆 is the constant scalar used to relate the intracellular and extracellular conductivity tensors, 𝛽 is the surface-to-volume ratio of the cell membrane, 𝐶_𝑚 is the membrane capacitance, 𝑉_th is the threshold potential, 𝑉_𝑝 is the plateau potential, 𝑐₁,𝑐₂,𝑐₃,and𝑐₄ are positive parameters, 𝑉(𝑥,𝑡) is the transmembrane potential, 𝑤(𝑥,𝑡) is the ionic current variable, 𝐼_𝑒(𝑥,𝑡) is the extracellular current density stimulus, 𝐼_ion(𝑉,𝑤) is the current density flowing through the ionic channels, and 𝑓(𝑉,𝑤) is the prescribed vector-value function.

Note that (2.1) is the cost functional, in which the transmembrane potential 𝑉 and the ionic current variable 𝑤 act as the state variables and the extracellular current 𝐼_𝑒 acts as the control variable. Also, the control variable 𝐼_𝑒 is restricted to be nontrivial only on the control domain Ω_𝑐. The monodomain model, as defined in (2.2), appears as the PDE constraints in the optimal control problem of monodomain model. Next, the boundary and initial conditions for this optimal control problem are given in (2.3), and, lastly, (2.4) is taken from the Rogers-modified FitzHugh-Nagumo model [16].

3. Discretization of the Optimal Control Problem

The optimal control problem of monodomain model is discretized using the optimize-then-discretize approach, where the infinite dimensional first-order optimality system is derived first and the resulting optimality system is then discretized. Let 𝑝(𝑥,𝑡) and 𝑞(𝑥,𝑡) denote the adjoint variables used to enforce the constraints. Define the Lagrange functional ℒ as 1ℒ=2^𝑇_0Ω𝑜|𝑉|²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐+𝑑𝑡_𝑇0Ω𝜆𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒+𝑝𝑑Ω𝑑𝑡_𝑇0Ω𝜕𝑤()𝜕𝑡−𝑓𝑉,𝑤𝑞𝑑Ω𝑑𝑡. (3.1) The first-order optimality system for this optimal control problem is given by: ℒ_𝑝:𝜆𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒ℒ=0,_𝑞:𝜕𝑤()ℒ𝜕𝑡−𝑓𝑉,𝑤=0,_𝑉:𝑉|_𝑜+𝜆𝐷1+𝜆∇⋅_𝑖∇𝑝+𝛽𝐶_𝑚𝜕𝑝𝐼𝜕𝑡−𝛽_ion𝑉𝑝−𝑓_𝑉ℒ𝑞=0,_𝑤:−𝛽𝐼_ion𝑤𝑝−𝜕𝑞−𝑓𝜕𝑡^𝑤ℒ𝑞=0,_𝐼𝑒:𝛼𝐼_𝑒−11+𝜆𝑝=0, (3.2)(3.3)(3.4) where [⋅]_∗ denotes the partial derivative with respect to ∗ and 𝑉|_𝑜 denotes the transmembrane potential in Ω_𝑜. Here, (3.2) is known as the state equations, (3.3) is known as the adjoint equations and (3.4) is known as the optimality condition. We further obtain the boundary and terminal conditions as follows: 𝐷_𝑖𝑝()()∇𝑝⋅𝜂=0,on𝜕𝐻,𝑥,𝑇=0,𝑞𝑥,𝑇=0,onΩ. (3.5)

Next, the state system can be formed using the boundary and initial conditions in (2.3) as well as the state equations in (3.2). As a result, the state system is given by 𝛽𝐶_𝑚𝜕𝑉=𝜆𝜕𝑡𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒,in𝐻,𝜕𝑤()𝐷𝜕𝑡=𝑓𝑉,𝑤,in𝐻,_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻,𝑥,0=𝑉⁰,𝑤(𝑥,0)=𝑤⁰,onΩ. (3.6) On the other hand, the adjoint system can be formed using the adjoint equations in (3.3) as well as the boundary and terminal conditions in (3.5). Thus, the adjoint system is given by 𝛽𝐶_𝑚𝜕𝑝𝜆𝜕𝑡=−𝐷1+𝜆∇⋅_𝑖∇𝑝+𝛽𝐼_ion𝑉𝑝+𝑓_𝑉𝑞−𝑉|_𝑜,in𝐻,𝜕𝑞𝐼𝜕𝑡=−𝛽_ion𝑤𝑝−𝑓_𝑤𝐷𝑞,in𝐻,_𝑖𝑝()()∇𝑝⋅𝜂=0,on𝜕𝐻,𝑥,𝑇=0,𝑞𝑥,𝑇=0,onΩ. (3.7)

Kunisch and Wagner [17] proved that the control-to-state mapping is well defined for the optimal control problem of monodomain model, that is, 𝐶∋𝐼_𝑒↦(𝑉(𝐼_𝑒),𝑤(𝐼_𝑒)). Thus, the cost functional in (2.1) may be rewritten as min𝐽̂𝐼_𝑒=𝐽𝑉𝐼_𝑒,𝑤𝐼_𝑒,𝐼_𝑒=12^𝑇_0Ω𝑜𝑉𝐼_𝑒²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐𝑑𝑡. (3.8) From the previous analysis, one can shows that the gradient of 𝐽(̂𝐼_𝑒) is given by ∇𝐽̂𝐼_𝑒=𝛼𝐼_𝑒−11+𝜆𝑝. (3.9) Clearly, the optimality condition in (3.4) serves as the reduced gradient in (3.9) for our computation. As a result, the first-order optimality system consists of the state system in (3.6), the adjoint system in (3.7), and the reduced gradient in (3.9).

Notice that the state and adjoint systems in (3.6) and (3.7) consist of a parabolic PDE coupled to a system of nonlinear ODEs. Consequently, these systems are computationally demanding. In order to overcome this difficulty, the operator splitting technique [18] is applied to (3.6) and (3.7) for decomposing the complex systems into subsystems. Eventually, the state system becomes linearPDE:𝛽𝐶_𝑚𝜕𝑉=𝜆𝜕𝑡𝐷1+𝜆∇⋅_𝑖∇𝑉,in𝐻,nonlinearODE:𝛽𝐶_𝑚𝜕𝑉𝜕𝑡=−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒,in𝐻,nonlinearODE:𝜕𝑤()𝐷𝜕𝑡=𝑓𝑉,𝑤,in𝐻,_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻,𝑥,0=𝑉⁰𝑤(),onΩ,𝑥,0=𝑤⁰,onΩ, (3.10) and the adjoint system becomes linearPDE:𝛽𝐶_𝑚𝜕𝑝𝜆𝜕𝑡=−𝐷1+𝜆∇⋅_𝑖∇𝑝,in𝐻,nonlinearODE:𝛽𝐶_𝑚𝜕𝑝𝐼𝜕𝑡=𝛽_ion𝑉𝑝+𝑓_𝑉𝑞−𝑉|_𝑜,in𝐻,nonlinearODE:𝜕𝑞𝐼𝜕𝑡=−𝛽_ion𝑤𝑝−𝑓_𝑤𝐷𝑞,in𝐻,_𝑖𝑝()𝑞()∇𝑝⋅𝜂=0,on𝜕𝐻,𝑥,𝑇=0,onΩ,𝑥,𝑇=0,onΩ. (3.11)

In order to complete the optimize-then-discretize approach, the first-order optimality system needs to be discretized. The linear PDEs in (3.10) and (3.11) are discretized using Galerkin finite element method in space and Crank-Nicolson method in time. For time discretization of the nonlinear ODEs in (3.10) and (3.11), a forward Euler method is used. The discretized state system is, therefore, given by 𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛+1=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛,𝐕^𝑛+1=𝐕^𝑛+Δ𝑡₂−𝐈^𝑛_𝐢𝐨𝐧𝐶_𝑚−𝐈^𝑛_𝐞𝛽𝐶_𝑚(),𝐰1+𝜆_𝑛+1=𝐰^𝑛+Δ𝑡₂𝐟^𝑛,𝐕(𝑥,0)=𝐕^𝟎,𝐰(𝑥,0)=𝐰^𝟎, (3.12) and the discretized adjoint system is given by 𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛+1,𝐩^𝑛=𝐩^𝑛+1+Δ𝑡₂𝐕^𝑛+1_𝑜𝛽𝐶_𝑚−𝐈^𝑛+1_{𝐢𝐨𝐧𝑉}𝐶_𝑚𝐩^𝑛+1−𝐟^𝑛+1_𝑉𝛽𝐶_𝑚𝐪^𝑛+1,𝐪^𝑛=𝐪^𝑛+1+Δ𝑡₂𝛽𝐈^𝑛+1_{𝐢𝐨𝐧𝑤}𝐩^𝑛+1+𝐟^𝑛+1_𝑤𝐪^𝑛+1,𝐩(𝑥,𝑇)=𝟎,𝐪(𝑥,𝑇)=𝟎. (3.13) Here, 𝑛 is the time index, Δ𝑡₁ and Δ𝑡₂ are the local time steps, 𝐌 is the mass matrix, 𝐊 is the stiffness matrix, and 𝐕, 𝐕|_𝑜, 𝐩, 𝐪,𝐰,  𝐟, 𝐈_𝐢𝐨𝐧, and  𝐈_𝐞 are the vectors of the same size.

4. Optimization Procedure

In this paper, we adopt the modified nonlinear conjugate gradient method, namely, the modified Dai-Yuan (MDY) method as proposed by Zhang [19] for solving the discretized optimal control problem of monodomain model. Starting from an initial guess 𝐈⁰_𝐞, the control variable is updated using the following recurrence: 𝐈_𝐞^𝑘+1=𝐈^𝑘_𝐞+𝛿^𝑘𝐝^𝑘,𝑘=0,1,…, (4.1) where 𝛿^𝑘>0 is the step length computed by the Armijo line search and 𝐝^𝑘 is the search direction generated by 𝐝^𝑘=−∇𝐽̂𝐈^𝑘_𝐞𝐽̂𝐈,if𝑘=0,−∇^𝑘_𝐞+𝜃^𝑘_MDY𝐝^𝑘−1,if𝑘>0, (4.2) where 𝜃^𝑘_MDY∈ℜ is the conjugate gradient update parameter given by 𝜃^𝑘_MDY=∇𝐽̂𝐈^𝑘_𝐞²𝐝^𝑘−1𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1+𝛾^𝑘𝛿^𝑘−1𝐝^𝑘−1, (4.3) where 𝛾^𝑘=10⁻³∇𝐽̂𝐈^𝑘_𝐞²𝛿+max0,−^𝑘−1𝐝^𝑘−1𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1𝛿^𝑘−1𝐝^𝑘−12. (4.4)

Recall that the step length 𝛿^𝑘 is computed by the Armijo line search. Given an initial step length 𝛿>0 and 𝜇,𝜌∈(0,1), the Armijo line search choose, 𝛿^𝑘=max{𝛿,𝛿𝜌,𝛿𝜌²,…} such that 𝐽̂𝐈^𝑘_𝐞+𝛿^𝑘𝐝^𝑘≤𝐽̂𝐈^𝑘_𝐞+𝜇𝛿^𝑘∇𝐽̂𝐈^𝑘_𝐞^𝑇𝐝^𝑘. (4.5) Consequently, the optimization algorithm for the MDY method can be summarized as follows.

4.1. Optimization Algorithm

Step 1. Provide an initial guess 𝐈⁰_𝐞 and set 𝑘=0.

Step 2. Solve the discretized state system in (3.12) for time interval [0,𝑇] as follows.(a) Solve the linear PDE for 𝑡∈[𝑡,𝑡+Δ𝑡/2]𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛+1=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛. (4.6)(b) Solve the nonlinear ODEs for 𝑡∈[𝑡,𝑡+Δ𝑡]𝐕^𝑛+1=𝐕^𝑛+Δ𝑡₂−𝐈^𝑛_𝐢𝐨𝐧𝐶_𝑚−𝐈^𝑛_𝐞𝛽𝐶_𝑚(),𝐰1+𝜆_𝑛+1=𝐰^𝑛+Δ𝑡₂𝐟^𝑛. (4.7)(c) Solve the linear PDE for 𝑡∈[𝑡+Δ𝑡/2,𝑡+Δ𝑡]𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛+1=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛, (4.8) where Δ𝑡 is the global time step.

Step 3. Evaluate the reduced cost functional 𝐽(̂𝐈^𝑘_𝐞).

Step 4. Solve the discretized adjoint system in (3.13) for time interval [𝑇,0] as follows.(a) Solve the linear PDE for 𝑡∈[𝑇,𝑇−Δ𝑡/2]𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛+1. (4.9)(b) Solve the nonlinear ODEs for 𝑡∈[𝑇,𝑇−Δ𝑡]𝐩^𝑛=𝐩^𝑛+1+Δ𝑡₂𝐕^𝑛+1_𝑜𝛽𝐶_𝑚−𝐈^𝑛+1_{𝐢𝐨𝐧𝑉}𝐶_𝑚𝐩^𝑛+1−𝐟^𝑛+1_𝑉𝛽𝐶_𝑚𝐪^𝑛+1𝐪^𝑛=𝐪^𝑛+1+Δ𝑡₂𝛽𝐈^𝑛+1_{𝐢𝐨𝐧𝑤}𝐩^𝑛+1+𝐟^𝑛+1_𝑤𝐪^𝑛+1. (4.10)(c) Solve the linear PDE for 𝑡∈[𝑇−Δ𝑡/2,𝑇−Δ𝑡]𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛+1. (4.11)

Step 5. Update the reduced gradient ∇𝐽(̂𝐈^𝑘_𝐞) using (3.9).

Step 6. For 𝑘≥1, check the following stopping criteria: 𝐽̂𝐈^𝑘_𝐞−𝐽̂𝐈_𝐞^𝑘−1≤10⁻³,∇𝐽̂𝐈^𝑘_𝐞≤10⁻³1+𝐽̂𝐈^𝑘_𝐞. (4.12)
If one of them is met, stop.

Step 7. Compute the conjugate gradient update parameter 𝜃^𝑘_MDY using (4.3).

Step 8. Compute the search direction 𝐝^𝑘 using (4.2).

Step 9. Compute the step length 𝛿^𝑘 that satisfies the condition in (4.5).

Step 10. Update the control variable 𝐈_𝐞^𝑘+1 using (4.1). Set 𝑘=𝑘+1 and go to Step 2.

5. Numerical Experiments

5.1. Experiment Setup

A two-dimensional computational domain of size 1×1 cm², that is, Ω=[0,1]×[0,1] is considered in our numerical experiments, and the final simulation time is set to be 𝑇=2 ms. Figure 1 displays the positions of the subdomains in the computational domain Ω for Test Case 1 and Test Case 2. As shown in the figures, Ω_𝑐1 and Ω_𝑐2 are the control domains, Ω_𝑐1 and Ω_𝑐2 are the neighborhoods of the control domains, Ω_𝑜=Ω⧵(Ω_𝑐1∪Ω_𝑐2) is the observation domain, and Ω_exi⊂Ω_𝑜 is the excitation domain where the cardiac arrhythmia first occurs.

(a)

(b)

Figure 1 

Computational domain Ω and its subdomains for (a) Test Case 1 and (b) Test Case 2.

Also, observe that the control domain for Test Case 1 is close to the excitation domain while the control domain for Test Case 2 is a little bit far away from the excitation domain. These positions of the control domain are purposely set for studying the effects of the control domain position on the optimal control problem of monodomain model. In our numerical experiments, the control domain might corresponds to the electrodes of the implantable cardioverter defibrillator (ICD) implanted in the chest of patients who are at significant risk of sudden cardiac death. More specifically, we are trying to locate a better position for implanting the electrodes of the ICD to ensure a successful and efficient defibrillation.

Table 1 lists the parameters that are used in our numerical experiments, with some of them adopted from Franzone et al. [20]. Moreover, the initial values for 𝑉, 𝑤, 𝐼_𝑒 are given as 𝑉(𝑥,0)=𝑉⁰=105mV,𝑥∈Ω_exi,𝑤()0mV,otherwise,𝑥,0=𝑤⁰𝐼=0,𝑥∈Ω,_𝑒(𝑥,0)=𝐼⁰_𝑒=0mAcm⁻³,𝑥∈Ω_𝑐,0mAcm⁻³,otherwise. (5.1)(5.2)(5.3) As shown in (5.3), the initial value for the control variable is set to be zero in the control domain. Meaning that no extracellular current is applied to the patient at the beginning of the numerical experiments.

Table 1 

Parameters used in the numerical experiment.

5.2. Experiment Results

The minimum values of the reduced cost functional 𝐽(̂𝐈^𝑘_𝐞) and the corresponding norm of reduced gradient ‖∇𝐽(̂𝐈^𝑘_𝐞)‖ along the optimization process for Test Case 1 and Test Case 2 are depicted in Figure 2. As shown in the figures, Test Case 1 and Test Case 2 successfully converge to the optimal solution by taking 18 and 27 optimization iterations. Also, observe that the minimum value for Test Case 1 is lower than that one for Test Case 2. This result implies that when the control domain is relocated to a position which is far away from the excitation domain, more iteration is needed for the MDY method to converge to the optimal solution.

(a)

(b)

Figure 2 

(a) Minimum values of reduced cost functional 𝐽̂ and (b) norm of reduced gradient ‖∇𝐽‖̂ for Test Case 1 and Test Case 2.

As illustrated in Figure 2, the reduced gradient for Test Case 1 decreased sharply at iteration 3, followed by a smooth decrease till the end of optimization iterations. On the other hand, the reduced gradient for Test Case 2 decreases less sharply than that of Test Case 1 at iteration 3 and slowly approaches zero with a smooth decrease.

Figure 3 illustrates the uncontrolled solutions for the optimal control problem of monodomain model at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms. Notice that the uncontrolled solutions are obtained by setting the control variable equal to zero at the control domain as well as at the outside of the control domain, that is, 𝐼_𝑒=0 for 𝑥∈Ω. As shown in the figures, the excitation wavefront spreads from the excitation domain to the nearby computational domain during the time interval from 0 ms to 2 ms. This result implies that if no extracellular current (𝐼_𝑒) is applied to the control domain, the excitation wavefront will continue to spread to the entire computational domain.

(a)

(b)

(c)

(d)

Figure 3 

The uncontrolled solutions (𝐕) at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.

Figure 4 illustrates the optimally controlled solutions for Test Case 1 at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms. As shown in the figures, the excitation wavefront is successfully dampened out by the optimal extracellular current 𝐈^opt_𝐞. Note that the excitation wavefront is almost completely dampened out at time 1.5 ms. Next, Figure 5 depicts the optimally controlled solutions for Test Case 2 at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms. Again, the excitation wavefront is successfully dampened out by 𝐈^opt_𝐞. Observe that at time 1.5 ms, the excitation wavefront is still not completely dampened out. Thus, by comparison to Test Case 1, longer time is needed for Test Case 2 to dampen out the excitation wavefront.

(a)

(b)

(c)

(d)

Figure 4 

The controlled solutions (𝐕opt) for Test Case 1 at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.

(a)

(b)

(c)

(d)

Figure 5 

The controlled solutions (𝐕opt) for Test Case 2 at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.

Figure 6 displays the optimal extracellular current 𝐈^opt_𝐞 at three fixed points, that is, (𝑥₁,𝑦₁)=(0.5625,0.5), (𝑥₂,𝑦₂)=(0.5625,0.5156), and (𝑥₃,𝑦₃)=(0.5625,0.5313) over the time evolution. As shown in the figures, the highest extracellular current required by Test Case 1 in the defibrillation process is around 1200mAcm⁻³. On the other hand, the highest extracellular current required by Test Case 2 in the defibrillation process is twice as large as Test Case 1, that is, approximately 3000mAcm⁻³. This result implies that higher extracellular current is necessary if the control domain is relocated to a position which is far away from the excitation domain.

(a)

(b)

Figure 6 

Time evolution of the optimal extracellular current 𝐈opt𝐞 at three fixed points for (a) Test Case 1 and (b) Test Case 2.

6. Conclusions

This work studies the effects of the control domain position on the optimal control problem of cardiac arrhythmia by using two test cases. The control domain for Test Case 1 is close to the excitation domain, while the control domain for Test Case 2 is a bit further away from it. Our numerical results show that when the control domain is relocated to a position which is far from the excitation domain (as demonstrated by Test Case 2), more iteration is needed for the MDY method to converge to the optimal solution. Numerical results also show that longer time is needed by Test Case 2 for dampening out the excitation wavefront. Moreover, if the control domain is far from the excitation domain, higher extracellular current is required during the defibrillation process in order to dampen out the excitation wavefront. Thus, if we know which part of the heart of the patient is frequently suffering from cardiac arrhythmia, the electrodes of the implantable cardioverter defibrillator (ICD) should be implanted close to that part for improving the performance of the defibrillation process.

Acknowledgments

The work is financed by Zamalah Scholarship and Research University Grant (GUP) Vote 06J26 provided by Universiti Teknologi Malaysia and the Ministry of Higher Education of Malaysia.