Cell Systems

Predicting single-cell transcriptional responses to perturbations is central to building the virtual cell, yet recent benchmarks show that simple baseline methods often outperform complex models, and model comparisons depend on the evaluation metric. Most studies assume that preprocessed RNA sequencing data are reliable ground truth for both training and evaluation. Here, we test this assumption by measuring the reliability of perturbations and their alignment with shared perturbation responses, classifying each perturbation as specific, shared, or unreliable. Among 7,170 perturbations from 29 datasets, 65% are unreliable, 11% shared, and 24% specific. Applying these quality labels to published benchmarks shows that model comparisons depend on perturbation quality. Training with reliable perturbations alone matches or outperforms full-data performance while using 55% of all training perturbations. Our framework also enables prospective experimental design: for most perturbations, a 28-cell pilot experiment accurately predicts how many cells a full screen needs to be reliable.

A central goal of conducting omics measurements is to understand how molecular features inform higher-order cell- and tissue-level phenotypes. In particular, multi-omics offers insights into how information encoded by the genome is coordinated through biological layers, resulting in functional outputs1. Due to myriad post-transcriptional regulatory processes, the coordination between mRNA and protein cannot be simply reduced to gene-wise correlation. Yet, both modalities have been shown to serve as representations of biological state, and multi-omics integration has been used to improve these representations. Multi-omics approaches typically do not focus on how mRNA and protein features coordinate, but rather use the additional information for improved prediction or feature selection. Here, instead, we showed that standard linear machine learning models provide an understanding of transcriptomic and proteomic coordination in the context of a biological phenotype of interest, in this case cancer metastasis. We find that, in the context of metastasis, a select subset of proteomic features--reflecting a more concentrated signal relative to the broadly distributed transcriptomic signal--offers additional information to that encoded by transcriptomics, as demonstrated by improved model performance when integrating the two modalities and the relative feature importance of proteomics. Top features show a depletion of gene-product overlap across modalities, indicating that the model primarily leverages instances in which the two modalities are providing complementary information with respect to phenotype. However, in instances when both modalities are selected for a given gene product, there is high information consistency that synergistically bolsters phenotype prediction. Altogether, by using model fits that relate both modalities to phenotype, we observe a nuanced coordination of protein and mRNA, in which both modalities tend to provide consistent information about phenotype, yet benefits remain to incorporating a combination of both complementary and reinforcing signals across modalities.

In recent years, high-throughput perturbational datasets have become an important tool for rapidly characterizing the function of large collections of chemical compounds. To overcome the biological and technical noise in these experiments, researchers have used consensus signatures - averages of multiple experiments - to summarize the effects of perturbations. In this work, we demonstrate that consensus signatures on the L1000 Connectivity Map show a pervasive similarity bias: as more signatures are averaged, the resulting consensus signatures are increasingly similar to each other, regardless of whether the signatures are related. We show that the distribution of Pearsons correlation changes as a function of the number of signatures averaged. The artifactual similarity bias is caused by skewness in the data and a consequence of using median normalization on non-normal distributions. Furthermore, we show that mean normalization can partly remedy this similarity bias and improve power to identify associations. The similarity bias introduced by consensus signatures is an important potential confounder of analysis of perturbational datasets, and our practical solution could easily be applied by practitioners in the field to improve the analysis of the L1000 Connectivity Map.

Predictive logical models have made substantial contributions to biomedical research by facilitating the identification of regulatory mechanisms driving biological processes. Yet, their construction remains a major bottleneck, as it typically relies on extensive prior knowledge at the genome scale. The emergence of high-throughput omics technologies has opened new opportunities for data-driven logical model construction by providing high-resolution observations of biological states. In parallel, formal constraint-based approaches have emerged, reformulating executable model inference as a constraint-satisfaction problem. However, transforming omics observations into biologically meaningful Boolean state abstractions and dynamical constraints suitable for logical synthesis engines remains challenging, particularly for poorly understood biological systems and multi-condition experimental settings. To address these challenges, we introduce scBOLT, a semi-automated framework for Boolean constraint engineering and logical model inference from multi-condition single-cell transcriptomes. Beyond these methodological contributions, scBOLT substantially improves computational tractability, enabling the inference of Boolean networks comprising several hundred components within practical runtimes.

