GillesPy2: A Biochemical Modeling Framework for Simulation Driven Biological Discovery.

Stochastic modeling has become an essential tool for studying biochemical reaction networks. There is a growing need for user-friendly and feature-complete software for model design and simulation. To address this need, we present GillesPy2, an open-source framework for building and simulating mathematical and biochemical models. GillesPy2, a major upgrade from the original GillesPy package, is now a stand-alone Python 3 package. GillesPy2 offers an intuitive interface for robust and reproducible model creation, facilitating rapid and iterative development. In addition to expediting the model creation process, GillesPy2 offers efficient algorithms to simulate stochastic, deterministic, and hybrid stochastic-deterministic models.

Interpretable polynomial neural ordinary differential equations.

Neural networks have the ability to serve as universal function approximators, but they are not interpretable and do not generalize well outside of their training region. Both of these issues are problematic when trying to apply standard neural ordinary differential equations (ODEs) to dynamical systems. We introduce the polynomial neural ODE, which is a deep polynomial neural network inside of the neural ODE framework. We demonstrate the capability of polynomial neural ODEs to predict outside of the training region, as well as to perform direct symbolic regression without using additional tools such as SINDy.

Decision heuristics in contexts integrating action selection and execution.

Heuristics can inform human decision making in complex environments through a reduction of computational requirements (accuracy-resource trade-off) and a robustness to overparameterisation (less-is-more). However, tasks capturing the efficiency of heuristics typically ignore action proficiency in determining rewards. The requisite movement parameterisation in sensorimotor control questions whether heuristics preserve efficiency when actions are nontrivial. We developed a novel action selection-execution task requiring joint optimisation of action selection and spatio-temporal skillful execution. State-appropriate choices could be determined by a simple spatial heuristic, or by more complex planning. Computational models of action selection parsimoniously distinguished human participants who adopted the heuristic from those using a more complex planning strategy. Broader comparative analyses then revealed that participants using the heuristic showed combined decisional (selection) and skill (execution) advantages, consistent with a less-is-more framework. In addition, the skill advantage of the heuristic group was predominantly in the core spatial features that also shaped their decision policy, evidence that the dimensions of information guiding action selection might be yoked to salient features in skill learning.

Protosequences in human cortical organoids model intrinsic states in the developing cortex.

Neuronal firing sequences are thought to be the basic building blocks of neural coding and information broadcasting within the brain. However, when sequences emerge during neurodevelopment remains unknown. We demonstrate that structured firing sequences are present in spontaneous activity of human brain organoids and ex vivo neonatal brain slices from the murine somatosensory cortex. We observed a balance between temporally rigid and flexible firing patterns that are emergent phenomena in human brain organoids and early postnatal murine somatosensory cortex, but not in primary dissociated cortical cultures. Our findings suggest that temporal sequences do not arise in an experience-dependent manner, but are rather constrained by an innate preconfigured architecture established during neurogenesis. These findings highlight the potential for brain organoids to further explore how exogenous inputs can be used to refine neuronal circuits and enable new studies into the genetic mechanisms that govern assembly of functional circuitry during early human brain development.

Preliminary study: quantification of chronic pain from physiological data.

Introduction: It is unknown if physiological changes associated with chronic pain could be measured with inexpensive physiological sensors. Recently, acute pain and laboratory-induced pain have been quantified with physiological sensors. Objectives: To investigate the extent to which chronic pain can be quantified with physiological sensors. Methods: Data were collected from chronic pain sufferers who subjectively rated their pain on a 0 to 10 visual analogue scale, using our recently developed pain meter. Physiological variables, including pulse, temperature, and motion signals, were measured at head, neck, wrist, and finger with multiple sensors. To quantify pain, features were first extracted from 10-second windows. Linear models with recursive feature elimination were fit for each subject. A random forest regression model was used for pain score prediction for the population-level model. Results: Predictive performance was assessed using leave-one-recording-out cross-validation and nonparametric permutation testing. For individual-level models, 5 of 12 subjects yielded intraclass correlation coefficients between actual and predicted pain scores of 0.46 to 0.75. For the population-level model, the random forest method yielded an intraclass correlation coefficient of 0.58. Bland-Altman analysis shows that our model tends to overestimate the lower end of the pain scores and underestimate the higher end. Conclusion: This is the first demonstration that physiological data can be correlated with chronic pain, both for individuals and populations. Further research and more extensive data will be required to assess whether this approach could be used as a "chronic pain meter" to assess the level of chronic pain in patients.

