Divergent Learning-Related Transcriptional States of Cortical Glutamatergic Neurons

Experience-dependent gene expression reshapes neural circuits, permitting the learning of knowledge and skills. Most learning involves repetitive experiences during which neurons undergo multiple stages of functional and structural plasticity. Currently, the diversity of transcriptional responses underlying dynamic plasticity during repetition-based learning is poorly understood. To close this gap, we analyzed single-nucleus transcriptomes of L2/3 glutamatergic neurons of the primary motor cortex after 3 d motor skill training or home cage control in water-restricted male mice. "Train" and "control" neurons could be discriminated with high accuracy based on expression patterns of many genes, indicating that recent experience leaves a widespread transcriptional signature across L2/3 neurons. These discriminating genes exhibited divergent modes of coregulation, differentiating neurons into discrete clusters of transcriptional states. Several states showed gene expressions associated with activity-dependent plasticity. Some of these states were also prominent in the previously published reference, suggesting that they represent both spontaneous and task-related plasticity events. Markedly, however, two states were unique to our dataset. The first state, further enriched by motor training, showed gene expression suggestive of late-stage plasticity with repeated activation, which is suitable for expected emergent neuronal ensembles that stably retain motor learning. The second state, equally found in both train and control mice, showed elevated levels of metabolic pathways and norepinephrine sensitivity, suggesting a response to common experiences specific to our experimental conditions, such as water restriction or circadian rhythm. Together, we uncovered divergent transcriptional responses across L2/3 neurons, each potentially linked with distinct features of repetition-based motor learning such as plasticity, memory, and motivation.

Cell and tissue-specific glycosylation pathways informed by single-cell transcriptomics

While single cell studies have made significant impacts in various subfields of biology, they lag in the Glycosciences. To address this gap, we analyzed single-cell glycogene expressions in the Tabula Sapiens dataset of human tissues and cell types using a recent glycosylation-specific gene ontology (GlycoEnzOnto). At the median sequencing (count) depth, ~40-50 out of 400 glycogenes were detected in individual cells. Upon increasing the sequencing depth, the number of detectable glycogenes saturates at ~200 glycogenes, suggesting that the average human cell expresses about half of the glycogene repertoire. Hierarchies in glycogene and glycopathway expressions emerged from our analysis: nucleotide-sugar synthesis and transport exhibited the highest gene expressions, followed by genes for core enzymes, glycan modification and extensions, and finally terminal modifications. Interestingly, the same cell types showed variable glycopathway expressions based on their organ or tissue origin, suggesting nuanced cell- and tissue-specific glycosylation patterns. Probing deeper into the transcription factors (TFs) of glycogenes, we identified distinct groupings of TFs controlling different aspects of glycosylation: core biosynthesis, terminal modifications, etc. We present webtools to explore the interconnections across glycogenes, glycopathways, and TFs regulating glycosylation in human cell/tissue types. Overall, the study presents an overview of glycosylation across multiple human organ systems.

Single-cell transcriptional uncertainty landscape of cell differentiation

Background: Single-cell studies have demonstrated the presence of significant cell-to-cell heterogeneity in gene expression. Whether such heterogeneity is only a bystander or has a functional role in the cell differentiation process is still hotly debated. Methods: In this study, we quantified and followed single-cell transcriptional uncertainty - a measure of gene transcriptional stochasticity in single cells - in 10 cell differentiation systems of varying cell lineage progressions, from single to multi-branching trajectories, using the stochastic two-state gene transcription model. Results: By visualizing the transcriptional uncertainty as a landscape over a two-dimensional representation of the single-cell gene expression data, we observed universal features in the cell differentiation trajectories that include: (i) a peak in single-cell uncertainty during transition states, and in systems with bifurcating differentiation trajectories, each branching point represents a state of high transcriptional uncertainty; (ii) a positive correlation of transcriptional uncertainty with transcriptional burst size and frequency; (iii) an increase in RNA velocity preceding the increase in the cell transcriptional uncertainty. Conclusions: Our findings suggest a possible universal mechanism during the cell differentiation process, in which stem cells engage stochastic exploratory dynamics of gene expression at the start of the cell differentiation by increasing gene transcriptional bursts, and disengage such dynamics once cells have decided on a particular terminal cell identity. Notably, the peak of single-cell transcriptional uncertainty signifies the decision-making point in the cell differentiation process.

GlycoEnzOnto: a GlycoEnzyme pathway and molecular function ontology

MOTIVATION

The 'glycoEnzymes' include a set of proteins having related enzymatic, metabolic, transport, structural and cofactor functions. Currently, there is no established ontology to describe glycoEnzyme properties and to relate them to glycan biosynthesis pathways.

RESULTS

We present GlycoEnzOnto, an ontology describing 403 human glycoEnzymes curated along 139 glycosylation pathways, 134 molecular functions and 22 cellular compartments. The pathways described regulate nucleotide-sugar metabolism, glycosyl-substrate/donor transport, glycan biosynthesis and degradation. The role of each enzyme in the glycosylation initiation, elongation/branching and capping/termination phases is described. IUPAC linear strings present systematic human/machine-readable descriptions of individual reaction steps and enable automated knowledge-based curation of biochemical networks. All GlycoEnzOnto knowledge is integrated with the Gene Ontology biological processes. GlycoEnzOnto enables improved transcript overrepresentation analyses and glycosylation pathway identification compared to other available schema, e.g. KEGG and Reactome. Overall, GlycoEnzOnto represents a holistic glycoinformatics resource for systems-level analyses.

AVAILABILITY AND IMPLEMENTATION

https://github.com/neel-lab/GlycoEnzOnto.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Global genome decompaction leads to stochastic activation of gene expression as a first step toward fate commitment in human hematopoietic cells

When human cord blood-derived CD34+ cells are induced to differentiate, they undergo rapid and dynamic morphological and molecular transformations that are critical for fate commitment. In particular, the cells pass through a transitory phase known as "multilineage-primed" state. These cells are characterized by a mixed gene expression profile, different in each cell, with the coexpression of many genes characteristic for concurrent cell lineages. The aim of our study is to understand the mechanisms of the establishment and the exit from this transitory state. We investigated this issue using single-cell RNA sequencing and ATAC-seq. Two phases were detected. The first phase is a rapid and global chromatin decompaction that makes most of the gene promoters in the genome accessible for transcription. It results 24 h later in enhanced and pervasive transcription of the genome leading to the concomitant increase in the cell-to-cell variability of transcriptional profiles. The second phase is the exit from the multilineage-primed phase marked by a slow chromatin closure and a subsequent overall down-regulation of gene transcription. This process is selective and results in the emergence of coherent expression profiles corresponding to distinct cell subpopulations. The typical time scale of these events spans 48 to 72 h. These observations suggest that the nonspecificity of genome decompaction is the condition for the generation of a highly variable multilineage expression profile. The nonspecific phase is followed by specific regulatory actions that stabilize and maintain the activity of key genes, while the rest of the genome becomes repressed again by the chromatin recompaction. Thus, the initiation of differentiation is reminiscent of a constrained optimization process that associates the spontaneous generation of gene expression diversity to subsequent regulatory actions that maintain the activity of some genes, while the rest of the genome sinks back to the repressive closed chromatin state.

Deep learning for neural decoding in motor cortex

