A coarse-grained resource allocation model of carbon and nitrogen metabolism in unicellular microbes

Coarse-grained resource allocation models (C-GRAMs) are simple mathematical models of cell physiology, where large components of the macromolecular composition are abstracted into single entities. The dynamics and steady-state behaviour of such models provides insights on optimal allocation of cellular resources and have explained experimentally observed cellular growth laws, but current models do not account for the uptake of compound sources of carbon and nitrogen. Here, we formulate a C-GRAM with nitrogen and carbon pathways converging on biomass production, with parametrizations accounting for respirofermentative and purely respiratory growth. The model describes the effects of the uptake of sugars, ammonium and/or compound nutrients such as amino acids on the translational resource allocation towards proteome sectors that maximized the growth rate. It robustly recovers cellular growth laws including the Monod law and the ribosomal growth law. Furthermore, we show how the growth-maximizing balance between carbon uptake, recycling, and excretion depends on the nutrient environment. Lastly, we find a robust linear correlation between the ribosome fraction and the abundance of amino acid equivalents in the optimal cell, which supports the view that simple regulation of translational gene expression can enable cells to achieve an approximately optimal growth state.

Modelling capture efficiency of single-cell RNA-sequencing data improves inference of transcriptome-wide burst kinetics

MOTIVATION

Gene expression is characterized by stochastic bursts of transcription that occur at brief and random periods of promoter activity. The kinetics of gene expression burstiness differs across the genome and is dependent on the promoter sequence, among other factors. Single-cell RNA sequencing (scRNA-seq) has made it possible to quantify the cell-to-cell variability in transcription at a global genome-wide level. However, scRNA-seq data are prone to technical variability, including low and variable capture efficiency of transcripts from individual cells.

RESULTS

Here, we propose a novel mathematical theory for the observed variability in scRNA-seq data. Our method captures burst kinetics and variability in both the cell size and capture efficiency, which allows us to propose several likelihood-based and simulation-based methods for the inference of burst kinetics from scRNA-seq data. Using both synthetic and real data, we show that the simulation-based methods provide an accurate, robust and flexible tool for inferring burst kinetics from scRNA-seq data. In particular, in a supervised manner, a simulation-based inference method based on neural networks proves to be accurate and useful when applied to both allele and nonallele-specific scRNA-seq data.

AVAILABILITY AND IMPLEMENTATION

The code for Neural Network and Approximate Bayesian Computation inference is available at https://github.com/WT215/nnRNA and https://github.com/WT215/Julia_ABC, respectively.

Data-driven spatio-temporal modelling of glioblastoma

Mathematical oncology provides unique and invaluable insights into tumour growth on both the microscopic and macroscopic levels. This review presents state-of-the-art modelling techniques and focuses on their role in understanding glioblastoma, a malignant form of brain cancer. For each approach, we summarize the scope, drawbacks and assets. We highlight the potential clinical applications of each modelling technique and discuss the connections between the mathematical models and the molecular and imaging data used to inform them. By doing so, we aim to prime cancer researchers with current and emerging computational tools for understanding tumour progression. By providing an in-depth picture of the different modelling techniques, we also aim to assist researchers who seek to build and develop their own models and the associated inference frameworks. Our article thus strikes a unique balance. On the one hand, we provide a comprehensive overview of the available modelling techniques and their applications, including key mathematical expressions. On the other hand, the content is accessible to mathematicians and biomedical scientists alike to accommodate the interdisciplinary nature of cancer research.

Selective advantage of epigenetically disrupted cancer cells via phenotypic inertia

The evolution of established cancers is driven by selection of cells with enhanced fitness. Subclonal mutations in numerous epigenetic regulator genes are common across cancer types, yet their functional impact has been unclear. Here, we show that disruption of the epigenetic regulatory network increases the tolerance of cancer cells to unfavorable environments experienced within growing tumors by promoting the emergence of stress-resistant subpopulations. Disruption of epigenetic control does not promote selection of genetically defined subclones or favor a phenotypic switch in response to environmental changes. Instead, it prevents cells from mounting an efficient stress response via modulation of global transcriptional activity. This "transcriptional numbness" lowers the probability of cell death at early stages, increasing the chance of long-term adaptation at the population level. Our findings provide a mechanistic explanation for the widespread selection of subclonal epigenetic-related mutations in cancer and uncover phenotypic inertia as a cellular trait that drives subclone expansion.

Fission yeast obeys a linear size law under nutrient titration

Steady-state cell size and geometry depend on growth conditions. Here, we use an experimental setup based on continuous culture and single-cell imaging to study how cell volume, length, width and surface-to-volume ratio vary across a range of growth conditions including nitrogen and carbon titration, the choice of nitrogen source, and translation inhibition. Overall, we find cell geometry is not fully determined by growth rate and depends on the specific mode of growth rate modulation. However, under nitrogen and carbon titrations, we observe that the cell volume and the growth rate follow the same linear scaling.

Efficient Bayesian inference for stochastic agent-based models

The modelling of many real-world problems relies on computationally heavy simulations of randomly interacting individuals or agents. However, the values of the parameters that underlie the interactions between agents are typically poorly known, and hence they need to be inferred from macroscopic observations of the system. Since statistical inference rests on repeated simulations to sample the parameter space, the high computational expense of these simulations can become a stumbling block. In this paper, we compare two ways to mitigate this issue in a Bayesian setting through the use of machine learning methods: One approach is to construct lightweight surrogate models to substitute the simulations used in inference. Alternatively, one might altogether circumvent the need for Bayesian sampling schemes and directly estimate the posterior distribution. We focus on stochastic simulations that track autonomous agents and present two case studies: tumour growths and the spread of infectious diseases. We demonstrate that good accuracy in inference can be achieved with a relatively small number of simulations, making our machine learning approaches orders of magnitude faster than classical simulation-based methods that rely on sampling the parameter space. However, we find that while some methods generally produce more robust results than others, no algorithm offers a one-size-fits-all solution when attempting to infer model parameters from observations. Instead, one must choose the inference technique with the specific real-world application in mind. The stochastic nature of the considered real-world phenomena poses an additional challenge that can become insurmountable for some approaches. Overall, we find machine learning approaches that create direct inference machines to be promising for real-world applications. We present our findings as general guidelines for modelling practitioners.

