Tight basis cycle representatives for persistent homology of large biological data sets.

Persistent homology (PH) is a popular tool for topological data analysis that has found applications across diverse areas of research. It provides a rigorous method to compute robust topological features in discrete experimental observations that often contain various sources of uncertainties. Although powerful in theory, PH suffers from high computation cost that precludes its application to large data sets. Additionally, most analyses using PH are limited to computing the existence of nontrivial features. Precise localization of these features is not generally attempted because, by definition, localized representations are not unique and because of even higher computation cost. Such a precise location is a sine qua non for determining functional significance, especially in biological applications. Here, we provide a strategy and algorithms to compute tight representative boundaries around nontrivial robust features in large data sets. To showcase the efficiency of our algorithms and the precision of computed boundaries, we analyze the human genome and protein crystal structures. In the human genome, we found a surprising effect of the impairment of chromatin loop formation on loops through chromosome 13 and the sex chromosomes. We also found loops with long-range interactions between functionally related genes. In protein homologs with significantly different topology, we found voids attributable to ligand-interaction, mutation, and differences between species.

Geometric and topological characterization of the cytoarchitecture of islets of Langerhans.

The islets of Langerhans are critical endocrine micro-organs that secrete hormones regulating energy metabolism in animals. Insulin and glucagon, secreted by beta and alpha cells, respectively, are responsible for metabolic switching between fat and glucose utilization. Dysfunction in their secretion and/or counter-regulatory influence leads to diabetes. Debate in the field centers on the cytoarchitecture of islets, as the signaling that governs hormonal secretion depends on structural and functional factors, including electrical connectivity, innervation, vascularization, and physical proximity. Much effort has therefore been devoted to elucidating which architectural features are significant for function and how derangements in these features are correlated or causative for dysfunction, especially using quantitative network science or graph theory characterizations. Here, we ask if there are non-local features in islet cytoarchitecture, going beyond standard network statistics, that are relevant to islet function. An example is ring structures, or cycles, of α and δ cells surrounding ß cell clusters or the opposite, ß cells surrounding α and δ cells. These could appear in two-dimensional islet section images if a sphere consisting of one cell type surrounds a cluster of another cell type. To address these issues, we developed two independent computational approaches, geometric and topological, for such characterizations. For the latter, we introduce an application of topological data analysis to determine locations of topological features that are biologically significant. We show that both approaches, applied to a large collection of islet sections, are in complete agreement in the context both of developmental and diabetes-related changes in islet characteristics. The topological approach can be applied to three-dimensional imaging data for islets as well.

Cryptic genetic variation enhances primate L1 retrotransposon survival by enlarging the functional coiled coil sequence space of ORF1p.

Accounting for continual evolution of deleterious L1 retrotransposon families, which can contain hundreds to thousands of members remains a major issue in mammalian biology. L1 activity generated upwards of 40% of some mammalian genomes, including humans where they remain active, causing genetic defects and rearrangements. L1 encodes a coiled coil-containing protein that is essential for retrotransposition, and the emergence of novel primate L1 families has been correlated with episodes of extensive amino acid substitutions in the coiled coil. These results were interpreted as an adaptive response to maintain L1 activity, however its mechanism remained unknown. Although an adventitious mutation can inactivate coiled coil function, its effect could be buffered by epistatic interactions within the coiled coil, made more likely if the family contains a diverse set of coiled coil sequences-collectively referred to as the coiled coil sequence space. Amino acid substitutions that do not affect coiled coil function (i.e., its phenotype) could be "hidden" from (not subject to) purifying selection. The accumulation of such substitutions, often referred to as cryptic genetic variation, has been documented in various proteins. Here we report that this phenomenon was in effect during the latest episode of primate coiled coil evolution, which occurred 30-10 MYA during the emergence of primate L1Pa7-L1Pa3 families. First, we experimentally demonstrated that while coiled coil function (measured by retrotransposition) can be eliminated by single epistatic mutations, it nonetheless can also withstand extensive amino acid substitutions. Second, principal component and cluster analysis showed that the coiled coil sequence space of each of the L1Pa7-3 families was notably increased by the presence of distinct, coexisting coiled coil sequences. Thus, sampling related networks of functional sequences rather than traversing discrete adaptive states characterized the persistence L1 activity during this evolutionary event.

