Monodomain Research Articles

Moreau-Yosida regularization to optimal control of the monodomain model with pointwise control and state constraints in cardiac electrophysiology

In this work, we study the optimal control problem of a coupled reaction-diffusion system, which is a monodomain model in cardiac electrophysiology with pointwise bilateral control and state constraints. We adopt the Moreau-Yosida regularization as a penalization technique to deal with the state constraints. The regularized problem’s first-order optimality condition is derived. In addition, sufficient second-order optimality condition is derived for the regularized problem using the virtual control concept by proving equivalence between Moreau-Yosida regularization and the virtual control concept. The convergence of optimal controls of the regularized problems to the optimal control of the original problem is proved. Moreover, the semi-smooth Newton method for numerically finding the optimal solution to the regularization problem is presented. Finally, numerical experiments are conducted, and the results allow us to understand the extinction of the wave excitation in cardiac defibrillation in the presence of both control and state constraints.

Evaluating computational efforts and physiological resolution of mathematical models of cardiac tissue

Computational techniques have significantly advanced our understanding of cardiac electrophysiology, yet they have predominantly concentrated on averaged models that do not represent the intricate dynamics near individual cardiomyocytes. Recently, accurate models representing individual cells have gained popularity, enabling analysis of the electrophysiology at the micrometer level. Here, we evaluate five mathematical models to determine their computational efficiency and physiological fidelity. Our findings reveal that cell-based models introduced in recent literature offer both efficiency and precision for simulating small tissue samples (comprising thousands of cardiomyocytes). Conversely, the traditional bidomain model and its simplified counterpart, the monodomain model, are more appropriate for larger tissue masses (encompassing millions to billions of cardiomyocytes). For simulations requiring detailed parameter variations along individual cell membranes, the EMI model emerges as the only viable choice. This model distinctively accounts for the extracellular (E), membrane (M), and intracellular (I) spaces, providing a comprehensive framework for detailed studies. Nonetheless, the EMI model’s applicability to large-scale tissues is limited by its substantial computational demands for subcellular resolution.

Unsupervised stochastic learning and reduced order modeling for global sensitivity analysis in cardiac electrophysiology models

Background and Objective:Numerical simulations in electrocardiology are often affected by various uncertainties inherited from the lack of precise knowledge regarding input values including those related to the cardiac cell model, domain geometry, and boundary or initial conditions used in the mathematical modeling. Conventional techniques for uncertainty quantification in modeling electrical activities of the heart encounter significant challenges, primarily due to the high computational costs associated with fine temporal and spatial scales. Additionally, the need for numerous model evaluations to quantify ubiquitous uncertainties increases the computational challenges even further. Methods:In the present study, we propose a non-intrusive surrogate model to perform uncertainty quantification and global sensitivity analysis in cardiac electrophysiology models. The proposed method combines an unsupervised machine learning technique with the polynomial chaos expansion to reconstruct a surrogate model for the propagation and quantification of uncertainties in the electrical activity of the heart. The proposed methodology not only accurately quantifies uncertainties at a very low computational cost but more importantly, it captures the targeted quantity of interest as either the whole spatial field or the whole temporal period. In order to perform sensitivity analysis, aggregated Sobol indices are estimated directly from the spectral mode of the polynomial chaos expansion. Results:We conduct Uncertainty Quantification (UQ) and global Sensitivity Analysis (SA) considering both spatial and temporal variations, rather than limiting the analysis to specific Quantities of Interest (QoIs). To assess the comprehensive performance of our methodology in simulating cardiac electrical activity, we utilize the monodomain model. Additionally, sensitivity analysis is performed on the parameters of the Mitchell-Schaeffer cell model. Conclusions:Unlike conventional techniques for uncertainty quantification in modeling electrical activities, the proposed methodology performs at a low computational cost the sensitivity analysis on the cardiac electrical activity parameters. The results are fully reproducible and easily accessible, while the proposed reduced-order model represents a significant contribution to enhancing global sensitivity analysis in cardiac electrophysiology.

Optimality Conditions for Optimal Control of the Monodomain Model with Pointwise Control and State Constraints

Optimality Conditions for Optimal Control of the Monodomain Model with Pointwise Control and State Constraints