Discrimination of Escherichia coli, Shigella flexneri, and Shigella sonnei using lipid profiling by MALDI-TOF mass spectrometry paired with machine learning

Matrix‐assisted laser desorption/ionization‐time of flight mass spectrometry (MALDI‐TOF MS) has become a staple in clinical microbiology laboratories. Protein‐profiling of bacteria using this technique has accelerated the identification of pathogens in diagnostic workflows. Recently, lipid profiling has emerged as a way to complement bacterial identification where protein‐based methods fail to provide accurate results. This study aimed to address the challenge of rapid discrimination between Escherichia coli and Shigella spp. using MALDI‐TOF MS in the negative ion mode for lipid profiling coupled with machine learning. Both E. coli and Shigella species are closely related; they share high sequence homology, reported for 16S rRNA gene sequence similarities between E. coli and Shigella spp. exceeding 99%, and a similar protein expression pattern but are epidemiologically distinct. A bacterial collection of 45 E. coli, 48 Shigella flexneri, and 62 Shigella sonnei clinical isolates were submitted to lipid profiling in negative ion mode using the MALDI Biotyper Sirius® system after treatment with mild‐acid hydrolysis (acetic acid 1% v/v for 15 min at 98°C). Spectra were then analyzed using our in‐house machine learning algorithm and top‐ranked features used for the discrimination of the bacterial species. Here, as a proof‐of‐concept, we showed that lipid profiling might have the potential to differentiate E. coli from Shigella species using the analysis of the top five ranked features obtained by MALDI‐TOF MS in the negative ion mode of the MALDI Biotyper Sirius® system. Based on this new approach, MALDI‐TOF MS analysis of lipids might help pave the way toward these goals.

Benchmarking imputation methods for network inference using a novel method of synthetic scRNA-seq data generation

Background

Single cell RNA-sequencing (scRNA-seq) has very rapidly become the new workhorse of modern biology providing an unprecedented global view on cellular diversity and heterogeneity. In particular, the structure of gene-gene expression correlation contains information on the underlying gene regulatory networks. However, interpretation of scRNA-seq data is challenging due to specific experimental error and biases that are unique to this kind of data including drop-out (or technical zeros).

Methods

To deal with this problem several methods for imputation of zeros for scRNA-seq have been developed. However, it is not clear how these processing steps affect inference of genetic networks from single cell data. Here, we introduce Biomodelling.jl, a tool for generation of synthetic scRNA-seq data using multiscale modelling of stochastic gene regulatory networks in growing and dividing cells.

Results

Our tool produces realistic transcription data with a known ground truth network topology that can be used to benchmark different approaches for gene regulatory network inference. Using this tool we investigate the impact of different imputation methods on the performance of several network inference algorithms.

Conclusions

Biomodelling.jl provides a versatile and useful tool for future development and benchmarking of network inference approaches using scRNA-seq data.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-022-04778-9

Emergent expression of fitness-conferring genes by phenotypic selection

Genotypic and phenotypic adaptation is the consequence of ongoing natural selection in populations and is key to predicting and preventing drug resistance. Whereas classic antibiotic persistence is all-or-nothing, here we demonstrate that an antibiotic resistance gene displays linear dose-responsive selection for increased expression in proportion to rising antibiotic concentration in growing Escherichia coli populations. Furthermore, we report the potentially wide-spread nature of this form of emergent gene expression (EGE) by instantaneous phenotypic selection process under bactericidal and bacteriostatic antibiotic treatment, as well as an amino acid synthesis pathway enzyme under a range of auxotrophic conditions. We propose an analogy to Ohm's law in electricity (V = IR), where selection pressure acts similarly to voltage (V), gene expression to current (I), and resistance (R) to cellular machinery constraints and costs. Lastly, mathematical modeling using agent-based models of stochastic gene expression in growing populations and Bayesian model selection reveal that the EGE mechanism requires variability in gene expression within an isogenic population, and a cellular "memory" from positive feedbacks between growth and expression of any fitness-conferring gene. Finally, we discuss the connection of the observed phenomenon to a previously described general fluctuation-response relationship in biology.

Growth-rate-dependent and nutrient-specific gene expression resource allocation in fission yeast

Cellular resources are limited and their relative allocation to gene expression programmes determines physiological states and global properties such as the growth rate. Here, we determined the importance of the growth rate in explaining relative changes in protein and mRNA levels in the simple eukaryote Schizosaccharomyces pombe grown on non-limiting nitrogen sources. Although expression of half of fission yeast genes was significantly correlated with the growth rate, this came alongside wide-spread nutrient-specific regulation. Proteome and transcriptome often showed coordinated regulation but with notable exceptions, such as metabolic enzymes. Genes positively correlated with growth rate participated in every level of protein production apart from RNA polymerase II-dependent transcription. Negatively correlated genes belonged mainly to the environmental stress response programme. Critically, metabolic enzymes, which represent ∼55-70% of the proteome by mass, showed mostly condition-specific regulation. In summary, we provide a rich account of resource allocation to gene expression in a simple eukaryote, advancing our basic understanding of the interplay between growth-rate-dependent and nutrient-specific gene expression.

Coordination of gene expression noise with cell size: analytical results for agent-based models of growing cell populations

The chemical master equation and the Gillespie algorithm are widely used to model the reaction kinetics inside living cells. It is thereby assumed that cell growth and division can be modelled through effective dilution reactions and extrinsic noise sources. We here re-examine these paradigms through developing an analytical agent-based framework of growing and dividing cells accompanied by an exact simulation algorithm, which allows us to quantify the dynamics of virtually any intracellular reaction network affected by stochastic cell size control and division noise. We find that the solution of the chemical master equation-including static extrinsic noise-exactly agrees with the agent-based formulation when the network under study exhibits stochastic concentration homeostasis, a novel condition that generalizes concentration homeostasis in deterministic systems to higher order moments and distributions. We illustrate stochastic concentration homeostasis for a range of common gene expression networks. When this condition is not met, we demonstrate by extending the linear noise approximation to agent-based models that the dependence of gene expression noise on cell size can qualitatively deviate from the chemical master equation. Surprisingly, the total noise of the agent-based approach can still be well approximated by extrinsic noise models.

Phenotypic Robustness of Epidermal Stem Cell Number in C. elegans Is Modulated by the Activity of the Conserved N-acetyltransferase nath-10/NAT10