A mathematical model of calcium dynamics: Obesity and mitochondria-associated ER membranes.

Multiple cellular organelles tightly orchestrate intracellular calcium (Ca2+) dynamics to regulate cellular activities and maintain homeostasis. The interplay between the endoplasmic reticulum (ER), a major store of intracellular Ca2+, and mitochondria, an important source of adenosine triphosphate (ATP), has been the subject of much research, as their dysfunction has been linked with metabolic diseases. Interestingly, throughout the cell's cytosolic domain, these two organelles share common microdomains called mitochondria-associated ER membranes (MAMs), where their membranes are in close apposition. The role of MAMs is critical for intracellular Ca2+ dynamics as they provide hubs for direct Ca2+ exchange between the organelles. A recent experimental study reported correlation between obesity and MAM formation in mouse liver cells, and obesity-related cellular changes that are closely associated with the regulation of Ca2+ dynamics. We constructed a mathematical model to study the effects of MAM Ca2+ dynamics on global Ca2+ activities. Through a series of model simulations, we investigated cellular mechanisms underlying the altered Ca2+ dynamics in the cells under obesity. We predict that, as the dosage of stimulus gradually increases, liver cells from obese mice will reach the state of saturated cytosolic Ca2+ concentration at a lower stimulus concentration, compared to cells from healthy mice.

Block network mapping approach to quantitative trait locus analysis.

BACKGROUND: Advances in experimental biology have enabled the collection of enormous troves of data on genomic variation in living organisms. The interpretation of this data to extract actionable information is one of the keys to developing novel therapeutic strategies to treat complex diseases. Network organization of biological data overcomes measurement noise in several biological contexts. Does a network approach, combining information about the linear organization of genomic markers with correlative information on these markers in a Bayesian formulation, lead to an analytic method with higher power for detecting quantitative trait loci? RESULTS: Block Network Mapping, combining Similarity Network Fusion (Wang et al., NM 11:333-337, 2014) with a Bayesian locus likelihood evaluation, leads to large improvements in area under the receiver operating characteristic and power over interval mapping with expectation maximization. The method has a monotonically decreasing false discovery rate as a function of effect size, unlike interval mapping. CONCLUSIONS: Block Network Mapping provides an alternative data-driven approach to mapping quantitative trait loci that leverages correlations in the sampled genotypes. The evaluation methodology can be combined with existing approaches such as Interval Mapping. Python scripts are available at http://lbm.niddk.nih.gov/vipulp/ . Genotype data is available at http://churchill-lab.jax.org/website/GattiDOQTL .

Metabolic scaling predicts posthepatectomy liver regeneration after accounting for hepatocyte hypertrophy.

We adapted a mathematical model of posthepatectomy liver regeneration using data from a subset of patients in the Adult-to-Adult Living Donor Liver Transplantation Cohort Study. The original model addressed changes in the number of quiescent, primed, and proliferating cells. Our adapted model takes into account hypertrophy of primed and replicating cells, and it is better able to predict liver volume. In addition, by building off the hypothesis that cell cycle parameters are approximately the same across all mammals, we found that changing only a single parameter characterizing metabolic load could model liver regeneration in 5 species of mammals. In conclusion, we improved a mathematical model of liver regeneration, predicted mammalian liver regeneration based on metabolism, and found correlations between model parameters and physiological measurements from liver donors.

Adaptation of pancreatic islet cyto-architecture during development.

Plasma glucose in mammals is regulated by hormones secreted by the islets of Langerhans embedded in the exocrine pancreas. Islets consist of endocrine cells, primarily α, ß, and δ cells, which secrete glucagon, insulin, and somatostatin, respectively. ß cells form irregular locally connected clusters within islets that act in concert to secrete insulin upon glucose stimulation. Varying demands and available nutrients during development produce changes in the local connectivity of ß cells in an islet. We showed in earlier work that graph theory provides a framework for the quantification of the seemingly stochastic cyto-architecture of ß cells in an islet. To quantify the dynamics of endocrine connectivity during development requires a framework for characterizing changes in the probability distribution on the space of possible graphs, essentially a Fokker-Planck formalism on graphs. With large-scale imaging data for hundreds of thousands of islets containing millions of cells from human specimens, we show that this dynamics can be determined quantitatively. Requiring that rearrangement and cell addition processes match the observed dynamic developmental changes in quantitative topological graph characteristics strongly constrained possible processes. Our results suggest that there is a transient shift in preferred connectivity for ß cells between 1-35 weeks and 12-24 months.