Explicit stabilized multirate methods for the monodomain model in cardiac electrophysiology

Fully explicit stabilized multirate (mRKC) methods are well-suited for the numerical solution of large multiscale systems of stiff ordinary differential equations thanks to their improved stability properties. To demonstrate their efficiency for the numerical solution of stiff, multiscale, nonlinear parabolic PDE's, we apply mRKC methods to the monodomain equation from cardiac electrophysiology. In doing so, we propose an improved version, specifically tailored to the monodomain model, which leads to the explicit exponential multirate stabilized (emRKC) method. Several numerical experiments are conducted to evaluate the efficiency of both mRKC and emRKC, while taking into account different finite element meshes (structured and unstructured) and realistic ionic models. The new emRKC method typically outperforms a standard implicit-explicit baseline method for cardiac electrophysiology. Code profiling and strong scalability results further demonstrate that emRKC is faster and inherently parallel without sacrificing accuracy.

A fast time integral finite difference method for a space-time fractional FitzHugh-Nagumo monodomain model in irregular domains

This work aims at proposing a fast time integral (FTI) method for the space-time fractional Fitzhugh-Nagumo (FHN) monodomain model in irregular domains, which is commonly used to characterize the transmembrane potential of the heart. In order to reduce the storage and the algorithm complexity due to the geometric configuration of the heart, the derivation of the sum of the exponentials (cf. [1]) in FTI method is improved by the integral transformation, integral truncation and Gauss-Legendre quadrature. Such strategy is applied to approximate the Caputo fractional derivative. The number of the exponentials, reduced in the FTI method, is related to the calculation efficiency. The CPU time of the FHN monodomain model using FTI method is reduced by an order of magnitude. Moreover, a second-order discrete method to the Riesz space fractional derivative adopts for the spatial discretization, then an implicit-explicit scheme is derived for the nonlinear FHN model under the finite difference method. Numerical results are reported to demonstrate the convergence behaviors, robustness and high efficiency of the proposed method.

Cardiac Electrophysiology Meshfree Modeling through the Mixed Collocation Method

We present the meshfree mixed collocation method (MCM) for cardiac electrophysiology simulation. Capitalizing on the meshfree property of MCM, we introduce an immersed grid approach for automated generation of meshfree node grids from medical image data. This approach allows us to avoid the time-consuming mesh generation and processing that mesh-based methods like the finite element method (FEM) require. We employ the MCM to solve the cardiac monodomain model considering electrical propagation in 2D tissue sheets, 3D tissue slabs, and a realistic biventricular anatomy. We demonstrate that the solutions obtained by the MCM are in good agreement with the FEM, particularly when immersed grid is used. These findings confirm the suitability of the MCM for cardiac electrophysiology simulation and make the MCM a promising alternative to the FEM for cardiac electrical investigations.

The simplified Kirchhoff network model (SKNM): a cell-based reaction–diffusion model of excitable tissue

Cell-based models of excitable tissues offer the advantage of cell-level precision, which cannot be achieved using traditional homogenized electrophysiological models. However, this enhanced accuracy comes at the cost of increased computational demands, necessitating the development of efficient cell-based models. The widely-accepted bidomain model serves as the standard in computational cardiac electrophysiology, and under certain anisotropy ratio conditions, it is well known that it can be reduced to the simpler monodomain model. Recently, the Kirchhoff Network Model (KNM) was developed as a cell-based counterpart to the bidomain model. In this paper, we aim to demonstrate that KNM can be simplified using the same steps employed to derive the monodomain model from the bidomain model. We present the cell-based Simplified Kirchhoff Network Model (SKNM), which produces results closely aligned with those of KNM while requiring significantly less computational resources.

Bubble domain evolution in well-ordered BiFeO3 nanocapacitors