Individual cells and organisms experience perturbations from internal and external sources, yet manage to buffer these to produce consistent phenotypes, a property known as robustness. While phenotypic robustness has often been examined in unicellular organisms, it has not been sufficiently studied in multicellular animals. Here, we investigate phenotypic robustness in Caenorhabditis elegans seam cells. Seam cells are stem cell-like epithelial cells along the lateral edges of the animal, which go through asymmetric and symmetric divisions contributing cells to the hypodermis and neurons, while replenishing the stem cell reservoir. The terminal number of seam cells is almost invariant in the wild-type population, allowing the investigation of how developmental precision is achieved. We report here that a loss-of-function mutation in the highly conserved N-acetyltransferase nath-10/NAT10 increases seam cell number variance in the isogenic population. RNA-seq analysis revealed increased levels of mRNA transcript variability in nath-10 mutant populations, which may have an impact on the phenotypic variability observed. Furthermore, we found disruption of Wnt signaling upon perturbing nath-10 function, as evidenced by changes in POP-1/TCF nuclear distribution and ectopic activation of its GATA transcription factor target egl-18. These results highlight that NATH-10/NAT-10 can influence phenotypic variability partly through modulation of the Wnt signaling pathway.

Integrating multi-OMICS data through sparse canonical correlation analysis for the prediction of complex traits: a comparison study

Motivation

Recent developments in technology have enabled researchers to collect multiple OMICS datasets for the same individuals. The conventional approach for understanding the relationships between the collected datasets and the complex trait of interest would be through the analysis of each OMIC dataset separately from the rest, or to test for associations between the OMICS datasets. In this work we show that integrating multiple OMICS datasets together, instead of analysing them separately, improves our understanding of their in-between relationships as well as the predictive accuracy for the tested trait. Several approaches have been proposed for the integration of heterogeneous and high-dimensional (p≫n) data, such as OMICS. The sparse variant of canonical correlation analysis (CCA) approach is a promising one that seeks to penalize the canonical variables for producing sparse latent variables while achieving maximal correlation between the datasets. Over the last years, a number of approaches for implementing sparse CCA (sCCA) have been proposed, where they differ on their objective functions, iterative algorithm for obtaining the sparse latent variables and make different assumptions about the original datasets.

Results

Through a comparative study we have explored the performance of the conventional CCA proposed by Parkhomenko et al., penalized matrix decomposition CCA proposed by Witten and Tibshirani and its extension proposed by Suo et al. The aforementioned methods were modified to allow for different penalty functions. Although sCCA is an unsupervised learning approach for understanding of the in-between relationships, we have twisted the problem as a supervised learning one and investigated how the computed latent variables can be used for predicting complex traits. The approaches were extended to allow for multiple (more than two) datasets where the trait was included as one of the input datasets. Both ways have shown improvement over conventional predictive models that include one or multiple datasets.

Availability and implementation

https://github.com/theorod93/sCCA.

Supplementary information

Supplementary data are available at Bioinformatics online.

The impact of random microenvironmental fluctuations on tumor control probability

The tumor control probability (TCP) is a metric used to calculate the probability of controlling or eradicating tumors through radiotherapy. Cancer cells vary in their response to radiation, and although many factors are involved, the tumor microenvironment is a crucial one that determines radiation efficacy. The tumor microenvironment plays a significant role in cancer initiation and propagation, as well as in treatment outcome. We have developed stochastic formulations to study the impact of arbitrary microenvironmental fluctuations on TCP and extinction probability (EP), which is defined as the probability of cancer cells removal in the absence of treatment. Since the derivation of analytical solutions are not possible for complicated cases, we employ a modified Gillespie algorithm to analyze TCP and EP, considering the random variations in cellular proliferation and death rates. Our results show that increasing the standard deviation in kinetic rates initially enhances the probability of tumor eradication. However, if the EP does not reach a probability of 1, the increase in the standard deviation subsequently has a negative impact on probability of cancer cells removal, decreasing the EP over time. The greatest effect on EP has been observed when both birth and death rates are being randomly modified and are anticorrelated. In addition, similar results are observed for TCP, where radiotherapy is included, indicating that increasing the standard deviation in kinetic rates at first enhances the probability of tumor eradication. But, it has a negative impact on treatment effectiveness if the TCP does not reach a probability of 1.

bayNorm: Bayesian gene expression recovery, imputation and normalization for single-cell RNA-sequencing data

Motivation

Normalization of single-cell RNA-sequencing (scRNA-seq) data is a prerequisite to their interpretation. The marked technical variability, high amounts of missing observations and batch effect typical of scRNA-seq datasets make this task particularly challenging. There is a need for an efficient and unified approach for normalization, imputation and batch effect correction.

Results

Here, we introduce bayNorm, a novel Bayesian approach for scaling and inference of scRNA-seq counts. The method’s likelihood function follows a binomial model of mRNA capture, while priors are estimated from expression values across cells using an empirical Bayes approach. We first validate our assumptions by showing this model can reproduce different statistics observed in real scRNA-seq data. We demonstrate using publicly available scRNA-seq datasets and simulated expression data that bayNorm allows robust imputation of missing values generating realistic transcript distributions that match single molecule fluorescence in situ hybridization measurements. Moreover, by using priors informed by dataset structures, bayNorm improves accuracy and sensitivity of differential expression analysis and reduces batch effect compared with other existing methods. Altogether, bayNorm provides an efficient, integrated solution for global scaling normalization, imputation and true count recovery of gene expression measurements from scRNA-seq data.

Availability and implementation

The R package ‘bayNorm’ is publishd on bioconductor at https://bioconductor.org/packages/release/bioc/html/bayNorm.html. The code for analyzing data in this article is available at https://github.com/WT215/bayNorm_papercode.

Supplementary information

Supplementary data are available at Bioinformatics online.

Epimutations driven by small RNAs arise frequently but most have limited duration in Caenorhabditis elegans

