Efficient computation of adjoint sensitivities at steady-state in ODE models of biochemical reaction networks.

Dynamical models in the form of systems of ordinary differential equations have become a standard tool in systems biology. Many parameters of such models are usually unknown and have to be inferred from experimental data. Gradient-based optimization has proven to be effective for parameter estimation. However, computing gradients becomes increasingly costly for larger models, which are required for capturing the complex interactions of multiple biochemical pathways. Adjoint sensitivity analysis has been pivotal for working with such large models, but methods tailored for steady-state data are currently not available. We propose a new adjoint method for computing gradients, which is applicable if the experimental data include steady-state measurements. The method is based on a reformulation of the backward integration problem to a system of linear algebraic equations. The evaluation of the proposed method using real-world problems shows a speedup of total simulation time by a factor of up to 4.4. Our results demonstrate that the proposed approach can achieve a substantial improvement in computation time, in particular for large-scale models, where computational efficiency is critical.

NEIL3-deficient bone marrow displays decreased hematopoietic capacity and reduced telomere length.

Deficiency of NEIL3, a DNA repair enzyme, has significant impact on mouse physiology, including vascular biology and gut health, processes related to aging. Leukocyte telomere length (LTL) is suggested as a marker of biological aging, and shortened LTL is associated with increased risk of cardiovascular disease. NEIL3 has been shown to repair DNA damage in telomere regions in vitro. Herein, we explored the role of NEIL3 in telomere maintenance in vivo by studying bone marrow cells from atherosclerosis-prone NEIL3-deficient mice. We found shortened telomeres and decreased activity of the telomerase enzyme in bone marrow cells derived from Apoe -/- Neil3 -/- as compared to Apoe -/- mice. Furthermore, Apoe -/- Neil3 -/- mice had decreased leukocyte levels as compared to Apoe -/- mice, both in bone marrow and in peripheral blood. Finally, RNA sequencing of bone marrow cells from Apoe -/- Neil3 -/- and Apoe -/- mice revealed different expression levels of genes involved in cell cycle regulation, cellular senescence and telomere protection. This study points to NEIL3 as a telomere-protecting protein in murine bone marrow in vivo.

AMICI: high-performance sensitivity analysis for large ordinary differential equation models.

SUMMARY: Ordinary differential equation models facilitate the understanding of cellular signal transduction and other biological processes. However, for large and comprehensive models, the computational cost of simulating or calibrating can be limiting. AMICI is a modular toolbox implemented in C++/Python/MATLAB that provides efficient simulation and sensitivity analysis routines tailored for scalable, gradient-based parameter estimation and uncertainty quantification. AVAILABILITYAND IMPLEMENTATION: AMICI is published under the permissive BSD-3-Clause license with source code publicly available on https://github.com/AMICI-dev/AMICI. Citeable releases are archived on Zenodo. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Ryanodine receptor dispersion disrupts Ca2+ release in failing cardiac myocytes.

Reduced cardiac contractility during heart failure (HF) is linked to impaired Ca2+ release from Ryanodine Receptors (RyRs). We investigated whether this deficit can be traced to nanoscale RyR reorganization. Using super-resolution imaging, we observed dispersion of RyR clusters in cardiomyocytes from post-infarction HF rats, resulting in more numerous, smaller clusters. Functional groupings of RyR clusters which produce Ca2+ sparks (Ca2+ release units, CRUs) also became less solid. An increased fraction of small CRUs in HF was linked to augmented 'silent' Ca2+ leak, not visible as sparks. Larger multi-cluster CRUs common in HF also exhibited low fidelity spark generation. When successfully triggered, sparks in failing cells displayed slow kinetics as Ca2+ spread across dispersed CRUs. During the action potential, these slow sparks protracted and desynchronized the overall Ca2+ transient. Thus, nanoscale RyR reorganization during HF augments Ca2+ leak and slows Ca2+ release kinetics, leading to weakened contraction in this disease.

