Balanced implicit Patankar-Euler methods for positive solutions of stochastic differential equations of biological regulatory systems.

Stochastic differential equations (SDEs) are a powerful tool to model fluctuations and uncertainty in complex systems. Although numerical methods have been designed to simulate SDEs effectively, it is still problematic when numerical solutions may be negative, but application problems require positive simulations. To address this issue, we propose balanced implicit Patankar-Euler methods to ensure positive simulations of SDEs. Instead of considering the addition of balanced terms to explicit methods in existing balanced methods, we attempt the deletion of possible negative terms from the explicit methods to maintain positivity of numerical simulations. The designed balanced terms include negative-valued drift terms and potential negative diffusion terms. The proposed method successfully addresses the issue of divisions with very small denominators in our recently designed stochastic Patankar method. Stability analysis shows that the balanced implicit Patankar-Euler method has much better stability properties than our recently designed composite Patankar-Euler method. Four SDE systems are used to examine the effectiveness, accuracy, and convergence properties of balanced implicit Patankar-Euler methods. Numerical results suggest that the proposed balanced implicit Patankar-Euler method is an effective and efficient approach to ensure positive simulations when any appropriate stepsize is used in simulating SDEs of biological regulatory systems.

Composite Patankar-Euler methods for positive simulations of stochastic differential equation models for biological regulatory systems.

Stochastic differential equations (SDE) are a powerful tool to model biological regulatory processes with intrinsic and extrinsic noise. However, numerical simulations of SDE models may be problematic if the values of noise terms are negative and large, which is not realistic for biological systems since the molecular copy numbers or protein concentrations should be non-negative. To address this issue, we propose the composite Patankar-Euler methods to obtain positive simulations of SDE models. A SDE model is separated into three parts, namely, the positive-valued drift terms, negative-valued drift terms, and diffusion terms. We first propose the deterministic Patankar-Euler method to avoid negative solutions generated from the negative-valued drift terms. The stochastic Patankar-Euler method is designed to avoid negative solutions generated from both the negative-valued drift terms and diffusion terms. These Patankar-Euler methods have the strong convergence order of a half. The composite Patankar-Euler methods are the combinations of the explicit Euler method, deterministic Patankar-Euler method, and stochastic Patankar-Euler method. Three SDE system models are used to examine the effectiveness, accuracy, and convergence properties of the composite Patankar-Euler methods. Numerical results suggest that the composite Patankar-Euler methods are effective methods to ensure positive simulations when any appropriate stepsize is used.

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.

Persistence as an Optimal Hedging Strategy.

Bacteria invest in a slow-growing subpopulation, called persisters, to ensure survival in the face of uncertainty. This hedging strategy is remarkably similar to financial hedging, where diversifying an investment portfolio protects against economic uncertainty. We provide a new, to our knowledge, theoretical foundation for understanding cellular hedging by unifying the study of biological population dynamics and the mathematics of financial risk management through optimal control theory. Motivated by the widely accepted role of volatility in the emergence of persistence, we consider several models of environmental volatility described by continuous-time stochastic processes. This allows us to study an emergent cellular hedging strategy that maximizes the expected per capita growth rate of the population. Analytical and simulation results probe the optimal persister strategy, revealing results that are consistent with experimental observations and suggest new opportunities for experimental investigation and design. Overall, we provide a new, to our knowledge, way of conceptualizing and modeling cellular decision making in volatile environments by explicitly unifying theory from mathematical biology and finance.

Effect of hydroxylysine-O-glycosylation on the structure of type I collagen molecule: A computational study.