Epigenetic regulation involves changes in gene expression independent of DNA sequence variation that are inherited through cell division. In addition to a fundamental role in cell differentiation, some epigenetic changes can also be transmitted transgenerationally through meiosis. Epigenetic alterations (epimutations) could thus contribute to heritable variation within populations and be subject to evolutionary processes such as natural selection and drift. However, the rate at which epimutations arise and their typical persistence are unknown, making it difficult to evaluate their potential for evolutionary adaptation. Here, we perform a genome-wide study of epimutations in a metazoan organism. We use experimental evolution to characterize the rate, spectrum and stability of epimutations driven by small silencing RNAs in the model nematode Caenorhabditis elegans. We show that epimutations arise spontaneously at a rate approximately 25 times greater than DNA sequence changes and typically have short half-lives of two to three generations. Nevertheless, some epimutations last at least ten generations. Epimutations mediated by small RNAs may thus contribute to evolutionary processes over a short timescale but are unlikely to bring about long-term divergence in the absence of selection.

Discrimination of bovine milk from non-dairy milk by lipids fingerprinting using routine matrix-assisted laser desorption ionization mass spectrometry

An important sustainable development goal for any country is to ensure food security by producing a sufficient and safe food supply. This is the case for bovine milk where addition of non-dairy milks such as vegetables (e.g., soya or coconut) has become a common source of adulteration and fraud. Conventionally, gas chromatography techniques are used to detect key lipids (e.g., triacylglycerols) has an effective read-out of assessing milks origins and to detect foreign milks in bovine milks. However, such approach requires several sample preparation steps and a dedicated laboratory environment, precluding a high throughput process. To cope with this need, here, we aimed to develop a novel and simple method without organic solvent extractions for the detection of bovine and non-dairy milks based on lipids fingerprint by routine MALDI-TOF mass spectrometry (MS). The optimized method relies on the simple dilution of milks in water followed by MALDI-TOF MS analyses in the positive linear ion mode and using a matrix consisting of a 9:1 mixture of 2,5-dihydroxybenzoic acid and 2-hydroxy-5-methoxybenzoic acid (super-DHB) solubilized at 10 mg/mL in 70% ethanol. This sensitive, inexpensive, and rapid method has potential for use in food authenticity applications.

LiBiNorm: an htseq-count analogue with improved normalisation of Smart-seq2 data and library preparation diagnostics

Protocols for preparing RNA sequencing (RNA-seq) libraries, most prominently “Smart-seq” variations, introduce global biases that can have a significant impact on the quantification of gene expression levels. This global bias can lead to drastic over- or under-representation of RNA in non-linear length-dependent fashion due to enzymatic reactions during cDNA production. It is currently not corrected by any RNA-seq software, which mostly focus on local bias in coverage along RNAs. This paper describes LiBiNorm, a simple command line program that mimics the popular htseq-count software and allows diagnostics, quantification, and global bias removal. LiBiNorm outputs gene expression data that has been normalized to correct for global bias introduced by the Smart-seq2 protocol. In addition, it produces data and several plots that allow insights into the experimental history underlying library preparation. The LiBiNorm package includes an R script that allows visualization of the main results. LiBiNorm is the first software application to correct for the global bias that is introduced by the Smart-seq2 protocol. It is freely downloadable at http://www2.warwick.ac.uk/fac/sci/lifesci/research/libinorm.

MALDI-TOF mass spectrometry on intact bacteria combined with a refined analysis framework allows accurate classification of MSSA and MRSA

Fast and reliable detection coupled with accurate data-processing and analysis of antibiotic-resistant bacteria is essential in clinical settings. In this study, we use MALDI-TOF on intact cells combined with a refined analysis framework to demonstrate discrimination between methicillin-susceptible (MSSA) and methicillin-resistant (MRSA) Staphylococcus aureus. By combining supervised and unsupervised machine learning methods, we firstly show that the mass spectroscopy data contains strong signal for the clustering of MSSA and MRSA. Then we concentrate on applying supervised learning to extract and verify the important features. A new workflow is proposed that allows for extracting a fixed set of reference peaks so that any new data can be aligned to it and hence consistent feature matrices can be obtained. Also note that by doing so we are able to examine the robustness of the important features that have been found. We also show that appropriate size of the benchmark data, appropriate alignment of the testing data and use of an optimal set of features via feature selection results in prediction accuracy over 90%. In summary, as proof-of-principle, our integrated experimental and bioinformatics study suggests a novel intact cell MALDI-TOF to be of great promise for fast and reliable detection of MRSA strains.

Analysis of ParAB dynamics in mycobacteria shows active movement of ParB and differential inheritance of ParA

Correct chromosomal segregation, coordinated with cell division, is crucial for bacterial survival, but despite extensive studies, the mechanisms underlying this remain incompletely understood in mycobacteria. We report a detailed investigation of the dynamic interactions between ParA and ParB partitioning proteins in Mycobacterium smegmatis using microfluidics and time-lapse fluorescence microscopy to observe both proteins simultaneously. During growth and division, ParB presents as a focused fluorescent spot that subsequently splits in two. One focus moves towards a higher concentration of ParA at the new pole, while the other moves towards the old pole. We show ParB movement is in part an active process that does not rely on passive movement associated with cell growth. In some cells, another round of ParB segregation starts before cell division is complete, consistent with initiation of a second round of chromosome replication. ParA fluorescence distribution correlates with cell size, and in sister cells, the larger cell inherits a local peak of concentrated ParA, while the smaller sister inherits more homogeneously distributed protein. Cells which inherit more ParA grow faster than their sister cell, raising the question of whether inheritance of a local concentration of ParA provides a growth advantage. Alterations in levels of ParA and ParB were also found to disturb cell growth.

Division rate, cell size and proteome allocation: impact on gene expression noise and implications for the dynamics of genetic circuits

The cell division rate, size and gene expression programmes change in response to external conditions. These global changes impact on average concentrations of biomolecule and their variability or noise. Gene expression is inherently stochastic, and noise levels of individual proteins depend on synthesis and degradation rates as well as on cell-cycle dynamics. We have modelled stochastic gene expression inside growing and dividing cells to study the effect of division rates on noise in mRNA and protein expression. We use assumptions and parameters relevant to Escherichia coli, for which abundant quantitative data are available. We find that coupling of transcription, but not translation rates to the rate of cell division can result in protein concentration and noise homeostasis across conditions. Interestingly, we find that the increased cell size at fast division rates, observed in E. coli and other unicellular organisms, buffers noise levels even for proteins with decreased expression at faster growth. We then investigate the functional importance of these regulations using gene regulatory networks that exhibit bi-stability and oscillations. We find that network topology affects robustness to changes in division rate in complex and unexpected ways. In particular, a simple model of persistence, based on global physiological feedback, predicts increased proportion of persister cells at slow division rates. Altogether, our study reveals how cell size regulation in response to cell division rate could help controlling gene expression noise. It also highlights that understanding circuits' robustness across growth conditions is key for the effective design of synthetic biological systems.