Ferroelectric nanocapacitors have attracted intensive research interest due to their novel functionalities and potential application in nanodevices. However, due to the lack of knowledge of domain evolution in isolated nanocapacitors, precise manipulation of topological domain switching in the nanocapacitor is still a challenge. Here, we report unique bubble and cylindrical domains in the well-ordered BiFeO3 nanocapacitor array. The transformation of bubble, cylindrical and mono domains in isolated ferroelectric nanocapacitor has been demonstrated via scanning probe microscopy (SPM). The bubble domain can be erased to mono domain or written to cylindrical domain and mono domain by positive and negative voltage, respectively. Additionally, the domain evolution rules, which are mainly affected by the depolarization field, have been observed in the nanocapacitors with different domain structures. This work will be helpful in understanding the domain evolution in ferroelectric nanocapacitors and providing guidance on the manipulation of nanoscale topological domains.

Performance of atrial conduction velocity algorithms: a comparative in-silico study

Abstract Funding Acknowledgements Type of funding sources: Foundation. Main funding source(s): British Heart Foundation Background Atrial conduction velocity is a key determinant of re-entry. Measurement of conduction velocity from electroanatomic mapping data is challenging due to sparsity, distribution, and uncertainty of measurements. Several mathematical methods have been proposed to address these challenges but their relative sensitivity to the typical uncertainty of clinical data is unknown. Furthermore, prior clinical studies utilising these algorithms provide conflicting results. Purpose To assess the performance of conduction velocity methods (triangulation, planar fitting, and radial basis function interpolation) for the calculation of atrial conduction velocity. Method Atrial activation was simulated with known conduction velocity using the monodomain model, Courtemanche cell model and atrial geometry derived from cardiac magnetic resonance imaging with mapped fibre orientations. The model was paced from the pulmonary veins, the mid coronary sinus, and the insertion of Bachmann's bundle. Reconstructed electroanatomic maps were used for conduction velocity calculation in OpenEP[1]. To test the robustness of algorithms to electrogram sampling density, simulated sampling densities were reduced from 16.45 to 8.43 points/cm2 (1950-1000 points/map). To test the robustness to uncertainty in local activation time assignment, random uniform noises with the maximum amplitude of ± 20ms was added to the local activation times from simulations. Mean squared velocity error was calculated using the differences between the ground truth and estimated conduction velocities. Result In electroanatomic maps recreated with high electrogram sampling density (16.45 points/cm2), planar fitting showed a slightly higher error compared to triangulation and the radial basis function methods (0.16 ± 0.03 m/s, 0.11 ± 0.01 m/s and 0.12 ± 0.03 m/s, respectively, P < 0.05). As mapping density was decreased, the relative performance of the planar fitting method increased and in low-density maps (13.5 to 8.43 points/cm2) the planar fitting method resulted in the lowest error. When local activation times were accurately assigned, all methods performed well. However, when uncertainty in local activation time assignment increased, conduction velocities calculated using planar fitting had a significantly lower error compared to those calculated with the triangulation and radial basis function methods (0.19 ± 0.03 m/s, 0.36 ± 0.04 and 0.29 ± 0.04 m/s, respectively, P < 0.05). Conclusion This study highlights the importance of spatial and temporal accuracy for the assessment of conduction velocity. At high electroanatomic mapping sampling densities (>16.9 points/cm2), all methods performed well. Below 13.5 points/cm2 (~1600 points per map), the planar fitting method outperformed the triangulation and radial basis function methods. The differential effects of clinical measurement uncertainty on the accuracy of each method may explain the conflicting results of prior studies.

Microscopic tridomain model of electrical activity in the heart with dynamical gap junctions. Part 2 – Derivation of the macroscopic tridomain model by unfolding homogenization method

We study the homogenization of a novel microscopic tridomain system, allowing for a more detailed analysis of the properties of cardiac conduction than the classical bidomain and monodomain models. In (Acta Appl.Math. 179 (2022) 1–35), we detail this model in which gap junctions are considered as the connections between adjacent cells in cardiac muscle and could serve as alternative or supporting pathways for cell-to-cell electrical signal propagation. Departing from this microscopic cellular model, we apply the periodic unfolding method to derive the macroscopic tridomain model. Several difficulties prevent the application of unfolding homogenization results, including the degenerate temporal structure of the tridomain equations and a nonlinear dynamic boundary condition on the cellular membrane. To prove the convergence of the nonlinear terms, especially those defined on the microscopic interface, we use the boundary unfolding operator and a Kolmogorov–Riesz compactness’s result.

Efficient time splitting schemes for the monodomain equation in cardiac electrophysiology.

