Research Article | Open Access
The Effects of Control Domain Position on Optimal Control of Cardiac Arrhythmia
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.
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.  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.  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.  for solving the optimal control problem of monodomain model, namely, the Dai-Yuan (DY) method , the Hager-Zhang (HZ) method , 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 . 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 , 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 and . The optimal control problem of monodomain model with Rogers-modified FitzHugh-Nagumo ion kinetic can be summarized as follows: where 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, is the threshold potential, is the plateau potential, are positive parameters, is the transmembrane potential, is the ionic current variable, is the extracellular current density stimulus, 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 .
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 The first-order optimality system for this optimal control problem is given by: 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:
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 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
Kunisch and Wagner  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 From the previous analysis, one can shows that the gradient of is given by 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  is applied to (3.6) and (3.7) for decomposing the complex systems into subsystems. Eventually, the state system becomes and the adjoint system becomes
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 and the discretized adjoint system is given by 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  for solving the discretized optimal control problem of monodomain model. Starting from an initial guess , the control variable is updated using the following recurrence: where is the step length computed by the Armijo line search and is the search direction generated by where is the conjugate gradient update parameter given by where
Recall that the step length is computed by the Armijo line search. Given an initial step length and , the Armijo line search choose, such that 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 .
Step 2. Solve the discretized state system in (3.12) for time interval as follows.(a) Solve the linear PDE for (b) Solve the nonlinear ODEs for (c) Solve the linear PDE for 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 as follows.(a) Solve the linear PDE for (b) Solve the nonlinear ODEs for (c) Solve the linear PDE for
Step 5. Update the reduced gradient using (3.9).
Step 6. For , check the following stopping criteria:
If one of them is met, stop.
Step 7. Compute the conjugate gradient update parameter 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).
5. Numerical Experiments
5.1. Experiment Setup
A two-dimensional computational domain of size cm2, that is, is considered in our numerical experiments, and the final simulation time is set to be 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, and are the control domains, and are the neighborhoods of the control domains, is the observation domain, and is the excitation domain where the cardiac arrhythmia first occurs.
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. . Moreover, the initial values for , , are given as 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.
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.
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, 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.
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 . 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 . 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.
Figure 6 displays the optimal extracellular current at three fixed points, that is, , , and over the time evolution. As shown in the figures, the highest extracellular current required by Test Case 1 in the defibrillation process is around . 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 . 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.
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.
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.
- A. Amann, R. Tratnig, and K. Unterkofler, “A new ventricular fibrillation detection algorithm for automated external defibrillators,” Computers in Cardiology, vol. 32, pp. 559–562, 2005.
- D. J. Dosdall, V. G. Fast, and R. E. Ideker, “Mechanisms of defibrillation,” Annual Review of Biomedical Engineering, vol. 12, pp. 233–258, 2010.
- R. C. Klein, M. H. Raitt, B. L. Wilkoff et al., “Analysis of implantable cardioverter defibrillator therapy in the antiarrhythmics versus implantable defibrillators (AVID) trial,” Journal of Cardiovascular Electrophysiology, vol. 14, no. 9, pp. 940–948, 2003.
- C. Nagaiah, K. Kunisch, and G. Plank, “Numerical solution for optimal control of the reaction-diffusion equations in cardiac electrophysiology,” Computational Optimization and Applications, vol. 49, no. 1, pp. 149–178, 2011.
- S. M. Shuaiby, M. A. Hassan, and M. El-Melegy, “Modeling and simulation of the action potential in human cardiac tissues using finite element method,” Journal of Communications and Computer Engineering, vol. 2, no. 3, pp. 21–27, 2012.
- Y. Belhamadia, A. Fortin, and Y. Bourgault, “Towards accurate numerical method for monodomain models using a realistic heart geometry,” Mathematical Biosciences, vol. 220, no. 2, pp. 89–101, 2009.
- R. H. Clayton and A. V. Panfilov, “A guide to modelling cardiac electrical activity in anatomically detailed ventricles,” Progress in Biophysics and Molecular Biology, vol. 96, no. 1–3, pp. 19–43, 2008.
- Y. H. Dai and Y. Yuan, “A nonlinear conjugate gradient method with a strong global convergence property,” SIAM Journal on Optimization, vol. 10, no. 1, pp. 177–182, 1999.
- W. W. Hager and H. Zhang, “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM Journal on Optimization, vol. 16, no. 1, pp. 170–192, 2005.
- E. Polak and G. Ribière, “Note sur la convergence de méthodes de directions conjuguées,” Revue Française d'Informatique et de Recherche Opérationnelle, vol. 16, pp. 35–43, 1969.
- B. T. Polyak, “The conjugate gradient method in extremal problems,” USSR Computational Mathematics and Mathematical Physics, vol. 9, no. 4, pp. 94–112, 1969.
- C. Nagaiah and K. Kunisch, “Higher order optimization and adaptive numerical solution for optimal control of monodomain equations in cardiac electrophysiology,” Applied Numerical Mathematics, vol. 61, no. 1, pp. 53–65, 2011.
- K. W. Ng and A. Rohanin, “Numerical solution for PDE-constrained optimization problem in cardiac electrophysiology,” International Journal of Computer Applications, vol. 44, no. 12, pp. 11–15, 2012.
- K. W. Ng and A. Rohanin, “Modified Fletcher-Reeves and Dai-Yuan conjugate gradient methods for solving optimal control problem of monodomain model,” Applied Mathematics, vol. 3, no. 8, pp. 864–872, 2012.
- K. W. Ng and A. Rohanin, “The effects of control domain size on optimal control problem of monodomain model,” International Journal of Computer Applications, vol. 47, no. 10, pp. 6–11, 2012.
- J. M. Rogers and A. D. McCulloch, “A collocation-Galerkin finite element model of cardiac action potential propagation,” IEEE Transactions on Biomedical Engineering, vol. 41, no. 8, pp. 743–757, 1994.
- K. Kunisch and M. Wagner, “Optimal control of the bidomain system (I): the monodomain approximation with the Rogers-McCulloch model,” Nonlinear Analysis: Real World Applications, vol. 13, no. 4, pp. 1525–1550, 2012.
- Z. Qu and A. Garfinkel, “An advanced algorithm for solving partial differential equation in cardiac conduction,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 9, pp. 1166–1168, 1999.
- L. Zhang, “Two modified Dai-Yuan nonlinear conjugate gradient methods,” Numerical Algorithms, vol. 50, no. 1, pp. 1–16, 2009.
- P. C. Franzone, P. Deuflhard, B. Erdmann, J. Lang, and L. F. Pavarino, “Adaptivity in space and time for reaction-diffusion systems in electrocardiology,” SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 942–962, 2006.
Copyright © 2012 Kin Wei Ng and Ahmad Rohanin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.