Abstract Plasticity enables cells to change their gene expression state in the absence of a genetic change. At the single-cell level, these gene expression states can persist for different lengths of time which is a quantitative measurement referred to as gene expression memory. Because plasticity is not encoded by genetic changes, these cell states can be reversible, and therefore, are amenable to modulation by disrupting gene expression memory. However, we currently do not have robust methods to find the regulators of memory or to track state switching in plastic cell populations. Here, we developed a lineage tracing-based technique to quantify gene expression memory and to identify single cells as they undergo cell state transitions. Applied to human melanoma cells, we quantified long-lived fluctuations in gene expression that underlie resistance to targeted therapy. Further, we identified the PI3K and TGF-β pathways as modulators of these state dynamics. Applying the gene expression signatures derived from this technique, we find that these expression states are generalizable to in vivo models and present in scRNA-seq from patient tumors. Leveraging the PI3K and TGF-β pathways as dials on memory between plastic states, we propose a “ pretreatment” model in which we first use a PI3K inhibitor to modulate the expression states of the cell population and then apply targeted therapy. This plasticity informed dosing scheme ultimately yields fewer resistant colonies than targeted therapy alone. Taken together, we describe a technique to find modulators of gene expression memory and then apply this knowledge to alter plastic cell states and their connected cell fates.

Abstract In the noisy cellular environment, RNAs and proteins are subject to considerable stochastic fluctuations in copy numbers over time. As a consequence, single cells within the same isoclonal population can differ in their expression profile and reside in different phenotypic states. The dynamic nature of this intercellular variation, where individual cells can transition between different states over time makes it a particularly hard phenomenon to characterize. Here we propose a novel fluctuation-test approach to infer the kinetics of transitions between cell states. More specifically, single cells are randomly drawn from the population and grown into cell colonies. After growth for a fixed number of generations, the number of cells residing in different states is assayed for each colony. In a simple system with reversible switching between two cell states, our analysis shows that the extent of colony-to-colony fluctuations in the fraction of cells in a given state is monotonically related to the switching kinetics. Several closed-form formulas for inferring the switching rates from experimentally quantified fluctuations are presented. We further extend this approach to multiple cell states where harnessing fluctuation signatures can reveal both the topology and the rates of cell-state switching. In summary, our analysis provides a powerful approach for dissecting cell-state transitions based on a single time point measurement. This is especially important for scenarios where a measurement involves killing the cell (for example, performing single-cell RNA-seq or assaying whether a microbial/cancer cell is in a drug-sensitive or drug-tolerant state), and hence the state of the same cell cannot be measured at different time points.

Abstract Inference of gene regulatory networks from single-cell expression data, such as single-cell RNA sequencing, is a popular problem in computational biology. Despite diverse methods spanning information theory, machine learning, and statistics, it is unsolved. This shortcoming can be attributed to measurement errors, lack of perturbation data, or difficulty in causal inference. Yet, it is not known if kinetic properties of gene expression also cause an issue. We show how the relative stability of mRNA and protein hampers inference. Available inference methods perform benchmarking on synthetic data lacking protein species, which is biologically incorrect. We use a simple model of gene expression, incorporating both mRNA and protein, to show that a more stable protein than mRNA can cause loss in correlation between the mRNA of a transcription factor and its target gene. This can also happen when mRNA and protein are on the same timescale. The relative difference in timescales affects true interactions more strongly than false positives, which may not be suppressed. Besides correlation, we find that information-theoretic nonlinear measures are also prone to this problem. Finally, we demonstrate these principles in real single-cell RNA sequencing data for over 1700 yeast genes.

Although the challenge of gene regulatory network inference has been studied for more than a decade, it is still unclear how well network inference methods work when applied to real data. Attempts to benchmark these methods on experimental data have yielded mixed results, in which sometimes even the best methods fail to outperform random guessing, and in other cases they perform reasonably well. So, one of the most valuable contributions one can currently make to the field of network inference is to benchmark methods on experimental data for which the true underlying network is already known, and report the results so that we can get a clearer picture of their efficacy. In this paper, we report results from the first, to our knowledge, benchmarking of network inference methods on single cell E. coli transcriptomic data. We report a moderate level of accuracy for the methods, better than random chance but still far from perfect. We also find that some methods that were quite strong and accurate on microarray and bulk RNA-seq data did not perform as well on the single cell data. Additionally, we benchmark a simple network inference method (Pearson correlation), on data generated through computer simulations in order to draw conclusions about general best practices in network inference studies. We predict that network inference would be more accurate using proteomic data rather than transcriptomic data, which could become relevant if high-throughput proteomic experimental methods are developed in the future. We also show through simulations that using a simplified model of gene expression that skips the mRNA step tends to substantially overestimate the accuracy of network inference methods, and advise against using this model for future in silico benchmarking studies.

Abstract One of the most difficult and pressing problems in computational cell biology is the inference of gene regulatory network structure from transcriptomic data. Benchmarking network inference methods on model organism datasets has yielded mixed results, in which the methods sometimes perform reasonably well and other times fail to outperform random guessing. In this paper, we analyze the feasibility of network inference under different noise conditions using stochastic simulations. We show that gene regulatory interactions with extrinsic noise appear to be more amenable to inference than those with only intrinsic noise, especially when the extrinsic noise causes the system to switch between distinct expression states. Furthermore, we analyze the problem of false positives between genes that have no direct interaction but share a common upstream regulator, and explore a strategy for distinguishing between these false positives and true interactions based on noise profiles of mRNA expression levels. Lastly, we derive mathematical formulas for the mRNA noise levels and correlation using moment analysis techniques, and show how these levels change as the mean mRNA expression level changes.