The Beta Cell in Its Cluster: Stochastic Graphs of Beta Cell Connectivity in the Islets of Langerhans.

Pancreatic islets of Langerhans consist of endocrine cells, primarily α, ß and δ cells, which secrete glucagon, insulin, and somatostatin, respectively, to regulate plasma glucose. ß cells form irregular locally connected clusters within islets that act in concert to secrete insulin upon glucose stimulation. Due to the central functional significance of this local connectivity in the placement of ß cells in an islet, it is important to characterize it quantitatively. However, quantification of the seemingly stochastic cytoarchitecture of ß cells in an islet requires mathematical methods that can capture topological connectivity in the entire ß-cell population in an islet. Graph theory provides such a framework. Using large-scale imaging data for thousands of islets containing hundreds of thousands of cells in human organ donor pancreata, we show that quantitative graph characteristics differ between control and type 2 diabetic islets. Further insight into the processes that shape and maintain this architecture is obtained by formulating a stochastic theory of ß-cell rearrangement in whole islets, just as the normal equilibrium distribution of the Ornstein-Uhlenbeck process can be viewed as the result of the interplay between a random walk and a linear restoring force. Requiring that rearrangements maintain the observed quantitative topological graph characteristics strongly constrained possible processes. Our results suggest that ß-cell rearrangement is dependent on its connectivity in order to maintain an optimal cluster size in both normal and T2D islets.

Modeling glucose and free fatty acid kinetics in glucose and meal tolerance test.

BACKGROUND: Quantitative evaluation of insulin regulation on plasma glucose and free fatty acid (FFA) in response to external glucose challenge is clinically important to assess the development of insulin resistance (World J Diabetes 1:36-47, 2010). Mathematical minimal models (MMs) based on insulin modified frequently-sampled intravenous glucose tolerance tests (IM-FSIGT) are widely applied to ascertain an insulin sensitivity index (IEEE Rev Biomed Eng 2:54-96, 2009). Furthermore, it is important to investigate insulin regulation on glucose and FFA in postprandial state as a normal physiological condition. A simple way to calculate the appearance rate (Ra) of glucose and FFA would be especially helpful to evaluate glucose and FFA kinetics for clinical applications. METHODS: A new MM is developed to simulate the insulin modulation of plasma glucose and FFA, combining IM-FSIGT with a mixed meal tolerance test (MT). A novel simple functional form for the appearance rate (Ra) of glucose or FFA in the MT is developed. Model results are compared with two other models for data obtained from 28 non-diabetic women (13 African American, 15 white). RESULTS: The new functional form for Ra of glucose is an acceptable empirical approximation to the experimental Ra for a subset of individuals. When both glucose and FFA are included in FSIGT and MT, the new model is preferred using the Bayes Information Criterion (BIC). CONCLUSIONS: Model simulations show that the new MM allows consistent application to both IM-FSIGT and MT data, balancing model complexity and data fitting. While the appearance of glucose in the circulation has an important effect on FFA kinetics in MT, the rate of appearance of FFA can be neglected for the time-period modeled.

Biochemical and clinical deficiency is uncommon in African immigrants despite a high prevalence of low vitamin D: the Africans in America study.