Division rate, cell size and proteome allocation: impact on gene expression noise and implications for the dynamics of genetic circuits

The cell division rate, size and gene expression programmes change in response to external conditions. These global changes impact on average concentrations of biomolecule and their variability or noise. Gene expression is inherently stochastic, and noise levels of individual proteins depend on synthesis and degradation rates as well as on cell-cycle dynamics. We have modelled stochastic gene expression inside growing and dividing cells to study the effect of division rates on noise in mRNA and protein expression. We use assumptions and parameters relevant to Escherichia coli, for which abundant quantitative data are available. We find that coupling of transcription, but not translation rates to the rate of cell division can result in protein concentration and noise homeostasis across conditions. Interestingly, we find that the increased cell size at fast division rates, observed in E. coli and other unicellular organisms, buffers noise levels even for proteins with decreased expression at faster growth. We then investigate the functional importance of these regulations using gene regulatory networks that exhibit bi-stability and oscillations. We find that network topology affects robustness to changes in division rate in complex and unexpected ways. In particular, a simple model of persistence, based on global physiological feedback, predicts increased proportion of persister cells at slow division rates. Altogether, our study reveals how cell size regulation in response to cell division rate could help controlling gene expression noise. It also highlights that understanding circuits' robustness across growth conditions is key for the effective design of synthetic biological systems.

Mycobacteria Modify Their Cell Size Control under Sub-Optimal Carbon Sources

The decision to divide is the most important one that any cell must make. Recent single cell studies suggest that most bacteria follow an “adder” model of cell size control, incorporating a fixed amount of cell wall material before dividing. Mycobacteria, including the causative agent of tuberculosis Mycobacterium tuberculosis, are known to divide asymmetrically resulting in heterogeneity in growth rate, doubling time, and other growth characteristics in daughter cells. The interplay between asymmetric cell division and adder size control has not been extensively investigated. Moreover, the impact of changes in the environment on growth rate and cell size control have not been addressed for mycobacteria. Here, we utilize time-lapse microscopy coupled with microfluidics to track live Mycobacterium smegmatis cells as they grow and divide over multiple generations, under a variety of growth conditions. We demonstrate that, under optimal conditions, M. smegmatis cells robustly follow the adder principle, with constant added length per generation independent of birth size, growth rate, and inherited pole age. However, the nature of the carbon source induces deviations from the adder model in a manner that is dependent on pole age. Understanding how mycobacteria maintain cell size homoeostasis may provide crucial targets for the development of drugs for the treatment of tuberculosis, which remains a leading cause of global mortality.

Modeling Enzyme Processivity Reveals that RNA-Seq Libraries Are Biased in Characteristic and Correctable Ways

Experimental procedures for preparing RNA-seq and single-cell (sc) RNA-seq libraries are based on assumptions regarding their underlying enzymatic reactions. Here, we show that the fairness of these assumptions varies within libraries: coverage by sequencing reads along and between transcripts exhibits characteristic, protocol-dependent biases. To understand the mechanistic basis of this bias, we present an integrated modeling framework that infers the relationship between enzyme reactions during library preparation and the characteristic coverage patterns observed for different protocols. Analysis of new and existing (sc)RNA-seq data from six different library preparation protocols reveals that polymerase processivity is the mechanistic origin of coverage biases. We apply our framework to demonstrate that lowering incubation temperature increases processivity, yield, and (sc)RNA-seq sensitivity in all protocols. We also provide correction factors based on our model for increasing accuracy of transcript quantification in existing samples prepared at standard temperatures. In total, our findings improve our ability to accurately reflect in vivo transcript abundances in (sc)RNA-seq libraries.

•

Characterization of global RNA-seq biases specific to library preparation protocols

•

Mathematical framework to reverse engineer enzyme reactions that cause bias

•

Insights from reverse engineering allow optimization of RNA-seq protocols

•

Lowered incubation temperatures during library preparation improve sensitivity

Biosensor architectures for high-fidelity reporting of cellular signaling

Understanding mechanisms of information processing in cellular signaling networks requires quantitative measurements of protein activities in living cells. Biosensors are molecular probes that have been developed to directly track the activity of specific signaling proteins and their use is revolutionizing our understanding of signal transduction. The use of biosensors relies on the assumption that their activity is linearly proportional to the activity of the signaling protein they have been engineered to track. We use mechanistic mathematical models of common biosensor architectures (single-chain FRET-based biosensors), which include both intramolecular and intermolecular reactions, to study the validity of the linearity assumption. As a result of the classic mechanism of zero-order ultrasensitivity, we find that biosensor activity can be highly nonlinear so that small changes in signaling protein activity can give rise to large changes in biosensor activity and vice versa. This nonlinearity is abolished in architectures that favor the formation of biosensor oligomers, but oligomeric biosensors produce complicated FRET states. Based on this finding, we show that high-fidelity reporting is possible when a single-chain intermolecular biosensor is used that cannot undergo intramolecular reactions and is restricted to forming dimers. We provide phase diagrams that compare various trade-offs, including observer effects, which further highlight the utility of biosensor architectures that favor intermolecular over intramolecular binding. We discuss challenges in calibrating and constructing biosensors and highlight the utility of mathematical models in designing novel probes for cellular signaling.

Connecting growth with gene expression: of noise and numbers

Growth is a dynamic process whereby cells accumulate mass. Growth rates of single cells are connected to RNA and protein synthesis rates, and therefore with biomolecule numbers. Noise in gene expression depends on these numbers, and is thus linked with cellular growth. Whether these global attributes of the cell participate in gene regulation is still largely unexplored. New experimental and modelling studies suggest that systemic variations in biomolecule numbers can coordinate cellular processes, including growth itself, through global regulatory feedback that acts in addition to genetic regulatory networks. Here, we review these findings and speculate on possible implications of this less appreciated layer of gene regulation for cellular physiology and adaptation to changing environments.

General transient solution of the one-step master equation in one dimension