Simple T-Wave Metrics May Better Predict Early Ischemia as Compared to ST Segment.

There is pressing clinical need to identify developing heart attack (infarction) in patients as early as possible. However, current state-of-the-art tools in clinical practice, underpinned by the evaluation of elevation of the ST segment of the 12-lead electrocardiogram (ECG), do not identify all patients suffering from lack of blood flow to the heart muscle (cardiac ischemia), worsening the risk for further adverse events and patient outcome overall. In this study, we aimed to explore and compare the portions of cardiac repolarization in the ECG that best capture the electrophysiological changes associated with ischemia. We developed three-dimensional electrophysiological models of the human ventricles and torso, incorporating biophysically-based membrane kinetics and realistic activation sequence, to compute simulated ECGs and their alteration with the application of simulated ischemia of differing severity in diverse regions of the heart. Results suggest that metrics based on the T-wave in addition to the ST segment may be more sensitive to detecting ischemia than those using the ST segment alone. Further research into how such simulation-aided risk assessment methods may aid workflows in extant clinical practice, with the ultimate goal of multimodality clinical support, is warranted.

Slow Calcium-Depolarization-Calcium waves may initiate fast local depolarization waves in ventricular tissue.

Intercellular calcium waves in cardiac myocytes are a well-recognized, if incompletely understood, phenomenon. In a variety of preparations, investigators have reported multi-cellular calcium waves or triggered propagated contractions, but the mechanisms of propagation and pathological importance of these events remain unclear. Here, we review existing experimental data and present a computational approach to investigate the mechanisms of multi-cellular calcium wave propagation. Over the past 50 years, the standard modeling paradigm for excitable cardiac tissue has seen increasingly detailed models of the dynamics of individual cells coupled in tissue solely by intercellular and interstitial current flow. Although very successful, this modeling regime has been unable to capture two important phenomena: 1) the slow intercellular calcium waves observed experimentally, and 2) how intercellular calcium events resulting in delayed after depolarizations at the cellular level could overcome a source-sink mismatch to initiate depolarization waves in tissue. In this paper, we introduce a mathematical model with subcellular spatial resolution, in which we allow both inter- and intracellular current flow and calcium diffusion. In simulations of coupled cells employing this model, we observe: a) slow inter-cellular calcium waves propagating at about 0.1 mm/s, b) faster Calcium-Depolarization-Calcium (CDC) waves, traveling at about 1 mm/s, and c) CDC-waves that can set off fast depolarization-waves (50 cm/s) in tissue with varying gap-junction conductivity.

Instabilities of the resting state in a mathematical model of calcium handling in cardiac myocytes.

We analyze a recently published model of calcium handling in cardiac myocytes in order to find conditions for the presence of instabilities in the resting state of the model. Such instabilities can create calcium waves which in turn may be able to initiate cardiac arrhythmias. The model was developed by Swietach, Spitzer and Vaughan-Jones in order to study the effect, on calcium waves, of varying ryanodine receptor (RyR)-permeability, sarco/endoplasmic reticulum calcium ATPase (SERCA) and calcium diffusion. We study the model using the extracellular calcium concentration c(e) and the maximal velocity of the SERCA-pump v(SERCA) as control parameters. In the (c(e),v(SERCA))-domain we derive an explicit function v∗=v∗(c(e)), and we claim that any resting state based on parameters that lie above the curve (i.e. any pair (c(e),v(SERCA)) such that with v(SERCA) > v∗(c(e))) is unstable in the sense that small perturbations will grow and can eventually turn into a calcium wave. And conversely; any pair (c(e),v(SERCA)) below the curve is stable in the sense that small perturbations to the resting state will decay to rest. This claim is supported by analyzing the stability of the system in terms of computing the eigenmodes of the linearized model. Furthermore, the claim is supported by direct simulations based on the non-linear model. Since the curve separating stable from unstable states is given as an explicit function, we can show how stability depends on other parameters of the model.

Note on a possible proarrhythmic property of antiarrhythmic drugs aimed at improving gap-junction coupling.