Collagen undergoes many types of post-translational modifications (PTMs), including intracellular modifications and extracellular modifications. Among these PTMs, glycosylation of hydroxylysine (Hyl) is the most complicated. Experimental studies demonstrated that this PTM ceases once the collagen triple helix is formed and that Hyl-O-glycosylation modulates collagen fibrillogenesis. However, the underlying atomic-level mechanisms of these phenomena remain unclear. In this study, we first adapted the force field parameters for O-linkages between Hyl and carbohydrates and then investigated the influence of Hyl-O-glycosylation on the structure of type I collagen molecule, by performing comprehensive molecular dynamic simulations in explicit solvent of collagen molecule segment with and without the glycosylation of Hyl. Data analysis demonstrated that (i) collagen triple helices remain in a triple-helical structure upon glycosylation of Hyl; (ii) glycosylation of Hyl modulates the peptide backbone conformation and their solvation environment in the vicinity and (iii) the attached sugars are arranged such that their hydrophilic faces are well exposed to the solvent, while their hydrophobic faces point towards the hydrophobic portions of collagen. The adapted force field parameters for O-linkages between Hyl and carbohydrates will aid future computational studies on proteins with Hyl-O-glycosylation. In addition, this work, for the first time, presents the detailed effect of Hyl-O-glycosylation on the structure of human type I collagen at the atomic level, which may provide insights into the design and manufacture of collagenous biomaterials and the development of biomedical therapies for collagen-related diseases.

Designing combination therapies using multiple optimal controls.

Strategic management of populations of interacting biological species routinely requires interventions combining multiple treatments or therapies. This is important in key research areas such as ecology, epidemiology, wound healing and oncology. Despite the well developed theory and techniques for determining single optimal controls, there is limited practical guidance supporting implementation of combination therapies. In this work we use optimal control theory to calculate optimal strategies for applying combination therapies to a model of acute myeloid leukaemia. We present a versatile framework to systematically explore the trade-offs that arise in designing combination therapy protocols using optimal control. We consider various combinations of continuous and bang-bang (discrete) controls, and we investigate how the control dynamics interact and respond to changes in the weighting and form of the pay-off characterising optimality. We demonstrate that the optimal controls respond non-linearly to treatment strength and control parameters, due to the interactions between species. We discuss challenges in appropriately characterising optimality in a multiple control setting and provide practical guidance for applying multiple optimal controls. Code used in this work to implement multiple optimal controls is available on GitHub.

A population of bang-bang switches of defective interfering particles makes within-host dynamics of dengue virus controllable.

The titre of virus in a dengue patient and the duration of this viraemia has a profound effect on whether or not a mosquito will become infected when it feeds on the patient and this, in turn, is a key driver of the magnitude of a dengue outbreak. The assessment of the heterogeneity of viral dynamics in dengue-infected patients and its precise treatment are still uncertain. Infection onset, patient physiology and immune response are thought to play major roles in the development of the viral load. Research has explored the interference and spontaneous generation of defective virus particles, but have not examined both the antibody and defective particles during natural infection. We explore the intrinsic variability in the within-host dynamics of viraemias for a population of patients using the method of population of models (POMs). A dataset from 208 patients is used to initially calibrate 20,000 models for the infection kinetics for each of the four dengue virus serotypes. The calibrated POMs suggests that naturally generated defective particles may interfere with the viraemia, but the generated defective virus particles are not adequate to reduce high fever and viraemia duration. The effect of adding excess defective dengue virus interfering particles to patients as a therapeutic is evaluated using the calibrated POMs in a bang-bang (on-off or two-step) optimal control setting. Bang-bang control is a class of binary feedback control that turns either 'ON' or 'OFF' at different time points, determined by the system feedback. Here, the bang-bang control estimates the mathematically optimal dose and duration of the intervention for each model in the POM set.

Variability in electrophysiological properties and conducting obstacles controls re-entry risk in heterogeneous ischaemic tissue.