Exact analytical solutions of the master equation are limited to special cases and exact numerical methods are inefficient. Even the generic one-dimensional, one-step master equation has evaded exact solution, aside from the steady-state case. This type of master equation describes the dynamics of a continuous-time Markov process whose range consists of positive integers and whose transitions are allowed only between adjacent sites. The solution of any master equation can be written as the exponential of a (typically huge) matrix, which requires the calculation of the eigenvalues and eigenvectors of the matrix. Here we propose a linear algebraic method for simplifying this exponential for the general one-dimensional, one-step process. In particular, we prove that the calculation of the eigenvectors is actually not necessary for the computation of exponential, thereby we dramatically cut the time of this calculation. We apply our new methodology to examples from birth-death processes and biochemical networks. We show that the computational time is significantly reduced compared to existing methods.

Mycobacterial Growth

In this work, we review progress made in understanding the molecular underpinnings of growth and division in mycobacteria, concentrating on work published since the last comprehensive review ( Hett and Rubin 2008). We have focused on exciting work making use of new time-lapse imaging technologies coupled with reporter-gene fusions and antimicrobial treatment to generate insights into how mycobacteria grow and divide in a heterogeneous manner. We try to reconcile the different observations reported, providing a model of how they might fit together. We also review the topic of mycobacterial spores, which has generated considerable discussion during the last few years. Resuscitation promoting factors, and regulation of growth and division, have also been actively researched, and we summarize progress in these areas.

Stochastic cellular fate decision making by multiple infecting lambda phage

Bacteriophage lambda is a classic system for the study of cellular decision making. Both experiments and mathematical models have demonstrated the importance of viral concentration in the lysis-lysogeny decision outcome in lambda phage. However, a recent experimental study using single cell and single phage resolution reported that cells with the same viral concentrations but different numbers of infecting phage (multiplicity of infection) can have markedly different rates of lysogeny. Thus the decision depends on not only viral concentration, but also directly on the number of infecting phage. Here, we attempt to provide a mechanistic explanation of these results using a simple stochastic model of the lambda phage genetic network. Several potential factors including intrinsic gene expression noise, spatial dynamics and cell-cycle effects are investigated. We find that interplay between the level of intrinsic noise and viral protein decision threshold is a major factor that produces dependence on multiplicity of infection. However, simulations suggest spatial segregation of phage particles does not play a significant role. Cellular image processing is used to re-analyse the original time-lapse movies from the recent study and it is found that higher numbers of infecting phage reduce the cell elongation rate. This could also contribute to the observed phenomena as cellular growth rate can affect transcription rates. Our model further predicts that rate of lysogeny is dependent on bacterial growth rate, which can be experimentally tested. Our study provides new insight on the mechanisms of individual phage decision making. More generally, our results are relevant for the understanding of gene-dosage compensation in cellular systems.

Uncaging calcium in neurons

Changes in intracellular free calcium concentration (Δ[Ca(2+)]i) driving physiological events such as neurotransmitter release or Ca(2+)-dependent currents can be monitored using Ca(2+)-sensitive fluorescent dyes. Although these dyes can correlate Δ[Ca(2+)]i with a physiological event, they cannot directly test for causality between changes in [Ca(2+)]i and that event. Photolabile Ca(2+) chelators are Ca(2+)-binding molecules that can alter and, to a certain extent, control [Ca(2+)]i in an inducible manner and with temporal and spatial resolution that surpasses microinjection or ionophore application. Here we discuss the properties of caged Ca(2+) compounds as well as some practical considerations for their use in neuronal cells, where they have proven particularly effective.

Cell division site placement and asymmetric growth in mycobacteria

Mycobacteria are members of the actinomycetes that grow by tip extension and lack apparent homologues of the known cell division regulators found in other rod-shaped bacteria. Previous work using static microscopy on dividing mycobacteria led to the hypothesis that these cells can grow and divide asymmetrically, and at a wide range of sizes, in contrast to the cell growth and division patterns observed in the model rod-shaped organisms. In this study, we test this hypothesis using live-cell time-lapse imaging of dividing Mycobacterium smegmatis labelled with fluorescent PBP1a, to probe peptidoglycan synthesis and label the cell septum. We demonstrate that the new septum is placed accurately at mid-cell, and that the asymmetric division observed is a result of differential growth from the cell tips, with a more than 2-fold difference in growth rate between fast and slow growing poles. We also show that the division site is not selected at a characteristic cell length, suggesting this is not an important cue during the mycobacterial cell cycle.

Bell-shaped and ultrasensitive dose-response in phosphorylation-dephosphorylation cycles: the role of kinase-phosphatase complex formation

Background

Phosphorylation-dephosphorylation cycles (PDCs) mediated by kinases and phosphatases are common in cellular signalling. Kinetic modelling of PDCs has shown that these systems can exhibit a variety of input-output (dose-response) behaviors including graded response, ultrasensitivity and bistability. In addition to proteins, there are a class of lipids known as phosphoinositides (PIs) that can be phosphorylated. Experimental studies have revealed the formation of an antagonistic kinase-phosphatase complex in regulation of phosphorylation of PIs. However, the functional significance of this type of complex formation is not clear.

Results

We first revisit the basic PDC and show that partial asymptotic phosphorylation of substrate limits ultrasensitivity. Also, substrate levels are changed one can obtain non-monotonic bell-shaped dose-response curves over a narrow range of parameters. Then we extend the PDC to include kinase-phosphatase complex formation. We report the possibility of robust bell-shaped dose-response for a specific class of the model with complex formation. Also, we show that complex formation can produce ultrasensitivity outside the Goldbeter-Koshland zero-order ultrasensitivity regime through a mechanism similar to competitive inhibition between an enzyme and its inhibitor.

Conclusions

We conclude that the novel PDC module studied here exhibits new dose-response behaviour. In particular, we show that the bell-shaped response could result in transient phosphorylation of substrate. We discuss the relevance of this result in the context of experimental observations on PI regulation in endosomal trafficking.

Modelling the influence of foot-and-mouth disease vaccine antigen stability and dose on the bovine immune response