Approximating the fast dynamics of depolarization waves in the human heart described by the monodomain model is numerically challenging. Splitting methods for the PDE-ODE coupling enable the computation with very fine space and time discretizations. Here, we compare different splitting approaches regarding convergence, accuracy, and efficiency. Simulations were performed for a benchmark problem with the Beeler-Reuter cell model on a truncated ellipsoid approximating the left ventricle including a localized stimulation. For this configuration, we provide a reference solution for the transmembrane potential. We found a semi-implicit approach with state variable interpolation to be the most efficient scheme. The results are transferred to a more physiological setup using a bi-ventricular domain with a complex external stimulation pattern to evaluate the accuracy of the activation time for different resolutions in space and time.

Abstract 9947: Prospective Prediction of Left Atrial Flutter From Computer Simulation

Introduction: Atypical left atrial flutter (ALAF) is often the result of ablation treatment for AF. Many patients report a reduction in quality of life from AF to ALAF. If pharmacological rate control fails additional ablation can provide solutions. However, success rates are low due to difficulties in mapping all potential flutters accurately. Hypothesis: LGE-MRI can inform computation simulation to predict all possible ALAF pathways and increase the information available to the electrophysiologist. Methods: Informed and consented patients scheduled for ALAF ablation received an MRA and LGE-MRI. Fig1. From the MRA the left atrial wall was segmented and registered to the LGE-MRI. Scar and gray areas were estimated by threshold segmentation. The threshold was calculated with the intensity value ratio. The information was integrated into a monodomain model. Ionic currents were modeled according to Courtemanche with a reduction of some conductance in fibrotic regions. To induce ALAF S1S2 pacing was done with an S1 cycle length of 600ms and S2 decreased from 250 to 170ms. Induction was tested at 10 different locations. Results were provided to the electrophysiologist. Results: Prospective prediction of the ALAF pathway for six patients was successful and the clinically observed flutter was predicted. Fig. 2) shows the value of the prediction, electroanatomic mapping was inconclusive even unable to exclude the right side. Ablation on the right side did not alter the flutter, but ablation near the mitral valve. Conclusions: We successfully predicted clinically observed ALAF pathways. The next step is to find ablation points for all predicted pathways.

A meshless fragile points method for the solution of the monodomain model for cardiac electrophysiology simulation

Meshless methods for in silico modeling and simulation of cardiac electrophysiology are gaining more and more popularity. These methods do not require a mesh and are more suitable than the Finite Element Method (FEM) to simulate the activity of complex geometrical structures like the human heart. However, challenges such as numerical integration accuracy and time efficiency remain, which limit their applicability. Recently, the Fragile Points Method (FPM) has been introduced in the meshless methods family. It uses local, simple, polynomial, discontinuous functions to construct trial and test functions in the Galerkin weak form. This allows for accurate integration and improved efficiency while enabling the imposition of essential and natural boundary conditions as in FEM. In this work, we consider the application of FPM for cardiac electrophysiology simulation. We derive the cardiac monodomain model using the FPM formulation and we solve several benchmark problems in 2D and 3D. We show that FPM leads to solutions of accuracy and efficiency similar to FEM while alleviating the need for a mesh. Additionally, FPM demonstrates comparable convergence to FEM in the considered benchmarks.

Indentation of Anisotropic Tissue Using a Three-Dimensional Mechanical Bidomain Model

Computation-based mathematical models of tissue indentation are capable of predicting the distribution of forces and mechanical properties of soft tissues. This paper presents a three-dimensional mathematical model of anisotropic tissue indentation developed using the mechanical bidomain model. The mechanical bidomain model hypothesizes that the relative displacement between intra- and extracellular spaces triggers a force on the mechanosensitive proteins in the membrane: integrins. Some soft tissues, such as cardiac muscle, are anisotropic, a property which arises from the fibrous structure of the tissue. The degree of anisotropy in intra- and extracellular spaces can be different. Tissue indentation for different anisotropy ratios that indicate isotropy, equal anisotropy and unequal anisotropy, were tested using the model. Results of the tissue indentation analysis compared the spatial distribution of the magnitude of bidomain displacement for different anisotropy conditions between monodomain and bidomain models. The proposed mathematical model predicted unexpected spatial patterns of cardiac mechanotransduction for unequal anisotropy ratios of mechanical modulus.