Ischaemia, in which inadequate blood supply compromises and eventually kills regions of cardiac tissue, can cause many types of arrhythmia, some life-threatening. A significant component of this is the effects of the resulting hypoxia, and concomitant hyperklaemia and acidosis, on the electrophysiological properties of myocytes. Clinical and experimental data have also shown that regions of structural heterogeneity (fibrosis, necrosis, fibro-fatty infiltration) can act as triggers for arrhythmias under acute ischaemic conditions. Mechanistic models have successfully captured these effects in silico. However, the relative significance of these separate facets of the condition, and how sensitive arrhythmic risk is to the extents of each, is far less explored. In this work, we use partitioned Gaussian process emulation and new metrics for source-sink mismatch that rely on simulations of bifurcating cardiac fibres to interrogate a model of heterogeneous ischaemic tissue. Re-entries were most sensitive to the level of hypoxia and the fraction of non-excitable tissue. In addition, our results reveal both protective and pro-arrhythmic effects of hyperklaemia, and present the levels of hyperklaemia, hypoxia and percentage of non-excitable tissue that pose the highest arrhythmic risks. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Adsorption of Collagen-like Peptides onto Gold Nanosurfaces.

The molecular behavior of proteins in the presence of inorganic surfaces is of fundamental biological significance. Examples include extracellular matrix proteins interacting with gold nanoparticles and metallic implant biomaterials, such as titanium and stainless steels. Uncharged inorganic surfaces that interact strongly with the solution phase (hydrophilic surfaces) have been commonly used in disease treatments. A deep understanding of the molecular behavior of body proteins in the presence of hydrophilic surfaces is important in terms of clinical applications. However, the adsorption mechanism of proteins onto hydrophilic surfaces remains not fully understood. Here, comprehensive molecular dynamics simulations are carried out to study the molecular response of a human collagen molecule segment (CMS) to the presence of a planar gold surface (AuNS) in explicit solvent, aiming to unravel the adsorption mechanism of proteins onto hydrophilic surfaces. The results demonstrate that in the presence of AuNS, the CMS first biasedly diffuses toward AuNS, followed by anchoring to the gold surface, and finally adsorbs stepwise onto AuNS, where the protein adjusts its structure to maximize the interaction with AuNS. We conclude that adsorption of proteins onto hydrophilic surfaces adheres to three steps, namely, biased diffusion, anchoring, and stepwise adsorption accompanied by structural adaptation. The obtained adsorption mechanism provides insights into the development of inorganic surfaces for biomedical and therapeutic applications.

Optimal control of acute myeloid leukaemia.

Acute myeloid leukaemia (AML) is a blood cancer affecting haematopoietic stem cells. AML is routinely treated with chemotherapy, and so it is of great interest to develop optimal chemotherapy treatment strategies. In this work, we incorporate an immune response into a stem cell model of AML, since we find that previous models lacking an immune response are inappropriate for deriving optimal control strategies. Using optimal control theory, we produce continuous controls and bang-bang controls, corresponding to a range of objectives and parameter choices. Through example calculations, we provide a practical approach to applying optimal control using Pontryagin's Maximum Principle. In particular, we describe and explore factors that have a profound influence on numerical convergence. We find that the convergence behaviour is sensitive to the method of control updating, the nature of the control, and to the relative weighting of terms in the objective function. All codes we use to implement optimal control are made available.

Interaction of gold nanosurfaces/nanoparticles with collagen-like peptides.

Nanotechnology has quickly emerged as a promising research field with potential effects in disease treatments. For example, gold nanoparticles (AuNPs) have been extensively used in diagnostics and therapeutics. When administrated into human tissues, AuNPs first encounter extracellular matrix (ECM) molecules. Amongst all the ECM components, collagen is the main tension-resisting constituent, whose biofunctional and mechanical properties are strongly dependent on its hierarchical structure. Therefore, an in-depth understanding of the structural response of collagen to the presence of gold nanosurfaces (AuNS) and AuNPs is crucial in terms of clinical applications of AuNPs. However, detailed understanding of the molecular-level and atomic-level interaction between AuNS/AuNPs and collagen in the ECM is elusive. In this study, comprehensive molecular dynamics (MD) simulations have been performed to investigate the molecular behaviour of a collagen molecule segment (CMS) in the presence of AuNS/AuNPs in explicit water, aiming to explore the interaction of AuNS/AuNPs with collagen triple helices at the molecular and atomic levels. The results show that the CMS forms a rapid association with AuNS/AuNPs and undergoes a severe unfolding upon adsorption on AuNS/AuNPs, indicating an unfolding propensity of gold surfaces. We conclude that collagen triple helices unfold readily on AuNS and bare AuNPs, due to the interaction of gold surfaces with the protein backbone. The revealed clear unfolding nature and the unravelled atomic-level unfolding mechanism of collagen triple helices onto AuNPs contribute to the development of AuNPs for biomedical and therapeutic applications, and the design of gold-binding proteins.