Functional neuronal circuitry and oscillatory dynamics in human brain organoids.

Human brain organoids replicate much of the cellular diversity and developmental anatomy of the human brain. However, the physiology of neuronal circuits within organoids remains under-explored. With high-density CMOS microelectrode arrays and shank electrodes, we captured spontaneous extracellular activity from brain organoids derived from human induced pluripotent stem cells. We inferred functional connectivity from spike timing, revealing a large number of weak connections within a skeleton of significantly fewer strong connections. A benzodiazepine increased the uniformity of firing patterns and decreased the relative fraction of weakly connected edges. Our analysis of the local field potential demonstrate that brain organoids contain neuronal assemblies of sufficient size and functional connectivity to co-activate and generate field potentials from their collective transmembrane currents that phase-lock to spiking activity. These results point to the potential of brain organoids for the study of neuropsychiatric diseases, drug action, and the effects of external stimuli upon neuronal networks.

Dynamic Phase Extraction: Applications in Pulse Rate Variability.

Pulse rate variability is a physiological parameter that has been extensively studied and correlated with many physical ailments. However, the phase relationship between inter-beat interval, IBI, and breathing has very rarely been studied. Develop a technique by which the phase relationship between IBI and breathing can be accurately and efficiently extracted from photoplethysmography (PPG) data. A program based on Lock-in Amplifier technology was written in Python to implement a novel technique, Dynamic Phase Extraction. It was tested using a breath pacer and a PPG sensor on 6 subjects who followed a breath pacer at varied breathing rates. The data were then analyzed using both traditional methods and the novel technique (Dynamic Phase Extraction) utilizing a breath pacer. Pulse data was extracted using a PPG sensor. Dynamic Phase Extraction (DPE) gave the magnitudes of the variation in IBI associated with breathing [Formula: see text] measured with photoplethysmography during paced breathing (with premature ventricular contractions, abnormal arrhythmias, and other artifacts edited out). [Formula: see text] correlated well with two standard measures of pulse rate variability: the Standard Deviation of the inter-beat interval (SDNN) (ρ = 0.911) and with the integrated value of the Power Spectral Density between 0.04 and 0.15 Hz (Low Frequency Power or LF Power) (ρ = 0.885). These correlations were comparable to the correlation between the SDNN and the LF Power (ρ = 0.877). In addition to the magnitude [Formula: see text], Dynamic Phase Extraction also gave the phase between the breath pacer and the changes in the inter-beat interval (IBI) due to respiratory sinus arrythmia (RSA), and correlated well with the phase extracted using a Fourier transform (ρ = 0.857). Dynamic Phase Extraction can extract both the phase between the breath pacer and the changes in IBI due to the respiratory sinus arrhythmia component of pulse rate variability ([Formula: see text], but is limited by needing a breath pacer.

Noise resistant synchronization and collective rhythm switching in a model of animal group locomotion.

Biology is suffused with rhythmic behaviour, and interacting biological oscillators often synchronize their rhythms with one another. Colonies of some ant species are able to synchronize their activity to fall into coherent bursts, but models of this phenomenon have neglected the potential effects of intrinsic noise and interspecific differences in individual-level behaviour. We investigated the individual and collective activity patterns of two Leptothorax ant species. We show that in one species (Leptothorax sp. W), ants converge onto rhythmic cycles of synchronized collective activity with a period of about 20 min. A second species (Leptothorax crassipilis) exhibits more complex collective dynamics, where dominant collective cycle periods range from 16 min to 2.8 h. Recordings that last 35 h reveal that, in both species, the same colony can exhibit multiple oscillation frequencies. We observe that workers of both species can be stimulated by nest-mates to become active after a refractory resting period, but the durations of refractory periods differ between the species and can be highly variable. We model the emergence of synchronized rhythms using an agent-based model informed by our empirical data. This simple model successfully generates synchronized group oscillations despite the addition of noise to ants' refractory periods. We also find that adding noise reduces the likelihood that the model will spontaneously switch between distinct collective cycle frequencies.

A positive feedback loop involving the Spa2 SHD domain contributes to focal polarization.

Focal polarization is necessary for finely arranged cell-cell interactions. The yeast mating projection, with its punctate polarisome, is a good model system for this process. We explored the critical role of the polarisome scaffold protein Spa2 during yeast mating with a hypothesis motivated by mathematical modeling and tested by in vivo experiments. Our simulations predicted that two positive feedback loops generate focal polarization, including a novel feedback pathway involving the N-terminal domain of Spa2. We characterized the latter using loss-of-function and gain-of-function mutants. The N-terminal region contains a Spa2 Homology Domain (SHD) which is conserved from yeast to humans, and when mutated largely reproduced the spa2Δ phenotype. Our work together with published data show that the SHD domain recruits Msb3/4 that stimulates Sec4-mediated transport of Bud6 to the polarisome. There, Bud6 activates Bni1-catalyzed actin cable formation, recruiting more Spa2 and completing the positive feedback loop. We demonstrate that disrupting this loop at any point results in morphological defects. Gain-of-function perturbations partially restored focal polarization in a spa2 loss-of-function mutant without restoring localization of upstream components, thus supporting the pathway order. Thus, we have collected data consistent with a novel positive feedback loop that contributes to focal polarization during pheromone-induced polarization in yeast.

Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods.

Identifying the reactions that govern a dynamical biological system is a crucial but challenging task in systems biology. In this work, we present a data-driven method to infer the underlying biochemical reaction system governing a set of observed species concentrations over time. We formulate the problem as a regression over a large, but limited, mass-action constrained reaction space and utilize sparse Bayesian inference via the regularized horseshoe prior to produce robust, interpretable biochemical reaction networks, along with uncertainty estimates of parameters. The resulting systems of chemical reactions and posteriors inform the biologist of potentially several reaction systems that can be further investigated. We demonstrate the method on two examples of recovering the dynamics of an unknown reaction system, to illustrate the benefits of improved accuracy and information obtained.

Extracellular detection of neuronal coupling.

We developed a method to non-invasively detect synaptic relationships among neurons from in vitro networks. Our method uses microelectrode arrays on which neurons are cultured and from which propagation of extracellular action potentials (eAPs) in single axons are recorded at multiple electrodes. Detecting eAP propagation bypasses ambiguity introduced by spike sorting. Our methods identify short latency spiking relationships between neurons with properties expected of synaptically coupled neurons, namely they were recapitulated by direct stimulation and were sensitive to changing the number of active synaptic sites. Our methods enabled us to assemble a functional subset of neuronal connectivity in our cultures.

Multiscale modeling meets machine learning: What can we learn?

Machine learning is increasingly recognized as a promising technology in the biological, biomedical, and behavioral sciences. There can be no argument that this technique is incredibly successful in image recognition with immediate applications in diagnostics including electrophysiology, radiology, or pathology, where we have access to massive amounts of annotated data. However, machine learning often performs poorly in prognosis, especially when dealing with sparse data. This is a field where classical physics-based simulation seems to remain irreplaceable. In this review, we identify areas in the biomedical sciences where machine learning and multiscale modeling can mutually benefit from one another: Machine learning can integrate physics-based knowledge in the form of governing equations, boundary conditions, or constraints to manage ill-posted problems and robustly handle sparse and noisy data; multiscale modeling can integrate machine learning to create surrogate models, identify system dynamics and parameters, analyze sensitivities, and quantify uncertainty to bridge the scales and understand the emergence of function. With a view towards applications in the life sciences, we discuss the state of the art of combining machine learning and multiscale modeling, identify applications and opportunities, raise open questions, and address potential challenges and limitations. We anticipate that it will stimulate discussion within the community of computational mechanics and reach out to other disciplines including mathematics, statistics, computer science, artificial intelligence, biomedicine, systems biology, and precision medicine to join forces towards creating robust and efficient models for biological systems.

Accelerated regression-based summary statistics for discrete stochastic systems via approximate simulators.

BACKGROUND: Approximate Bayesian Computation (ABC) has become a key tool for calibrating the parameters of discrete stochastic biochemical models. For higher dimensional models and data, its performance is strongly dependent on having a representative set of summary statistics. While regression-based methods have been demonstrated to allow for the automatic construction of effective summary statistics, their reliance on first simulating a large training set creates a significant overhead when applying these methods to discrete stochastic models for which simulation is relatively expensive. In this τ work, we present a method to reduce this computational burden by leveraging approximate simulators of these systems, such as ordinary differential equations and τ-Leaping approximations. RESULTS: We have developed an algorithm to accelerate the construction of regression-based summary statistics for Approximate Bayesian Computation by selectively using the faster approximate algorithms for simulations. By posing the problem as one of ratio estimation, we use state-of-the-art methods in machine learning to show that, in many cases, our algorithm can significantly reduce the number of simulations from the full resolution model at a minimal cost to accuracy and little additional tuning from the user. We demonstrate the usefulness and robustness of our method with four different experiments. CONCLUSIONS: We provide a novel algorithm for accelerating the construction of summary statistics for stochastic biochemical systems. Compared to the standard practice of exclusively training from exact simulator samples, our method is able to dramatically reduce the number of required calls to the stochastic simulator at a minimal loss in accuracy. This can immediately be implemented to increase the overall speed of the ABC workflow for estimating parameters in complex systems.

An arbitrary Lagrangian Eulerian smoothed particle hydrodynamics (ALE-SPH) method with a boundary volume fraction formulation for fluid-structure interaction.

We present a new weakly-compressible smoothed particle hydrodynamics (SPH) method capable of modeling non-slip fixed and moving wall boundary conditions. The formulation combines a boundary volume fraction (BVF) wall approach with the transport-velocity SPH method. The resulting method, named SPH-BVF, offers detection of arbitrarily shaped solid walls on-the-fly, with small computational overhead due to its local formulation. This simple framework is capable of solving problems that are difficult or infeasible for standard SPH, namely flows subject to large shear stresses or at moderate Reynolds numbers, and mass transfer in deformable boundaries. In addition, the method extends the transport-velocity formulation to reaction-diffusion transport of mass in Newtonian fluids and linear elastic solids, which is common in biological structures. Taken together, the SPH-BVF method provides a good balance of simplicity and versatility, while avoiding some of the standard obstacles associated with SPH: particle penetration at the boundaries, tension instabilities and anisotropic particle alignments, that hamper SPH from being applied to complex problems such as fluid-structure interaction in a biological system.

Associations of longitudinal D-Dimer and Factor II on early trauma survival risk.

BACKGROUND: Trauma-induced coagulopathy (TIC) is a disorder that occurs in one-third of severely injured trauma patients, manifesting as increased bleeding and a 4X risk of mortality. Understanding the mechanisms driving TIC, clinical risk factors are essential to mitigating this coagulopathic bleeding and is therefore essential for saving lives. In this retrospective, single hospital study of 891 trauma patients, we investigate and quantify how two prominently described phenotypes of TIC, consumptive coagulopathy and hyperfibrinolysis, affect survival odds in the first 25 h, when deaths from TIC are most prevalent. METHODS: We employ a joint survival model to estimate the longitudinal trajectories of the protein Factor II (% activity) and the log of the protein fragment D-Dimer ([Formula: see text]g/ml), representative biomarkers of consumptive coagulopathy and hyperfibrinolysis respectively, and tie them together with patient outcomes. Joint models have recently gained popularity in medical studies due to the necessity to simultaneously track continuously measured biomarkers as a disease evolves, as well as to associate them with patient outcomes. In this work, we estimate and analyze our joint model using Bayesian methods to obtain uncertainties and distributions over associations and trajectories. RESULTS: We find that a unit increase in log D-Dimer increases the risk of mortality by 2.22 [1.57, 3.28] fold while a unit increase in Factor II only marginally decreases the risk of mortality by 0.94 [0.91,0.96] fold. This suggests that, while managing consumptive coagulopathy and hyperfibrinolysis both seem to affect survival odds, the effect of hyperfibrinolysis is much greater and more sensitive. Furthermore, we find that the longitudinal trajectories, controlling for many fixed covariates, trend differently for different patients. Thus, a more personalized approach is necessary when considering treatment and risk prediction under these phenotypes. CONCLUSION: This study reinforces the finding that hyperfibrinolysis is linked with poor patient outcomes regardless of factor consumption levels. Furthermore, it quantifies the degree to which measured D-Dimer levels correlate with increased risk. The single hospital, retrospective nature can be understood to specify the results to this particular hospital's patients and protocol in treating trauma patients. Expanding to a multi-hospital setting would result in better estimates about the underlying nature of consumptive coagulopathy and hyperfibrinolysis with survival, regardless of protocol. Individual trajectories obtained with these estimates can be used to provide personalized dynamic risk prediction when making decisions regarding management of blood factors.

Experimentally Validated Reconstruction and Analysis of a Genome-Scale Metabolic Model of an Anaerobic Neocallimastigomycota Fungus.

Anaerobic gut fungi in the phylum Neocallimastigomycota typically inhabit the digestive tracts of large mammalian herbivores, where they play an integral role in the decomposition of raw lignocellulose into its constitutive sugar monomers. However, quantitative tools to study their physiology are lacking, partially due to their complex and unresolved metabolism that includes the largely uncharacterized fungal hydrogenosome. Modern omics approaches combined with metabolic modeling can be used to establish an understanding of gut fungal metabolism and develop targeted engineering strategies to harness their degradation capabilities for lignocellulosic bioprocessing. Here, we introduce a high-quality genome of the anaerobic fungus Neocallimastix lanati from which we constructed the first genome-scale metabolic model of an anaerobic fungus. Relative to its size (200 Mbp, sequenced at 62× depth), it is the least fragmented publicly available gut fungal genome to date. Of the 1,788 lignocellulolytic enzymes annotated in the genome, 585 are associated with the fungal cellulosome, underscoring the powerful lignocellulolytic potential of N. lanati The genome-scale metabolic model captures the primary metabolism of N. lanati and accurately predicts experimentally validated substrate utilization requirements. Additionally, metabolic flux predictions are verified by 13C metabolic flux analysis, demonstrating that the model faithfully describes the underlying fungal metabolism. Furthermore, the model clarifies key aspects of the hydrogenosomal metabolism and can be used as a platform to quantitatively study these biotechnologically important yet poorly understood early-branching fungi.IMPORTANCE Recent genomic analyses have revealed that anaerobic gut fungi possess both the largest number and highest diversity of lignocellulolytic enzymes of all sequenced fungi, explaining their ability to decompose lignocellulosic substrates, e.g., agricultural waste, into fermentable sugars. Despite their potential, the development of engineering methods for these organisms has been slow due to their complex life cycle, understudied metabolism, and challenging anaerobic culture requirements. Currently, there is no framework that can be used to combine multi-omic data sets to understand their physiology. Here, we introduce a high-quality PacBio-sequenced genome of the anaerobic gut fungus Neocallimastix lanati Beyond identifying a trove of lignocellulolytic enzymes, we use this genome to construct the first genome-scale metabolic model of an anaerobic gut fungus. The model is experimentally validated and sheds light on unresolved metabolic features common to gut fungi. Model-guided analysis will pave the way for deepening our understanding of anaerobic gut fungi and provides a systematic framework to guide strain engineering efforts of these organisms for biotechnological use.

Epidemiological modeling in StochSS Live!

SUMMARY: We present StochSS Live!, a web-based service for modeling, simulation and analysis of a wide range of mathematical, biological and biochemical systems. Using an epidemiological model of COVID-19, we demonstrate the power of StochSS Live! to enable researchers to quickly develop a deterministic or a discrete stochastic model, infer its parameters and analyze the results. AVAILABILITY AND IMPLEMENTATION: StochSS Live! is freely available at https://live.stochss.org/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Coordinating cell polarization and morphogenesis through mechanical feedback.

Many cellular processes require cell polarization to be maintained as the cell changes shape, grows or moves. Without feedback mechanisms relaying information about cell shape to the polarity molecular machinery, the coordination between cell polarization and morphogenesis, movement or growth would not be possible. Here we theoretically and computationally study the role of a genetically-encoded mechanical feedback (in the Cell Wall Integrity pathway) as a potential coordination mechanism between cell morphogenesis and polarity during budding yeast mating projection growth. We developed a coarse-grained continuum description of the coupled dynamics of cell polarization and morphogenesis as well as 3D stochastic simulations of the molecular polarization machinery in the evolving cell shape. Both theoretical approaches show that in the absence of mechanical feedback (or in the presence of weak feedback), cell polarity cannot be maintained at the projection tip during growth, with the polarization cap wandering off the projection tip, arresting morphogenesis. In contrast, for mechanical feedback strengths above a threshold, cells can robustly maintain cell polarization at the tip and simultaneously sustain mating projection growth. These results indicate that the mechanical feedback encoded in the Cell Wall Integrity pathway can provide important positional information to the molecular machinery in the cell, thereby enabling the coordination of cell polarization and morphogenesis.

Motor adaptation via distributional learning.

Objective. Both artificial and biological controllers experience errors during learning that are probabilistically distributed. We develop a framework for modeling distributions of errors and relating deviations in these distributions to neural activity.Approach. The biological system we consider is a task where human subjects are required to learn to minimize the roll of an inverted T-shaped object with an unbalanced weight (i.e. one side of the object is heavier than the other side) during lift. We also collect BOLD activity during this process. For our experimental setup, we define the state of the system to be the maximum magnitude roll of the object after lift onset and give subjects the goal of achieving the zero state.Main Results. We derive a model for this problem from a variant of Temporal Difference Learning. We then combine this model with Distributional Reinforcement Learning (DRL), a framework that involves defining a value distribution by treating the reward as stochastic. This model transforms the goal of the controller from achieving a target state, to achieving a distribution over distances from the target state. We call it a Distributional Temporal Difference Model (DTDM). The DTDM allows us to model errors in unsuccessfully minimizing object roll using deviations in the value distribution when the center of mass of the unbalanced object is changed. We compute deviations in global neural activity and show that they vary continuously with deviations in the value distribution. Different aspects might contribute to this global shift or signal difference, including a difference in grasp and lift force at lift onset, as well as sensory feedback of error/roll after lift onset. We predict that there exists a coordinated, global response to errors that incorporates all of this information, which is encoding the DTDM objective and used on subsequent trials enabling success. We validate the utility of the DTDM as a model for biological adaptation by using it to engineer a robotic controller to solve a similar problem.Significance. We develop a novel theoretical framework and show that it can be used to model a non-trivial motor learning task. Because this theoretical framework is consistent with state-of-the-art reinforcement learning, we can also use it to program a robot to perform a similar task. These results suggest a way to model the multiple subsystems composing global neural activity in a way that transfers well to engineering artificial intelligence.

Computational model of tranexamic acid on urokinase mediated fibrinolysis.

Understanding the coagulation process is critical to developing treatments for trauma and coagulopathies. Clinical studies on tranexamic acid (TXA) have resulted in mixed reports on its efficacy in improving outcomes in trauma patients. The largest study, CRASH-2, reported that TXA improved outcomes in patients who received treatment prior to 3 hours after the injury, but worsened outcomes in patients who received treatment after 3 hours. No consensus has been reached about the mechanism behind the duality of these results. In this paper we use a computational model for coagulation and fibrinolysis to propose that deficiencies or depletions of key anti-fibrinolytic proteins, specifically antiplasmin, a1-antitrypsin and a2-macroglobulin, can lead to worsened outcomes through urokinase-mediated hyperfibrinolysis.