Foot and mouth disease virus causes a livestock disease of significant global socio-economic importance. Advances in its control and eradication depend critically on improvements in vaccine efficacy, which can be best achieved by better understanding the complex within-host immunodynamic response to inoculation. We present a detailed and empirically parametrised dynamical mathematical model of the hypothesised immune response in cattle, and explore its behaviour with reference to a variety of experimental observations relating to foot and mouth immunology. The model system is able to qualitatively account for the observed responses during in-vivo experiments, and we use it to gain insight into the incompletely understood effect of single and repeat inoculations of differing dosage using vaccine formulations of different structural stability.

Ultrasensitivity in multisite phosphorylation of membrane-anchored proteins

Cellular signaling is initially confined to the plasma membrane, where the cytoplasmic tails of surface receptors and other membrane-anchored proteins are phosphorylated in response to ligand binding. These proteins often contain multiple phosphorylation sites that are regulated by membrane-confined enzymes. Phosphorylation of these proteins is thought to be tightly regulated, because they initiate and regulate signaling cascades leading to cellular activation, yet how their phosphorylation is regulated is poorly understood. Ultrasensitive or switchlike responses in their phosphorylation state are not expected because the modifying enzymes are in excess. Here, we describe a novel mechanism of ultrasensitivity exhibited by multisite membrane-anchored proteins, but not cytosolic proteins, even when enzymes are in excess. The mechanism underlying this concentration-independent ultrasensitivity is the local saturation of a single enzyme by multiple sites on the substrate. Local saturation is a passive process arising from slow membrane diffusion, steric hindrances, and multiple sites, and therefore may be widely applicable. Critical to this ultrasensitivity is the brief enzymatic inactivation that follows substrate modification. Computations are presented using ordinary differential equations and stochastic spatial simulations. We propose a new role, to our knowledge, for multisite membrane-anchored proteins, discuss experiments that can be used to probe the model, and relate our findings to previous theoretical work.

Scalable rule-based modelling of allosteric proteins and biochemical networks

Much of the complexity of biochemical networks comes from the information-processing abilities of allosteric proteins, be they receptors, ion-channels, signalling molecules or transcription factors. An allosteric protein can be uniquely regulated by each combination of input molecules that it binds. This “regulatory complexity” causes a combinatorial increase in the number of parameters required to fit experimental data as the number of protein interactions increases. It therefore challenges the creation, updating, and re-use of biochemical models. Here, we propose a rule-based modelling framework that exploits the intrinsic modularity of protein structure to address regulatory complexity. Rather than treating proteins as “black boxes”, we model their hierarchical structure and, as conformational changes, internal dynamics. By modelling the regulation of allosteric proteins through these conformational changes, we often decrease the number of parameters required to fit data, and so reduce over-fitting and improve the predictive power of a model. Our method is thermodynamically grounded, imposes detailed balance, and also includes molecular cross-talk and the background activity of enzymes. We use our Allosteric Network Compiler to examine how allostery can facilitate macromolecular assembly and how competitive ligands can change the observed cooperativity of an allosteric protein. We also develop a parsimonious model of G protein-coupled receptors that explains functional selectivity and can predict the rank order of potency of agonists acting through a receptor. Our methodology should provide a basis for scalable, modular and executable modelling of biochemical networks in systems and synthetic biology.

The complexity of biochemical networks challenges our ability to create quantitative and predictive models of cellular responses to extracellular changes. In these networks, the regulation of allosteric receptors and proteins by multiple drugs or endogenous ligands introduces “regulatory complexity” because a large number of parameters is required to describe such interactions. Protein interactions also give rise to “combinatorial complexity” by generating large numbers of protein complexes and covalent modification states. To address these twin problems, we propose a modelling framework that combines a modular description of protein structure and function with a rule-based description of protein interactions. We define the input-output function of an allosteric protein through its thermodynamic properties and structural components. We show that our “biomolecule-centric” methodology, in contrast to ad hoc approaches that emphasize the regulatory logic of interactions, can reduce the number of parameters required to model experimental observations. We also demonstrate how the application of our framework gives insights into the assembly of macromolecular complexes and increases the predictive power of a standard model of G protein-coupled receptors. These benefits are possible in many systems, given the ubiquity of allostery in biochemical networks. Our research delineates a fundamental relationship between allostery, modularity, and complexity in biochemical networks.

The stochastic nature of biochemical networks

Cell behaviour and the cellular environment are stochastic. Phenotypes vary across isogenic populations and in individual cells over time. Here we will argue that to understand the abilities of cells we need to understand their stochastic nature. New experimental techniques allow gene expression to be followed in single cells over time and reveal stochastic bursts of both mRNA and protein synthesis in many different types of organisms. Stochasticity has been shown to be exploited by bacteria and viruses to decide between different behaviours. In fluctuating environments, cells that respond stochastically can out-compete those that sense environmental changes, and stochasticity may even have contributed to chromosomal gene order. We will focus on advances in modelling stochasticity, in understanding its effects on evolution and cellular design, and on means by which it may be exploited in biotechnology and medicine.

Colored extrinsic fluctuations and stochastic gene expression

Stochasticity is both exploited and controlled by cells. Although the intrinsic stochasticity inherent in biochemistry is relatively well understood, cellular variation, or ‘noise', is predominantly generated by interactions of the system of interest with other stochastic systems in the cell or its environment. Such extrinsic fluctuations are nonspecific, affecting many system components, and have a substantial lifetime, comparable to the cell cycle (they are ‘colored'). Here, we extend the standard stochastic simulation algorithm to include extrinsic fluctuations. We show that these fluctuations affect mean protein numbers and intrinsic noise, can speed up typical network response times, and can explain trends in high-throughput measurements of variation. If extrinsic fluctuations in two components of the network are correlated, they may combine constructively (amplifying each other) or destructively (attenuating each other). Consequently, we predict that incoherent feedforward loops attenuate stochasticity, while coherent feedforwards amplify it. Our results demonstrate that both the timescales of extrinsic fluctuations and their nonspecificity substantially affect the function and performance of biochemical networks.

Ca2+ from one or two channels controls fusion of a single vesicle at the frog neuromuscular junction

