Regularization of Ill-Posed Point Neuron Models.

Point neuron models with a Heaviside firing rate function can be ill-posed. That is, the initial-condition-to-solution map might become discontinuous in finite time. If a Lipschitz continuous but steep firing rate function is employed, then standard ODE theory implies that such models are well-posed and can thus, approximately, be solved with finite precision arithmetic. We investigate whether the solution of this well-posed model converges to a solution of the ill-posed limit problem as the steepness parameter of the firing rate function tends to infinity. Our argument employs the Arzelà-Ascoli theorem and also yields the existence of a solution of the limit problem. However, we only obtain convergence of a subsequence of the regularized solutions. This is consistent with the fact that models with a Heaviside firing rate function can have several solutions, as we show. Our analysis assumes that the vector-valued limit function v, provided by the Arzelà-Ascoli theorem, is threshold simple: That is, the set containing the times when one or more of the component functions of v equal the threshold value for firing, has zero Lebesgue measure. If this assumption does not hold, we argue that the regularized solutions may not converge to a solution of the limit problem with a Heaviside firing function.

Ill-Posed Point Neuron Models.

We show that point-neuron models with a Heaviside firing rate function can be ill posed. More specifically, the initial-condition-to-solution map might become discontinuous in finite time. Consequently, if finite precision arithmetic is used, then it is virtually impossible to guarantee the accurate numerical solution of such models. If a smooth firing rate function is employed, then standard ODE theory implies that point-neuron models are well posed. Nevertheless, in the steep firing rate regime, the problem may become close to ill posed, and the error amplification, in finite time, can be very large. This observation is illuminated by numerical experiments. We conclude that, if a steep firing rate function is employed, then minor round-off errors can have a devastating effect on simulations, unless proper error-control schemes are used.

Computing ischemic regions in the heart with the bidomain model--first steps towards validation.

We investigate whether it is possible to use the bidomain model and body surface potential maps (BSPMs) to compute the size and position of ischemic regions in the human heart. This leads to a severely ill posed inverse problem for a potential equation. We do not use the classical inverse problems of electrocardiography, in which the unknown sources are the epicardial potential distribution or the activation sequence. Instead we employ the bidomain theory to obtain a model that also enables identification of ischemic regions transmurally. This approach makes it possible to distinguish between subendocardial and transmural cases, only using the BSPM data. The main focus is on testing a previously published algorithm on clinical data, and the results are compared with images taken with perfusion scintigraphy. For the four patients involved in this study, the two modalities produce results that are rather similar: The relative differences between the center of mass and the size of the ischemic regions, suggested by the two modalities, are 10.8% ± 4.4% and 7.1% ± 4.6%, respectively. We also present some simulations which indicate that the methodology is robust with respect to uncertainties in important model parameters. However, in contrast to what has been observed in investigations only involving synthetic data, inequality constraints are needed to obtain sound results.

Towards a computational method for imaging the extracellular potassium concentration during regional ischemia.

We investigate the possibility of using body surface potential maps to image the extracellular potassium concentration during regional ischemia. The problem is formulated as an inverse problem based on a linear approximation of the bidomain model, where we minimize the difference between the results of the model and observations of body surface potentials. The minimization problem is solved by a one-shot technique, where the original PDE system, an adjoint problem, and the relation describing the minimum, are solved simultaneously. This formulation of the problem requires the solution of a 5 x 5 system of linear partial differential equations. The performance of the model is investigated by performing tests based on synthetic data. We find that the model will in many cases detect the correct position and approximate size of the ischemic regions, while some cases are more difficult to locate. It is observed that a simple post-processing of the results produces images that are qualitatively very similar to the true solution.

A computationally efficient method for determining the size and location of myocardial ischemia.