Objective. Neural decoding is an important tool in neural engineering and neural data analysis. Of various machine learning algorithms adopted for neural decoding, the recently introduced deep learning is promising to excel. Therefore, we sought to apply deep learning to decode movement trajectories from the activity of motor cortical neurons.Approach. In this paper, we assessed the performance of deep learning methods in three different decoding schemes, concurrent, time-delay, and spatiotemporal. In the concurrent decoding scheme where the input to the network is the neural activity coincidental to the movement, deep learning networks including artificial neural network (ANN) and long-short term memory (LSTM) were applied to decode movement and compared with traditional machine learning algorithms. Both ANN and LSTM were further evaluated in the time-delay decoding scheme in which temporal delays are allowed between neural signals and movements. Lastly, in the spatiotemporal decoding scheme, we trained convolutional neural network (CNN) to extract movement information from images representing the spatial arrangement of neurons, their activity, and connectomes (i.e. the relative strengths of connectivity between neurons) and combined CNN and ANN to develop a hybrid spatiotemporal network. To reveal the input features of the CNN in the hybrid network that deep learning discovered for movement decoding, we performed a sensitivity analysis and identified specific regions in the spatial domain.Main results. Deep learning networks (ANN and LSTM) outperformed traditional machine learning algorithms in the concurrent decoding scheme. The results of ANN and LSTM in the time-delay decoding scheme showed that including neural data from time points preceding movement enabled decoders to perform more robustly when the temporal relationship between the neural activity and movement dynamically changes over time. In the spatiotemporal decoding scheme, the hybrid spatiotemporal network containing the concurrent ANN decoder outperformed single-network concurrent decoders.Significance. Taken together, our study demonstrates that deep learning could become a robust and effective method for the neural decoding of behavior.

Metabolic adaptation of ovarian tumors in patients treated with an IDO1 inhibitor constrains antitumor immune responses

To uncover underlying mechanisms associated with failure of indoleamine 2,3-dioxygenase 1 (IDO1) blockade in clinical trials, we conducted a pilot, window-of-opportunity clinical study in 17 patients with newly diagnosed advanced high-grade serous ovarian cancer before their standard tumor debulking surgery. Patients were treated with the IDO1 inhibitor epacadostat, and immunologic, transcriptomic, and metabolomic characterization of the tumor microenvironment was undertaken in baseline and posttreatment tumor biopsies. IDO1 inhibition resulted in efficient blockade of the kynurenine pathway of tryptophan degradation and was accompanied by a metabolic adaptation that shunted tryptophan catabolism toward the serotonin pathway. This resulted in elevated nicotinamide adenine dinucleotide (NAD⁺), which reduced T cell proliferation and function. Because NAD⁺ metabolites could be ligands for purinergic receptors, we investigated the impact of blocking purinergic receptors in the presence or absence of NAD⁺ on T cell proliferation and function in our mouse model. We demonstrated that A2a and A2b purinergic receptor antagonists, SCH58261 or PSB1115, respectively, rescued NAD⁺-mediated suppression of T cell proliferation and function. Combining IDO1 inhibition and A2a/A2b receptor blockade improved survival and boosted the antitumor immune signature in mice with IDO1 overexpressing ovarian cancer. These findings elucidate the downstream adaptive metabolic consequences of IDO1 blockade in ovarian cancers that may undermine antitumor T cell responses in the tumor microenvironment.

Integrating image analysis with single cell RNA-seq data to study podocyte-specific changes in diabetic kidney disease

Podocyte injury plays a crucial role in the progression of diabetic kidney disease (DKD). Injured podocytes demonstrate variations in nuclear shape and chromatin distribution. These morphometric changes have not yet been quantified in podocytes. Furthermore, the molecular mechanisms underlying these variations are poorly understood. Recent advances in omics have shed new lights into the biological mechanisms behind podocyte injury. However, there currently exists no study analyzing the biological mechanisms underlying podocyte morphometric variations during DKD. First, to study the importance of nuclear morphometrics, we performed morphometric quantification of podocyte nuclei from whole slide images of renal tissue sections obtained from murine models of DKD. Our results indicated that podocyte nuclear textural features demonstrate statistically significant difference in diabetic podocytes when compared to control. Additionally, the morphometric features demonstrated the existence of multiple subpopulations of podocytes suggesting a potential cause for their varying response to injury. Second, to study the underlying pathophysiology, we employed single cell RNA sequencing data from the murine models. Our results again indicated five subpopulations of podocytes in control and diabetic mouse models, validating the morphometrics-based results. Additionally, gene set enrichment analysis revealed epithelial to mesenchymal transition and apoptotic pathways in a subgroup of podocytes exclusive to diabetic mice, suggesting the molecular mechanism behind injury. Lastly, our results highlighted two distinct lineages of podocytes in control and diabetic cases suggesting a phenotypical change in podocytes during DKD. These results suggest that textural variations in podocyte nuclei may be key to understanding the pathophysiology behind podocyte injury.

FARCI: Fast and Robust Connectome Inference

The inference of neuronal connectome from large-scale neuronal activity recordings, such as two-photon Calcium imaging, represents an active area of research in computational neuroscience. In this work, we developed FARCI (Fast and Robust Connectome Inference), a MATLAB package for neuronal connectome inference from high-dimensional two-photon Calcium fluorescence data. We employed partial correlations as a measure of the functional association strength between pairs of neurons to reconstruct a neuronal connectome. We demonstrated using in silico datasets from the Neural Connectomics Challenge (NCC) and those generated using the state-of-the-art simulator of Neural Anatomy and Optimal Microscopy (NAOMi) that FARCI provides an accurate connectome and its performance is robust to network sizes, missing neurons, and noise levels. Moreover, FARCI is computationally efficient and highly scalable to large networks. In comparison with the best performing connectome inference algorithm in the NCC, Generalized Transfer Entropy (GTE), and Fluorescence Single Neuron and Network Analysis Package (FluoroSNNAP), FARCI produces more accurate networks over different network sizes, while providing significantly better computational speed and scaling.

ΔFBA-Predicting metabolic flux alterations using genome-scale metabolic models and differential transcriptomic data

Genome-scale metabolic models (GEMs) provide a powerful framework for simulating the entire set of biochemical reactions in a cell using a constraint-based modeling strategy called flux balance analysis (FBA). FBA relies on an assumed metabolic objective for generating metabolic fluxes using GEMs. But, the most appropriate metabolic objective is not always obvious for a given condition and is likely context-specific, which often complicate the estimation of metabolic flux alterations between conditions. Here, we propose a new method, called ΔFBA (deltaFBA), that integrates differential gene expression data to evaluate directly metabolic flux differences between two conditions. Notably, ΔFBA does not require specifying the cellular objective. Rather, ΔFBA seeks to maximize the consistency and minimize inconsistency between the predicted flux differences and differential gene expression. We showcased the performance of ΔFBA through several case studies involving the prediction of metabolic alterations caused by genetic and environmental perturbations in Escherichia coli and caused by Type-2 diabetes in human muscle. Importantly, in comparison to existing methods, ΔFBA gives a more accurate prediction of flux differences.

Metabolic alterations are often used as hallmarks of observable phenotypes. In this regard, reconstructed genome-scale metabolic models (GEMs) provide a rich and computable representation of the entire set of biochemical reactions in a cell. However, the performance of analytical tools for predicting metabolic reaction rates or fluxes using GEMs is sensitive to the assumed metabolic objective that is often unknown and likely context-specific. Here, we propose a novel method called ΔFBA that combines differential gene expression data and GEMs to evaluate differences in the metabolic fluxes between two conditions (perturbation vs. control) without the need for specifying a metabolic objective. In our demonstration, ΔFBA outperformed other existing methods in predicting metabolic flux alterations.