Homogenisation for the monodomain model in the presence of microscopic fibrotic structures

Computational models in cardiac electrophysiology are notorious for long runtimes, restricting the numbers of nodes and mesh elements in the numerical discretisations used for their solution. This makes it particularly challenging to incorporate structural heterogeneities on small spatial scales, preventing a full understanding of the critical arrhythmogenic effects of conditions such as cardiac fibrosis. In this work, we explore the technique of homogenisation by volume averaging for the inclusion of non-conductive micro-structures into larger-scale cardiac meshes with minor computational overhead. Importantly, our approach is not restricted to periodic patterns, enabling homogenised models to represent, for example, the intricate patterns of collagen deposition present in different types of fibrosis. We first highlight the importance of appropriate boundary condition choice for the closure problems that define the parameters of homogenised models. Then, we demonstrate the technique’s ability to correctly upscale the effects of fibrotic patterns with a spatial resolution of 10µm into much larger numerical mesh sizes of 100-250µm. The homogenised models using these coarser meshes correctly predict critical pro-arrhythmic effects of fibrosis, including slowed conduction, source/sink mismatch, and stabilisation of re-entrant activation patterns. As such, this approach to homogenisation represents a significant step towards whole organ simulations that unravel the effects of microscopic cardiac tissue heterogeneities.

A fast cardiac electromechanics model coupling the Eikonal and the nonlinear mechanics equations

We present a new model of human cardiac electromechanics for the left ventricle where electrophysiology is described by a Reaction–Eikonal model and which enables an off-line resolution of the reaction model, thus entailing a big saving of computational time. Subcellular dynamics is coupled with a model of tissue mechanics, which is in turn coupled with a Windkessel model for blood circulation. Our numerical results show that the proposed model is able to provide a physiological response to changes in certain variables (end-diastolic volume, total peripheral resistance, contractility). We also show that our model is able to reproduce with high accuracy and with a considerably lower computational time the results that we would obtain if the monodomain model should be used in place of the Eikonal model.

Centres of spiral waves can be detected as phase defect lines in optical voltage mapping data and numerical simulations

Abstract Funding Acknowledgements Type of funding sources: Public grant(s) – National budget only. Main funding source(s): KU Leuven and FWO-Flanders Introduction To study the excitation patterns of spiral waves in two spatial dimensions over time, it is useful to abstract them to much simpler mathematical structures. Previously, those waves have been abstracted to phase singularities: single points around which spiral moves [1,2]. Recently, it was proposed to upgrade these points to phase defect lines [3,4]. The dynamics of either of these structures can then be studied in place of the spiral waves. Purpose To work with phase defects, firstly, we need a reliable method of locating them. Methods We recorded the transmembrane voltage over space and time in a numerical simulation of the mono-domain model as well as by using optical voltage mapping data of 10 cm² monolayers of conditionally immortalised human atrial myocytes [5]. Optical voltage mapping was performed with a sampling time of 6 ms and 100×100 square pixels of 0.25 mm width. The recorded intensity were smoothed and rescaled to get an estimate of the transmembrane voltages. Detailed protocols can be found in [5]. For the numerical simulations, we have implemented the linear-core tissue model by Bueno-Orovio, Cherry, and Fenton (2008) [6] with finite differences on a 1200×1200 rectangular grid with a grid size of 0.25 mm and a time step of 9 μs. We normalised the transmembrane voltage u such that zero corresponds to the resting state and one to the plateau of the action potential. A phase φ is an angle in radians which monotonously increases from zero to 2π during one action potential. Phases are defined in such a manner that, even for the sharp wave front of a cardiac action potential, this increase happens in a slow, gradual way. A phase can for instance be defined by using the local arrival time and the characteristic action potential duration [4]. At the centre of spiral waves, neighbouring parts of tissue have vastly different phase values φ. A phase singularity is a point where all phases meet. Phase defects are large jumps across space in phase. Phase defect lines can then be detected as lines in the phase defect field ρ with a value much larger than zero [7]. Results For the simulation, we can clearly see that the linear core of the Bueno-Orovio, Cherry, and Fenton (2008) [6] model is identified as a phase defect line of up to 100 mm in length (Fig.1). The spiral wave moves along it and turns at the end points. The same behaviour can be seen in vitro in the optical voltage mapping data (Fig.2). Multiple spirals can be observed, all of which revolve around phase defect lines of lengths of up to 4 mm. Conclusion In this study, we have shown that phase defect lines can be detected both in silico and in vitro. Their properties such as length, rotational frequency, and their dynamics can be used to gain additional insight into spiral waves in arrhythmias. We have introduced multiple methods to localise phase defects for which the code is publicly available.