From ionic to cellular variability in human atrial myocytes: an integrative computational and experimental study.

Variability refers to differences in physiological function between individuals, which may translate into different disease susceptibility and treatment efficacy. Experiments in human cardiomyocytes face wide variability and restricted tissue access; under these conditions, computational models are a useful complementary tool. We conducted a computational and experimental investigation in cardiomyocytes isolated from samples of the right atrial appendage of patients undergoing cardiac surgery to evaluate the impact of variability in action potentials (APs) and subcellular ionic densities on Ca2+ transient dynamics. Results showed that 1) variability in APs and ionic densities is large, even within an apparently homogenous patient cohort, and translates into ±100% variation in ionic conductances; 2) experimentally calibrated populations of models with wide variations in ionic densities yield APs overlapping with those obtained experimentally, even if AP characteristics of the original generic model differed significantly from experimental APs; 3) model calibration with AP recordings restricts the variability in ionic densities affecting upstroke and resting potential, but redundancy in repolarization currents admits substantial variability in ionic densities; and 4) model populations constrained with experimental APs and ionic densities exhibit three Ca2+ transient phenotypes, differing in intracellular Ca2+ handling and Na+/Ca2+ membrane extrusion. These findings advance our understanding of the impact of variability in human atrial electrophysiology. NEW & NOTEWORTHY Variability in human atrial electrophysiology is investigated by integrating for the first time cellular-level and ion channel recordings in computational electrophysiological models. Ion channel calibration restricts current densities but not cellular phenotypic variability. Reduced Na+/Ca2+ exchanger is identified as a primary mechanism underlying diastolic Ca2+ fluctuations in human atrial myocytes.

In Vivo and In Silico Investigation Into Mechanisms of Frequency Dependence of Repolarization Alternans in Human Ventricular Cardiomyocytes.

RATIONALE: Repolarization alternans (RA) are associated with arrhythmogenesis. Animal studies have revealed potential mechanisms, but human-focused studies are needed. RA generation and frequency dependence may be determined by cell-to-cell variability in protein expression, which is regulated by genetic and external factors. OBJECTIVE: To characterize in vivo RA in human and to investigate in silico using human models, the ionic mechanisms underlying the frequency-dependent differences in RA behavior identified in vivo. METHODS AND RESULTS: In vivo electrograms were acquired at 240 sites covering the epicardium of 41 patients at 6 cycle lengths (600-350 ms). In silico investigations were conducted using a population of biophysically detailed human models incorporating variability in protein expression and calibrated using in vivo recordings. Both in silico and in vivo, 2 types of RA were identified, with Fork- and Eye-type restitution curves, based on RA persistence or disappearance, respectively, at fast pacing rates. In silico simulations show that RA are strongly correlated with fluctuations in sarcoplasmic reticulum calcium, because of strong release and weak reuptake. Large L-type calcium current conductance is responsible for RA disappearance at fast frequencies in Eye-type (30% larger in Eye-type versus Fork-type; P<0.01), because of sarcoplasmic reticulum Ca(2+) ATPase pump potentiation caused by frequency-induced increase in intracellular calcium. Large Na(+)/Ca(2+) exchanger current is the main driver in translating Ca(2+) fluctuations into RA. CONCLUSIONS: In human in vivo and in silico, 2 types of RA are identified, with RA persistence/disappearance as frequency increases. In silico, L-type calcium current and Na(+)/Ca(2+) exchanger current determine RA human cell-to-cell differences through intracellular and sarcoplasmic reticulum calcium regulation.

Modelling biochemical reaction systems by stochastic differential equations with reflection.