Genetic perturbations can reshape cell populations by altering the relative abundance of specific cell types and states within the profiled population, including increases, decreases and states that become detectable after perturbation. Pooled single-cell screens, such as Perturb-seq, measure such responses at scale. However, only a small fraction of possible perturbations can be tested experimentally. A central challenge is therefore to predict compositional shifts induced by unseen perturbations. Many perturbation-prediction methods do not directly optimize for this outcome; instead, they predict gene-expression responses and infer cell-type and cell-state composition downstream. Surprisingly, we find that even models that accurately predict perturbation-induced changes in average gene expression perform poorly at forecasting these compositional shifts. To address this gap, we present Pop-Corn, a method that directly predicts how a perturbation reshapes cell-type composition without reconstructing gene expression. In the primary T-cell benchmark, Pop-Corn predicted the overall cell-state composition of held-out perturbations more accurately than the evaluated expression-prediction pipelines, while better preserving the diversity of observed cell states. We further extend Pop-Corn to intact tissue, where it predicts perturbation-induced cell-type proportion changes in local cellular neighborhoods and uses attention patterns to generate hypotheses about context-dependent cellular interactions. Retrospective virtual screens support the use of Pop-Corn to prioritize perturbations for experimental follow-up according to their predicted effects on cell-state composition.

Machine learning models in single-cell biology increasingly forecast differentiation, reprogramming and therapeutic response from early transcriptomic profiles. Testing whether a model has learned real biology requires held-out cells. Single-cell data, however, are grouped: cells from the same clone, patient or batch share the same label. A random split therefore places relatives of each test cell, carrying its label, in the training set, and a model can score well by memorizing a relative instead of learning a transferable rule. Grouped validation removes this leakage but leaves far fewer independent units behind each error bar. Here we show how to estimate this leakage before training any model, from two properties of the data: exposure, the fraction of test cells with relatives in training, and retrievability, how often a nearest-neighbor search returns such a relative rather than an unrelated cell. Across lineage-barcoded and patient data, exposure determines whether a random split opens a leakage channel, and retrievability determines how much it can inflate the score. The inflation is negligible where cell state has decoupled from ancestry, much larger where clonal sisters remain close in expression space, and in a patient cohort large enough to overturn a clinical conclusion. We also provide eakcheck, which computes both properties in seconds, before the outcome model is fitted.

Biological systems use protein circuits to organize cellular activities in space and time, but engineering synthetic dynamics is challenging due to stochastic effects of genetic and biochemical variation on circuit behavior. Genetically encoded oscillators (GEOs) built from bacterial MinDE-family ATPase and Activator modules generate fast orthogonal protein waves in eukaryotic cells, providing an experimental model system for genetic and biochemical coordination of synthetic protein dynamics. Here, we use budding yeast to experimentally define and model phase portraits that reveal how the breadth of frequencies and amplitudes available to a GEO are genetically controlled by ATPase and Activator expression levels and noise. GEO amplitude is encoded by ATPase absolute abundance, making it sensitive to extrinsic noise on a population level. In contrast, GEO frequency is remarkably stable because it is controlled by the Activator:ATPase ratio and thus affected primarily by intrinsic noise. These features facilitate noise-guided design of different expression strategies that act as filters on GEO waveform, enabling us to construct clonal populations that oscillate at different frequencies as well as independently tune frequency and amplitude variation within a single population. By characterizing 169 biochemically distinct GEOs, we provide a rich assortment of phase portraits as starting points for application of our waveform engineering approach. Our findings suggest noise-guided design may be a valuable strategy for achieving precision control over dynamic protein circuits.