A systems-based framework to computationally describe putative transcription factors and signaling pathways regulating glycan biosynthesis

Glycosylation is a common posttranslational modification, and glycan biosynthesis is regulated by a set of glycogenes. The role of transcription factors (TFs) in regulating the glycogenes and related glycosylation pathways is largely unknown. In this work, we performed data mining of TF–glycogene relationships from the Cistrome Cancer database (DB), which integrates chromatin immunoprecipitation sequencing (ChIP-Seq) and RNA-Seq data to constitute regulatory relationships. In total, we observed 22,654 potentially significant TF–glycogene relationships, which include interactions involving 526 unique TFs and 341 glycogenes that span 29 the Cancer Genome Atlas (TCGA) cancer types. Here, TF–glycogene interactions appeared in clusters or so-called communities, suggesting that changes in single TF expression during both health and disease may affect multiple carbohydrate structures. Upon applying the Fisher’s exact test along with glycogene pathway classification, we identified TFs that may specifically regulate the biosynthesis of individual glycan types. Integration with Reactome DB knowledge provided an avenue to relate cell-signaling pathways to TFs and cellular glycosylation state. Whereas analysis results are presented for all 29 cancer types, specific focus is placed on human luminal and basal breast cancer disease progression. Overall, the article presents a computational approach to describe TF–glycogene relationships, the starting point for experimental system-wide validation.

Stochastic and Age-Dependent Proteostasis Decline Underlies Heterogeneity in Heat-Shock Response Dynamics

Significant non-genetic stochastic factors affect aging, causing lifespan differences among individuals, even those sharing the same genetic and environmental background. In Caenorhabditis elegans, differences in heat-shock response (HSR) are predictive of lifespan. However, factors contributing to the heterogeneity of HSR are still not fully elucidated. Here, the authors characterized HSR dynamics in isogenic C. elegans expressing GFP reporter for hsp-16.2 for identifying the key contributors of HSR heterogeneity. Specifically, microfluidic devices that enable cross-sectional and longitudinal measurements of HSR dynamics in C. elegans at different scales are developed: in populations, within individuals, and in embryos. The authors adapted a mathematical model of HSR to single C. elegans and identified model parameters associated with proteostasis-maintenance of protein homeostasis-more specifically, protein turnover, as the major drivers of heterogeneity in HSR dynamics. It is verified that individuals with enhanced proteostasis fidelity in early adulthood live longer. The model-based comparative analysis of protein turnover in day-1 and day-2 adult C. elegans revealed a stochastic-onset of age-related proteostasis decline that increases the heterogeneity of HSR capacity. Finally, the analysis of C. elegans embryos showed higher HSR and proteostasis capacity than young adults and established transgenerational contribution to HSR heterogeneity that depends on maternal age.

Thermodynamic analysis of DNA hybridization signatures near mitochondrial DNA deletion breakpoints

Broad evidence in the literature supports double-strand breaks (DSBs) as initiators of mitochondrial DNA (mtDNA) deletion mutations. While DNA misalignment during DSB repair is commonly proposed as the mechanism by which DSBs cause deletion mutations, details such as the specific DNA repair errors are still lacking. Here, we used DNA hybridization thermodynamics to infer the sequence lengths of mtDNA misalignments that are associated with mtDNA deletions. We gathered and analyzed 9,921 previously reported mtDNA deletion breakpoints in human, rhesus monkey, mouse, rat, and Caenorhabditis elegans. Our analysis shows that a large fraction of mtDNA breakpoint positions can be explained by the thermodynamics of short ≤ 5-nt misalignments. The significance of short DNA misalignments supports an important role for erroneous non-homologous and micro-homology-dependent DSB repair in mtDNA deletion formation. The consistency of the results of our analysis across species further suggests a shared mode of mtDNA deletion mutagenesis.

DeltaNeTS+: elucidating the mechanism of drugs and diseases using gene expression and transcriptional regulatory networks

BACKGROUND

Knowledge on the molecular targets of diseases and drugs is crucial for elucidating disease pathogenesis and mechanism of action of drugs, and for driving drug discovery and treatment formulation. In this regard, high-throughput gene transcriptional profiling has become a leading technology, generating whole-genome data on the transcriptional alterations caused by diseases or drug compounds. However, identifying direct gene targets, especially in the background of indirect (downstream) effects, based on differential gene expressions is difficult due to the complexity of gene regulatory network governing the gene transcriptional processes.

RESULTS

In this work, we developed a network analysis method, called DeltaNeTS+, for inferring direct gene targets of drugs and diseases from gene transcriptional profiles. DeltaNeTS+ uses a gene regulatory network model to identify direct perturbations to the transcription of genes using gene expression data. Importantly, DeltaNeTS+ is able to combine both steady-state and time-course expression profiles, as well as leverage information on the gene network structure. We demonstrated the power of DeltaNeTS+ in predicting gene targets using gene expression data in complex organisms, including Caenorhabditis elegans and human cell lines (T-cell and Calu-3). More specifically, in an application to time-course gene expression profiles of influenza A H1N1 (swine flu) and H5N1 (avian flu) infection, DeltaNeTS+ shed light on the key differences of dynamic cellular perturbations caused by the two influenza strains.

CONCLUSION

DeltaNeTS+ is a powerful network analysis tool for inferring gene targets from gene expression profiles. As demonstrated in the case studies, by incorporating available information on gene network structure, DeltaNeTS+ produces accurate predictions of direct gene targets from a small sample size (~ 10 s). Integrating static and dynamic expression data with transcriptional network structure extracted from genomic information, as enabled by DeltaNeTS+, is crucial toward personalized medicine, where treatments can be tailored to individual patients. DeltaNeTS+ can be freely downloaded from http://www.github.com/cabsel/deltanetsplus .

CALISTA: Clustering and LINEAGE Inference in Single-Cell Transcriptional Analysis

We present Clustering and Lineage Inference in Single-Cell Transcriptional Analysis (CALISTA), a numerically efficient and highly scalable toolbox for an end-to-end analysis of single-cell transcriptomic profiles. CALISTA includes four essential single-cell analyses for cell differentiation studies, including single-cell clustering, reconstruction of cell lineage specification, transition gene identification, and cell pseudotime ordering, which can be applied individually or in a pipeline. In these analyses, we employ a likelihood-based approach where single-cell mRNA counts are described by a probabilistic distribution function associated with stochastic gene transcriptional bursts and random technical dropout events. We illustrate the efficacy of CALISTA using single-cell gene expression datasets from different single-cell transcriptional profiling technologies and from a few hundreds to tens of thousands of cells. CALISTA is freely available on https://www.cabselab.com/calista.

Metabolic stress is a primary pathogenic event in transgenic Caenorhabditis elegans expressing pan-neuronal human amyloid beta