Reduced conduction velocity (CV) in the myocardium is well known to increase the probability of arrhythmia and can be caused by structural changes, reduced excitability of individual myocytes, or decreased electrical coupling in the tissue. Recently, investigators have developed antiarrhythmic drugs that target the connections between individual myocytes with the goal of restoring tissue CV, specifically through increasing gap-junction coupling. In a simple but qualitatively relevant mathematical model, we show here that the introduction of a drug that improves intercellular conductance will indeed increase the CV. However, conditions that would require such a drug, such as fibrotic remodeling, may also increase the load of fibroblasts. Fibroblasts may couple to myocytes in much the same way as myocytes couple to each other, and therefore the use of such an agent may also improve coupling between myocytes and fibroblasts. We present numerical examples illustrating that when the load of coupled fibroblasts on myocytes is low or nonexistent, the drug works as expected, i.e., the drug increases CV. On the other hand, when the fibroblast load is high, changes in CV are nonmonotonic, i.e., the CV first increases and then decreases with an increase in dosage. The existence of coupled fibroblasts may therefore impair the effect of the drug, and under unfortunate conditions may be proarrhythmic.

On the frequency of automaticity during ischemia in simulations based on stochastic perturbations of the Luo-Rudy 1 model.

AIMS: To compute the effects of parameter perturbations for single ischemic cardiac cells, and to determine how perturbations influenced the tendency for the cells to undergo spontaneous depolarization (automaticity) during 20 min of acute ischemia. METHODS: A modified Luo-Rudy 1 cell model was used. Since the range of biological variation and measurement errors is largely unknown, we conducted our study of the consequences of perturbations under the assumption that cell model parameters have a normal distribution with a 10% standard deviation. A total of 10000 random cell realizations were tested while varying important Luo-Rudy cell model parameters. Ischemia was modelled by deterministic functions chosen for the expected values of crucial ion concentrations and gating parameters as they developed with time, while realizing the respective parameter values from static normal distributions with a 10% standard deviation. RESULTS AND CONCLUSION: It was found that the tendency towards automaticity did increase as the stochastic parameters were varied. In particular, cells with standard Luo-Rudy parameter values did not become automatic during ischemia, whereas a significant portion of the cells with randomized parameter values did. The relative importance of model parameter variations was also determined and a sodium m-gate activation parameter was identified as the most critical parameter. The frequency of arrhythmic events during acute ischemia is known to be bell-shaped, with a peak at around 7-8 min after the onset of ischemia. Our simulations display a similar peak in the frequency of automaticity.

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.

Simulation of ST segment changes during subendocardial ischemia using a realistic 3-D cardiac geometry.

The mechanisms underlying the ST segment shifts associated with subendocardial ischemia remain unclear. The aim of this paper is to shed further light on the subject through numerical simulations of these shifts. A realistic three-dimensional model of the ventricles, including fiber rotation and anisotropy, is embedded in a nonhomogeneous torso model. A simplification of the bidomain model is used to calculate only the ST segment shift, assuming known values of the transmembrane potential during the plateau and rest phases. A similar simulation is performed in two dimensions. The simulation results suggest that subendocardial ischemia can be located by ST segment shift on the epicardial and torso surfaces. It is shown that ST elevation is associated with the transmural ischemic boundary, while ST depression is associated with the lateral ischemic boundaries.

An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso.

In this paper we present a numerical method for the bidomain model, which describes the electrical activity in the heart. The model consists of two partial differential equations (PDEs), which are coupled to systems of ordinary differential equations (ODEs) describing electrochemical reactions in the cardiac cells. Many applications require coupling these equations to a third PDE, describing the electrical fields in the torso surrounding the heart. The resulting system is challenging to solve numerically, because of its complexity and very strict resolution requirements in time and space. We propose a method based on operator splitting and a fully coupled discretization of the three PDEs. Numerical experiments show that for simple simulation cases and fine discretizations, the algorithm is second-order accurate in space and time.