The precise delineation of cell types is fundamental to single-cell transcriptomics, yet current clustering pipelines often violate an axiomatic principle: hierarchical consistency. Existing methods measure cell-to-cell distances within a fixed global feature space, disregarding the fact that biological distinctions are inherently context-dependent lineage separation requires different gene programs than subtype resolution. Mathematically, this implies that the similarity metric itself should not be a static functional, but a pair-dependent energy functional evaluated within a specific Hilbert subspace determined by the biological comparison at hand. The challenge lies in the fact that allowing pair-dependent metrics typically destroys the global geometric consistency required for downstream analysis, unless the family of Hilbert subspaces is given strong biological structure. To resolve this geometric dilemma, we introduce GeCCo (Gene Co-expression Constructed identity), which constructs identities by projecting cells onto a rigorously derived hierarchy of gene programs. To construct this hierarchy, GeCCo first quantifies Boolean regulatory logic via the{phi} coefficient, and subsequently employs a greedy topological inference to organize genes based on their synergistic and antagonistic relationships. Benchmarking on human immune atlases demonstrates that GeCCo achieves superior hierarchical consistency, ensuring that globally inferred cell identities rigorously match locally refined subtypes. Furthermore, in pancreatic endocrine progenitors, GeCCo resolves a hidden mitotic bridge state, suggesting a concentrated division phase prior to differentiation. Ultimately, GeCCo shifts the paradigm from ad hoc clustering to programmatic cell typing, offering a mathematically grounded framework for scalable atlases of cellular discovery.

It is widely accepted in genomics that coexpression of RNA transcripts suggests a commonality of function. This intuition is explicitly leveraged in machine learning methods that predict gene function, where it is often combined with other features such as protein interactions and sequence similarity. For example, including coexpression data from human tissue expression boosts performance for predicting Gene Ontology annotations. However, the biological underpinnings of this observation have not been well-investigated. Building on earlier results from our group, in this work we show that gene function is predictable from coexpression substantially because it reflects differences in expression between cell types, and these differences are also intrinsic to the ground truth labels. Using simulations and analyses of real data, we show that variance in the cellular composition of bulk samples impacts function learnability and attribute this to cell type marker gene content in the GO terms. We further show that cell type profiles, where the relationship between gene expression and cell type is made transparent, are effective for predicting gene function while increasing interpretability. These results indicate that function prediction models trained on bulk coexpression are largely limited to cell-type-level resolution rather than fine-grained biochemical function, with direct consequences for how such predictions should be interpreted.

Predicting how a genetic perturbation reshapes a cells transcriptome is a central goal of computational biology. Previous studies report that the mean response across training perturbations rivals specialized models on standard accuracy metrics, even though it cannot distinguish which perturbation occurred. Across 2,270 CRISPRi perturbations measured in each of six cell lines, we show that this apparent paradox reflects a conserved organization of perturbation responses. Perturbations span a continuum from responses strongly aligned with the mean to more targeted responses that depart from it. Crucially, a perturbations position along this continuum is conserved across cell lines (Kendalls W = 0.59) and predictable from STRING protein-interaction embeddings (R2 = 0.35). We formalize this structure with COMPASS, an interpretable linear model that decomposes each response into shared and gene-specific components and estimates them separately. The shared-response component is modeled as a cell-line-wide response scaled by a perturbation-specific coefficient. This coefficient is strongly conserved across cell lines. The residual gene-specific component--which is moderately conserved across cell lines--recovers pathway-level programs. COMPASS outperforms scGPT, CPA, GEARS, GenePert, and SO_SCPLOWTATEC_SCPLOW in both response accuracy (de-biased Pearson delta 0.34 vs. [≤] 0.32) and perturbation discrimination (cosine PDS gain 0.23 vs. [≤] 0.08). These results recast perturbation prediction across cellular contexts as component-wise inference, with each component estimated from the evidence best suited to it.