Alzheimer’s disease (AD) is the most common neurodegenerative disease affecting the elderly worldwide. Mitochondrial dysfunction has been proposed as a key event in the etiology of AD. We have previously modeled amyloid-beta (Aβ)-induced mitochondrial dysfunction in a transgenic Caenorhabditis elegans strain by expressing human Aβ peptide specifically in neurons (GRU102). Here, we focus on the deeper metabolic changes associated with this Aβ-induced mitochondrial dysfunction. Integrating metabolomics, transcriptomics and computational modeling, we identify alterations in Tricarboxylic Acid (TCA) cycle metabolism following even low-level Aβ expression. In particular, GRU102 showed reduced activity of a rate-limiting TCA cycle enzyme, alpha-ketoglutarate dehydrogenase. These defects were associated with elevation of protein carbonyl content specifically in mitochondria. Importantly, metabolic failure occurred before any significant increase in global protein aggregate was detectable. Treatment with an anti-diabetes drug, Metformin, reversed Aβ-induced metabolic defects, reduced protein aggregation and normalized lifespan of GRU102. Our results point to metabolic dysfunction as an early and causative event in Aβ-induced pathology and a promising target for intervention.

Network perturbation analysis of gene transcriptional profiles reveals protein targets and mechanism of action of drugs and influenza A viral infection

Genome-wide transcriptional profiling provides a global view of cellular state and how this state changes under different treatments (e.g. drugs) or conditions (e.g. healthy and diseased). Here, we present ProTINA (Protein Target Inference by Network Analysis), a network perturbation analysis method for inferring protein targets of compounds from gene transcriptional profiles. ProTINA uses a dynamic model of the cell-type specific protein-gene transcriptional regulation to infer network perturbations from steady state and time-series differential gene expression profiles. A candidate protein target is scored based on the gene network's dysregulation, including enhancement and attenuation of transcriptional regulatory activity of the protein on its downstream genes, caused by drug treatments. For benchmark datasets from three drug treatment studies, ProTINA was able to provide highly accurate protein target predictions and to reveal the mechanism of action of compounds with high sensitivity and specificity. Further, an application of ProTINA to gene expression profiles of influenza A viral infection led to new insights of the early events in the infection.

Constraints on Human CD34+ Cell Fate due to Lentiviral Vectors Can Be Relieved by Valproic Acid

The initial stages following the in vitro cytokine stimulation of human cord blood CD34+ cells overlap with the period when lentiviral gene transfer is typically performed. Single-cell transcriptional profiling and time-lapse microscopy were used to investigate how the vector-cell crosstalk impacts on the fate decision process. The single-cell transcription profiles were analyzed using a new algorithm, and it is shown that lentiviral transduction during the early stages of stimulation modifies the dynamics of the fate choice process of the CD34+ cells. The cells transduced with a lentiviral vector are biased toward the common myeloid progenitor lineage. Valproic acid, a histone deacetylase inhibitor known to increase the grafting potential of the CD34+ cells, improves the transduction efficiency to almost 100%. The cells transduced in the presence of valproic acid can subsequently undergo normal fate commitment. The higher gene transfer efficiency did not alter the genomic integration profile of the vector. These observations open the way to substantially improving lentiviral gene transfer protocols.

Clonal expansion of mitochondrial DNA deletions is a private mechanism of aging in long-lived animals

Disruption of mitochondrial metabolism and loss of mitochondrial DNA (mtDNA) integrity are widely considered as evolutionarily conserved (public) mechanisms of aging (López-Otín et al., Cell, 153, 2013 and 1194). Human aging is associated with loss in skeletal muscle mass and function (Sarcopenia), contributing significantly to morbidity and mortality. Muscle aging is associated with loss of mtDNA integrity. In humans, clonally expanded mtDNA deletions colocalize with sites of fiber breakage and atrophy in skeletal muscle. mtDNA deletions may therefore play an important, possibly causal role in sarcopenia. The nematode Caenorhabditis elegans also exhibits age-dependent decline in mitochondrial function and a form of sarcopenia. However, it is unclear if mtDNA deletions play a role in C. elegans aging. Here, we report identification of 266 novel mtDNA deletions in aging nematodes. Analysis of the mtDNA mutation spectrum and quantification of mutation burden indicates that (a) mtDNA deletions in nematode are extremely rare, (b) there is no significant age-dependent increase in mtDNA deletions, and (c) there is little evidence for clonal expansion driving mtDNA deletion dynamics. Thus, mtDNA deletions are unlikely to drive the age-dependent functional decline commonly observed in C. elegans. Computational modeling of mtDNA dynamics in C. elegans indicates that the lifespan of short-lived animals such as C. elegans is likely too short to allow for significant clonal expansion of mtDNA deletions. Together, these findings suggest that clonal expansion of mtDNA deletions is likely a private mechanism of aging predominantly relevant in long-lived animals such as humans and rhesus monkey and possibly in rodents.

SINCERITIES: inferring gene regulatory networks from time-stamped single cell transcriptional expression profiles

Motivation

Single cell transcriptional profiling opens up a new avenue in studying the functional role of cell-to-cell variability in physiological processes. The analysis of single cell expression profiles creates new challenges due to the distributive nature of the data and the stochastic dynamics of gene transcription process. The reconstruction of gene regulatory networks (GRNs) using single cell transcriptional profiles is particularly challenging, especially when directed gene-gene relationships are desired.

Results

We developed SINCERITIES (SINgle CEll Regularized Inference using TIme-stamped Expression profileS) for the inference of GRNs from single cell transcriptional profiles. We focused on time-stamped cross-sectional expression data, commonly generated from transcriptional profiling of single cells collected at multiple time points after cell stimulation. SINCERITIES recovers directed regulatory relationships among genes by employing regularized linear regression (ridge regression), using temporal changes in the distributions of gene expressions. Meanwhile, the modes of the gene regulations (activation and repression) come from partial correlation analyses between pairs of genes. We demonstrated the efficacy of SINCERITIES in inferring GRNs using in silico time-stamped single cell expression data and single cell transcriptional profiles of THP-1 monocytic human leukemia cells. The case studies showed that SINCERITIES could provide accurate GRN predictions, significantly better than other GRN inference algorithms such as TSNI, GENIE3 and JUMP3. Moreover, SINCERITIES has a low computational complexity and is amenable to problems of extremely large dimensionality. Finally, an application of SINCERITIES to single cell expression data of T2EC chicken erythrocytes pointed to BATF as a candidate novel regulator of erythroid development.

Availability and implementation

MATLAB and R version of SINCERITIES are freely available from the following websites: http://www.cabsel.ethz.ch/tools/sincerities.html and https://github.com/CABSEL/SINCERITIES. The single cell THP-1 and T2EC transcriptional profiles are available from the original publications (Kouno et al., 2013; Richard et al., 2016). The in silico single cell data are available on SINCERITIES websites.

Supplementary information

Supplementary data are available at Bioinformatics online.

Assessing and Resolving Model Misspecifications in Metabolic Flux Analysis

Metabolic flux analysis (MFA) is an indispensable tool in metabolic engineering. The simplest variant of MFA relies on an overdetermined stoichiometric model of the cell’s metabolism under the pseudo-steady state assumption to evaluate the intracellular flux distribution. Despite its long history, the issue of model error in overdetermined MFA, particularly misspecifications of the stoichiometric matrix, has not received much attention. We evaluated the performance of statistical tests from linear least square regressions, namely Ramsey’s Regression Equation Specification Error Test (RESET), the F-test, and the Lagrange multiplier test, in detecting model misspecifications in the overdetermined MFA, particularly missing reactions. We further proposed an iterative procedure using the F-test to correct such an issue. Using Chinese hamster ovary and random metabolic networks, we demonstrated that: (1) a statistically significant regression does not guarantee high accuracy of the flux estimates; (2) the removal of a reaction with a low flux magnitude can cause disproportionately large biases in the flux estimates; (3) the F-test could efficiently detect missing reactions; and (4) the proposed iterative procedure could robustly resolve the omission of reactions. Our work demonstrated that statistical analysis and tests could be used to systematically assess, detect, and resolve model misspecifications in the overdetermined MFA.