Simulation of atrial fibrillation in a non-ohmic propagation model with dynamic gap junctions.

Gap junctions exhibit nonlinear electrical properties that have been hypothesized to be relevant to arrhythmogenicity in a structurally remodeled tissue. Large-scale implementation of gap junction dynamics in 3D propagation models remains challenging. We aim to quantify the impact of nonlinear diffusion during episodes of arrhythmias simulated in a left atrial model. Homogenization of conduction properties in the presence of nonlinear gap junctions was performed by generalizing a previously developed mathematical framework. A monodomain model was solved in which conductivities were time-varying and depended on transjunctional potentials. Gap junction conductances were derived from a simplified Vogel-Weingart model with first-order gating and adjustable time constant. A bilayer interconnected cable model of the left atrium with 100 μm resolution was used. The diffusion matrix was recomputed at each time step according to the state of the gap junctions. Sinus rhythm and atrial fibrillation episodes were simulated in remodeled tissue substrates. Slow conduction was induced by reduced coupling and by diffuse or stringy fibrosis. Simulations starting from the same initial conditions were repeated with linear and nonlinear gap junctions. The discrepancy in activation times between the linear and nonlinear diffusion models was quantified. The results largely validated the linear approximation for conduction velocities >20 cm/s. In very slow conduction substrates, the discrepancy accumulated over time during atrial fibrillation, eventually leading to qualitative differences in propagation patterns, while keeping the descriptive statistics, such as cycle lengths, unchanged. The discrepancy growth rate was increased by impaired conduction, fibrosis, conduction heterogeneity, lateral uncoupling, fast gap junction time constant, and steeper action potential duration restitution.

Role of Scar and Border Zone Geometry on the Genesis and Maintenance of Re-Entrant Ventricular Tachycardia in Patients With Previous Myocardial Infarction.

In patients with healed myocardial infarction, the left ventricular ejection fraction is characterized by low sensitivity and specificity in the prediction of future malignant arrhythmias. Thus, there is the need for new parameters in daily practice to perform arrhythmic risk stratification. The aim of this study is to identify some features of proarrhythmic geometric configurations of scars and border zones (BZ), by means of numerical simulations based on left ventricular models derived from post myocardial infarction patients. Two patients with similar clinical characteristics were included in this study. Both patients exhibited left ventricular scars characterized by subendo- and subepicardial BZ and a transmural BZ isthmus. The scar of patient #1 was significantly larger than that of patient #2, whereas the transmural BZ isthmus and the subdendo- and subepicardial BZs of patient #2 were thicker than those of patient #1. Patient #1 was positive at electrophysiologic testing, whereas patient #2 was negative. Based on the cardiac magnetic resonance (CMR) data, we developed a geometric model of the left ventricles of the two patients, taking into account the position, extent, and topological features of scars and BZ. The numerical simulations were based on the anisotropic monodomain model of electrocardiology. In the model of patient #1, sustained ventricular tachycardia (VT) was inducible by an S2 stimulus delivered at any of the six stimulation sites considered, while in the model of patient #2 we were not able to induce sustained VT. In the model of patient #1, making the subendo- and subepicardial BZs as thick as those of patient #2 did not affect the inducibility and maintenance of VT. On the other hand, in the model of patient #2, making the subendo- and subepicardial BZs as thin as those of patient #1 yielded sustained VT. In conclusion, the results show that the numerical simulations have an effective predictive capability in discriminating patients at high arrhythmic risk. The extent of the infarct scar and the presence of transmural BZ isthmuses and thin subendo- and subepicardial BZs promote sustained VT.