Combination therapies are central to cancer treatment, but exhaustive screening is impractical. We introduce PHAROS, a framework that turns a pretrained single-cell perturbation model into a target-directed search engine for drug combinations. PHAROS predicts how a cell population changes under a drug, one drug at a time, then chains these predictions together to simulate drug combinations. It scores each simulated outcome against the desired target state and uses a search algorithm to find the most promising combinations, all without retraining the underlying model. Across two independent combinatorial perturbation datasets, PHAROS recovered exact or mechanism-matched two-drug responses in cell lines, both seen and unseen during model training. Its rankings were specific to the requested conversion and were not explained by single-drug effects, additive effects, or shared mechanism of action. In exploratory analyses of patient-derived metastatic HR+/HER2$-$ breast tumors and basal cell carcinoma (BCC), PHAROS prioritized FDA-approved regimens, distinguished combinations by their predicted tumor-versus-immune objective profiles, and nominated pathway-level hypotheses, while explicitly identifying both tumor cohorts as outside the model's supported distribution. PHAROS provides a modular route from pretrained virtual-cell models to inverse, single-cell combination-screening platforms.

Mechanistic models of biological networks predict but cannot be calibrated at scale, while learned representations describe but expose no dynamical operator. Both lack an inverse step mapping data onto the few dynamical quantities the data determine. We develop Unified Network Spectroscopy (UNS) as that operator and establish when it is well posed. Reading fluctuations about a physiological setpoint through the Lyapunov relation of the linear noise approximation, UNS returns dominant timescales and frequencies, the subspace on which they act, and the projected noise. The inversion itself is established; our contribution is architectural. Three results govern its success. First, biochemical noise is generically anisotropic, and extrinsic modulation makes the dynamics non-normal. Assuming isotropy can invert mode rankings and destroy precision-matrix topology recovery; we give the exact, checkable condition under which recovery survives. Second, equal-time covariance identifies neither the dynamics nor their spectrum: the residual freedom is precisely the probability current that non-equilibrium living systems generically carry. Third, population snapshots replace time series only under ergodicity; heterogeneity biases inferred relaxation rates downward, making systems appear closer to bifurcation, whatever the sample size. Lagged covariance or per-cell velocity restores identifiability and turns noise structure into a recovered, falsifiable quantity. We provide a three-regime ergodicity taxonomy with diagnostics and show that the variational operator remains well posed when samples are fewer than variables. The output is a typed interface contract for hierarchical model composition, carrying its uncertainty, attainment and validity domain, with controllability inheritance as a falsification test.

Spatial transcriptomics enables inference of cell-cell communication through ligand-receptor (LR) interactions, but current prioritization strategies often rely on expression strength or interface-associated enrichment without explicitly modeling tissue geometry. As a result, interactions associated with population interfaces are frequently interpreted as spatially localized even when their underlying expression is broadly distributed. Here, we present a geometry-aware framework for LR prioritization that explicitly separates interface structure from spatial localization within a locked and reproducible analysis pipeline. We quantify interface-associated communication using a distance-weighted boundary score defined on a spatial neighbor graph, evaluate interface specificity using a label-permutation null model that preserves spatial geometry, and compute an LR-specific localization score that captures the proximity of ligand and receptor expression to the corresponding interface. This framework distinguishes interface-associated compatibility from interaction-level spatial concentration. Across spatial transcriptomics datasets from breast cancer, colorectal cancer, melanoma, and pancreatic ductal adenocarcinoma, interface-aware ranking consistently recovers pathway families associated with extracellular matrix, adhesion, inflammatory, and immune-related processes. However, interface enrichment frequently shows limited separation from the null model, indicating that interface structure alone does not establish spatial specificity. Incorporating geometric localization substantially alters LR prioritization, distinguishing interactions that are concentrated near interfaces from those that are more diffusely distributed. Under a fixed, deterministic pipeline applied identically across datasets without parameter tuning, discrete spatial communication regimes were not reproducibly recovered. Instead, variation across samples is more consistently captured as continuous differences in geometry-aware attenuation, reflecting the degree to which inferred interactions are spatially constrained by tissue architecture. Together, these results demonstrate that interface-associated enrichment and spatial localization are distinct properties of inferred LR interactions, and that accurate interpretation of spatial communication requires explicit modeling of tissue geometry. Under this framework, tumor communication is more consistently described as a continuum of spatial constraint.