Bioprocess optimization under uncertainty using ensemble modeling

The performance of model-based bioprocess optimizations depends on the accuracy of the mathematical model. However, models of bioprocesses often have large uncertainty due to the lack of model identifiability. In the presence of such uncertainty, process optimizations that rely on the predictions of a single "best fit" model, e.g. the model resulting from a maximum likelihood parameter estimation using the available process data, may perform poorly in real life. In this study, we employed ensemble modeling to account for model uncertainty in bioprocess optimization. More specifically, we adopted a Bayesian approach to define the posterior distribution of the model parameters, based on which we generated an ensemble of model parameters using a uniformly distributed sampling of the parameter confidence region. The ensemble-based process optimization involved maximizing the lower confidence bound of the desired bioprocess objective (e.g. yield or product titer), using a mean-standard deviation utility function. We demonstrated the performance and robustness of the proposed strategy in an application to a monoclonal antibody batch production by mammalian hybridoma cell culture.

Single-Cell-Based Analysis Highlights a Surge in Cell-to-Cell Molecular Variability Preceding Irreversible Commitment in a Differentiation Process

In some recent studies, a view emerged that stochastic dynamics governing the switching of cells from one differentiation state to another could be characterized by a peak in gene expression variability at the point of fate commitment. We have tested this hypothesis at the single-cell level by analyzing primary chicken erythroid progenitors through their differentiation process and measuring the expression of selected genes at six sequential time-points after induction of differentiation. In contrast to population-based expression data, single-cell gene expression data revealed a high cell-to-cell variability, which was masked by averaging. We were able to show that the correlation network was a very dynamical entity and that a subgroup of genes tend to follow the predictions from the dynamical network biomarker (DNB) theory. In addition, we also identified a small group of functionally related genes encoding proteins involved in sterol synthesis that could act as the initial drivers of the differentiation. In order to assess quantitatively the cell-to-cell variability in gene expression and its evolution in time, we used Shannon entropy as a measure of the heterogeneity. Entropy values showed a significant increase in the first 8 h of the differentiation process, reaching a peak between 8 and 24 h, before decreasing to significantly lower values. Moreover, we observed that the previous point of maximum entropy precedes two paramount key points: an irreversible commitment to differentiation between 24 and 48 h followed by a significant increase in cell size variability at 48 h. In conclusion, when analyzed at the single cell level, the differentiation process looks very different from its classical population average view. New observables (like entropy) can be computed, the behavior of which is fully compatible with the idea that differentiation is not a “simple” program that all cells execute identically but results from the dynamical behavior of the underlying molecular network.

A single-cell transcriptomics analysis offers a new dynamical view of the differentiation process, involving an increase in between-cell variability prior to commitment.

The differentiation process has classically been seen as a stereotyped program leading from one progenitor toward a functional cell. This vision was based upon cell population-based analyses averaged over millions of cells. However, new methods have recently emerged that allow interrogation of the molecular content at the single-cell level, challenging this view with a new model suggesting that cell-to-cell gene expression stochasticity could play a key role in differentiation. We took advantage of a physiologically relevant avian cellular model to analyze the expression level of 92 genes in individual cells collected at several time-points during differentiation. We first observed that the process analyzed at the single-cell level is very different and much less well ordered than the population-based average view. Furthermore, we showed that cell-to-cell variability in gene expression peaks transiently before strongly decreasing. This rise in variability precedes two key events: an irreversible commitment to differentiation, followed by a significant increase in cell size variability. Altogether, our results support the idea that differentiation is not a “simple” series of well-ordered molecular events executed identically by all cells in a population but likely results from dynamical behavior of the underlying molecular network.

Phase Sensitivity Analysis of Circadian Rhythm Entrainment

As a biological clock, circadian rhythms evolve to accomplish a stable (robust) entrainment to environmental cycles, of which light is the most obvious. The mechanism of photic entrainment is not known, but two models of entrainment have been proposed based on whether light has a continuous (parametric) or discrete (nonparametric) effect on the circadian pacemaker. A novel sensitivity analysis is developed to study the circadian entrainment in silico based on a limit cycle approach and applied to a model of Drosophila circadian rhythm. The comparative analyses of complete and skeleton photoperiods suggest a trade-off between the contribution of period modulation (parametric effect) and phase shift (nonparametric effect) in Drosophila circadian entrainment. The results also give suggestions for an experimental study to (in)validate the two models of entrainment.

TRaCE+: Ensemble inference of gene regulatory networks from transcriptional expression profiles of gene knock-out experiments

Background

The inference of gene regulatory networks (GRNs) from transcriptional expression profiles is challenging, predominantly due to its underdetermined nature. One important consequence of underdetermination is the existence of many possible solutions to this inference. Our previously proposed ensemble inference algorithm TRaCE addressed this issue by inferring an ensemble of network directed graphs (digraphs) using differential gene expressions from gene knock-out (KO) experiments. However, TRaCE could not deal with the mode of the transcriptional regulations (activation or repression), an important feature of GRNs.

Results

In this work, we developed a new algorithm called TRaCE+ for the inference of an ensemble of signed GRN digraphs from transcriptional expression data of gene KO experiments. The sign of the edges indicates whether the regulation is an activation (positive) or a repression (negative). TRaCE+ generates the upper and lower bounds of the ensemble, which define uncertain regulatory interactions that could not be verified by the data. As demonstrated in the case studies using Escherichia coli GRN and 100-gene gold-standard GRNs from DREAM 4 network inference challenge, by accounting for regulatory signs, TRaCE+ could extract more information from the KO data than TRaCE, leading to fewer uncertain edges. Importantly, iterating TRaCE+ with an optimal design of gene KOs could resolve the underdetermined issue of GRN inference in much fewer KO experiments than using TRaCE.

Conclusions

TRaCE+ expands the applications of ensemble GRN inference strategy by accounting for the mode of the gene regulatory interactions. In comparison to TRaCE, TRaCE+ enables a better utilization of gene KO data, thereby reducing the cost of tackling underdetermined GRN inference. TRaCE+ subroutines for MATLAB are freely available at the following website: http://www.cabsel.ethz.ch/tools/trace.html.

Inferring gene targets of drugs and chemical compounds from gene expression profiles

Motivation: Finding genes which are directly perturbed or targeted by drugs is of great interest and importance in drug discovery. Several network filtering methods have been created to predict the gene targets of drugs from gene expression data based on an ordinary differential equation model of the gene regulatory network (GRN). A critical step in these methods involves inferring the GRN from the expression data, which is a very challenging problem on its own. In addition, existing network filtering methods require computationally intensive parameter tuning or expression data from experiments with known genetic perturbations or both.

Results: We developed a method called DeltaNet for the identification of drug targets from gene expression data. Here, the gene target predictions were directly inferred from the data without a separate step of GRN inference. DeltaNet formulation led to solving an underdetermined linear regression problem, for which we employed least angle regression (DeltaNet-LAR) or LASSO regularization (DeltaNet-LASSO). The predictions using DeltaNet for expression data of Escherichia coli, yeast, fruit fly and human were significantly more accurate than those using network filtering methods, namely mode of action by network identification (MNI) and sparse simultaneous equation model (SSEM). Furthermore, DeltaNet using LAR did not require any parameter tuning and could provide computational speed-up over existing methods.