UNLABELLED: African ancestry is associated with low vitamin D levels but high bone density. Fifty percent of African immigrants had low vitamin D levels, but <10 % had evidence of deficiency. The value of providing vitamin D supplementation to African immigrants without evidence of deficiency needs to be determined. INTRODUCTION: The Endocrine Society and Institute of Medicine (IOM) have concluded from studies in largely white populations that 25(OH)D is necessary for bone health. However, their definition of vitamin D insufficiency differs. The Endocrine Society recommends a 25(OH)D threshold of <30 ng/mL. The IOM uses a lower threshold of 25(OH)D of <20 ng/mL. As African ancestry is associated with decreased 25(OH)D but increased bone mineral density (BMD), the applicability of these thresholds to Africans is unknown. Therefore, we examined in African immigrants the relationship of 25(OH)D to parathyroid hormone (PTH) and BMD. METHODS: One hundred eighty-six African immigrants(69 % male, age 38 ± 10 (mean ± SD), range 20-64 years) living in metropolitan Washington, DC, were enrolled. BMD was determined from whole-body dual-energy X-ray absorptiometry (DXA) scans. Decreased BMD required T-scores ≤-1.0. The threshold for low 25(OH)D was the concentration of 25(OH)D at which PTH became suppressed. This is known as the inflection point. Biochemical deficiency required low 25(OH)D and PTH of >65 pg/mL. Clinical deficiency required low 25(OH)D and T-scores ≤-1.0. RESULTS: 25(OH)D <30 and <20 ng/mL occurred in 83 and 46 % of African immigrants, respectively. PTH inversely correlated with 25(OH)D (r = -0.31, P = 0.002). The inflection point occurred at a 25(OH)D concentration of 20 ng/mL. Biochemical and clinical deficiency occurred in only 8 and 3 % of immigrants, respectively. CONCLUSION: As PTH became suppressed at 25(OH)D of 20 ng/mL, the 25(OH)D <20 ng/mL threshold for insufficiency may apply to African immigrants. However, ~50 % of African immigrants have 25(OH)D <20 ng/mL, but only <10 % had evidence of deficiency. The value of providing vitamin D supplementation to the large number of African immigrants with 25(OH)D <20 ng/mL and no detectable evidence of deficiency needs to be determined.

The jigsaw puzzle of sequence phenotype inference: Piecing together Shannon entropy, importance sampling, and Empirical Bayes.

A nucleotide sequence 35 base pairs long can take 1,180,591,620,717,411,303,424 possible values. An example of systems biology datasets, protein binding microarrays, contain activity data from about 40,000 such sequences. The discrepancy between the number of possible configurations and the available activities is enormous. Thus, albeit that systems biology datasets are large in absolute terms, they oftentimes require methods developed for rare events due to the combinatorial increase in the number of possible configurations of biological systems. A plethora of techniques for handling large datasets, such as Empirical Bayes, or rare events, such as importance sampling, have been developed in the literature, but these cannot always be simultaneously utilized. Here we introduce a principled approach to Empirical Bayes based on importance sampling, information theory, and theoretical physics in the general context of sequence phenotype model induction. We present the analytical calculations that underlie our approach. We demonstrate the computational efficiency of the approach on concrete examples, and demonstrate its efficacy by applying the theory to publicly available protein binding microarray transcription factor datasets and to data on synthetic cAMP-regulated enhancer sequences. As further demonstrations, we find transcription factor binding motifs, predict the activity of new sequences and extract the locations of transcription factor binding sites. In summary, we present a novel method that is efficient (requiring minimal computational time and reasonable amounts of memory), has high predictive power that is comparable with that of models with hundreds of parameters, and has a limited number of optimized parameters, proportional to the sequence length.

A network characteristic that correlates environmental and genetic robustness.

As scientific advances in perturbing biological systems and technological advances in data acquisition allow the large-scale quantitative analysis of biological function, the robustness of organisms to both transient environmental stresses and inter-generational genetic changes is a fundamental impediment to the identifiability of mathematical models of these functions. An approach to overcoming this impediment is to reduce the space of possible models to take into account both types of robustness. However, the relationship between the two is still controversial. This work uncovers a network characteristic, transient responsiveness, for a specific function that correlates environmental imperturbability and genetic robustness. We test this characteristic extensively for dynamic networks of ordinary differential equations ranging up to 30 interacting nodes and find that there is a power-law relating environmental imperturbability and genetic robustness that tends to linearity as the number of nodes increases. Using our methods, we refine the classification of known 3-node motifs in terms of their environmental and genetic robustness. We demonstrate our approach by applying it to the chemotaxis signaling network. In particular, we investigate plausible models for the role of CheV protein in biochemical adaptation via a phosphorylation pathway, testing modifications that could improve the robustness of the system to environmental and/or genetic perturbation.

Dory: Computation of persistence diagrams up to dimension two for Vietoris-Rips filtrations of large data sets.

