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

Authors: Kin Wei Ng; Ahmad Rohanin
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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 Section2, 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

LetT be the final simulation time, Ω⊂ℜ2 the computational domain with Lipschitz boundary ∂Ω, Ωc⊂Ω the control domain, and Ωo⊂Ω the observation domain. We further set H=Ω×[0,T] and ∂H=∂Ω×[0,T]. The optimal control problem of monodomain model with Rogers-modified FitzHugh-Nagumo ion kinetic can be summarized as follows:
(2.1)minJ(V,w,Ie)=12∫0T(∫Ωo|V|2dΩo+α∫Ωc|Ie|2dΩc)dt,(2.2)s.t.λ1+λ∇·(Di∇V)-βCm∂V∂t-βIion(V,w)-11+λIe=0,inH,∂w∂t-f(V,w)=0,inH,(2.3)(Di∇V)·η=0,on∂H,V(x,0)=V0,w(x,0)=w0,onΩ,
where
(2.4)Iion(V,w)=c1V(1-VVth)(1-VVp)+c2wV,f(V,w)=c3(VVp-c4w).
Here, α is the regularization parameter, η is the outer normal to Ω, Di 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, Cm is the membrane capacitance, Vth is the threshold potential, Vp is the plateau potential, c1,c2,c3,andc4 are positive parameters, V(x,t) is the transmembrane potential, w(x,t) is the ionic current variable, Ie(x,t) is the extracellular current density stimulus, Iion(V,w) is the current density flowing through the ionic channels, and f(V,w) is the prescribed vector-value function.Note that (2.1) is the cost functional, in which the transmembrane potential V and the ionic current variable w act as the state variables and the extracellular current Ie acts as the control variable. Also, the control variable Ie is restricted to be nontrivial only on the control domain Ωc. 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. Letp(x,t) and q(x,t) denote the adjoint variables used to enforce the constraints. Define the Lagrange functional ℒ as
(3.1)ℒ=12∫0T(∫Ωo|V|2dΩo+α∫Ωc|Ie|2dΩc)dt+∫0T∫Ω(λ1+λ∇·(Di∇V)-βCm∂V∂t-βIion(V,w)-11+λIe)pdΩdt+∫0T∫Ω(∂w∂t-f(V,w))qdΩdt.
The first-order optimality system for this optimal control problem is given by:
(3.2)ℒp:λ1+λ∇·(Di∇V)-βCm∂V∂t-βIion(V,w)-11+λIe=0,ℒq:∂w∂t-f(V,w)=0,(3.3)ℒV:V|o+λ1+λ∇·(Di∇p)+βCm∂p∂t-β[Iion]Vp-[f]Vq=0,ℒw:-β[Iion]wp-∂q∂t-[f]wq=0,(3.4)ℒIe:αIe-11+λp=0,
where [·]* denotes the partial derivative with respect to * and V|o denotes the transmembrane potential in Ωo. 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:
(3.5)(Di∇p)·η=0,on∂H,p(x,T)=0,q(x,T)=0,onΩ.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
(3.6)βCm∂V∂t=λ1+λ∇·(Di∇V)-βIion(V,w)-11+λIe,inH,∂w∂t=f(V,w),inH,(Di∇V)·η=0,on∂H,V(x,0)=V0,w(x,0)=w0,onΩ.
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
(3.7)βCm∂p∂t=-λ1+λ∇·(Di∇p)+β[Iion]Vp+[f]Vq-V|o,inH,∂q∂t=-β[Iion]wp-[f]wq,inH,(Di∇p)·η=0,on∂H,p(x,T)=0,q(x,T)=0,onΩ.Kunisch and Wagner [17] proved that the control-to-state mapping is well defined for the optimal control problem of monodomain model, that is, C∋Ie↦(V(Ie),w(Ie)). Thus, the cost functional in (2.1) may be rewritten as
(3.8)minJ^(Ie)=J(V(Ie),w(Ie),Ie)=12∫0T(∫Ωo|V(Ie)|2dΩo+α∫Ωc|Ie|2dΩc)dt.
From the previous analysis, one can shows that the gradient of J^(Ie) is given by
(3.9)∇J^(Ie)=αIe-11+λp.
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
(3.10)linearPDE:βCm∂V∂t=λ1+λ∇·(Di∇V),inH,nonlinearODE:βCm∂V∂t=-βIion(V,w)-11+λIe,inH,nonlinearODE:∂w∂t=f(V,w),inH,(Di∇V)·η=0,on∂H,V(x,0)=V0,onΩ,w(x,0)=w0,onΩ,
and the adjoint system becomes
(3.11)linearPDE:βCm∂p∂t=-λ1+λ∇·(Di∇p),inH,nonlinearODE:βCm∂p∂t=β[Iion]Vp+[f]Vq-V|o,inH,nonlinearODE:∂q∂t=-β[Iion]wp-[f]wq,inH,(Di∇p)·η=0,on∂H,p(x,T)=0,onΩ,q(x,T)=0,onΩ.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
(3.12)(βCmM+Δt12[λ1+λ]K)Vn+1=(βCmM-Δt12[λ1+λ]K)Vn,Vn+1=Vn+Δt2(-IionnCm-IenβCm(1+λ)),wn+1=wn+Δt2fn,V(x,0)=V0,w(x,0)=w0,
and the discretized adjoint system is given by
(3.13)(βCmM+Δt12[λ1+λ]K)pn=(βCmM-Δt12[λ1+λ]K)pn+1,pn=pn+1+Δt2(Vn+1|oβCm-[Iionn+1]VCmpn+1-[fn+1]VβCmqn+1),qn=qn+1+Δt2(β[Iionn+1]wpn+1+[fn+1]wqn+1),p(x,T)=0,q(x,T)=0.
Here, n is the time index, Δt1 and Δt2 are the local time steps, M is the mass matrix, K is the stiffness matrix, and V, V|o, p, q,w,  f, Iion, and  Ie 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 Ie0, the control variable is updated using the following recurrence:
(4.1)Iek+1=Iek+δkdk,k=0,1,…,
where δk>0 is the step length computed by the Armijo line search and dk is the search direction generated by
(4.2)dk={-∇J^(Iek),ifk=0,-∇J^(Iek)+θMDYkdk-1,ifk>0,
where θMDYk∈ℜ is the conjugate gradient update parameter given by
(4.3)θMDYk=∥∇J^(Iek)∥2(dk-1)T(∇J^(Iek)-∇J^(Iek-1)+γkδk-1dk-1),
where
(4.4)γk=10-3∥∇J^(Iek)∥2+max{0,-(δk-1dk-1)T(∇J^(Iek)-∇J^(Iek-1))∥δk-1dk-1∥2}.Recall that the step lengthδk is computed by the Armijo line search. Given an initial step length δ>0 and μ,ρ∈(0,1), the Armijo line search choose, δk=max{δ,δρ,δρ2,…} such that
(4.5)J^(Iek+δkdk)≤J^(Iek)+μδk∇J^(Iek)Tdk.
Consequently, the optimization algorithm for the MDY method can be summarized as follows.


### 4.1. Optimization Algorithm

Step 1.
Provide an initial guessIe0 and set k=0.Step 2.
Solve the discretized state system in (3.12) for time interval [0,T] as follows.(a)
Solve the linear PDE fort∈[t,t+Δt/2](4.6)(βCmM+Δt12[λ1+λ]K)Vn+1=(βCmM-Δt12[λ1+λ]K)Vn.(b)
Solve the nonlinear ODEs fort∈[t,t+Δt](4.7)Vn+1=Vn+Δt2(-IionnCm-IenβCm(1+λ)),wn+1=wn+Δt2fn.(c)
Solve the linear PDE fort∈[t+Δt/2,t+Δt](4.8)(βCmM+Δt12[λ1+λ]K)Vn+1=(βCmM-Δt12[λ1+λ]K)Vn,whereΔt is the global time step.Step 3.
Evaluate the reduced cost functionalJ^(Iek).Step 4.
Solve the discretized adjoint system in (3.13) for time interval [T,0] as follows.(a)
Solve the linear PDE fort∈[T,T-Δt/2](4.9)(βCmM+Δt12[λ1+λ]K)pn=(βCmM-Δt12[λ1+λ]K)pn+1.(b)
Solve the nonlinear ODEs fort∈[T,T-Δt](4.10)pn=pn+1+Δt2(Vn+1|oβCm-[Iionn+1]VCmpn+1-[fn+1]VβCmqn+1)qn=qn+1+Δt2(β[Iionn+1]wpn+1+[fn+1]wqn+1).(c)
Solve the linear PDE fort∈[T-Δt/2,T-Δt](4.11)(βCmM+Δt12[λ1+λ]K)pn=(βCmM-Δt12[λ1+λ]K)pn+1.Step 5.
Update the reduced gradient∇J^(Iek) using (3.9).Step 6.
Fork≥1, check the following stopping criteria:
(4.12)|J^(Iek)-J^(Iek-1)|≤10-3,∥∇J^(Iek)∥≤10-3(1+|J^(Iek)|).
If one of them is met, stop.Step 7.
Compute the conjugate gradient update parameterθMDYk using (4.3).Step 8.
Compute the search directiondk using (4.2).Step 9.
Compute the step lengthδk that satisfies the condition in (4.5).Step 10.
Update the control variableIek+1 using (4.1). Set k=k+1 and go to Step 2.


## 4.1. Optimization Algorithm

Step 1.
Provide an initial guessIe0 and set k=0.Step 2.
Solve the discretized state system in (3.12) for time interval [0,T] as follows.(a)
Solve the linear PDE fort∈[t,t+Δt/2](4.6)(βCmM+Δt12[λ1+λ]K)Vn+1=(βCmM-Δt12[λ1+λ]K)Vn.(b)
Solve the nonlinear ODEs fort∈[t,t+Δt](4.7)Vn+1=Vn+Δt2(-IionnCm-IenβCm(1+λ)),wn+1=wn+Δt2fn.(c)
Solve the linear PDE fort∈[t+Δt/2,t+Δt](4.8)(βCmM+Δt12[λ1+λ]K)Vn+1=(βCmM-Δt12[λ1+λ]K)Vn,whereΔt is the global time step.Step 3.
Evaluate the reduced cost functionalJ^(Iek).Step 4.
Solve the discretized adjoint system in (3.13) for time interval [T,0] as follows.(a)
Solve the linear PDE fort∈[T,T-Δt/2](4.9)(βCmM+Δt12[λ1+λ]K)pn=(βCmM-Δt12[λ1+λ]K)pn+1.(b)
Solve the nonlinear ODEs fort∈[T,T-Δt](4.10)pn=pn+1+Δt2(Vn+1|oβCm-[Iionn+1]VCmpn+1-[fn+1]VβCmqn+1)qn=qn+1+Δt2(β[Iionn+1]wpn+1+[fn+1]wqn+1).(c)
Solve the linear PDE fort∈[T-Δt/2,T-Δt](4.11)(βCmM+Δt12[λ1+λ]K)pn=(βCmM-Δt12[λ1+λ]K)pn+1.Step 5.
Update the reduced gradient∇J^(Iek) using (3.9).Step 6.
Fork≥1, check the following stopping criteria:
(4.12)|J^(Iek)-J^(Iek-1)|≤10-3,∥∇J^(Iek)∥≤10-3(1+|J^(Iek)|).
If one of them is met, stop.Step 7.
Compute the conjugate gradient update parameterθMDYk using (4.3).Step 8.
Compute the search directiondk using (4.2).Step 9.
Compute the step lengthδk that satisfies the condition in (4.5).Step 10.
Update the control variableIek+1 using (4.1). Set k=k+1 and go to Step 2.


## 5. Numerical Experiments

### 5.1. Experiment Setup

A two-dimensional computational domain of size1×1 cm2, that is, Ω=[0,1]×[0,1] is considered in our numerical experiments, and the final simulation time is set to be T=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, Ωc1 and Ωc2 are the control domains, Ω~c1 and Ω~c2 are the neighborhoods of the control domains, Ωo=Ω∖(Ω~c1∪Ω~c2) is the observation domain, and Ωexi⊂Ωo is the excitation domain where the cardiac arrhythmia first occurs.Computational domainΩ and its subdomains for (a) Test Case 1 and (b) Test Case 2.
(a)
(b)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.Table1 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 V, w, Ie are given as
(5.1)V(x,0)=V0={105mV,x∈Ωexi,0mV,otherwise,(5.2)w(x,0)=w0=0,x∈Ω,(5.3)Ie(x,0)=Ie0={0mAcm-3,x∈Ωc,0mAcm-3,otherwise.
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.
Parameter
Value
Units
β
1
0
3
c
m
-
1
C
m
1
0
-
3
mF c
m
-
2
D
i
l
3
×
1
0
-
3
S c
m
-
1
D
i
t
3.1525
×
1
0
-
4
S c
m
-
1
V
th
1.3
×
1
0
1
mV
V
p
1
0
2
mV
c
1
1.5
mS c
m
-
2
c
2
4.4
mS c
m
-
2
c
3
1.2
×
1
0
-
2
m
s
-
1
c
4
1
Dimensionless
α
1
0
-
4
Dimensionless
λ
7.062
×
1
0
-
1
Dimensionless
δ
1
Dimensionless
μ
1
0
-
4
Dimensionless
ρ
1
0
-
1
Dimensionless
Δ
t
0.1
ms
Δ
t
1
2.5
×
1
0
-
2
ms
Δ
t
2
1
0
-
3
ms


### 5.2. Experiment Results

The minimum values of the reduced cost functionalJ^(Iek) and the corresponding norm of reduced gradient ∥∇J^(Iek)∥ 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) Minimum values of reduced cost functionalJ^ and (b) norm of reduced gradient ∥∇J^∥ for Test Case 1 and Test Case 2.
(a)
(b)As illustrated in Figure2, 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.Figure3 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, Ie=0 for x∈Ω. 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 (Ie) is applied to the control domain, the excitation wavefront will continue to spread to the entire computational domain.The uncontrolled solutions(V) at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.
(a)
(b)
(c)
(d)Figure4 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 Ieopt. 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 Ieopt. 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.The controlled solutions(Vopt) 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)The controlled solutions(Vopt) for Test Case 2 at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.
(a)
(b)
(c)
(d)Figure6 displays the optimal extracellular current Ieopt  at three fixed points, that is, (x1,y1)=(0.5625,0.5), (x2,y2)=(0.5625,0.5156), and (x3,y3)=(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-3. 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-3. 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.Time evolution of the optimal extracellular currentIeopt at three fixed points for (a) Test Case 1 and (b) Test Case 2.
(a)
(b)


## 5.1. Experiment Setup

A two-dimensional computational domain of size1×1 cm2, that is, Ω=[0,1]×[0,1] is considered in our numerical experiments, and the final simulation time is set to be T=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, Ωc1 and Ωc2 are the control domains, Ω~c1 and Ω~c2 are the neighborhoods of the control domains, Ωo=Ω∖(Ω~c1∪Ω~c2) is the observation domain, and Ωexi⊂Ωo is the excitation domain where the cardiac arrhythmia first occurs.Computational domainΩ and its subdomains for (a) Test Case 1 and (b) Test Case 2.
(a)
(b)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.Table1 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 V, w, Ie are given as
(5.1)V(x,0)=V0={105mV,x∈Ωexi,0mV,otherwise,(5.2)w(x,0)=w0=0,x∈Ω,(5.3)Ie(x,0)=Ie0={0mAcm-3,x∈Ωc,0mAcm-3,otherwise.
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.
Parameter
Value
Units
β
1
0
3
c
m
-
1
C
m
1
0
-
3
mF c
m
-
2
D
i
l
3
×
1
0
-
3
S c
m
-
1
D
i
t
3.1525
×
1
0
-
4
S c
m
-
1
V
th
1.3
×
1
0
1
mV
V
p
1
0
2
mV
c
1
1.5
mS c
m
-
2
c
2
4.4
mS c
m
-
2
c
3
1.2
×
1
0
-
2
m
s
-
1
c
4
1
Dimensionless
α
1
0
-
4
Dimensionless
λ
7.062
×
1
0
-
1
Dimensionless
δ
1
Dimensionless
μ
1
0
-
4
Dimensionless
ρ
1
0
-
1
Dimensionless
Δ
t
0.1
ms
Δ
t
1
2.5
×
1
0
-
2
ms
Δ
t
2
1
0
-
3
ms


## 5.2. Experiment Results

The minimum values of the reduced cost functionalJ^(Iek) and the corresponding norm of reduced gradient ∥∇J^(Iek)∥ 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) Minimum values of reduced cost functionalJ^ and (b) norm of reduced gradient ∥∇J^∥ for Test Case 1 and Test Case 2.
(a)
(b)As illustrated in Figure2, 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.Figure3 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, Ie=0 for x∈Ω. 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 (Ie) is applied to the control domain, the excitation wavefront will continue to spread to the entire computational domain.The uncontrolled solutions(V) at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.
(a)
(b)
(c)
(d)Figure4 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 Ieopt. 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 Ieopt. 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.The controlled solutions(Vopt) 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)The controlled solutions(Vopt) for Test Case 2 at (a) 0.2 ms, (b) 0.8 ms, (c) 1.5 ms, and (d) 2.0 ms.
(a)
(b)
(c)
(d)Figure6 displays the optimal extracellular current Ieopt  at three fixed points, that is, (x1,y1)=(0.5625,0.5), (x2,y2)=(0.5625,0.5156), and (x3,y3)=(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-3. 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-3. 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.Time evolution of the optimal extracellular currentIeopt at three fixed points for (a) Test Case 1 and (b) Test Case 2.
(a)
(b)


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



---

*Source: 468189-2012-12-19.xml*