Conclusion: DeltaNet is a robust and numerically efficient tool for identifying gene perturbations from gene expression data. Importantly, the method requires little to no expert supervision, while providing accurate gene target predictions.

Availability and implementation: DeltaNet is available on http://www.cabsel.ethz.ch/tools/DeltaNet.

Contact:rudi.gunawan@chem.ethz.ch

Supplementary information:Supplementary data are available at Bioinformatics online.

Optimal design of gene knockout experiments for gene regulatory network inference

MOTIVATION

We addressed the problem of inferring gene regulatory network (GRN) from gene expression data of knockout (KO) experiments. This inference is known to be underdetermined and the GRN is not identifiable from data. Past studies have shown that suboptimal design of experiments (DOE) contributes significantly to the identifiability issue of biological networks, including GRNs. However, optimizing DOE has received much less attention than developing methods for GRN inference.

RESULTS

We developed REDuction of UnCertain Edges (REDUCE) algorithm for finding the optimal gene KO experiment for inferring directed graphs (digraphs) of GRNs. REDUCE employed ensemble inference to define uncertain gene interactions that could not be verified by prior data. The optimal experiment corresponds to the maximum number of uncertain interactions that could be verified by the resulting data. For this purpose, we introduced the concept of edge separatoid which gave a list of nodes (genes) that upon their removal would allow the verification of a particular gene interaction. Finally, we proposed a procedure that iterates over performing KO experiments, ensemble update and optimal DOE. The case studies including the inference of Escherichia coli GRN and DREAM 4 100-gene GRNs, demonstrated the efficacy of the iterative GRN inference. In comparison to systematic KOs, REDUCE could provide much higher information return per gene KO experiment and consequently more accurate GRN estimates.

CONCLUSIONS

REDUCE represents an enabling tool for tackling the underdetermined GRN inference. Along with advances in gene deletion and automation technology, the iterative procedure brings an efficient and fully automated GRN inference closer to reality.

AVAILABILITY AND IMPLEMENTATION

MATLAB and Python scripts of REDUCE are available on www.cabsel.ethz.ch/tools/REDUCE CONTACT: rudi.gunawan@chem.ethz.ch

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

REDEMPTION: reduced dimension ensemble modeling and parameter estimation

Summary: Here, we present REDEMPTION (Reduced Dimension Ensemble Modeling and Parameter estimation), a toolbox for parameter estimation and ensemble modeling of ordinary differential equations (ODEs) using time-series data. For models with more reactions than measured species, a common scenario in biological modeling, the parameter estimation is formulated as a nested optimization problem based on incremental parameter estimation strategy. REDEMPTION also includes a tool for the identification of an ensemble of parameter combinations that provide satisfactory goodness-of-fit to the data. The functionalities of REDEMPTION are accessible through a MATLAB user interface (UI), as well as through programming script. For computational speed-up, REDEMPTION provides a numerical parallelization option using MATLAB Parallel Computing toolbox.

Availability and implementation: REDEMPTION can be downloaded from http://www.cabsel.ethz.ch/tools/redemption.

Contact:rudi.gunawan@chem.ethz.ch

Are mutagenic non D-loop direct repeat motifs in mitochondrial DNA under a negative selection pressure?

Non D-loop direct repeats (DRs) in mitochondrial DNA (mtDNA) have been commonly implicated in the mutagenesis of mtDNA deletions associated with neuromuscular disease and ageing. Further, these DRs have been hypothesized to put a constraint on the lifespan of mammals and are under a negative selection pressure. Using a compendium of 294 mammalian mtDNA, we re-examined the relationship between species lifespan and the mutagenicity of such DRs. Contradicting the prevailing hypotheses, we found no significant evidence that long-lived mammals possess fewer mutagenic DRs than short-lived mammals. By comparing DR counts in human mtDNA with those in selectively randomized sequences, we also showed that the number of DRs in human mtDNA is primarily determined by global mtDNA properties, such as the bias in synonymous codon usage (SCU) and nucleotide composition. We found that SCU bias in mtDNA positively correlates with DR counts, where repeated usage of a subset of codons leads to more frequent DR occurrences. While bias in SCU and nucleotide composition has been attributed to nucleotide mutational bias, mammalian mtDNA still exhibit higher SCU bias and DR counts than expected from such mutational bias, suggesting a lack of negative selection against non D-loop DRs.

Ensemble inference and inferability of gene regulatory networks

The inference of gene regulatory network (GRN) from gene expression data is an unsolved problem of great importance. This inference has been stated, though not proven, to be underdetermined implying that there could be many equivalent (indistinguishable) solutions. Motivated by this fundamental limitation, we have developed new framework and algorithm, called TRaCE, for the ensemble inference of GRNs. The ensemble corresponds to the inherent uncertainty associated with discriminating direct and indirect gene regulations from steady-state data of gene knock-out (KO) experiments. We applied TRaCE to analyze the inferability of random GRNs and the GRNs of E. coli and yeast from single- and double-gene KO experiments. The results showed that, with the exception of networks with very few edges, GRNs are typically not inferable even when the data are ideal (unbiased and noise-free). Finally, we compared the performance of TRaCE with top performing methods of DREAM4 in silico network inference challenge.

Assessment of network inference methods: how to cope with an underdetermined problem

The inference of biological networks is an active research area in the field of systems biology. The number of network inference algorithms has grown tremendously in the last decade, underlining the importance of a fair assessment and comparison among these methods. Current assessments of the performance of an inference method typically involve the application of the algorithm to benchmark datasets and the comparison of the network predictions against the gold standard or reference networks. While the network inference problem is often deemed underdetermined, implying that the inference problem does not have a (unique) solution, the consequences of such an attribute have not been rigorously taken into consideration. Here, we propose a new procedure for assessing the performance of gene regulatory network (GRN) inference methods. The procedure takes into account the underdetermined nature of the inference problem, in which gene regulatory interactions that are inferable or non-inferable are determined based on causal inference. The assessment relies on a new definition of the confusion matrix, which excludes errors associated with non-inferable gene regulations. For demonstration purposes, the proposed assessment procedure is applied to the DREAM 4 In Silico Network Challenge. The results show a marked change in the ranking of participating methods when taking network inferability into account.

Effects of lithium on age-related decline in mitochondrial turnover and function in Caenorhabditis elegans

Aging has been associated with the accumulation of damages in molecules and organelles in cells, particularly mitochondria. The rate of damage accumulation is closely tied to the turnover of the affected cellular components. Perturbing mitochondrial turnover has been shown to significantly affect the rate of deterioration of mitochondrial function with age and to alter lifespan of model organisms. In this study, we investigated the effects of upregulating autophagy using lithium in Caenorhabditis elegans. We found that lithium treatment increased both the lifespan and healthspan of C. elegans without any significant change in the mortality rate and oxidative damages to proteins. The increase in healthspan was accompanied by improved mitochondrial energetic function. In contrast, mitochondrial DNA copy number decreased faster with age under lithium. To better understand the interactions among mitochondrial turnover, damage, and function, we created a mathematical model that described the dynamics of functional and dysfunctional mitochondria population. The combined analysis of model and experimental observations showed how preferential (selective) autophagy of dysfunctional mitochondria could lead to better mitochondrial functionality with age, despite a lower population size. However, the results of model analysis suggest that the benefit of increasing autophagy for mitochondrial function is expected to diminish at higher levels of upregulation due to a shrinking mitochondrial population.