In this paper, we gave a new framework for modelling and simulating biochemical reaction systems by stochastic differential equations with reflection not in a heuristic way but in a mathematical way. The model is computationally efficient compared with the discrete-state Markov chain approach, and it ensures that both analytic and numerical solutions remain in a biologically plausible region. Specifically, our model mathematically ensures that species numbers lie in the domain D, which is a physical constraint for biochemical reactions, in contrast to the previous models. The domain D is actually obtained according to the structure of the corresponding chemical Langevin equations, i.e., the boundary is inherent in the biochemical reaction system. A variant of projection method was employed to solve the reflected stochastic differential equation model, and it includes three simple steps, i.e., Euler-Maruyama method was applied to the equations first, and then check whether or not the point lies within the domain D, and if not perform an orthogonal projection. It is found that the projection onto the closure D¯ is the solution to a convex quadratic programming problem. Thus, existing methods for the convex quadratic programming problem can be employed for the orthogonal projection map. Numerical tests on several important problems in biological systems confirmed the efficiency and accuracy of this approach.

Application of stochastic phenomenological modelling to cell-to-cell and beat-to-beat electrophysiological variability in cardiac tissue.

Variability in the action potential of isolated myocytes and tissue samples is observed in experimental studies. Variability is manifested as both differences in the action potential (AP) morphology between cells (extrinsic variability), and also 'intrinsic' or beat-to-beat variability of repolarization (BVR) in the AP duration of each cell. We studied the relative contributions of experimentally recorded intrinsic and extrinsic variability to dispersion of repolarization in tissue. We developed four cell-specific parameterizations of a phenomenological stochastic differential equation AP model exhibiting intrinsic variability using APs recorded from isolated guinea pig ventricular myocytes exhibiting BVR. We performed simulations in tissue using the four different model parameterizations in the presence and the absence of both intrinsic and extrinsic variability. We altered the coupling of the tissue to determine how inter-cellular coupling affected the dispersion of the AP duration in tissue. Both intrinsic and extrinsic variability were gradually revealed by reduction of tissue coupling. However, the recorded extrinsic variability between individual myocytes produced a greater degree of dispersion in repolarization in tissue than the intrinsic variability of each myocyte.

Stochastic dynamics of interacting haematopoietic stem cell niche lineages.

Since we still know very little about stem cells in their natural environment, it is useful to explore their dynamics through modelling and simulation, as well as experimentally. Most models of stem cell systems are based on deterministic differential equations that ignore the natural heterogeneity of stem cell populations. This is not appropriate at the level of individual cells and niches, when randomness is more likely to affect dynamics. In this paper, we introduce a fast stochastic method for simulating a metapopulation of stem cell niche lineages, that is, many sub-populations that together form a heterogeneous metapopulation, over time. By selecting the common limiting timestep, our method ensures that the entire metapopulation is simulated synchronously. This is important, as it allows us to introduce interactions between separate niche lineages, which would otherwise be impossible. We expand our method to enable the coupling of many lineages into niche groups, where differentiated cells are pooled within each niche group. Using this method, we explore the dynamics of the haematopoietic system from a demand control system perspective. We find that coupling together niche lineages allows the organism to regulate blood cell numbers as closely as possible to the homeostatic optimum. Furthermore, coupled lineages respond better than uncoupled ones to random perturbations, here the loss of some myeloid cells. This could imply that it is advantageous for an organism to connect together its niche lineages into groups. Our results suggest that a potential fruitful empirical direction will be to understand how stem cell descendants communicate with the niche and how cancer may arise as a result of a failure of such communication.

Stochastic linear multistep methods for the simulation of chemical kinetics.

In this paper, we introduce the Stochastic Adams-Bashforth (SAB) and Stochastic Adams-Moulton (SAM) methods as an extension of the τ-leaping framework to past information. Using the Θ-trapezoidal τ-leap method of weak order two as a starting procedure, we show that the k-step SAB method with k ≥ 3 is order three in the mean and correlation, while a predictor-corrector implementation of the SAM method is weak order three in the mean but only order one in the correlation. These convergence results have been derived analytically for linear problems and successfully tested numerically for both linear and non-linear systems. A series of additional examples have been implemented in order to demonstrate the efficacy of this approach.