Persistent homology (PH) is an approach to topological data analysis (TDA) that computes multi-scale topologically invariant properties of high-dimensional data that are robust to noise. While PH has revealed useful patterns across various applications, computational requirements have limited applications to small data sets of a few thousand points. We present Dory, an efficient and scalable algorithm that can compute the persistent homology of sparse Vietoris-Rips complexes on larger data sets, up to and including dimension two and over the field Z2. As an application, we compute the PH of the human genome at high resolution as revealed by a genome-wide Hi-C data set containing approximately three million points. Extant algorithms were unable to process it, whereas Dory processed it within five minutes, using less than five GB of memory. Results show that the topology of the human genome changes significantly upon treatment with auxin, a molecule that degrades cohesin, corroborating the hypothesis that cohesin plays a crucial role in loop formation in DNA.

Annealing approach to root finding.

The Newton-Raphson method is a fundamental root-finding technique with numerous applications in physics. In this study, we propose a parameterized variant of the Newton-Raphson method, inspired by principles from physics. Through analytical and empirical validation, we demonstrate that this approach offers increased robustness and faster convergence during root-finding iterations. Furthermore, we establish connections to the Adomian series method and provide a natural interpretation within a series framework. Remarkably, the introduced parameter, akin to a temperature variable, enables an annealing approach. This advancement sets the stage for a fresh exploration of numerical iterative root-finding methodologies.

Annealing approach to root-finding.

The Newton-Raphson method is a fundamental root-finding technique with numerous applications in physics. In this study, we propose a parameterized variant of the Newton-Raphson method, inspired by principles from physics. Through analytical and empirical validation, we demonstrate that this novel approach offers increased robustness and faster convergence during root-finding iterations. Furthermore, we establish connections to the Adomian series method and provide a natural interpretation within a series framework. Remarkably, the introduced parameter, akin to a temperature variable, enables an annealing approach. This advancement sets the stage for a fresh exploration of numerical iterative root-finding methodologies.

Synergy in free radical generation is blunted by high-fat diet induced alterations in skeletal muscle mitochondrial metabolism.

Due to their role in cellular energetics and metabolism, skeletal muscle mitochondria appear to play a key role in the development of insulin resistance and type II diabetes. High-fat diet can induce higher levels of reactive oxygen species (ROS), evidenced by hydrogen peroxide (H2O2) emission from mitochondria, which may be causal for insulin resistance in skeletal muscle. The underlying mechanisms are unclear. Recent published data on single substrate (pyruvate, succinate, fat) metabolism in both normal diet (CON) and high-fat diet (HFD) states of skeletal muscle allowed us to develop an integrated mathematical model of skeletal muscle mitochondrial metabolism. Model simulations suggested that long-term HFD may affect specific metabolic reaction/pathways by altering enzyme activities. Our model allows us to predict oxygen consumption and ROS generation for any combination of substrates. In particular, we predict a synergy between (iso-membrane potential) combinations of pyruvate and fat in ROS production compared to the sum of ROS production with each substrate singly in both CON and HFD states. This synergy is blunted in the HFD state.

The fractal spatial distribution of pancreatic islets in three dimensions: a self-avoiding growth model.

The islets of Langerhans, responsible for controlling blood glucose levels, are dispersed within the pancreas. A universal power law governing the fractal spatial distribution of islets in two-dimensional pancreatic sections has been reported. However, the fractal geometry in the actual three-dimensional pancreas volume, and the developmental process that gives rise to such a self-similar structure, has not been investigated. Here, we examined the three-dimensional spatial distribution of islets in intact mouse pancreata using optical projection tomography and found a power law with a fractal dimension of 2.1. Furthermore, based on two-dimensional pancreatic sections of human autopsies, we found that the distribution of human islets also follows a universal power law with a fractal dimension of 1.5 in adult pancreata, which agrees with the value previously reported in smaller mammalian pancreas sections. Finally, we developed a self-avoiding growth model for the development of the islet distribution and found that the fractal nature of the spatial islet distribution may be associated with the self-avoidance in the branching process of vascularization in the pancreas.

Expression in aneuploid Drosophila S2 cells.