Mathematical modeling of the role of mitochondrial fusion and fission in mitochondrial DNA maintenance

Accumulation of mitochondrial DNA (mtDNA) mutations has been implicated in a wide range of human pathologies, including neurodegenerative diseases, sarcopenia, and the aging process itself. In cells, mtDNA molecules are constantly turned over (i.e. replicated and degraded) and are also exchanged among mitochondria during the fusion and fission of these organelles. While the expansion of a mutant mtDNA population is believed to occur by random segregation of these molecules during turnover, the role of mitochondrial fusion-fission in this context is currently not well understood. In this study, an in silico modeling approach is taken to investigate the effects of mitochondrial fusion and fission dynamics on mutant mtDNA accumulation. Here we report model simulations suggesting that when mitochondrial fusion-fission rate is low, the slow mtDNA mixing can lead to an uneven distribution of mutant mtDNA among mitochondria in between two mitochondrial autophagic events leading to more stochasticity in the outcomes from a single random autophagic event. Consequently, slower mitochondrial fusion-fission results in higher variability in the mtDNA mutation burden among cells in a tissue over time, and mtDNA mutations have a higher propensity to clonally expand due to the increased stochasticity. When these mutations affect cellular energetics, nuclear retrograde signalling can upregulate mtDNA replication, which is expected to slow clonal expansion of these mutant mtDNA. However, our simulations suggest that the protective ability of retrograde signalling depends on the efficiency of fusion-fission process. Our results thus shed light on the interplay between mitochondrial fusion-fission and mtDNA turnover and may explain the mechanism underlying the experimentally observed increase in the accumulation of mtDNA mutations when either mitochondrial fusion or fission is inhibited.

Incremental parameter estimation of kinetic metabolic network models

BACKGROUND

An efficient and reliable parameter estimation method is essential for the creation of biological models using ordinary differential equation (ODE). Most of the existing estimation methods involve finding the global minimum of data fitting residuals over the entire parameter space simultaneously. Unfortunately, the associated computational requirement often becomes prohibitively high due to the large number of parameters and the lack of complete parameter identifiability (i.e. not all parameters can be uniquely identified).

RESULTS

In this work, an incremental approach was applied to the parameter estimation of ODE models from concentration time profiles. Particularly, the method was developed to address a commonly encountered circumstance in the modeling of metabolic networks, where the number of metabolic fluxes (reaction rates) exceeds that of metabolites (chemical species). Here, the minimization of model residuals was performed over a subset of the parameter space that is associated with the degrees of freedom in the dynamic flux estimation from the concentration time-slopes. The efficacy of this method was demonstrated using two generalized mass action (GMA) models, where the method significantly outperformed single-step estimations. In addition, an extension of the estimation method to handle missing data is also presented.

CONCLUSIONS

The proposed incremental estimation method is able to tackle the issue on the lack of complete parameter identifiability and to significantly reduce the computational efforts in estimating model parameters, which will facilitate kinetic modeling of genome-scale cellular metabolism in the future.

Role of direct repeat and stem-loop motifs in mtDNA deletions: cause or coincidence?

Deletion mutations within mitochondrial DNA (mtDNA) have been implicated in degenerative and aging related conditions, such as sarcopenia and neuro-degeneration. While the precise molecular mechanism of deletion formation in mtDNA is still not completely understood, genome motifs such as direct repeat (DR) and stem-loop (SL) have been observed in the neighborhood of deletion breakpoints and thus have been postulated to take part in mutagenesis. In this study, we have analyzed the mitochondrial genomes from four different mammals: human, rhesus monkey, mouse and rat, and compared them to randomly generated sequences to further elucidate the role of direct repeat and stem-loop motifs in aging associated mtDNA deletions. Our analysis revealed that in the four species, DR and SL structures are abundant and that their distributions in mtDNA are not statistically different from randomized sequences. However, the average distance between the reported age associated mtDNA breakpoints and their respective nearest DR motifs is significantly shorter than what is expected of random chance in human (p<10⁻⁴) and rhesus monkey (p = 0.0034), but not in mouse (p = 0.0719) and rat (p = 0.0437), indicating the existence of species specific difference in the relationship between DR motifs and deletion breakpoints. In addition, the frequencies of large DRs (>10 bp) tend to decrease with increasing lifespan among the four mammals studied here, further suggesting an evolutionary selection against stable mtDNA misalignments associated with long DRs in long-living animals. In contrast to the results on DR, the probability of finding SL motifs near a deletion breakpoint does not differ from random in any of the four mtDNA sequences considered. Taken together, the findings in this study give support for the importance of stable mtDNA misalignments, aided by long DRs, as a major mechanism of deletion formation in long-living, but not in short-living mammals.

Maximizing signal-to-noise ratio in the random mutation capture assay

The ‘Random Mutation Capture’ assay allows for the sensitive quantitation of DNA mutations at extremely low mutation frequencies. This method is based on PCR detection of mutations that render the mutated target sequence resistant to restriction enzyme digestion. The original protocol prescribes an end-point dilution to about 0.1 mutant DNA molecules per PCR well, such that the mutation burden can be simply calculated by counting the number of amplified PCR wells. However, the statistical aspects associated with the single molecular nature of this protocol and several other molecular approaches relying on binary (on/off) output can significantly affect the quantification accuracy, and this issue has so far been ignored. The present work proposes a design of experiment (DoE) using statistical modeling and Monte Carlo simulations to obtain a statistically optimal sampling protocol, one that minimizes the coefficient of variance in the measurement estimates. Here, the DoE prescribed a dilution factor at about 1.6 mutant molecules per well. Theoretical results and experimental validation revealed an up to 10-fold improvement in the information obtained per PCR well, i.e. the optimal protocol achieves the same coefficient of variation using one-tenth the number of wells used in the original assay. Additionally, this optimization equally applies to any method that relies on binary detection of a small number of templates.

Parameter estimation of kinetic models from metabolic profiles: two-phase dynamic decoupling method

MOTIVATION

Time-series measurements of metabolite concentration have become increasingly more common, providing data for building kinetic models of metabolic networks using ordinary differential equations (ODEs). In practice, however, such time-course data are usually incomplete and noisy, and the estimation of kinetic parameters from these data is challenging. Practical limitations due to data and computational aspects, such as solving stiff ODEs and finding global optimal solution to the estimation problem, give motivations to develop a new estimation procedure that can circumvent some of these constraints.

RESULTS

In this work, an incremental and iterative parameter estimation method is proposed that combines and iterates between two estimation phases. One phase involves a decoupling method, in which a subset of model parameters that are associated with measured metabolites, are estimated using the minimization of slope errors. Another phase follows, in which the ODE model is solved one equation at a time and the remaining model parameters are obtained by minimizing concentration errors. The performance of this two-phase method was tested on a generic branched metabolic pathway and the glycolytic pathway of Lactococcus lactis. The results showed that the method is efficient in getting accurate parameter estimates, even when some information is missing.

Understanding dynamics using sensitivity analysis: caveat and solution

Background

