Exponential integrators for a Markov chain model of the fast sodium channel of cardiomyocytes

The modern Markov chain models of ionic channels in excitable membranes are numerically stiff. The popular numerical methods for these models require very small time steps to ensure stability. Our objective is to formulate and test two methods addressing this issue, so that the timestep can be chosen based on accuracy rather than stability. Both proposed methods extend Rush-Larsen technique, which was originally developed to Hogdkin-Huxley type gate models. One method, "Matrix Rush-Larsen" (MRL) uses a matrix reformulation of the Rush-Larsen scheme, where the matrix exponentials are calculated using precomputed tables of eigenvalues and eigenvectors. The other, "hybrid operator splitting" (HOS) method exploits asymptotic properties of a particular Markov chain model, allowing explicit analytical expressions for the substeps. We test both methods on the Clancy and Rudy (2002) INa Markov chain model. With precomputed tables for functions of the transmembrane voltage, both methods are comparable to the forward Euler method in accuracy and computational cost, but allow longer time steps without numerical instability. We conclude that both methods are of practical interest. MRL requires more computations than HOS, but is formulated in general terms which can be readily extended to other Markov Chain channel models, whereas the utility of HOS depends on the asymptotic properties of a particular model. The significance of the methods is that they allow a considerable speed-up of large-scale computations of cardiac excitation models by increasing the time step, while maintaining acceptable accuracy and preserving numerical stability.Comment: 9 pages, 5 figures main text + 14 pages, 1 figure appendix, as submitted in final form to IEEE TBME 2014/11/11. Copyright IEEE (2014

Mathematical and Computational Study of Markovian Models of Ion Channels in Cardiac Excitation

This thesis studies numerical methods for integrating the master equations describing Markov chain models of cardiac ion channels. Such models describe the time evolution of the probability that ion channels are in a particular state. Numerical simulations of such models are often computationally demanding because many solvers require relatively small time steps to ensure numerical stability. The aim of this project is to analyse selected Markov chains and develop more efficient and accurate solvers. We separate a Markov chain model into fast and slow time-scales based on the speed of transitions between states. Eliminating the fast transitions, we find an asymptotic reduction of zeroth-order and first-order in a small parameter describing the time-scales separation. We apply the theory to a Markov chain model of the fast sodium channel INa. We consider several variants for classifying some transitions as fast in order to find reduced systems that yield a good accuracy. However, the time step size is still restricted by numerical instabilities. We adapt the Rush-Larsen technique originally developed for gate models. Assuming that a transition matrix can be considered constant during each time step, we solve the Markov chain model analytically. The solution provides a recipe for a stable exponential solver, which we call "Matrix Rush-Larsen" (MRL). Using operator splitting we design an even more flexible "hybrid" method that combines the MRL with other solvers. The resulting improvement in stability allows a large increase in the time step size. In some models, we obtain reasonably accurate results 27 times faster using a hybrid method than with the forward Euler method, even with the maximal time step allowed by the stability constraint. Finally, we extend the cardiac simulation package BeatBox by the developed exponential solvers. We upgrade a format of "ionic" modules which describe a cardiac cell, in order to allow for a specific definition of Markov chain models. We also modify a particular integrator for ionic modules to include the MRL and the hybrid method. To test the functionality of the code, we have converted a number of cellular models into the ionic format. The documented code is available in the official BeatBox package distribution

Modelling and analysis of neurons coupled by electrical synapses

The objective of this thesis is to analyze the role of the intrinsic properties of neurons in the communication through electrical synapses. Mesencephalic trigeminal neurons constitute an excellent experimental model to study the communication between neurons, because of its easy experimental access experimental and simple to model and analyze a biological system. Among the contributions of this thesis are: the complete modeling of the sodium currents and other ionic current (and its modulation); the explanation preference subthreshold frequency transfer between neuronfor example and its coupling. Some preliminary results of this work have been presented at international conferences.morphology. However, the analysis of real neurons is limited by experimental constraints that do not allow to explore all aspects of the model. Within the context of this thesis, a mathematical model is built, based on electrophysiological recordings made by Sebastián Curti at the School of Medicine of Universidad de la República. The model consists of a set of differential equations, which can be represented by a nonlinear electrical circuit. Some of the differential equations are obtained from literature and only some minor parameters’ adjustments are made. Moreover, during the thesis we have found that more data was needed in order to explain some of the most important features of the behavior of neurons, such as the duration of the action potential. Therefore, more experimental recordings were made, allowing to refine the model. The model allows to evaluate the response of the neuron to different stimuli (currents or voltages imposed by an electrode), making possible to make new “experiments” that are not possible in a laboratory. Alternatives models are analyzed (varying ionic currents and morphology) using experimental information to validate them. Then the model is used to understand some unusual features of the communication between neurons. First, it is studied the subthreshold transfer function (i.e. without action potentials) between neurons coupled by electrical synapses. A reduced model is used and then linearized, in order to derive an analytical expression of the transfer function, whose behaviour is consistent with experimental results. Moreover, numerical simulations are performed to analyze the rol of the intrinsic properties of neurons in their synchronization. It is shown that the same properties that determine the subthreshold behavior are relevant to improve synchronization between neurons too. Finally, this thesis contributes not only with new models and answers, but with new questions, which should be studied using experimental models as well. This thesis applies several tools used for electrical engineering (frequency response of systems, cable equation, Markov chains, evolutionary algorithms, etc.) to model and analyze a biological system. Among the contributions of this thesis are: the complete modeling of the sodium currents and other ionic current (and its modulation); the explanation preference subthreshold frequency transfer between neuronfor example and its coupling. Some preliminary results of this work have been presented at international conferences.El objetivo de esta tesis es analizar el rol de las propiedades intrínsecas de las neuronas en la comunicación a través de sinapsis eléctricas. Las neuronas del nervio trigeminal del mesencéfalo constituyen un excelente modelo experimental para estudiar la comunicación entre neuronas, debido a su fácil acceso experimental y su sencilla morfología. Sin embargo, el análisis de neuronas reales está limitado por restricciones experimentales que impiden explorar todos los aspectos del modelo. En el marco de esta tesis, se construye un modelo matemático basado en registros electrofisiológicos realizados por Sebastián Curti en la Facultad de Medicina de la Universidad de la República. El modelo consiste en un sistema de ecuaciones diferenciales, que puede ser representado por un circuito eléctrico con componentes no lineales. Algunas de las ecuaciones diferenciales son obtenidas de bibliografía y se realizan algunos ajustes menores de parámetros. Por otro lado, durante la tesis evaluamos que se necesitaba más información para reproducir algunas de las características más importantes del comportamientos de las neuronas, como la duración del potencial de acción. Por eso, se debieron realizar nuevos registros experimentales, que permitieron refinar el modelo. El modelo permite evaluar la respuesta de la neurona ante diferentes estímulos (corrientes o voltajes impuestos por un electrodo), posibilitando nuevos “experimentos” que no son posibles en un laboratorio. Se analizan diversas alternativas de modelado (variando corrientes iónicas y morfología) usando información experimental para validarlos. Luego, el modelo es utilizado para entender algunas características inusuales de la comunicación entre neuronas. En primer lugar, se estudia la transferencia subumbral (i.e.: sin potenciales de acción) entre neuronas acopladas por sinapsis eléctricas. Se utiliza un modelo reducido, que es linealizado para obtener una expresión analítica de la transferencia, cuyo comportamiento es coherente con los resultados experimentales. Asimismo, se realizan simulaciones numéricas para analizar el rol en la sincronización de las propiedades intrínsecas de las neuronas. Se muestra que las mismas propiedades que determinan el comportamiento subumbral son relevantes para mejorar la sincronización entre neuronas. Finalmente, esta tesis no sólo contribuye con nuevos modelos y respuestas, sino con nuevas preguntas, que deberán ser estudiadas usando modelos experimentales también. Esta tesis hace uso de diversas herramientas utilizadas por la ingeniería eléctrica (comportamiento en frecuencia de sistemas, ecuación del cable, cadenas de Markov, algoritmos evolutivos, etc) para modelar y analizar un sistema biológico. Se realizan diversos aportes, por ejemplo: modelado completo de las corrientes de sodio, así como de la modulación de otra corriente; explicación de la preferencia en frecuencia de la transferencia subumbral entre neuronas; estudio de la sincronización en función de las propiedades de los osciladores y de su acople. Algunos resultados preliminares de este trabajo han sido presentados en congresos internacionales

A TIME-AND-SPACE PARALLELIZED ALGORITHM FOR THE CABLE EQUATION

Electrical propagation in excitable tissue, such as nerve fibers and heart muscle, is described by a nonlinear diffusion-reaction parabolic partial differential equation for the transmembrane voltage $V(x,t)$, known as the cable equation. This equation involves a highly nonlinear source term, representing the total ionic current across the membrane, governed by a Hodgkin-Huxley type ionic model, and requires the solution of a system of ordinary differential equations. Thus, the model consists of a PDE (in 1-, 2- or 3-dimensions) coupled to a system of ODEs, and it is very expensive to solve, especially in 2 and 3 dimensions. In order to solve this equation numerically, we develop an algorithm, extended from the Parareal Algorithm, to efficiently incorporate space-parallelized solvers into the framework of the Parareal algorithm, to achieve time-and-space parallelization. Numerical results and comparison of the performance of several serial, space-parallelized and time-and-space-parallelized time-stepping numerical schemes in one-dimension and in two-dimensions are also presented

Approaches for unveiling the kinetic mechanisms of voltage gated ion channels in neurons

Includes vitaIonic currents drive cellular function within all living cells to perform highly specific tasks. For excitable cells, such as muscle and neurons, voltage-gated ion channels have finely tuned kinetics that allow the transduction of Action potentials to other cells. Voltage-gated ion channels are molecular machines that open and close depending on electrical potential. Neuronal firing rates are largely determined by the overall availability of voltage-gated Na+ and K+ currents.This work describes new approaches for collecting and analyzing experimental data that can be used to streamline experiments. Ion channels are composed of multimeric complexes regulated by intracellular factors producing complex kinetics. The stochastic behavior of thousands of individual ion hannels coordinates to produce cellular activity. To describe their activity and test hypotheses about the channel, experimenters often fit Markov models to a set of experimental data. Markov models are defined by a set of states, whose transitions described by rate constants. To improve the modeling process, we have developed computational approaches to introduce kinetic constraints that reduces the parameter search space. This work describes the implementation and mathematical transformations required to describe linear and non-linear parameter constraints that govern rate constants. Not all ion channel behaviors can easily be described by rate constants. Therefore, we developed and implemented a penalty-based mechanism that can be used to guide the optimization engine to produce a model with a desired behavior, such as single-channel open probability and use dependent effects. To streamline data collection for experiments in brain slice preparations, we developed a 3D virtual software environment that incorporates data from micro-positioning motors and scientific cameras in real-time. This environment provides positional feedback to the investigator and allows for the creation of data maps including both images and electrical recordings. We have also produced semi-automatic targeting procedures that simplifies the overall experimental experience. Experimentally, this work also examines how the kinetic mechanism of voltage gated Na channels regulates the neuronal firing of brainstem respiratory neurons. These raphe neurons are intrinsic pacemakers that do not rely on synaptic connections to elicit activity. I explored how intracellular calcium regulates the kinetics of TTX-sensitive Na+ currents using whole-cell patch clamp electrophysiology. Established with intracellular Ca2+ buffers, high [Ca2+] levels greater than ~7 [micro]M did not change the voltage dependence of steady-state activation and inactivation, but slightly slowed inactivation time course. However, the recovery from inactivation and use dependence inactivation is slowed by high intracellular [Ca2+]. Overall, these approaches described in this work have improved data acquisition and data analysis to create better ion channel models and enhance the electrophysiology experience.Includes bibliographical reference

Enhancing multi-scale cardiac simulations by coupling electrophysiology and mechanics: a flexible high performance approach to cardiac electromechanics

This work focuses on the development of computational methods for the simulation of the propagation of the electrical potential in the heart and of the resulting mechanical contraction. The interaction of these two physical phenomena is described by an electromechanical model which consists of the monodomain system, which describes the propagation of the action potential in the cardiac tissue, and the equations of incompressible elasticity, which describe its mechanical response. In fully-coupled electromechanical simulations, two main computational challenges are usually identified in literature: the time integration of the monodomain system and the efficient solution of the equations of incompressible elasticity. These two are the actual bottlenecks in the realization of accurate and efficient fully-coupled electromechanical simulations. The first computational challenge arises from the discretization in time of the equations that describe the electrical activation of cardiac tissue. The monodomain system should be discretized employing both fine spatial grids and small time-steps, to capture the spatial steep gradients typical of the action potential and the behavior of the stiff gating variables, respectively. To obtain an accurate and computationally-cheap numerical solution, we propose a novel method based on coupling high-order backward differentiation formulae with high-order exponential time stepping schemes for the time integration of the monodomain system. We propose a novel quasi-Newton approach for the implicit discretization of the monodomain equation. We also compare this latter approach against a complex step differentiation-based approach. As a result, we show by means of numerical tests the accuracy of the developed strategies and how the use of high-order time integration schemes affects the simulation of post- processed quantities of clinical relevance such as the conduction velocity. The second computational challenge is due to the structure the discretization of the equations of incompressible elasticity. Due to the incompressibility constraint, the arising linear system has a saddle point structure for which standard solution methods such as multigrid or domain de- composition do not provide optimal convergence if not properly adapted. In order to overcome this problematic, we propose a segregated multigrid preconditioned solution method. The segregated approach allows to recast the saddle-point problem into two elliptic problems for which multigrid methods are shown to provide optimal convergence. The electromechanical model is employed to evaluate the effects of geometrical changes due to the contraction of the heart on simulated electrocardiograms. Finally, the effect of different electrical activations on the resulting pressure-volume loops is investigated by coupling the electromechanical model with a lumped model of the circulatory system

Structural mechanisms of gating at the selectivity filter of the human cardiac ryanodine receptor (hRyR2) channel

The cardiac ryanodine receptor (RyR2) contains structural elements within the channel pore that function as gates to regulate the release of intracellular calcium, initiating cardiac muscle contraction. The precise regulation of these gates is critical in maintaining normal cardiac function, and channel dysfunction, resulting in altered calcium handling, underlies the mechanisms of arrhythmia and sudden cardiac death. The enormous size of RyR2 has impeded the gathering of detailed structural information, hence the structural determinants for channel gating remain unknown. Structural modelling studies have revealed similarities between the RyR2 pore and the K+ channel, KcsA, providing a framework in which to test channel gating mechanisms. A region termed the selectivity filter is a gating component in K+ channels, involved in inactivation and flicker closings, and its conformation is maintained by a transient hydrogen-bonding network. This project examined the role of the RyR2 selectivity filter in channel gating by generating mutants at analogous positions to KcsA that either disrupted (Y4813A, D4829A and Y4839A) or maintained (Y4813W, D4829E and Y4839W) a proposed hydrogen-bonding network, and assessed their intracellular Ca2+ release, ryanodine modification and biophysical properties. Y4813A and D4829A had drastic effects on channel function, whereas retaining physicochemical properties of conservative mutations, Y4813W and Y4839W, maintained the functional characteristics of WT RyR2. Flicker closings were affected by Y4839A mutation however, in general, single-channel gating for Y4813W, Y4839A and Y4839W was comparable to WT RyR2. Interestingly, monitoring single-channels for prolonged periods revealed novel insights into channel behaviour, characterised by inherent fluctuations in channel activity under steady-state conditions. This thesis reveals that the selectivity filter region is an important component for RyR2 channel function. However, it remains unclear whether the selectivity filter regulates channel gating, as the proposed hydrogen-bonding network would not be possible due to altered residue distances revealed from recent high-resolution RyR structural models