Digital twinning of the human ventricular activation sequence to Clinical 12-lead ECGs and magnetic resonance imaging using realistic Purkinje networks for in silico clinical trials.

Cardiac in silico clinical trials can virtually assess the safety and efficacy of therapies using human-based modelling and simulation. These technologies can provide mechanistic explanations for clinically observed pathological behaviour. Designing virtual cohorts for in silico trials requires exploiting clinical data to capture the physiological variability in the human population. The clinical characterisation of ventricular activation and the Purkinje network is challenging, especially non-invasively. Our study aims to present a novel digital twinning pipeline that can efficiently generate and integrate Purkinje networks into human multiscale biventricular models based on subject-specific clinical 12-lead electrocardiogram and magnetic resonance recordings. Essential novel features of the pipeline are the human-based Purkinje network generation method, personalisation considering ECG R wave progression as well as QRS morphology, and translation from reduced-order Eikonal models to equivalent biophysically-detailed monodomain ones. We demonstrate ECG simulations in line with clinical data with clinical image-based multiscale models with Purkinje in four control subjects and two hypertrophic cardiomyopathy patients (simulated and clinical QRS complexes with Pearson's correlation coefficients > 0.7). Our methods also considered possible differences in the density of Purkinje myocardial junctions in the Eikonal-based inference as regional conduction velocities. These differences translated into regional coupling effects between Purkinje and myocardial models in the monodomain formulation. In summary, we demonstrate a digital twin pipeline enabling simulations yielding clinically consistent ECGs with clinical CMR image-based biventricular multiscale models, including personalised Purkinje in healthy and cardiac disease conditions.

Recurrence quantification analysis for fine-scale characterisation of arrhythmic patterns in cardiac tissue.

This paper uses recurrence quantification analysis (RQA) combined with entropy measures and organization indices to characterize arrhythmic patterns and dynamics in computer simulations of cardiac tissue. We performed different simulations of cardiac tissues of sizes comparable to the human heart atrium. In these simulations, we observed four classic arrhythmic patterns: a spiral wave anchored to a highly fibrotic region resulting in sustained re-entry, a meandering spiral wave, fibrillation, and a spiral wave anchored to a scar region that breaks up into wavelets away from the main rotor. A detailed analysis revealed that, within the same simulation, maps of RQA metrics could differentiate regions with regular AP propagation from ones with chaotic activity. In particular, the combination of two RQA metrics, the length of the longest diagonal string of recurrence points and the mean length of diagonal lines, was able to identify the location of rotor tips, which are the active elements that maintain spiral waves and fibrillation. By proposing low-dimensional models based on the mean value and spatial correlation of metrics calculated from membrane potential time series, we identify RQA-based metrics that successfully separate the four different types of cardiac arrhythmia into distinct regions of the feature space, and thus might be used for automatic classification, in particular distinguishing between fibrillation driven by self-sustaining chaos and that created by a persistent rotor and wavebreak. We also discuss the practical applicability of such an approach.

Implications of different membrane compartmentalization models in particle-based in silico studies.

Studying membrane dynamics is important to understand the cellular response to environmental stimuli. A decisive spatial characteristic of the plasma membrane is its compartmental structure created by the actin-based membrane-skeleton (fences) and anchored transmembrane proteins (pickets). Particle-based reaction-diffusion simulation of the membrane offers a suitable temporal and spatial resolution to analyse its spatially heterogeneous and stochastic dynamics. Fences have been modelled via hop probabilities, potentials or explicit picket fences. Our study analyses the different approaches' constraints and their impact on simulation results and performance. Each of the methods comes with its own constraints; the picket fences require small timesteps, potential fences might induce a bias in diffusion in crowded systems, and probabilistic fences, in addition to carefully scaling the probability with the timesteps, induce higher computational costs for each propagation step.