Abstract In isogenic cell populations, cells can switch back and forth between different gene expression states. These expression states can be biologically relevant. For example, a certain expression state may cause a tumor cell to be resistant to treatment, while another state may leave it vulnerable to treatment. However, estimating the rates of state-switching can be difficult, because experimentally measuring a cell’s transcriptome often involves destroying the cell, so it can only be measured once. In this paper, we propose a computational method to estimate the rate of switching between expression states, given data from a Luria-Delbrück style fluctuation test that is experimentally simple and feasible. We then benchmark this method using simulated data to test its efficacy, with varying assumptions made about cell cycle timing distribution in the simulations.

ABSTRACT Mechanistic models of how single cells respond to different perturbagens can help integrate disparate big data sets or predict response to varied drug combinations. However, the construction and simulation of such models have proved challenging. Our lab previously constructed one of the largest mechanistic models for single mammalian cell regulation of proliferation and death (774 species, 141 genes, 8 ligands, 2400 reactions). However, this, as many other large-scale models, was written using licensed software (MATLAB) with intricate programming structure, impeding alteration, expansion, and sharing. Here, we generated a new foundation for this model, which includes a python-based creation and simulation pipeline converting a few structured text files into an SBML-compatible format. This new open-source model (named SPARCED) is high-performance- and cloud-computing compatible and enables the study of virtual cell population responses at the single-cell level. We applied this new model to a subset of the LINCS MCF10A Data Cube, which observed that IFNγ acts as an anti-proliferative factor, but the reasons why were unknown. After expanding the SPARCED model with an IFNγ signaling module (to 950 species, 150 genes, 9 ligands, 2500 reactions), we ran stochastic single-cell simulations for two different putative crosstalk mechanisms and looked at the number of cycling cells in each case. Our model-based analysis suggested, and experiments support that these observations are better explained by IFNγ-induced SOCS1 expression sequestering activated EGF receptors, thereby downregulating AKT activity, as opposed to direct IFNγ-induced upregulation of p21 expression. This work forms a foundation for increased mechanistic model-based data integration on a single-cell level, an important building block for clinically predictive mechanistic models.

Abstract One of the most difficult challenges in cancer therapy is the emergence of drug resistance within tumors. Sometimes drug resistance can emerge as the result of mutations and Darwinian selection. However, recently another phenomenon has been discovered, in which tumor cells switch back and forth between drug-sensitive and pre-resistant states. Upon exposure to the drug, sensitive cells die off, and pre-resistant cells become locked in to a state of permanent drug resistance. In this paper, we explore the implications of this transient state switching for therapy scheduling. We propose a model to describe the phenomenon and estimate parameters from experimental melanoma data. We then compare the performance of continuous and alternating drug schedules, and use sensitivity analysis to explore how different conditions affect the efficacy of each schedule. We find that for our estimated parameters, a continuous therapy schedule is optimal. However we also find that an alternating schedule can be optimal for other, hypothetical parameter sets, depending on the difference in growth rate between pre-drug and post-drug cells, the delay between exposure to the drug and emergence of resistance, and the rate at which pre-resistant cells become resistant relative to the rate at which they switch back to the sensitive state.

SUMMARY Non-genetic transcriptional variability at the single-cell level is a potential mechanism for therapy resistance in melanoma. Specifically, rare subpopulations of melanoma cells occupy a transient pre-resistant state characterized by coordinated high expression of several genes. Importantly, these rare cells are able to survive drug treatment and develop resistance. How might these extremely rare states arise and disappear within the population? It is unclear whether the canonical stochastic models of probabilistic transcriptional pulsing can explain this behavior, or if it requires special, hitherto unidentified molecular mechanisms. Here we use mathematical modeling to show that a minimal network comprising of transcriptional bursting and interactions between genes can give rise to rare coordinated high states. We next show that although these states occur across networks of different sizes, they depend strongly on three (out of seven) model parameters and require network connectivity to be ≤ 6. Interestingly, we find that while entry into the rare coordinated high state is initiated by a long transcriptional burst that also triggers entry of other genes, the exit from it occurs through the independent inactivation of individual genes. Finally, our model predicts that increased network connectivity can lead to transcriptionally stable states, which we verify using network inference analysis of experimental data. In sum, we demonstrate that established principles of gene regulation are sufficient to describe this new class of rare cell variability and argue for its general existence in other biological contexts.

Abstract During development, cells not only adopt specialized identities but also maintain those identities. Endoreduplication is thought to maintain cell identity. High concentrations of ARABIDOPSIS THALIANA MERISTEM LAYER1 (ATML1) specify giant cell identity and induce endoreduplication in sepals. How different concentrations of ATML1 can specify different identities remains unclear. Here, we show that high concentrations of ATML1 induce the biosynthesis of both long-chain and very long-chain fatty acids (LCFAs/VLCFAs), and these fatty acids are required for the maintenance of giant cell identity. Inhibition of VLCFA biosynthesis causes endoreduplicated giant cells to resume division and lose their identity, indicating that endoreduplication is not sufficient to maintain cell identity. Structural predictions suggest that LCFA-containing lipids bind to the START domain 2 of ATML1, causing ATML1 dimerization and its auto-activation. Our data and modeling imply that ATML1 induces biosynthesis of its own lipid ligands in a positive feedback loop, shedding light on the intricate network dynamics that specify and maintain giant cell identity. Teaser: Endoreduplicated cells in Arabidopsis thaliana sepals divide and de-differentiate in the absence of VLCFA biosynthesis.