Neurotransmitter release is triggered by the cooperative action of approximately five Ca2+ ions entering the presynaptic terminal through Ca2+ channels. Depending on the organization of the active zone (AZ), influx through one or many channels may be needed to cause fusion of a vesicle. Using a combination of experiments and modeling, we examined the number of channels that contribute Ca2+ for fusion of a single vesicle in a frog neuromuscular AZ. We compared Ca2+ influx to neurotransmitter release by measuring presynaptic action potential-evoked (AP-evoked) Ca2+ transients simultaneously with postsynaptic potentials. Ca2+ influx was manipulated by changing extracellular [Ca2+] (Ca(ext)) to alter the flux per channel or by reducing the number of open Ca2+ channels with omega-conotoxin GVIA (omega-CTX). When Ca(ext) was reduced, the exponent of the power relationship relating release to Ca2+ influx was 4.16 +/- 0.62 (SD; n = 4), consistent with a biochemical cooperativity of approximately 5. In contrast, reducing influx with omega-CTX yielded a power relationship of 1.7 +/- 0.44 (n = 5) for Ca(ext) of 1.8 mM and 2.12 +/- 0.44 for Ca(ext) of 0.45 mM (n = 5). Using geometrically realistic Monte Carlo simulations, we tracked Ca2+ ions as they entered through each channel and diffused in the terminal. Experimental and modeling data were consistent with two to six channel openings per AZ per AP; the Ca2+ that causes fusion of a single vesicle originates from one or two channels. Channel cooperativity depends mainly on the physical relationship between channels and vesicles and is insensitive to changes in the non-geometrical parameters of our model.

Brevity of the Ca2+ Microdomain and Active Zone Geometry Prevent Ca2+-Sensor Saturation for Neurotransmitter Release

The brief time course of the calcium (Ca2+) channel opening combined with the molecular-level colocalization of Ca2+ channels and synaptic vesicles in presynaptic terminals predict sub-millisecond calcium concentration ([Ca2+]) transients of ≥100 μM in the immediate vicinity of the vesicle. This [Ca2+] is much higher than some of the recent estimates for the equilibrium dissociation constant of the Ca2+ sensor(s) that control neurotransmitter release, suggesting release should be close to saturation, yet it is well known that release is highly sensitive to changes in Ca2+ influx. We show that due to the brevity of the Ca2+ influx the binding kinetics of the Ca2+ sensor rather than its equilibrium affinity determine receptor occupancy. For physiologically relevant Ca2+ currents and forward Ca2+ binding rates, the effective affinity of the Ca2+ sensor can be several-fold lower than the equilibrium affinity. Using simple models, we show redundant copies of the binding sites increase effective affinity of the Ca2+ sensor for release. Our results predict that different levels of expression of Ca2+ binding sites could account for apparent differences in Ca2+ sensor affinities between synapses. Using Monte Carlo simulations of Ca2+ dynamics with nanometer resolution, we demonstrate that these kinetic constraints combined with vesicles acting as diffusion barriers can prevent saturation of the Ca2+-sensor(s) for neurotransmitter release. We further show the random positioning of the Ca2+-sensor molecules around the vesicle can result in the emergence of two distinct populations of the vesicles with low and high release probability. These considerations allow experimental evidence for the Ca2+ channel-vesicle colocalization to be reconciled with a high equilibrium affinity for the Ca2+ sensor of the release machinery.

Consequences of molecular-level Ca2+ channel and synaptic vesicle colocalization for the Ca2+ microdomain and neurotransmitter exocytosis: a monte carlo study

Morphological and biochemical studies indicate association between voltage-gated Ca2+ channels and the vesicle docking complex at vertebrate presynaptic active zones, which constrain the separation between some Ca2+ channels and vesicles to 20 nm or less. To address the effect of the precise geometrical relationship among the vesicles, the Ca2+ channel, and the proteins of the release machinery on neurotransmitter release, we developed a Monte Carlo simulation of Ca2+ diffusion and buffering with nanometer resolution. We find that the presence of a vesicle as a diffusion barrier alters the shape of the Ca2+ microdomain of a single Ca2+ channel around the vesicle. This effect is maximal in the vicinity of the vesicle and depends critically on the vesicle's distance from the plasmalemma. Ca2+-sensor(s) for release would be exposed to markedly different [Ca2+], varying by up to 13-fold, depending on their position around the vesicle. As a result, the precise position of Ca2+-sensor(s) with respect to the vesicle and the channel can be critical to determining the release probability. Variation in the position of Ca2+-sensor molecule(s) and their accessibility could be an important source of heterogeneity in vesicle release probability.

Competition between phasic and asynchronous release for recovered synaptic vesicles at developing hippocampal autaptic synapses

Developing hippocampal neurons in microisland culture undergo rapid and extensive transmitter release-dependent depression of evoked (phasic) excitatory synaptic activity in response to 1 sec trains of 20 Hz stimulation. Although evoked phasic release was attenuated by repeated stimuli, asynchronous (miniature like) release continued at a high rate equivalent to approximately 2.8 readily releasable pools (RRPs) of quanta/sec. Asynchronous release reflected the recovery and immediate release of quanta because it was resistant to sucrose-induced depletion of the RRP. Asynchronous and phasic release appeared to compete for a common limited supply of release-ready quanta because agents that block asynchronous release, such as EGTA-AM, led to enhanced steady-state phasic release, whereas prolongation of the asynchronous release time course by LiCl delayed recovery of phasic release from depression. Modeling suggested that the resistance of asynchronous release to depression was associated with its ability to out-compete phasic release for recovered quanta attributable to its relatively low release rate (up to 0.04/msec per vesicle) stimulated by bulk intracellular Ca2+ concentration ([Ca2+]i) that could function over prolonged intervals between successive stimuli. Although phasic release was associated with a considerably higher peak rate of release (0.4/msec per vesicle), the [Ca2+]i microdomains that trigger it are brief (1 msec), and with asynchronous release present, relatively few quanta can accumulate within the RRP to be available for phasic release. We conclude that despite depression of phasic release during train stimulation, transmission can be maintained at a near-maximal rate by switching to an asynchronous mode that takes advantage of a bulk presynaptic [Ca2+]i.

Protein ground state candidates in a simple model: an enumeration study

The concept of the reduced set of contact maps is introduced. Using this concept we find the ground state candidates for a hydrophobic-polar lattice model on a two-dimensional square lattice. Using these results we exactly enumerate the native states of all proteins for a wide range of energy parameters. In this way, we show that there are some sequences which have an absolute native state. Moreover, we study the scale dependence of the number of members of the reduced set, the number of ground state candidates, and the number of perfectly stable sequences by comparing the results for sequences with lengths of 6 up to 20.