Transforming growth factor-{beta} (TGF{beta}) ligands and receptors interact with overlapping selectivity to form different signaling complexes. This many-to-many wiring provides versatility in how ligands are perceived and cell types become activated, but the systems-level impact of coreceptors, which bind TGF{beta} ligands yet do not signal, remains unresolved. We examined the role of canonical surface-bound coreceptors by numerically simulating TGF{beta} ligand-receptor systems of increasing complexity and pairing results with single-cell measurements in acutely stimulated cells. Using sampled combinations of biologically plausible rate parameters and initial conditions, we find that coreceptors inhibit downstream signaling six times more often than they promote it, and the models reconcile this inhibitory bias with the prevailing view that coreceptors either promote or inhibit signaling. In multi-ligand systems, coreceptor inhibition causes stimuli to be perceived differently by altering competition, and this threshold is confirmed in cells engineered to induce a single coreceptor heterogeneously. Coreceptors are among the most variably expressed components of TGF{beta} signaling, suggesting that cells exploit coreceptors as a way to switch active ligand-receptor landscapes with a single gene.

Discovering biochemical circuits that exhibit a desired behavior is an outstanding problem in biological engineering. The traditional approach of enumerating every possible circuit topology becomes intractable for circuits with more than four components due to combinatorial scaling of the search space. Here, we use Monte Carlo Tree Search (MCTS), a reinforcement learning (RL) algorithm, to optimize circuit topology for a target phenotype by approaching circuit design as a sequence of assembly decisions. Our RL-based design framework, which we call CircuiTree, efficiently and comprehensively finds robust designs for three-component oscillators by prioritizing sparsity. CircuiTree can also infer candidate network motifs from its search results, producing similar results to enumeration. Using parallel MCTS, we scale this workflow up to five components and find that highly fault-tolerant designs use a novel strategy, which we call motif multiplexing. Multiplexed circuits contain many overlapping network motifs that each activate in different mutational scenarios. The evolutionary robustness of multiplexing may explain the ubiquity of multiple sub-oscillators in circadian clock circuits. Overall, CircuiTree provides the first scalable computational platform for designing biochemical circuits.

Microbiome measurements are commonly judged by how faithfully they reconstruct an underlying biological profile, using taxonomic similarity, prediction accuracy, or compositional distance. Yet faithfulness and usefulness are not the same: a distorted measurement may still support the correct intervention, while a small error can be decisive when it falls near the boundary between two treatment choices. We therefore evaluate measurement quality by whether information loss changes the decisions and outcomes that remain attainable. We treat the intervention-relevant state of the gut community as a hidden quantity seen only through an imperfect, information losing measurement, mapped by a decision rule to a graded intervention. Each decision rule is evaluated by net host benefit and ecological margin, generating a feasible set of attainable outcome pairs. We define the measurement-induced decision loss (MIDL) as the maximum net host benefit lost under imperfect measurement relative to perfect observation. For scalar affine information loss, increasing degradation is a Blackwell garbling; under post-processing closure, more informative measurement weakly expands the feasible set and cannot reduce optimal host benefit. In the binary-state model, MIDL increases linearly while an observation-dependent decision rule remains optimal and plateaus once the optimal policy becomes signal-independent. We extend the framework to finite states, state-specific loss, multi-segment stool aggregation, and ecological-margin constraints. The framework distinguishes measurement faithfulness from decision usefulness and provides a principled basis for comparing stool and richer measurements by the interventions and trade-offs they support.

Generative probabilistic models of biological sequences have widespread existing and potential applications in analyzing, predicting and designing proteins, RNA and genomes. To test the predictions of such a model experimentally, the standard approach is to draw samples, and then synthesize each sample individually in the laboratory. However, often orders of magnitude more sequences can be experimentally assayed than can affordably be synthesized individually. In this article, we propose instead to use stochastic synthesis methods, such as mixed nucleotides or trimers. We describe a black-box algorithm for optimizing stochastic synthesis protocols to produce approximate samples from any target generative model. We establish theoretical bounds on the methods performance, and validate it in simulation using held-out sequence-to-function predictors trained on real experimental data. We show that using optimized stochastic synthesis protocols in place of individual synthesis can increase the number of hits in protein engineering efforts by orders of magnitude, e.g. from zero to a thousand.