Parametric sensitivity analysis (PSA) has become one of the most commonly used tools in computational systems biology, in which the sensitivity coefficients are used to study the parametric dependence of biological models. As many of these models describe dynamical behaviour of biological systems, the PSA has subsequently been used to elucidate important cellular processes that regulate this dynamics. However, in this paper, we show that the PSA coefficients are not suitable in inferring the mechanisms by which dynamical behaviour arises and in fact it can even lead to incorrect conclusions.

Results

A careful interpretation of parametric perturbations used in the PSA is presented here to explain the issue of using this analysis in inferring dynamics. In short, the PSA coefficients quantify the integrated change in the system behaviour due to persistent parametric perturbations, and thus the dynamical information of when a parameter perturbation matters is lost. To get around this issue, we present a new sensitivity analysis based on impulse perturbations on system parameters, which is named impulse parametric sensitivity analysis (iPSA). The inability of PSA and the efficacy of iPSA in revealing mechanistic information of a dynamical system are illustrated using two examples involving switch activation.

Conclusions

The interpretation of the PSA coefficients of dynamical systems should take into account the persistent nature of parametric perturbations involved in the derivation of this analysis. The application of PSA to identify the controlling mechanism of dynamical behaviour can be misleading. By using impulse perturbations, introduced at different times, the iPSA provides the necessary information to understand how dynamics is achieved, i.e. which parameters are essential and when they become important.

Elucidating cytochrome C release from mitochondria: insights from an in silico three-dimensional model

Mitochondrial regulation of apoptosis depends on the programmed release of proapoptotic proteins such as cytochrome c (Cyt c) through the outer mitochondrial membrane (OMM). Although a few key processes involved in this release have been identified, including the liberation of inner membrane-bound Cyt c and formation of diffusible pores on the OMM, other details like the transport of Cyt c within complex mitochondrial compartments, e.g., the cristae and crista junctions, are not yet fully understood (to our knowledge). In particular, a remodeling of the inner mitochondrial membrane accompanying apoptosis seen in a few studies, in which crista junctions widen, has been hypothesized to be a necessary step in the Cyt c release. Using a three-dimensional spatial modeling of mitochondrial crista and the crista junction, model simulations and analysis illustrated how the interplay among solubilization of Cyt c, fast diffusion of Cyt c, and OMM permeabilization gives rise to the observed experimental release profile. Importantly, the widening of the crista junction was found to have a negligible effect on the transport of free Cyt c from cristae. Finally, model simulations showed that increasing the fraction of free/loosely-bound Cyt c can sensitize the cell to apoptotic stimuli in a threshold manner, which may explain increased sensitivity to cell death associated with aging.

Impact of repeated nicotine and alcohol coexposure on in vitro and in vivo chlorpyrifos dosimetry and cholinesterase inhibition

Chlorpyrifos (CPF) is an organophosphorus insecticide, and neurotoxicity results from inhibition of acetylcholinesterase (AChE) by its metabolite, chlorpyrifos-oxon. Routine consumption of alcohol and tobacco modifies metabolic and physiological processes impacting the metabolism and pharmacokinetics of other xenobiotics, including pesticides. This study evaluated the influence of repeated ethanol and nicotine coexposure on in vivo CPF dosimetry and cholinesterase (ChE) response (ChE- includes AChE and/or butyrylcholinesterase (BuChE)). Hepatic microsomes were prepared from groups of naive, ethanol-only (1 g/kg/d, 7 d, po), and ethanol + nicotine (1 mg/kg/d 7 d, sc)-treated rats, and the in vitro metabolism of CPF was evaluated. For in vivo studies, rats were treated with saline or ethanol (1 g/kg/d, po) + nicotine (1 mg/kg/d, sc) in addition to CPF (1 or 5 mg/kg/d, po) for 7 d. The major CPF metabolite, 3,5,6-trichloro-2-pyridinol (TCPy), in blood and urine and the plasma ChE and brain acetylcholinesterase (AChE) activities were measured in rats. There were differences in pharmacokinetics, with higher TCPy peak concentrations and increased blood TCPy AUC in ethanol + nicotine groups compared to CPF only (approximately 1.8- and 3.8-fold at 1 and 5 mg CPF doses, respectively). Brain AChE activities after ethanol + nicotine treatments showed significantly less inhibition following repeated 5 mg CPF/kg dosing compared to CPF only (96 ± 13 and 66 ± 7% of naive at 4 h post last CPF dosing, respectively). Although brain AChE activity was minimal inhibited for the 1-mg CPF/kg/d groups, the ethanol + nicotine pretreatment resulted in a similar trend (i.e., slightly less inhibition). No marked differences were observed in plasma ChE activities due to the alcohol + nicotine treatments. In vitro, CPF metabolism was not markedly affected by repeated ethanol or both ethanol + nicotine exposures. Compared with a previous study of nicotine and CPF exposure, there were no apparent additional exacerbating effects due to ethanol coexposure.

Global parameter estimation methods for stochastic biochemical systems

Background

The importance of stochasticity in cellular processes having low number of molecules has resulted in the development of stochastic models such as chemical master equation. As in other modelling frameworks, the accompanying rate constants are important for the end-applications like analyzing system properties (e.g. robustness) or predicting the effects of genetic perturbations. Prior knowledge of kinetic constants is usually limited and the model identification routine typically includes parameter estimation from experimental data. Although the subject of parameter estimation is well-established for deterministic models, it is not yet routine for the chemical master equation. In addition, recent advances in measurement technology have made the quantification of genetic substrates possible to single molecular levels. Thus, the purpose of this work is to develop practical and effective methods for estimating kinetic model parameters in the chemical master equation and other stochastic models from single cell and cell population experimental data.

Results

Three parameter estimation methods are proposed based on the maximum likelihood and density function distance, including probability and cumulative density functions. Since stochastic models such as chemical master equations are typically solved using a Monte Carlo approach in which only a finite number of Monte Carlo realizations are computationally practical, specific considerations are given to account for the effect of finite sampling in the histogram binning of the state density functions. Applications to three practical case studies showed that while maximum likelihood method can effectively handle low replicate measurements, the density function distance methods, particularly the cumulative density function distance estimation, are more robust in estimating the parameters with consistently higher accuracy, even for systems showing multimodality.

Conclusions

The parameter estimation methodologies described in this work have provided an effective and practical approach in the estimation of kinetic parameters of stochastic systems from either sparse or dense cell population data. Nevertheless, similar to kinetic parameter estimation in other modelling frameworks, not all parameters can be estimated accurately, which is a common problem arising from the lack of complete parameter identifiability from the available data.

Parameter identifiability of power-law biochemical system models

Mathematical modeling has become an integral component in biotechnology, in which these models are frequently used to design and optimize bioprocesses. Canonical models, like power-laws within the Biochemical Systems Theory, offer numerous mathematical and numerical advantages, including built-in flexibility to simulate general nonlinear behavior. The construction of such models relies on the estimation of unknown case-specific model parameters by way of experimental data fitting, also known as inverse modeling. Despite the large number of publications on this topic, this task remains the bottleneck in canonical modeling of biochemical systems. The focus of this paper concerns with the question of identifiability of power-law models from dynamic data, that is, whether the parameter values can be uniquely and accurately identified from time-series data. Existing and newly developed parameter identifiability methods were applied to two power-law models of biochemical systems, and the results pointed to the lack of parametric identifiability as the root cause of the difficulty faced in the inverse modeling. Despite the focus on power-law models, the analyses and conclusions are extendable to other canonical models, and the issue of parameter identifiability is expected to be a common problem in biochemical system modeling.