The purpose of this paper is to introduce a new method for solving the inverse problem of locating ischemic regions in the heart. The electrical activity in the human heart is modeled by the bidomain equations, which can be expanded to compute the potentials on the body surface. The associated inverse problem is to use ECG recordings to gain information about ischemias. We propose an algorithm for doing this, combining the level set method with a simpler minimization problem. Instead of trying to determine the shape, as in the level set method, we simply make the approximation that the ischemia is spherical. The effects of ischemia on the electrical attributes of heart tissue are examined. The new method is tested with computer simulations on synthetic body surface potential maps (BSPMs) in a realistic geometry, using realistic values for the parameters. It is shown to be, in some respects, superior to the level set approach and may be a step toward a practical algorithm useful in medical diagnostics.

On the possibility for computing the transmembrane potential in the heart with a one shot method: an inverse problem.

We analyze the possibility for using body surface potential maps (BSPMs), a priori information about the voltage distribution in the heart and the bidomain equations to compute the transmembrane potential throughout the myocardium. Our approach is defined in terms of an inverse problem for elliptic partial differential equations (PDEs). More precisely, we formulate it in terms of an output least squares framework in which a goal functional is minimized subject to suitable PDE constraints. The problem is highly unstable and, even under optimal recording conditions, it does not have a unique solution. We propose a methodology for stabilizing and enforcing uniqueness for this inverse problem. Moreover, a fully implicit method for solving the involved minimization problem is presented. In other words, we show how one may solve it in terms of a system consisting of three linear elliptic PDEs, i.e. we derive a so-called one shot method (also commonly referred to as an all-at-once method). Finally, our theoretical findings are illuminated by a series of numerical experiments. These examples indicate that, in the presence of regional ischemia, it might be possible to approximately recover the transmembrane potential during the resting and plateau phases of the heart cycle. This is probably due to the fact that rather accurate a priori information is available during these time intervals. The problem of computing the transmembrane potential at an arbitrary time instance during a heart beat is still an open problem.

A linear system of partial differential equations modeling the resting potential of a heart with regional ischemia.

Ischemic ST-segment shift has been modeled using scalar stationary approximations of the bidomain model. Here, we propose an alternative simplification of the bidomain equations: a linear system modeling the resting potential, to be used in determining ischemic TP shift. Results of 2D and 3D simulations show that the linear system model is much more accurate than the scalar model. This improved accuracy is important if the model is to be used for solving the inverse problem of determining the size and location of an ischemic region. Furthermore, the model can provide insight into how the resting potential depends on the variations in the extracellular potassium concentration that characterize ischemic regions.

Computing the size and location of myocardial ischemia using measurements of ST-segment shift.

It is well known that the presence of myocardial ischemia can be observed as a shift in the ST segment of an electrocardiogram (ECG) recording. The question we address in this paper is whether or not ST shift can be used to compute approximations of the size and location of the ischemic region. We begin by investigating a cost functional (measuring the difference between synthetic recorded data and simulated values of ST shift) for a parameter identification problem to locate the ischemic region. We then formulate a more flexible representation of the ischemia using a level set framework and solve the associated minimization problem for the size and position of the ischemia. We apply this framework to a set of ECG data generated by the Bidomain model using the cell model of Winslow et al. Based on this data, we show that values of ST shift recorded at the body surface are capable of identifying the position and (roughly) the size of the ischemia.

On the computational complexity of the bidomain and the monodomain models of electrophysiology.

The bidomain model, coupled with accurate models of cell membrane kinetics, is generally believed to provide a reasonable basis for numerical simulations of cardiac electrophysiology. Because of changes occurring in very short time intervals and over small spatial domains, discretized versions of these models must be solved on fine computational grids, and small time-steps must be applied. This leads to huge computational challenges that have been addressed by several authors. One popular way of reducing the CPU demands is to approximate the bidomain model by the monodomain model, and thus reducing a two by two set of partial differential equations to one scalar partial differential equation; both of which are coupled to a set of ordinary differential equations modeling the cell membrane kinetics. A reduction in CPU time of two orders of magnitude has been reported. It is the purpose of the present paper to provide arguments that such a reduction is not present when order-optimal numerical methods are applied. Theoretical considerations and numerical experiments indicate that the reduction factor of the CPU requirements from bidomain to monodomain computations, using order-optimal methods, typically is about 10 for simple cell models and less than two for more complex cell models.