Multi-omic signatures are widely used in biomarker discovery, precision oncology, and systems biology, yet they are typically treated as vectors or composite scores that collapse intrinsically multidimensional biological organization into one-dimensional summaries. As a result, their internal structure, contextual dependencies, and mechanistic coherence remain largely inaccessible. Here, we introduce a geometric framework that reconceptualizes omic signatures as multidimensional informational entities whose biological meaning arises from structural organization rather than molecular membership alone. Each signature is embedded in a shared latent space integrating regulatory, phenotypic, microenvironmental, immune, and clinical constraints, and represented as a convex polytope. This representation preserves internal organization and enables intrinsic geometric measurements--including barycenter distance, volume, anisotropy, and asymmetry--that quantify concordance, divergence, and latent complexity. We apply this framework to 24,796 metabolic regulatory circuitries reconstructed across 32 TCGA cancer types, encoded as paired regulatory and metabolic signatures in an 18-dimensional latent space. Geometric analysis shows that discordance predominates: most circuitries occupy strong or extreme discordance regimes and display high-dimensional, frequently asymmetric geometries, whereas fully concordant circuitries are rare and structurally constrained. These geometric phenotypes stratify metabolic pathways and superfamilies in reproducible, non-uniform patterns that are not detectable with vector- or network-based representations. By transforming omic signatures into measurable geometric objects, this framework enables principled comparison, de-redundancy, and mechanistic interpretation of multi-omic biomarkers, providing a scalable approach for analyzing complex regulatory systems across cancer and beyond. All geometric representations and derived descriptors are available through the SigPolytope Shiny application (https://sigpolytope.shinyapps.io/geometricatlas/).

Pooled Perturb-seq screens are often interpreted through a linear steady-state model in which each targeted gene supplies an additive perturbation input. Although this formulation identifies an operator on the observed response subspace in closed form, identifiability alone does not establish that inferred operators generalize to new perturbation targets. We evaluated fixed program-space additive inputs in three Perturb-seq screens -- Replogle K562 essential, Replogle RPE1 essential, and Jost 2020 -- using target-grouped nested cross-validation and matched signal-to-noise linear-truth controls. For inverse operator prediction, real held-out error remained at the predict-zero baseline (0.96, 1.00, and 1.00), whereas matched linear controls achieved 0.18, 0.53, and 0.72. This discrepancy persisted across random K562 target panels, an overdetermined Jost dimension sweep, direction-only fits, ridge and truncated-SVD estimators, and five linear ground-truth ensembles. Forward target-held-out prediction was more nuanced: fixed, footprint, and learned linear encodings modestly beat a predict-training-mean baseline in K562 but not in RPE1 or Jost, and the K562 gains recovered only about 7% of the gap to the matched linear control. On the tested screens, fitted fixed program-space operators should not be interpreted as validated target-generalizing regulatory operators. We release anchor-op, a reproducible toolkit for target-grouped evaluation, identifiability checks, and matched-geometry positive controls for perturbation-response models.

The true accuracy of a machine learning model is a population-level statistic that cannot be observed directly. In practice, predictor performance is estimated against one or more test datasets, and the accuracy of this estimate strongly depends on how well the test sets represent all possible unseen datasets. Here we present paired evaluation, a simple approach for increasing the robustness of performance evaluation by systematic pairing of test samples, and use it to evaluate predictors of drug response in breast cancer cell lines and of disease severity in patients with Alzheimers Disease. Our results demonstrate that the choice of test data can cause estimates of performance to vary by as much as 30%, and that paired evaluation makes it possible to identify outliers, improve the accuracy of performance estimates in the presence of known confounders, and assign statistical significance when comparing machine learning models.