Extensive departures from balanced gene dose in aneuploids are highly deleterious. However, we know very little about the relationship between gene copy number and expression in aneuploid cells. We determined copy number and transcript abundance (expression) genome-wide in Drosophila S2 cells by DNA-Seq and RNA-Seq. We found that S2 cells are aneuploid for >43 Mb of the genome, primarily in the range of one to five copies, and show a male genotype ( approximately two X chromosomes and four sets of autosomes, or 2X;4A). Both X chromosomes and autosomes showed expression dosage compensation. X chromosome expression was elevated in a fixed-fold manner regardless of actual gene dose. In engineering terms, the system "anticipates" the perturbation caused by X dose, rather than responding to an error caused by the perturbation. This feed-forward regulation resulted in precise dosage compensation only when X dose was half of the autosome dose. Insufficient compensation occurred at lower X chromosome dose and excessive expression occurred at higher doses. RNAi knockdown of the Male Specific Lethal complex abolished feed-forward regulation. Both autosome and X chromosome genes show Male Specific Lethal-independent compensation that fits a first order dose-response curve. Our data indicate that expression dosage compensation dampens the effect of altered DNA copy number genome-wide. For the X chromosome, compensation includes fixed and dose-dependent components.

Longitudinal multi-omics analyses of the gut-liver axis reveals metabolic dysregulation in hepatitis C infection and cirrhosis.

The gut and liver are connected via the portal vein, and this relationship, which includes the gut microbiome, is described as the gut-liver axis. Hepatitis C virus (HCV) can infect the liver and cause fibrosis with chronic infection. HCV has been associated with an altered gut microbiome; however, how these changes impact metabolism across the gut-liver axis and how this varies with disease severity and time is unclear. Here we used multi-omics analysis of portal and peripheral blood, faeces and liver tissue to characterize the gut-liver axis of patients with HCV across a fibrosis severity gradient before (n = 29) and 6 months after (n = 23) sustained virologic response, that is, no detection of the virus. Fatty acids were the major metabolites perturbed across the liver, portal vein and gut microbiome in HCV, especially in patients with cirrhosis. Decreased fatty acid degradation by hepatic peroxisomes and mitochondria was coupled with increased free fatty acid (FFA) influx to the liver via the portal vein. Metatranscriptomics indicated that Anaerostipes hadrus-mediated fatty acid synthesis influences portal FFAs. Both microbial fatty acid synthesis and portal FFAs were associated with enhanced hepatic fibrosis. Bacteroides vulgatus-mediated intestinal glycan breakdown was linked to portal glycan products, which in turn correlated with enhanced portal inflammation in HCV. Paired comparison of patient samples at both timepoints showed that hepatic metabolism, especially in peroxisomes, is persistently dysregulated in cirrhosis independently of the virus. Sustained virologic response was associated with a potential beneficial role for Methanobrevibacter smithii, which correlated with liver disease severity markers. These results develop our understanding of the gut-liver axis in HCV and non-HCV liver disease aetiologies and provide a foundation for future therapies.

Bayesian functional integral method for inferring continuous data from discrete measurements.

Inference of the insulin secretion rate (ISR) from C-peptide measurements as a quantification of pancreatic ß-cell function is clinically important in diseases related to reduced insulin sensitivity and insulin action. ISR derived from C-peptide concentration is an example of nonparametric Bayesian model selection where a proposed ISR time-course is considered to be a "model". An inferred value of inaccessible continuous variables from discrete observable data is often problematic in biology and medicine, because it is a priori unclear how robust the inference is to the deletion of data points, and a closely related question, how much smoothness or continuity the data actually support. Predictions weighted by the posterior distribution can be cast as functional integrals as used in statistical field theory. Functional integrals are generally difficult to evaluate, especially for nonanalytic constraints such as positivity of the estimated parameters. We propose a computationally tractable method that uses the exact solution of an associated likelihood function as a prior probability distribution for a Markov-chain Monte Carlo evaluation of the posterior for the full model. As a concrete application of our method, we calculate the ISR from actual clinical C-peptide measurements in human subjects with varying degrees of insulin sensitivity. Our method demonstrates the feasibility of functional integral Bayesian model selection as a practical method for such data-driven inference, allowing the data to determine the smoothing timescale and the width of the prior probability distribution on the space of models. In particular, our model comparison method determines the discrete time-step for interpolation of the unobservable continuous variable that is supported by the data. Attempts to go to finer discrete time-steps lead to less likely models.