Fluctuation structure predicts genome-wide perturbation outcomes

Fluctuation structure predicts genome-wide perturbation outcomes | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Fluctuation structure predicts genome-wide perturbation outcomes View ORCID Profile Benjamin Kuznets-Speck , Leon Schwartz , Hanxiao Sun , View ORCID Profile Madeline E. Melzer , Nitu Kumari , Benjamin Haley , View ORCID Profile Ekta Prashnani , Suriyanarayanan Vaikuntanathan , View ORCID Profile Yogesh Goyal doi: https://doi.org/10.1101/2025.06.27.661814 Benjamin Kuznets-Speck 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Benjamin Kuznets-Speck For correspondence: yogesh.goyal{at}northwestern.edu benks{at}northwestern.edu Leon Schwartz 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA 9 Data Science of Systems Biology, Technical University of Munich , Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Hanxiao Sun 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Madeline E. Melzer 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Madeline E. Melzer Nitu Kumari 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Benjamin Haley 5 Department of Ophthalmology, Université de Montréal , Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ekta Prashnani 6 NVIDIA , Santa Clara CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Ekta Prashnani Suriyanarayanan Vaikuntanathan 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA 7 Department of Chemistry, University of Chicago , Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Yogesh Goyal 1 Department of Cell and Developmental Biology, Feinberg School of Medicine, Northwestern University , Chicago IL, USA 2 Center for Synthetic Biology, Northwestern University , Chicago IL, USA 3 Robert H. Lurie Comprehensive Cancer Center, Northwestern University, Feinberg School of Medicine , Chicago IL, USA 4 NSF-Simons National Institute for Theory and Mathematics in Biology , Chicago IL, USA 8 Chan-Zuckerberg Biohub Chicago , LLC, Chicago IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yogesh Goyal For correspondence: yogesh.goyal{at}northwestern.edu benks{at}northwestern.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Pooled single-cell perturbation screens represent powerful experimental platforms for functional genomics, yet interpreting these rich datasets for meaningful biological conclusions remains challenging. Most current methods fall at one of two extremes: either opaque deep learning models that obscure biological meaning, or simplified frameworks that treat genes as isolated units. As such, these approaches overlook a crucial insight: gene co-fluctuations in unperturbed cellular states can be harnessed to model perturbation responses. Here we present CIPHER (Covariance Inference for Perturbation and High-dimensional Expression Response), a conceptual framework leveraging linear response theory from statistical physics to predict transcriptome-wide perturbation outcomes using gene co-fluctuations in unperturbed cells. We validated CIPHER on synthetic regulatory networks before applying it to 11 large-scale single-cell perturbation datasets covering 4,234 perturbations and over 1.36M cells. CIPHER robustly recapitulated genome-wide responses to single and double perturbations by exploiting baseline gene covariance structure. Importantly, eliminating gene-gene covariances, while retaining gene-intrinsic variances, reduced model performance by 11-fold, demonstrating the rich information stored within baseline fluctuation structures. Moreover, gene-gene correlations transferred successfully across independent experiments of the same cell type, revealing stereotypic fluctuation structures. Furthermore, CIPHER outperformed conventional differential expression metrics in identifying true perturbations while providing uncertainty-aware effect size estimates through Bayesian inference. Finally, most genome-wide responses propagated through the covariance matrix along approximately three independent and global gene modules. CIPHER underscores the importance of theoretically-grounded models in capturing complex biological responses, highlighting fundamental design principles encoded in cellular fluctuation patterns. INTRODUCTION Single-cell profiling technologies enable measurements of cell and tissue molecular makeup at an unprecedented resolution and scale 1 – 3 . Reference collections of the molecular profiles of individual cells across tissues, organisms, and conditions (“atlases”) comprise multimodal measurements, documenting variability within and across cell types through global consortia. Moving beyond the baseline states catalogued by these atlases, recent breakthroughs in experimental platforms, such as Perturb-seq and CROP-seq 4 – 11 , have enabled the investigation of genetically-perturbed cells by integrating single-cell RNA sequencing with CRISPR screens. Such approaches are transforming functional genomics by providing a systematic, scalable, and high-resolution view of how genetic perturbations drive phenotypic changes. These large-scale perturbation datasets lay the groundwork for mechanistic models that move beyond correlation, toward truly causal maps of gene regulation in single cells 12 , 13 . While powerful, high-content pooled single-cell perturbation screens face significant challenges, including few cells per perturbation in addition to the noise and sparsity inherent in single-cell data 2 , 12 , 14 , 15 . These experimental limitations impede reliable detection of genome-wide perturbation effects, especially when using conventional differential expression analyses and significance testing 12 , 16 . A suite of computational methods have been developed to tackle these limitations 12 . These approaches broadly fall into two categories: interpretable statistical frameworks that analyze existing perturbation data, and predictive models that forecast cellular responses to unseen perturbations. Statistical frameworks include MIMOSCA 6 and scMAGeCK-LR 17 , which estimate the average impact of each perturbation by fitting linear regression models. Other interpretable methods use matrix factorization (e.g., SVD, NMF, SMAF) to perform dimensionality reduction and identify gene programs or modules that respond coherently to perturbations 18 , 19 . A notable recent interpretable statistical model is TRADE 16 , which models true differential expression distributions while accounting for noise, using a “transcriptome-wide impact” metric to quantify perturbation effects including weak signals missed by conventional significance testing. At the same time, traditional linear regression models have some limitations: they typically treat genes as independent, cannot explicitly identify the genes driving the response (or the degree of impact per driver gene), and, with the exception of TRADE 16 , produce only point estimates of the regression coefficients and are not uncertainty-aware. Predictive models such as scGen 20 , GEARS 21 , and scGPT 22 , on the other hand, employ deep learning frameworks to model perturbation responses in latent spaces and can predict outcomes of unseen perturbations. These methods often present “black-box” solutions with complex non-linear mappings, limiting interpretability 23 of the underlying biological mechanisms, and typically require substantial amounts of training data. Despite progress, we still lack interpretable, physically-grounded approaches with the ability to model genome-wide responses primarily based on the structure of the unperturbed system 24 – 27 . Nevertheless, techniques from statistical physics help provide a universal framework to predict the effect of perturbations from baseline fluctuations in other contexts 24 – 27 . If similar physically-grounded models could explain genome-wide perturbation responses, what might that imply about the underlying design principles of biological responses? Here we present CIPHER ( C ovariance I nference for P erturbation and H igh-dimensional E xpression R esponse), a framework leveraging linear response theory, which has been foundational in describing complex physical phenomena ( Fig. 1 ) 28 – 30 . Linear response theory postulates that response to small perturbations is encoded in correlations present within the unforced/unperturbed system’s components. Without necessarily requiring a predefined notion of what constitutes the underlying experimental complexity, linear response theory non-trivially recapitulates experimental findings in a variety of contexts, including the Marcus theory of electron transport and Onsager’s reciprocal relations 28 – 30 . Here we use CIPHER to identify the determinants of genome-wide responses to genetic perturbations. Download figure Open in new tab Fig. 1: Schematic of the CIPHER workflow and applications to synthetic and real experimental datasets. (column 1) The gene-gene covariance matrix informs on changes to gene expression upon perturbation. (column 2) Perturb-seq interrogates thousands of perturbations by combining RNA-seq with CRISPR screens. (column 3) Linear response can apply to progressively complex synthetic regulatory networks. (column 4) Given covariance and expression changes, CIPHER has three complementary modalities. The forward problem predicts changes in expression, the reverse problem predicts the driving perturbation(s) and the framework is fully interpretable telling us 1) how each genes change is made up of contributions from both itself and every other gene and 2) how response is propagated along dominant or subdominant modes of the covariance. We first implemented CIPHER for proof-of-concept recapitulation of responses to perturbations on three classes of synthetic genetic networks where ground truth is known, including stochastic linear and non-linear models of gene regulation. We next applied CIPHER to ten large single-cell perturbation sequencing datasets 16 , 31 – 35 and one imaging Perturb-FISH dataset 10 , including 9 CRISPR interference (CRISPRi) gene knockdown and 2 CRISPR gene activation (CRISPRa) datasets, covering 4,288 perturbations, 1,347,778 cells, and 9 cell types ( Table 1 ). CIPHER robustly recapitulated transcriptome-wide response to single perturbations by using only the information on gene fluctuation covariances in unperturbed (control) single cells. View this table: View inline View popup Download powerpoint Table 1: Specifications for single-cell perturbation datasets used in this study. Importantly, eliminating either gene-gene covariances (while retaining individual gene count distributions across cells) in the unperturbed datasets dramatically decreased the model performance. These correlation structures also transferred robustly across independent experiments of the same cell type to predict responses. Together, these findings demonstrate the rich information contained within baseline gene-gene fluctuation architectures. Moreover, CIPHER outperformed (AUROC across all perturbations = 0.92) existing metrics for predicting true perturbations given the genome-wide change upon perturbation by leveraging Bayesian inference to measure effect size distributions for all possible perturbations from the correlation matrix. Furthermore, CIPHER enabled decomposition of genome-wide responses into latent regulatory gene programs (“soft modes”), enabling interpretability of our findings. Together, our results on CIPHER highlight the power of simple theoretical models inspired by statistical physics in probing and dissecting biological complexity. RESULTS Overview of CIPHER Gene correlations have long served as proxies for gene-gene interactions 36 – 38 . Recent work reveals these correlations contain meaningful information about gene regulatory network dynamics, particularly during bifurcations 39 . CIPHER builds on these insights, using correlations and a variant of the linear response theory to predict perturbation-induced expression changes. Linear response theory is a classical result in (quasi)equilibrium statistical mechanics which states that the response of a physical system to relatively small external forces is encoded in the correlations of the unforced system. This theory has successfully explained various physical phenomena, including electron transfer kinetics, NMR spectroscopy, dielectric response in polar liquids, and protein conformations and allostery 24 – 27 . We reasoned that high-content pooled single-cell perturbation screens naturally offer within- and across-gene fluctuations in the unperturbed populations. Moreover, individual gene perturbations can be thought of as “small” given the scale of the system (usually several thousand expressed genes at a given time). Furthermore, as we illustrate below, CIPHER’s success at recapitulating experimental cellular responses suggests a general separation of timescales between fast molecular transcription events and slower state-fate transitions that drive phenotypic changes. This rapid equilibration appears to be a generic feature of gene regulatory and mesoscopic biophysical systems 40 , 41 . According to a general formulation of linear response theory, population averages before and after a perturbation vector u is applied are coupled by linear combinations of pre-perturbation gene-gene correlations: where Δ X = u − 0 and the elements of the covariance matrix Σ are Σ ij = for δ X i = X i − 0 . We denote averages of gene expression X over cells as and distinguish between the perturbed and control ensembles with subscript u and 0 respectively. In a variety of synthetic gene expression networks and experimental Perturb-seq datasets ( Fig. 1 ), we investigated the utility of this framework using a three-pronged approach. First, we considered the forward problem of estimating perturbation responses given true single or double gene perturbations and fluctuation structure Σ. Next, we tackled the u = Σ −1 Δ X reverse problem of predicting perturbations from Σ and Δ X . We subsequently examined eigengenes—groups of co-expressed genes that captured the dominant expression patterns—to interpret response as propagating broadly on the gene level but localized to a handful of underlying modes or regulatory modules. Detailed description of CIPHER as well as corresponding derivations and applications are provided in Methods 1 . CIPHER on synthetic regulatory networks We initially validated CIPHER through proof-of-concept testing on a suite of synthetic gene networks (see Methods 2 ). These toy networks provided controllable architectures and interaction parameters with known ground truth for benchmarking. We began with linear systems of random gene regulatory networks (dx/dt = Ax + noise, where A is a matrix of gene-gene interactions and x is a vector of gene expression). We reasoned that CIPHER should accurately recapitulate perturbation responses within the linear regime (Case 1, see Methods 2.1 ). As Robert May demonstrated in their seminal study 42 , such networks remain stable for N genes when coupling strength, beyond which they bifurcate to become globally unstable. Consistently, increasing coupling strength ( Fig. 2A ) caused gene-gene correlations to become modular and fractal-like ( Fig. 2B ). The transition from the subcritical to supercritical regime is marked by a striking increase in the percent variance explained by the first principal component. Accordingly, a collapse of the dynamics onto a low-dimensional subspace (with low participation ratio) where “team” 43 , 44 structure emerges in the covariance as measured by variance explained and participation ratio—a measure of how many components actively contribute to the response ( Fig. 2C,D ). Biologically, this corresponds to the dynamics passing through an unstable fixed point, as would happen, for example, during a differentiation trajectory 45 . Importantly, CIPHER’s predicted responses to genetic perturbations exhibited relatively low coefficient of variation within the linear regime, but diverged at and beyond the critical coupling threshold ( Fig. 2E ). Download figure Open in new tab Fig. 2: Predicting genome-wide response to perturbations in synthetic regulatory networks. A) Random linear regulatory networks of N genes– nodes represent genes and edge thickness represent the (absolute) strength of interactions between genes. Subcritical, critical and supercritical networks are shown. B) Steady-state gene-gene correlations for the networks in A). C) Percent variance explained by the first principal component as a function of the gene-gene interaction parameter. D) Participation ratio (effective dimension) as a function of the scaled interaction parameter for different sized networks. E) Points are the typical error in linear response, ||Δ X − Σ u ||/||Δ X ||, across all genes in the three-regimes, averaged over 500 trajectories. Standard error bars shown. (One-sided Wilcoxon signed-rank tests: subcritical error < critical error, p-value = 7.6×10 −7 ; critical error < supercritical error, p-value = 3.0×10 −51 ). F) Non-linear networks with activating Hill function interactions, n = 10. G) Error in linear response as a function of effective interaction strength, G eff . Points represent individual simulation runs and lines guide the eye. H) Error in linear response as a function of perturbation magnitude for several different Hill coefficients (n). I) prototypical ‘teams’ network wherein groups of genes mutually activate within teams and inhibit across them, each gene has its own promoter whose bursting activity is modulated by nonlinear Hill-type reactions with all other genes. J) Time-series of the average expression for the two teams. K) Steady state correlations over the teams time-series. L) Empirical expression change and corresponding linear response predictions. Each point represents an individual gene in the network. We next tested CIPHER on a prototypical nonlinear dynamical gene expression system consisting of a network of Hill functions (Case 2, Fig. 2F , Methods 2.2 ). By perturbing a single gene in this nonlinear network, we could directly compare the observed response with CIPHER’s predictions. First, we fixed the Hill-function parameters and varied the gene-gene interaction parameter G, finding a sharp transition at G eff = G/G* = 1 where G* is the critical interaction above which the trivial homogeneous state (zero expression for all genes) becomes unstable (see Methods 2.2 for derivation and Fig. S1A,B for additional examples) ( Fig. 2G ). Increasing the perturbation magnitude corresponded to a subsequent increase in the relative error of linear response (||Δ X − Σ u ||/||Δ X ||, Fig. 2G ). We next fixed G and varied the Hill coefficient, n. As expected, CIPHER performed optimally under linear conditions (Hill coefficient, n=0), yielding the lowest error between actual and predicted values ( Fig. 2H ). However, as nonlinearity increased (Hill coefficient >1), error rates peaked during transition periods but remained low otherwise—a pattern consistent with the characteristic response curves of Hill function dynamics ( Fig. 2H ). We hypothesized that CIPHER’s prediction accuracy should decline under more challenging perturbation conditions. Specifically, we tested whether increasing either the perturbation magnitude ( u i ) for single gene knockouts or the total number of simultaneously perturbed genes would degrade performance. Indeed, CIPHER performance systematically decreased as we increased both the magnitude ( u i ) and total number of simultaneous perturbations ( u i , where i ∈ {1,2,3,4,5}) ( Fig. 2F-H and S1C ), confirming that the method’s linear assumptions become limiting under extreme perturbation conditions. Lastly, we tested CIPHER on a stochastic non-linear network capturing the recently described “teams” topology 43 , 44 , where genes on the same team mutually activate each other and there is inhibition between genes on different teams (Case 3, Fig. 2I , see Methods 2.3 for details on the model construction). Teams can emerge within gene regulatory networks and have been implicated in driving phenotypic transitions (e.g., epithelial-mesenchymal transition) in multiple cancers and during development 43 , 44 . Interestingly, despite the heterogeneity in the interaction parameters and a sparse interaction network, the synthetic system functions as a collective bi-stable switch toggling stochastically between two fixed states (high A, low B, and conversely) ( Fig. 2I,J ). We found that CIPHER readily captured ( R 2 = 1 − ||Δ X − Σ u ||/||Δ X || = 1. 0) the responses to perturbations (‘knocking out’ a gene by setting its individual gene expression to zero) in this nonlinear stochastic model, which incorporates transcriptional bursting and a Hill function response ( Fig. 2K,L ). Collectively, our proof-of-concept implementation of CIPHER on multiple toy synthetic gene networks demonstrates its versatility across diverse dynamical systems under controlled conditions. CIPHER captures genome-wide perturbation responses in real datasets Complex gene regulatory networks can be viewed as out-of-equilibrium systems with underlying correlations between genes. We questioned whether the CIPHER framework might recapitulate genome-wide responses to perturbations, which are typically a single gene knockdown or activation at a time. To address this question, we first systematically collected and performed quality control on existing single-cell perturbation datasets (see Methods 3,6 ). In total, we analyzed ten high-throughput perturbation datasets, containing 8 CRISPR interference gene knockdown and 2 CRISPR gene activation datasets and covering 4,234 perturbations, 1,342,852 cells, and 8 cell types ( Fig. 3A , Table 1 ). We calculated covariance matrices for unperturbed single cells from each experimental gene count matrix (raw absolute counts for each gene) and from corresponding scrambled control matrices which constitute a null model ( Fig. 3B ). We then implemented CIPHER on the forward problem of predicting transcriptome-wide response, optimizing for a single u i * corresponding to the true perturbation i of interest and measuring the response by calculating an R 2 value ( Fig. 3C,D ). An R 2 = 1 implies that we could perfectly transport the initial vector of gene expression 0 to the final measured target u . Remarkably, we found that CIPHER alone was able to significantly recapitulate the transcriptome-wide responses as compared to the null model for both activating and interfering perturbations, with some R 2 values approaching 1 ( Fig. 3E,F ). We have illustrated gene-wise concordance between CIPHER predictions and experimental data for a few representative test cases in Fig. 3G-J and S2 . We explored the mechanistic underpinnings of variability in R 2 values across datasets later in the results section. We note, however, that the per-gene values of R 2 had little dependence on the control cell expression ( Fig. S3A ). Download figure Open in new tab Fig. 3: The forward problem: to what extent can covariance structure explain response to known gene perturbations from Perturb-seq. A) We consider 10 Perturb-seq datasets, comprising CRISPRi (8 datasets) and CRISPRa (2 datasets) gene perturbations (outermost ring) for 8 different cell types (middle ring). The innermost ring shows the log-number of control (grey) and perturbed (colored) samples for each perturbation in each dataset. B) Show control gene expression matrices before and after shuffling rows and columns, as well as corresponding covariance matrices C) Distributions of R 2 values over all CRISPRi perturbations, optimizing over the known perturbation with either the real or shuffled covariance matrix (KS test stat = 0.881, p-value = 5.444×10 −79 ). D) Distributions of R 2 values over all CRISPRa perturbations, with real or shuffled covariance (KS test stat = 0.876, p-value = 0). E) Average R 2 values over all CRISPRi perturbations (black points, black line) and per dataset (colored points, grey lines) F) and analogously for CRISPRa (p-value calculated using combined one-sided Wilcoxon signed-rank tests including all 10 datasets: p-value = 0.00098). G) The top scoring perturbation for Tian_19b is PPP2R1A. We show the empirical Δ X (lines) as well as that predicted by linear response (dots). H) The top scoring perturbation for Tian_21b is TUSC1. I) Predicted Δ X = Σ u vs experimental data for PPP2R1A, each point is a gene, J) and similarly for TUSC1. K) Real and ‘mean-field’ (shuffled over cells, not genes) gene expression count matrices, and corresponding covariance matrices. L) Distribution of R 2 values over all datasets for mean-field and real covariance matrices (KS test stat = 0.654, p-value = 0). M) The average R 2 values for CRISPRi/a datasets with the real Σ compared to mean field. Dark points are averages over all CRISPRi (green) and CRISPRa (yellow) perturbations, and light points are per dataset averages. (Combined one-sided Wilcoxon signed-rank tests including all 10 datasets: p-value = 0.00098) N) Histograms of R 2 values over all double perturbations to different genes A and B. We compare the full double perturbation response to that achieved for shuffled X 0 , perturbations A or B alone, and the additive solution. (Shuffled vs Single A: KS stat = 0.5250, p-value = 1.823×10 −15 ; Shuffled vs Single B: KS stat = 0.5167, p-value = 5.753×10 −15 ; Shuffled vs Additive: KS stat = 0.6333, p = 6.627×10 −23 ; Single A vs Single B: KS stat = 0.0833, p-value = 8.012×10 −1 ; Additive vs True Σ: KS stat = 0.3750, p-value = 6.651×10 −08 ). O) Mean R 2 across the conditions in N) (Shuffled vs Single A: p-value = 7.9×10 −6 ; Shuffled vs Single B: p = 0.00044; Shuffled vs Additive: p-value = 2.6×10 −6 ; Additive vs True Σ: p-value = 9.9×10 −22 ). Orange points are over Norman19 perturbations, blue points over Tian_19a perturbations and the black points are the averages over both the datasets. P) Responses in a host study from inter-study covariances either of the same or different cell type. Q) R 2 values across neuron datasets. R 2 values using the true covariance matrix from the host study compared to those from a different dataset but same cell type. Each point is a perturbed gene. R) Average R 2 values across all pairs of datasets and conditions, using a covariance matrix from either a shuffled X 0 (across datasets), the mean-field approximation (across datasets), the unshuffled X 0 (across datasets) and the unshuffled X 0 from the same dataset. Thin lines correspond to individual dataset pairs (averages over all perturbations in the host dataset with a Σ of the same cell type in red or different cell types in orange) and points are averages over the typical R 2 values for these pairs of datasets. Note that there are only three lines connecting the across and within dataset real covariance for the same cell type. This is because we are comparing the three datasets shown in Q, and that there is no true covariance from the same dataset and different cell type. (Same cell type Shuffled vs. mean-field: p = 0.0020; Same cell type mean-field vs cross-dataset covariance: p-value = 0.0020; Same cell type cross-dataset covariance vs same dataset covariance: p-value = 0.016. Different cell type Shuffled vs. mean-field: p-value = 4.8×10 −7 ; Different cell type mean-field vs true: p-value = 4.8×10 −7 ). We next reasoned that if the co-fluctuations between the genes in the unperturbed population were indeed critical to the CIPHER framework, simply shuffling the gene counts column-wise—which retains the original mean and distribution for each gene individually but scrambles the covariance—should result in reduced R 2 values ( Fig. 3K ). This kind of shuffling operation, equivalent to a widely used ‘mean-field approximation’ in physical systems, indeed resulted in reduced R 2 values ( Fig. 3L,M and S3B ), confirming the dominant role of off-diagonal elements of the covariance matrix in predominantly driving transcriptome-wide response to activating and repressive genetic perturbations. Even optimization for any random gene u * in this mean-field model did not improve the R 2 values ( Fig. S3C) , further consolidating the dominant role of the gene-gene covariances in driving transcriptome-wide outcomes. We also implemented CIPHER on recently proposed optical screens like Perturb-FISH 10 , with similar high-confidence results ( Fig. S3D ). CIPHER captures response to double perturbations Pooled single-cell CRISPR screens provide rich datasets for investigating epistasis and gene interactions in cases where the same cell receives guides perturbing multiple genes. We implemented CIPHER to evaluate the relationships between the genes perturbed and the respective outcomes. For each double perturbation (2 genes perturbed simultaneously), we performed a similar optimization, but in this case over u i and u j instead of just a single u i . This optimization contains a nonlinear cross term that couples u i and u j through the off-diagonal element of the covariance matrix Σ ij , so that epistasis manifests itself in this model directly through Σ ij . When tested on two datasets with double perturbations ( Table 1 ), we found that the full model outperformed the additive solution u i + u j in a vast majority of cases ( Fig. 3N,O , S3E ). Biologically, these results suggest that many true double perturbations are to some extent epistatic. Moreover, as was the case for synthetic networks, as the number of simultaneous perturbations per cell gets larger, the additive performance decreases ( Fig. S3F ). Our results also highlight how CIPHER can be leveraged to discriminate between additive or epistatic relationships. Predictive power of covariances are preserved across studies of same cell type We tested whether the gene-gene co-fluctuations in unperturbed cells from one dataset or experiment can be informative in predicting perturbation response in another independent experiment of the same cell type or cell line. If so, it would suggest a universal underlying structure of covariance that carries essential information for predicting responses to perturbations. Conversely, gene-gene covariances in each experiment—regardless of whether cell type is shared—are sensitive to specific experimental conditions. To test this idea, we focused on three unrelated datasets with a shared cell type (neuron) ( Table 1 , Fig. 3A ). We calculated covariances from two of the datasets and used CIPHER to predict responses in the third “host” dataset (Tian21_a, Fig. 3P ). Remarkably, we found strong correlation (R 2 =0.73) between transcriptome-wide responses calculated from cross-dataset covariance matrices of the same cell type ( Fig. 3Q ) and typically much poorer performance from cross-dataset of a different cell type ( Fig. S3G) . Similar analysis across distinct cell types/lines failed to recapitulate transcriptome-wide responses ( Fig. 3R ). However, real covariances of unrelated cell types outperformed their corresponding mean-field models (shuffling the gene counts column-wise, which retains the original mean and distribution for each gene individually), suggesting the presence of certain universal, albeit weakly-contributing, fluctuation structures that inform outcomes ( Fig. 3R ) . These findings are consistent with a recent study outlining the presence of shared and cell-type-specific correlative gene modules between the human osteosarcoma (U2OS) and human embryonic kidney (HEK293T) cell lines 46 . Collectively, these results suggest the presence of cell-type specific strongly-contributing and universal weakly-contributing latent fluctuation structures which can be harnessed to predict transcriptome-wide responses. CIPHER accurately predicts true perturbations Having tested CIPHER in the forward direction on capturing genome-wide response to single and double perturbations, we next investigated the inverse problem: given a covariance matrix and a transcriptome-wide change in gene expression, what gene (or genes) most likely induced the change ( Fig. 4A )? To address this question, we solved for a full combinatorial gene perturbation u i (a potentially nonzero element u for all genes) as Σ −1 Δ X . We argued that ordering the resulting u i by magnitude should, in principle, provide an effective ranking of gene i and a point estimate of whether they were likely knocked-down or activated. High-throughput Perturb-seq data provides an excellent platform to solve this inverse problem since the ground-truth perturbations are known a priori . To incorporate the possibility of error in the estimates of the averages, covariances, or deviations from linearity, we reframed the question as a Bayesian linear regression problem ( Fig. 4B ). Concretely, we sample from the posterior distribution p ( u i |Δ X , Σ) assuming some prior knowledge on sparsity (i.e., only a few gene perturbations should have significant nonzero effects). This is accomplished by introducing a prior distribution over the perturbations p ( u ) which has support over zero (no effect) and a mode away from zero occupied by true perturbation effects (see Methods 4) . This operation enabled us to put error bars on estimates of perturbation effects, and assign each gene a posterior inclusion probability (PIP) that it has true nonzero effect ( Fig. 4B , S4A ). Together, the Bayesian approach allowed CIPHER to provide uncertainty-aware estimates of the perturbation effect sizes that the naive point estimates could not. Download figure Open in new tab Fig. 4: The inverse problem: predicting causal drivers of perturbation response in transcriptome-wide screens. A) CIPHER ranks genes responsible for a measured perturbation using knowledge of the control-cell gene expression fluctuations. Each point corresponds to a transcriptome from a single cell. B) Left: Posterior probability of a nonzero perturbation effect across genes (points). True perturbed genes A and B have high posterior probabilities. Right: effect size distributions for different perturbations. Distributions of effect size for genes A and B as well as a composite distribution over other unperturbed genes. C) Receiver-operator characteristic curves for each CRISPRi dataset. Each curve corresponds to a dataset. D) Receiver-operator characteristic curves for each CRISPRa dataset. E) AUROC comparison across metrics: p-value vs. PIP. (One-sided Wilcoxon signed-rank test p-value = 0.042) F) AUROC comparison across metrics: log fold change vs maximum posterior effect size distribution spread vs. PIP. Each line connects AUROC scores for the same dataset over different conditions. (One-sided Wilcoxon signed-rank test p-value = 0.042) G) AUROC comparison across covariance conditions: shuffled, ZINB, mean-field and the true covariance. Note that ZINB only has 4 datasets associated with it, as reparameterization for the other datasets led to pathological results (negative mean, see Methods 4 ). (Shuffled Σ vs ZINB (n=4): one-sided p-value = 0.31; ZINB vs Meanfield (shuffled X 0 ) (n=4): one-sided p-value = 0.063; Meanfield (shuffled X 0 ) vs Real Σ (n=10): one-sided p-value = 0.042) H) Receiver-operator characteristic curves across datasets containing double (2-gene) perturbations. To make the problem tractable, we performed Markov-Chain Monte Carlo on the top 200 prospective perturbations according to Σ −1 Δ X , forcibly including the true perturbation if needed (see Methods 4 ). We found that often when the true perturbation is originally ranked low by the naive point estimate, it rises through the ranks during inference to become one of the, if not the most, highly probable genes with the largest magnitude effect size of the top k (median rank increase: 28.2, see Fig. S4B ). For each dataset in Table 1 , we compared the ranking metrics (PIP,max(+ / − Σ i )) with receiver-operator characteristic (ROC) curves, finding that ranking the perturbations by the maximum spread max(+ / − Σ i ) proved most informative of the true perturbation ( Fig. 4C-G ). We found that regardless of the dataset, cell type, or perturbation type (CRISPRi or CRISPRa), CIPHER achieved high values on the receiver-operator characteristic curves for maximum spread and PIP ( Fig. 4C-E ). Next, we wondered how CIPHER results compare to conventional metrics used in differential analysis of control and perturbation conditions, such as log fold-change (logFC) and associated p-values. We found that CIPHER outperformed these conventional metrics in identifying the true perturbation, particularly for datasets with relatively lower mean AUROC values ( Fig. 4E-G ) . CIPHER exhibited high performance for double perturbations datasets as well ( Fig. 4H ). Unlike conventional metrics, CIPHER calculates how much each gene contributes to other genes’ expression changes, powering it to reveal the mechanistic underpinnings of cellular response, which we cover in the next section. Soft modes drive genome-wide responses to perturbations Given CIPHER’s success in capturing transcriptome-wide responses, we next examined the underlying drivers of its high performance. CIPHER inherently quantifies how each gene influences every other gene’s expression, in principle allowing us to interpret both gene- and “regulatory module”-level contributions to responses. We first investigated the dimensionality of the response. That is, to what extent the transcriptome-wide response is propagated along specific directions or “soft modes” of the covariance matrix. (A soft mode is a direction of change that the system can move relatively easily along and mathematically corresponds to a dominant principal component.) We decomposed the covariance matrix into its eigen-representation, where eigengenes constitute effective regulatory modules of multiple co-expressed genes, each capturing a different amount of variance. Should the variance captured by a given eigengene dominate over others, it is considered a soft mode. We confirmed the presence of soft modes in the datasets, as revealed by the high magnitude of eigenvalues of the covariance matrix ( Fig. S5A ). Such soft modes have recently been argued to be ubiquitous in transducing response to perturbation in complex biophysical systems, illustrating the utility of adopting a linear response framework like the one employed in our study. 47 For each perturbation in each dataset, we additionally calculated a participation ratio. The participation ratio—a measure of how many eigengenes actively contribute to the response—of a vector is a measure of how many nonzero elements it has: it evaluates to 1 when there is one nonzero and to the length of the vector when all elements have equal size. Specifically, we evaluated the participation ratio of the squared overlap (dot product) of each eigengene with the total response vector Δ X , yielding an estimate of the effective dimensionality of response. We reasoned that a low participation ratio could contribute to higher CIPHER performance. Indeed, the participation ratio was generally low relative to the dimension of the datasets (mean = 3.07 eigengenes over all perturbations) with notable variation between datasets ( Fig. 5A and S5B-D ). Interestingly, a majority of responses in Tian19_a ( Table 1 )—which exhibited the highest R 2 values ( Fig. 3E )—could essentially be described with a single dimension. Conversely, datasets with relatively high participation ratios correspondingly exhibited lower R 2 values (correlation R 2 = 0.75; slope= -6.5) ( Fig. 5D ). We observed consistent results when comparing the normalized contribution of each eigengene for each perturbation, again illustrating that Tian19_a responses were typically well described by their first principal component, unlike the other datasets ( Fig. 5B ). Next, we checked whether certain gene programs underlie dominant eigenmodes by calculating the dominant genes from the eigengene loadings for each principal component, and performed a gene ontology 48 search on them. In line with a recent study on covariances in mammalian cells 46 , the majority of processes are related to translation and general cellular homeostasis ( Fig. 5C ), a trend that holds true over all datasets combined as well as within datasets from a single cell type (Neuron, Fig. S5F,G ). This structure reflects an intrinsic coordination essential for maintaining cellular homeostasis—specifically, at and near equilibrium states, cellular processes perhaps prioritize general growth pathways; these pathways may serve as the primary conduits through which perturbations are transduced in cells. Download figure Open in new tab Fig. 5: The effective dimensions of transcriptome-wide response. A) Distributions of participation ratio (effective dimension) across all perturbations for three Perturb-seq datasets. B) Clustered heatmaps of the fraction of the response that falls along each principal component for the datasets in A). C) Gene-ontology significance heatmaps (clustered) from gene sets with high-loadings in each principal component. D) Mean participation ratios vs. mean R 2 values over the CRISPRi/a Perturb-seq datasets considered. Each point corresponds to a dataset and we show the best fit line describing the linear trend. E) The effective number of genes driving global change,calculated from the entropy of the fractional contributions of each gene’s change from every other gene. F) The effective number of genes impacting the expression change of the true perturbed gene across datasets, calculated in the same way as E). CIPHER provides estimates on total genes affected during perturbations We next sought to calculate the effective number of differentially expressed genes upon perturbation. CIPHER provides the total normalized contribution to the expression change of gene i due to correlations with gene j , and the analogous contribution of each gene to the total response. We calculated the entropy of these two probability vectors and exponentiated it to obtain an effective number of genes responsible for two kinds of estimates: 1) the expression change in the true perturbed gene, and 2) the total transcriptomic change upon perturbation (see Methods 5.3 ). Entropy-based estimates suggest that the transcriptomic response is distributed across a large number of genes ( Fig. 5E , S5D ), consistent with prior findings 16 that essential gene perturbations trigger widespread but diffuse expression changes. Many perturbed genes showed high self-variance contributions to their expression changes, implying that the perturbations directly affected the gene being targeted, as expected ( Fig. S5E ). Furthermore, while specific gene interactions could significantly affect perturbation responses, most essential gene perturbation effects appeared to emerge from many smaller correlations between the perturbed gene and its partners ( Fig. 5F ). Collectively, our results highlight despite the presence of soft modes, each soft mode represents a complex, yet quantifiable milieu of gene interactions. DISCUSSION Here, we developed CIPHER, a framework combining linear response theory, Bayesian statistics, and regulatory networks to model rich single-cell perturbation experiment outcomes. We benchmarked CIPHER on synthetic networks, demonstrating its capability to predict systemic responses to perturbations. Applied to experimental single-cell perturbation datasets ( Table 1 ), CIPHER reliably captured transcriptome-wide responses to gene perturbations—an effect that disappeared when gene-gene fluctuations in unperturbed cells were eliminated. Moreover, gene-gene fluctuations from one dataset performed similarly well in other datasets of the same cell type, highlighting the presence of stereotypic genome-wide fluctuation structures within cell types. CIPHER also identified the true perturbation with high-fidelity given only the covariance matrix and change in gene expression upon perturbation, while also revealing that transcriptome-wide responses tend to align with soft eigenmodes. Each soft mode captures a quantifiable ensemble of weakly- and strongly-coupled gene interactions that confer predictability to responses. In sum, our results demonstrate that theoretical approaches drawn from statistical physics like CIPHER offer powerful tools for fundamentally understanding the responses to a large class of perturbations, despite the underlying biological complexity. CIPHER’s success in predicting cellular responses reflects the broader power of the linear response theory, which has proven remarkably effective across a diverse array of complex systems. In physical sciences, linear response theory underpins foundational theories like Marcus electron transport theory and Kubo’s Fluctuation-Dissipation theorem in non-equilibrium statistical mechanics, explaining phenomena ranging from NMR spectroscopy and dielectric response to protein allostery 24 – 27 . CIPHER’s ability to accurately predict complex, transcriptome-wide responses using a simple, interpretable linear model suggests that biological systems may operate in regimes where linear response theory remains surprisingly effective despite underlying nonlinear dynamics. The presence of coordinated fluctuations and their important contribution to responses to genetic and environmental perturbations has been reported for diverse microbial and mammalian behaviors across scales including transcriptional, protein, and metabolic cell states 46 , 49 – 51 . In fact, the authors of one particular experimental study found steady-state single-cell gene correlations could predict the effect of p53 perturbation more accurately than even chromatin immunoprecipitation studies 46 . Yet, these observations came from the perturbed gene correlations, and had not been explicitly rationalized or quantified within a theoretical framework. As such, our study provides a strong theoretical footing toward a fundamental biological principle: complexity can emerge from or be effectively captured by simpler mechanistic rules encoded in the fluctuations of unperturbed cellular states. Such fundamental insights offer a path toward mechanistic rather than purely correlative or black-box understanding cellular responses to perturbations or lack thereof 52 . Given our results within and across cell-types, in principle, CIPHER could be harnessed to reveal the hierarchical relationships between cell types on underlying genome-wide fluctuation structures. Going further, we imagine CIPHER could be employed across species to probe the evolutionary history of gene-gene correlation. In its present form, CIPHER takes each dataset as input separately. Future studies could couple the CIPHER framework with existing analytical and deep neural network approaches 53 – 59 for joint embeddings of disparate datasets. Furthermore, incorporating a CIPHER-based consistency loss that leverages pre-computed gene-gene covariance matrices could enhance the generative modeling of perturbation effects. By encouraging the model to preserve these covariance structures, this approach may improve generalization to unseen cell types and perturbations. CIPHER demonstrated that transcriptome-wide perturbation effects can emerge through different modes: concentrated effects from a small subset of highly responsive genes ( u i ) versus distributed effects where many genes contribute modestly but cumulatively significant changes. These findings, combined with effect size calculations and comparison from TRADE 16 , point to the potential application of CIPHER in designing efficient attention or, even, tokenization mechanisms 60 , 61 to enhance the predictive power of existing deep learning models or designing newer, more expressive, cellular foundation models. While CIPHER effectively captures responses in many cases, it sometimes fails due to nonlinearities and factors beyond direct gene-gene fluctuations, as well as higher response complexity in some datasets. Incorporating the recently proposed single-cell clone tracing technologies 2 , 62 – 64 with pooled CRISPR screening could enhance the accuracy of CIPHER by accounting for clone-specific fluctuations. Although beyond the scope of the current work, CIPHER could be harnessed to predict unknown perturbations driving transitions in complex disease traits such as cancer metastasis and drug resistance to targeted therapy 65 , genetic disorders with unknown causal genes 66 , and more, significantly shortening the list of possible causal drivers. In such cases for example, one would compute the covariance of the naive cancer cells before drug as well as the response to drug and apply the inverse problem as we have formulated it here. This could constitute a standalone procedure without necessarily the need for integration of prior perturbation data, as demonstrated by other recent approaches combining trait and perturbation datasets 66 . Scaling up perturbation experiments with AI and statistical physics could enable efficient virtual screens to map genotype to phenotype beyond current single-trait CRISPR studies. AUTHOR CONTRIBUTIONS YG and BKS conceived and designed the study. BKS designed and performed a majority of the computational experiments and analyses with input from YG. LS and MEM collected and harmonized the datasets. HS, BKS, and YG prepared the figures and tables. HS, NK, BH, and EK contributed to the analyses. SV assisted BKS and YG in formulating the theoretical framework. BKS and YG wrote the manuscript with input from all authors. CONFLICT OF INTEREST The authors declare no competing interests. DATA AND CODE AVAILABILITY This paper analyzes existing, publicly available datasets 10/11 of which are available on the scPerturb database. All 11 datasets can be accessed through google drive (see github README for link). All code for the analyses in this manuscript has been deposited at: https://github.com/GoyalLab/CIPHER . SUPPLEMENTARY CAPTIONS Fig. S1: Additional numerical examples on the Hill function regulatory network . A) Relative error in linear response as a function of Hill saturation parameter K for three perturbation magnitudes. B) Relative error in linear response as a function of decay rate γ for three perturbation magnitudes. C) Relative error in linear response with multiple genes perturbed. For increasing Hill coefficient (left to right) we show relative error as a function of perturbation magnitude for 1-5 genes perturbed with equal magnitude.Relative error in linear response (single gene perturbed) as a function of Hill saturation parameter K for three perturbation magnitudes. Each point represents a single simulation run and lines guide the eye. Fig. S2: Predicting perturbation response from covariance structure in single gene perturbation A) - F) Distributions of R 2 values over perturbations for additional datasets (labeled by dataset) for the full covariance, the mean field model and the shuffled covariance. G) - J) Transcriptome-wide changes in gene expression over 4 genes from different datasets (lines are experimental data and points are the linear response prediction). K) - N) Predicted vs. measured change in expression over genes (points) for the genes and datasets in G)-J). Fig. S3: Additional forward problem plots A) A scatter plot of per-gene R 2 as a function of mean control counts for a particular perturbation. All genes measured for all perturbations shown as points. B) Distributions of R 2 over perturbations across CRISPRi/a datasets: real covariance compared to the mean-field model. C) R 2 distributions across all perturbations and datasets. Shuffled and mean-field histograms contain data from all perturbations and every possible perturbed gene. The true covariance R 2 histogram only includes genes that were truly perturbed. C) R 2 histogram from linear response applied to perturb-fish imaging data: real compared to mean-field and shuffled null models. E) R 2 histograms for individual datasets over double perturbations. Density of R 2 values computed with the true Σ compared to that computed with a correlation matrix from a different cell type. F) Mean R 2 over all perturbations from datasets with combinatorial perturbations (2 or more genes perturbed). Results from the full, multi-gene optimization, and the additive approximation are shown. G) Density of R 2 over all perturbations using the true covariance vs. a covariance matrix from a different cell type. Fig. S4: Bayesian inference of true double perturbation effect sizes A) Top two rows: examples of effect size distributions for the constituent perturbations comprising double perturbations. Third row: posterior inclusion probabilities for the double perturbations shown above (points are individual genes; true perturbed genes shown in red). Bottom row: posterior means with standard deviation error bars. B) Distribution of rank differences of true perturbations (between the point estimate and Bayesian inference). Larger differences indicate a better final rank than initially calculated without inference. Fig. S5: Eigengene analysis within and across datasets . A) Ranked eigenvalues of covariance matrices from all 10 Perturb-seq datasets. B) Distributions of participation ratio for the remaining datasets not shown in the main figure 5A . C) Percent that each PC falls along the direction of the response as a function of perturbation. D) Distribution of participation ratio for all datasets and perturbations (left) and across cell type specific neuron datasets. E) Left: Distribution of self rank over all perturbed genes. Right: the effective number of genes needed to reconstruct 99% of the response. F) Clustered gene ontology significance heatmaps as a function of principal components over all datasets. G) Clustered gene ontology significance heatmaps as a function of principal components over all Neuron datasets. Methods for CIPHER Here we detail the methods and procedures for CIPHER. We start by deriving the results of linear response theory and motivating its use in describing transcriptional changes. Next, we discuss simulations for the synthetic regulatory networks considered in the main text. From here, we consider the ‘forward’ and ‘inverse’ problems in turn, and conclude by describing the analysis surrounding our discussion of soft modes in propagating the response of gene expression to cellular perturbation. 1 The case for transcriptome-wide linear response theory We now describe the theoretical foundations of CIPHER. Linear response is a classical result in (quasi)equilibrium statistical mechanics [ 1 , 2 ] which states that the typical response of a physical system to some relatively small external force is encoded in the correlations of the unforced system. Denoting the state of the unper-turbed system as x = ( x 1 , x 2 , …), the force as u and the covariance as Σ so that Σ ij = cov ( x i , x j ), the main result of static linear response is or This result naturally emerges from many realistic physical contexts. Namely, Eq. 2 holds: 1) for any quasi-stationary system where the perturbation u is conjugate to the expression state variable x 2) in a number of simple transcriptional bursting models, and 3) in systems with soft modes that dominate the spectrum of the covariance. We discuss each of these scenarios below. 1.1 Quasi-stationary systems with X conjugate to u We consider a system in a steady state described by a probability distribution p 0 ( x ), where x ∈ ℝ N is the state vector. The system is perturbed by a small external field u ∈ ℝ N , such that the perturbed distribution becomes: This form implies that x is the sufficient statistic and is conjugate to u in the exponential family sense. The normalization factor Z ( u ) is the moment generating function: The expectation of x under the perturbed distribution is which can be written as the gradient of log Z ( u ), Now taking the derivative of ⟨ x ⟩ u with respect to u to obtain the response Evaluating this at u = 0 , we recover the covariance matrix of the unperturbed system, Therefore, to leading order in u , the response of the mean is This is the static linear response formula for systems where the perturbation enters the distribution as , and is valid under the assumption of quasi-equilibrium (i.e., fast microscopic relaxation compared to the timescale of u ). 1.2 Heuristic examples showing the perturbation is conjugate to gene expression In the following examples we write the mean equations of motion and assume unless otherwise stated that gaussian white noise determines fluctuations around these mean equations, i.e. that the noise amplitude . 1.2.1 Example 1: Constant-rate transcription with linear perturbation Consider a gene with linear degradation and production perturbed by a small additive input u At steady state ( dx/dt = 0), the solution is The steady-state distribution is Gaussian, Under perturbation u , the mean shifts to x * = ( β + u ) /γ , which yields Thus, showing that x is conjugate to u in the exponential family sense. 1.2.2 Example 2: Bursting transcription from a two-state promoter We consider a minimal model of transcriptional bursting with promoter switching: Let s ( t ) ∈ {0, 1} denote the promoter state and x ( t ) ∈ ℕ 0 the mRNA count. In the bursting regime , we assume: Slow promoter switching: k on , k off ≪ γ Fast transcription when ON: β ≫ γ Under these conditions, transcription occurs in bursts when the promoter switches ON briefly. The steady-state mRNA distribution is approximately exponential: with mean This can be interpreted as a bursty birth-death process: Bursts arrive as a Poisson process at rate k on Each burst yields B ∼ Geom( q ) mRNAs with mean burst size mRNAs degrade independently at rate γ The corresponding steady-state distribution is a negative binomial: where with mean In the limit B → ∞, k on → 0, with k on B = const, the negative binomial converges to the exponential distribution above. Now consider a small perturbation u that increases k on → k on + u . The ON-probability becomes and the perturbed mean expression is Expanding to first order in u , we see that with response coefficient The perturbed distribution then becomes so to first order So x is conjugate to the perturbation u and with 1.2.3 Example 3: Hill-type regulation with small perturbation to activator Consider a gene x transcriptionally regulated by an activator R via a Hill function, with a small basal transcription rate β . The deterministic dynamics are Here α is the maximal inducible transcription rate, K is the Hill constant, n is the Hill coefficient, and γ is the degradation rate. At steady state, the system fluctuates around the fixed point We now consider a small perturbation u that increases the activator level: R → R + u . For small u , we expand the Hill function as The derivative is so the perturbed steady-state becomes Assuming the fluctuations around x * ( u ) are Gaussian with variance σ 2 = D/γ , the steady-state distribution is Substituting in the linear expansion for x * ( u ), we write: where So, finally, we have with 1.3 Perturbation-conjugate systems and the variational principle Many biological systems are robust to small perturbations: when subjected to weak inputs u , their steady-state distribution shifts minimally from its unperturbed form. This robustness can be formalized by posing a variational problem: among all distributions p ( x ) that produce a given shift in the mean ⟨ x ⟩ = ⟨ x ⟩ 0 + Δ x , the actual response minimizes the Kullback–Leibler (KL) divergence from the original distribution p 0 ( x ). This is equivalent to saying the system follows the path of least surprise or minimal information cost under perturbation. We define the optimization problem as follows, where KL divergence is To solve this constrained optimization, we introduce a Lagrange multiplier u ∈ ℝ N for the mean constraint and a multiplier λ for normalization. The Lagrangian is then defined as Taking the functional derivative with respect to p ( x ) and setting it to zero, Solving this yields with normalization constant Thus, the optimal distribution under the mean constraint is an exponential tilt of the original distribution. This derivation shows that in any system that responds minimally to perturbations in the information-theoretic sense, the perturbed distribution takes exponential form and x is conjugate to u , guaranteeing linear response in the small- u limit. 1.4 Soft modes re-enforce linear response We consider stochastic dynamics governed by the Langevin equation: where D is a constant, possibly anisotropic (i.e., non-diagonal) diffusion matrix. Linearizing around a stable trajectory x * ( t ), fluctuations δ x ( t ) = x ( t ) − x * ( t ) evolve as: Even if J is asymmetric and D is non-diagonal, the steady-state covariance matrix Σ = ⟨ δ x δ x ⊤ ⟩ satisfies the Lyapunov equation: Now consider adding a small constant force u to the drift: F ( x ) → F ( x ) + u . In the linearized system, this induces a steady-state shift: If the dynamics are dominated by a soft mode, i.e., a single eigenvalue λ d ≪ Re(λ i ) for i > 1, then the inverse Jacobian is approximately: where p d and q d are the right and left eigenvectors of J , normalized so that . In this same limit, the covariance matrix Σ is also approximately rank-1, and from Eq. (18) in [ 3 ], its leading behavior is: where is the diffusion matrix expressed in the eigenbasis of J , and is the soft-mode component. This shows that the leading contribution to Σ is a projection onto the soft eigendirection p d , with a scalar prefactor determined solely by λ d and even when D is anisotropic. Consequently, the linear response becomes: where the perturbation has been rescaled by a constant factor. Absorbing this factor, we see that holds approximately in the soft mode limit regardless of the symmetry of J or structure of D . It can be shown through similar arguments that when J is symmetric and D diagonal, we arrive at the same answer. Thus, the existence of a dominant soft mode ensures that the system exhibits effective equilibrium-like linear response behavior, even under highly non-equilibrium stochastic dynamics. 1.4.1 Corrections from higher modes While the leading soft-mode approximation captures the dominant structure of Σ in systems near bifurcation or criticality, we can make this approximation more precise by quantifying corrections from higher modes. Assume J is diagonalizable, with eigen-decomposition J = PΛP −1 , where Λ = diag(λ 1 , …, λ N ), and P contains the right eigenvectors . The Lyapunov equation can be solved by transforming into the eigenbasis of J . Again, define , and Then the Lyapunov equation becomes diagonal in this basis: where denotes a complex conjugate. To return to physical coordinates, the covariance is reconstructed as: Now assume the system has a dominant soft mode λ d → 0, with all other eigenvalues satisfying Re(λ k ) ≫ λ d for k ≠ d . Then the leading contribution is: The remaining contribution comes from all ( k, l )≠( d, d ), forming the correction: This yields the full expansion: The magnitude of the correction term Σ (1) is suppressed when there is a spectral gap between λ d and all other eigenvalues; the perturbation and noise project weakly onto the stiff modes , and the off-diagonal elements are small, i.e., noise is nearly aligned with the eigenbasis of J . Thus, although the soft mode approximation is leading-order accurate, this expansion systematically quantifies when and how it fails, that is, when λ d / λ 2 ≪ 1 begins to break down. 1.5 Fitting single and multi-gene response We can judge how well correlations explain the expression difference Δ X by computing the coefficient of determination As R 2 approaches 1, the correlations perfectly transport the initial vector of gene expression ⟨ x ⟩ 0 to the final measured target ⟨ x ⟩ u . When R 2 > 0.5, for example, the initial distance between the control and perturbed gene expression has been halved. If the cellular perturbations causing phenotypic change are known, as they are in Perturb-seq experiments, wherein it is typical for a single gene to be upregulated or knocked down per cell by CRISPRa or CRISPRi respectively, we can directly test the extent to which linear response to that single-gene perturbation effects the expected expression change. That is, for any known single-gene perturbation i , we can solve for the optimal u i with which to force the system towards the terminal state ⟨ x ⟩ u . Setting all other elements of u equal to zero, the u i that solves the optimization problem is where Σ i is the ith column of Σ. We use such optimal single-gene solutions in the ‘forward’ problem below and can derive it explicitly as follows. Assuming we know which gene is perturbed we set all other elements of u to zero. Now, Σ u = Σ 1 i u i + … + Σ Ni u i ≡ Σ i u i where Σ i is the i th column of Σ. This means that the optimization in Eq. 59 reduces to which implies Solving this equation for the scalar u i gives the desired result. The structure of steady-state linear response dictates that we can solve for independent perturbation components u * by ‘decorrelating’ the observed expression change, so that We interpret large elements u i as those driving the phenotypic change relatively independently above and beyond that encoded by the typical correlation. Relatively small elements of u i would then correspond to genes whose expression changes largely follow expected correlations, so that the perturbation response of that gene can largely be explained by the behavior of other genes. 2 Application to synthetic gene regulatory networks 2.1 Random linear network Having explained the theoretical foundations of the CIPHER framework, we now test it on a number of synthetic gene regulatory networks. We first consider the case of a linear system with random regulatory structure. That is, around a stationary fixed point centered at the origin, so that x * = 0 and δ x = x , the Jacobian becomes constant and we parameterize it as − J ≡ A = 1 + ϵ K , where K is an off-diagonal matrix of interactions with elements drawn from a standard normal distribution. As Robert May showed in his seminal paper on the stability of random ecosystems, A is almost surely stable when the dimension of the system N -genes satisfies the following criteria from random matrix theory: . As the network of genetic interactions becomes more and more dense, the coupling strength approaches its critical value and the smallest eigenvalue of A → 0. The system bifurcates at this critical coupling, becoming globally unstable. We simulate this system using an Euler-Maryama integrator with timestep dt = 0.1. For the correlations and linear response estimates in FIG 2 , we use N = 300 genes and integrate from the stationary state δ X = 0 for a total time of T = 300 with a constant diffusion coefficient D = 0.5 for all genes. For this system, the critical coupling is , and we simulate one system at ϵ = 0.95 ϵ c , one system at criticality ϵ = ϵ c and one above the threshold ϵ = 1.05 ϵ c . Linear response estimates are averaged over 500 independent simulation runs. The emergence of gene teams past criticality is a direct result of the dynamics begin projected (according to center manifold theory [ 4 ]) onto the null-space of A , i.e. the direction with zero or negative eigenvalue [ 3 ]. This effective dimensionality reduction is not a result of gene teams, but the approach to and past criticality. This insight sheds light on the natural emergence of gene teams in the vicinity of a phenotypic transition; as the Jacobian becomes unstable as it would when passing through a short lived unstable fixed point, or barrier, in phase-space, mutual activation within large groups of genes and inhibition between groups becomes inevitable. This along with the fact that such gene teams have recently been posited to induce low-dimensional phenotypic structure in the epithelial to mesenchymal transition is why we study an explicit teams system below. 2.2 Nonlinear Hill function network We now move on to test CIPHER on a prototypical nonlinear dynamical gene expression system comprised of a network of hill functions [ 5 ]. Specifically, the deterministic force on the i th gene x i is . We apply an additional constant external force to the first gene with magnitude u and compare the true average response (from solving the stationary fixed point equations F i (⟨ x ⟩ u ) + u = 0 and F i (⟨ x ⟩ 0 ) = 0 to that predicted by linear response (Eq. ?? ), where again J is evaluated at the unperturbed fixed point ⟨ x ⟩ 0 . For this system, we calculate the original fixed point (with u = 0) and the Jacobian evaluated at that fixed point, and relate it to the new driven fixed point after perturbation. We perform these calculations using deterministic gradient flow (towards the fixed point) for a maximum of 2000 iterations with time-step dt = 0.01, stopping if the maximum change in x over an iteration is less than a fixed tolerance of 10 −12 . We consider a relatively small system N = 10 to sweep over many parameters and fix K = 10, γ = 1, G ij ∼ N(1, 0.1), starting from an initial guess where all genes are at abundance𝒩⟨ G ij ⟩ /γ = 10. When sweeping over the other parameters, we fix the Hill coefficient to n = 4. For this system, can explain the sudden jumps in the linear response error as a function of the system parameters by examining its stability. We assume a homogeneous fixed point at zero and follow what happens when we perturb around it by an amount δx , so that, averaging over the noise and assuming fixed G ij = G , The system goes from stable to perturbation to unstable when , so that Solving for G * , γ * , or K * yield the critical parameter values past which the system has finite average gene expression and the linear response error estimates jump. 2.3 Teams network The last synthetic system we explore here is a stochastic hybrid model of team dynamics, where genes on the same team mutually activate each other and there is inhibition between genes on different teams. Each gene has a promoter that transitions from off s = 0 to on s = 1, and vice versa, with rates k on and k off . These rates take in input from all the other genes in the system which is modulated by a hill function. That is where B on / off is a matrix of activating/inhibitory interactions, respectively. The force on each gene is then F i ( x ) = βs i − γx i where β and γ are transcription and decay rates. This nonlinear model features collective transitions wherein all the genes on a team will cooperatively transition from off to on while the genes on the other team transition from on to off. Even in such a complex simulated system, gene expression changes upon knock-out of a single gene (setting β = 0 for some gene i ) are well captured by linear response. For our simulations, we fix , , K = 25 and n = 2. We set β = 50 and γ = 4 and simulate N = 20 genes (10 on each team) for a total time of T = 2 × 10 6 with timestep dt = 0.1. We set the sparsity to S = 0.7 by randomly setting elements of B on/off to zero with probability S . At the end of the long simulation, we calculate the gene-gene covariance and the optimal single-gene perturbation according to Eq. 60 after setting the production rate of the first gene to zero, knocking it out. 3 The forward problem: predicting responses from Perturb-seq We compute the optimal single gene perturbation from Eq. 60 for each single and double-gene perturbation in each dataset, after cleaning it (see Data Processing below). We then use to calculate an R 2 value for 1) the real Σ matrix, 2) the ‘mean-field’ matrix, shuffled only over cells, not genes, and 3) the fully-’shuffled’ matrix. First we regularize the denominator in both the R 2 expression and the expression for by adding a small number 10 −8 to avoid division by zero. In order to suppress spurious small and noisy correlations from artificially driving down R 2 , in the sum what uninteresting case that the expression of a given gene does not change (Δ x i = 0), we only compute R 2 over genes that have nonzero response to the perturbation. When comparing distributions and means of R 2 across datasets, we calculate p-values using Kolmogorov-Smirnoff tests and Wilcoxon signed rank t-tests (one-sided), as implemented in sciPy, respectively. 4 The inverse problem: predicting perturbations from correlations and response Next, we investigated how informative the linear response estimates of the full (polygenic) perturbations in Eq. 63 are at inferring the true perturbation. For each of the 10 datasets considered, we considered ranking the genes by 1) their u * magnitude, 2) log-fold change, and 3) minus the log-10 p-value between control and perturbation conditions. However, since Eq. 63 only provides a point estimate of the true perturbations driving the response, we wanted to venture further to quantify the distribution of effect sizes for each gene p ( u i ) as well as the posterior inclusion probability that the effect size is truly nonzero. To do this we re-framed the linear response problem as a Bayesian linear regression [ 6 ], where Δ X is the data, Σ u is the model with Σ fixed and ξ is a random Normally distributed error term, with mean zero and variance σ 2 . Here ξ can be thought of as modeling co-variates unaccounted for in linear response. Using Bayes’ rule, it can be shown that the posterior probability of observing u given the data is, up to a constant, where p ( u ) is a prior on the perturbation that we choose to be a so called Horseshoe prior [ 7 ] since it reflects the belief that the true vector of perturbations should be sparse. Since MCMC is extremely computationally intensive for large optimizations over many possible perturbations, we optimize log-likelihood ln p ( u |Δ X ) over a subset of the perturbations (and corresponding rows of Δ X and block of Σ) chosen as the top-k perturbations ranked by absolute u -value according to the point estimate. We then performed Markov-Chain Monte Carlo (MCMC). As a test, we force the inclusion of the true perturbation (even if it lies outside the top-k) to see if the Bayesian estimate places it as highly probable to have non-zero effect with posterior inclusion probability (PIP) close to 1. Indeed, true perturbations that were originally ranked low can rise through the ranks during inference to become one of if not the most highly probable gene with the largest magnitude effect size of the top-k. Even when perturbations have low PIPs, their effect size distributions often have heavy tails and variances much larger than the other genes considered in the optimization, which is one reason why our chosen metric (the maximum mean-shifted spread of p ( u i )) has AUROC scores close to 1. To maximize the posterior distribution of perturbation effects we employ the python MCMC package pyMC. We use the no u-turn sampler [ 8 ], taking 1000 MC steps after tuning for 1000 steps and accumulating samples from 8 such runs. We parameterize the sparse Horseshoe prior for each gene i as where λ i ∼ Half − Cauchy(0, 1) and ln τ ∼ 𝒩(−4, 1). We learn the noise parameter as ln σ ∼ 𝒩(−2, 2). We use a conservative target acceptance rate (0.95). For the ZINB control AUROC shown in Fig. 4G , we first fit each gene’s control expression counts to a zero-inflated negative binomial (ZINB) distribution, then randomly drew a new variance for each gene, resampled from the resulting distribution and calculated the resulting covariance. We only performed this procedure for 4 of the 10 Perturb-seq datasets considered as the others resulted in pathological individual gene distributions as they were originally too sparse to support an altered variance without pushing the mean below zero. 5 Eigengene analysis Let Σ = V Λ V ⊤ be the eigendecomposition [ 9 ] of the covariance matrix with V = [ V 1 , …, v G ]. Then, the projection coefficients which quantify the extent to which the response to perturbation j lies in the samedirection as the eigenvectors are and the normalized contribution of PC i to perturbation j can be estimated as: 5.1 Participation Ratio Another measure of effective dimensionality of response to perturbation j is the participation ratio: 5.2 GO enrichment analysis on principle components We perform gene ontology (GO) enrichment [ 10 ] on the dominant axes of transcriptional variation across multiple Perturb-seq datasets. For each dataset, the we begin by computing the top 30 eigenvectors of the control covariance ordered by explained variance, which represent the primary directions of variation in control gene expression. These eigenvectors are interpreted by identifying the top 200 genes with the highest absolute loadings for each. These top genes are then submitted to g:Profiler for functional enrichment analysis against GO Biological Process (BP), Molecular Function (MF), and Reactome pathway (REAC) databases. For each eigengene, the top 5 significantly enriched terms (based on p-value) are retained. The results are aggregated into a matrix where rows represent GO or Reactome terms and columns represent eigengenes, with cell values reflecting enrichment strength (log p-value). This matrix provides a compact, interpretable summary of which biological processes are most associated with the dominant axes of gene expression variance in each dataset. All outputs—including the eigenvectors, gene rankings, enrichment results, and GO enrichment matrices—are saved for downstream analysis and visualization. 5.3 Gene Contribution Fractions Define the absolute propagated contribution to gene i : which can in general be larger than the true expression change because of possible cancellation of positive and negative correlations. The normalized global contribution vector is then and the true perturbed target gene-specific contribution vector is 5.4 Effective Number of Contributing Genes (Entropy) For a normalized distribution f ∈ ℝ G , we define the entropy-based [ 11 ] effective number of genes as: 5.5 Self Rank We rank the perturbed gene’s index g in its own contribution vector f (target) , with genes sorted by descending f j . 6 Data pre-processing We apply a simple filtering procedure to the raw count matrices from the 10 studies considered here (accessed via scPerturb [ 12 ]). We begin by loading single-cell perturbation datasets in .h5ad format using Scanpy [ 13 ], computing the sparsity of the full expression matrix before filtering. We then standardize perturbation names by stripping replicate or guide suffixes (e.g., g1, g2) to extract a base perturbation label. To retain biologically meaningful genes and limit the amount of noise in the data due to extremely low expression counts, we apply an average expression filter, keeping only genes with a mean expression above a specified threshold (1 count) and ensuring that all measured perturbed genes are preserved regardless of their expression level. For the perturbations, we retain only those with at least a specified minimum number of cells (100). After filtering, we extract control cells and the selected perturbation cells. ACKNOWLEDGEMENTS We thank members of the Goyal labs, especially Emanuelle Grody, for helpful discussions and comments on the manuscript. We thank the Feinberg Information Technology, Northwestern University Information Technology, and Quest High Performance Computing Cluster at Northwestern University Feinberg School of Medicine for their services. We acknowledge support from the National Institute for Theory and Mathematics in Biology (NITMB) through the National Science Foundation (DMS-2235451) and the Simons Foundation (MPTMPS-00005320). BKS acknowledges support from NUCATS T32, NITMB, and startup funds to YG from Northwestern University. LS and HS were supported by funding to YG. MEM acknowledges support from Ryan Fellowship and the International Institute for Nanotechnology, NITMB, and Carcinogenesis NIH T32. NK acknowledges support from CSB postdoc fellowship and the HFSP long-term fellowship. B.H. was funded from IVADO and the Canada First Research Excellence Fund. S.V. was supported by the National Institute of General Medical Sciences of the NIH under Award No. R35GM147400. YG acknowledges support from the Pew-Stewart Scholars Program and Burroughs Wellcome Fund Career Awards at the Scientific Interface. YG is a CZ Biohub Investigator. Funder Information Declared National Science Foundation, https://ror.org/021nxhr62 , DMS-2235451 Simons Foundation, https://ror.org/01cmst727 , MPTMPS-00005320 NUCATS T32, Northwestern University Ryan Fellowship, Northwestern University CSB postdoc fellowship, Northwestern University IVADO Canada First Research Excellence Fund National Institute of General Medical Sciences of the NIH , R35GM147400 Pew-Stewart Scholars Program Burroughs Wellcome Fund Career Awards at the Scientific Interface Chan-Zuckerberg Biohub Chicago Footnotes ↵ † Lead contact https://github.com/GoyalLab/CIPHER References 1. ↵ Klein , A. M. et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells . Cell 161 , 1187 – 1201 ( 2015 ). OpenUrl CrossRef PubMed 2. ↵ Pillai , M. , Hojel , E. , Jolly , M. K. & Goyal , Y. Unraveling non-genetic heterogeneity in cancer with dynamical models and computational tools . Nature Computational Science 3 , 301 – 313 ( 2023 ). OpenUrl CrossRef PubMed 3. ↵ Macosko , E. Z. et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets . Cell 161 , 1202 – 1214 ( 2015 ). OpenUrl CrossRef PubMed 4. ↵ Adamson , B. et al. A Multiplexed Single-Cell CRISPR Screening Platform Enables Systematic Dissection of the Unfolded Protein Response . Cell 167 , 1867 – 1882 .e21 ( 2016 ). OpenUrl CrossRef PubMed 5. Datlinger , P. et al. Pooled CRISPR screening with single-cell transcriptome readout . Nat. Methods 14 , 297 – 301 ( 2017 ). OpenUrl CrossRef PubMed 6. ↵ Dixit , A. et al. Perturb-seq: Dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens . Cell 167 , 1853 – 1866 .e17 ( 2016 ). OpenUrl CrossRef PubMed 7. Jaitin , D. A. et al. Dissecting Immune Circuits by Linking CRISPR-Pooled Screens with Single-Cell RNA-Seq . Cell 167 , 1883 – 1896 .e15 ( 2016 ). OpenUrl CrossRef PubMed 8. Bock , C. et al. High-content CRISPR screening . Nat. Rev. Methods Primers 2 , 1 – 23 ( 2022 ). OpenUrl CrossRef 9. Baysoy , A. et al. Spatially resolved in vivo CRISPR screen sequencing via Perturb-DBiT . bioRxivorg 2024.11.18.624106 ( 2024 ) doi: 10.1101/2024.11.18.624106 . OpenUrl Abstract / FREE Full Text 10. ↵ Binan , L. et al. Simultaneous CRISPR screening and spatial transcriptomics reveal intracellular, intercellular, and functional transcriptional circuits . Cell 188 , 2141 – 2158 .e18 ( 2025 ). OpenUrl CrossRef PubMed 11. ↵ Binan , L. Investigating spatial gene circuits and gene-phenotype mechanisms with Perturb-FISH . Nat. Rev. Genet . ( 2025 ) doi: 10.1038/s41576-025-00857-8 . OpenUrl CrossRef 12. ↵ Rood , J. E. , Hupalowska , A. & Regev , A. Toward a foundation model of causal cell and tissue biology with a Perturbation Cell and Tissue Atlas . Cell 187 , 4520 – 4545 ( 2024 ). OpenUrl CrossRef PubMed 13. ↵ Bunne , C. et al. How to build the virtual cell with artificial intelligence: Priorities and opportunities . Cell 187 , 7045 – 7063 ( 2024 ). OpenUrl CrossRef PubMed 14. ↵ Jiang , L. et al. Systematic reconstruction of molecular pathway signatures using scalable single-cell perturbation screens . Nat. Cell Biol . 27 , 505 – 517 ( 2025 ). OpenUrl CrossRef PubMed 15. ↵ Song , B. et al. Decoding heterogeneous single-cell perturbation responses . Nat. Cell Biol . 27 , 493 – 504 ( 2025 ). OpenUrl CrossRef PubMed 16. ↵ Nadig , A. et al. Transcriptome-wide analysis of differential expression in perturbation atlases . Nat. Genet . 1 – 10 ( 2025 ). 17. ↵ Yang , L. et al. scMAGeCK links genotypes with multiple phenotypes in single-cell CRISPR screens . Genome Biol . 21 , 19 ( 2020 ). OpenUrl CrossRef PubMed 18. ↵ Cleary , B. , Cong , L. , Cheung , A. , Lander , E. S. & Regev , A. Efficient generation of transcriptomic profiles by random composite measurements . Cell 171 , 1424 – 1436 .e18 ( 2017 ). OpenUrl CrossRef PubMed 19. ↵ Yao , D. et al. Scalable genetic screening for regulatory circuits using compressed Perturb-seq . Nat. Biotechnol . 42 , 1282 – 1295 ( 2024 ). OpenUrl CrossRef PubMed 20. ↵ Lotfollahi , M. , Wolf , F. A. & Theis , F. J. scGen predicts single-cell perturbation responses . Nat. Methods 16 , 715 – 721 ( 2019 ). OpenUrl CrossRef PubMed 21. ↵ Roohani , Y. , Huang , K. & Leskovec , J. Predicting transcriptional outcomes of novel multigene perturbations with GEARS . Nat. Biotechnol . 42 , 927 – 935 ( 2024 ). OpenUrl CrossRef PubMed 22. ↵ Cui , H. et al. scGPT: toward building a foundation model for single-cell multi-omics using generative AI . Nat. Methods 21 , 1470 – 1480 ( 2024 ). OpenUrl CrossRef PubMed 23. ↵ Bereska , L. & Gavves , E. Mechanistic Interpretability for AI Safety — A Review . TMLR ( 2024 ). 24. ↵ Kubo , R. The fluctuation-dissipation theorem . Rep. Prog. Phys . 29 , 255 – 284 ( 1966 ). OpenUrl CrossRef Web of Science 25. Marcus , R. A. Electron transfer reactions in chemistry. Theory and experiment . Rev. Mod. Phys . 65 , 599 – 610 ( 1993 ). OpenUrl CrossRef Web of Science 26. Warshel , A. Dynamics of reactions in polar solvents. Semiclassical trajectory studies of electron-transfer and proton-transfer reactions . J. Phys. Chem . 86 , 2218 – 2224 ( 1982 ). OpenUrl CrossRef Web of Science 27. ↵ Changeux , J.-P. & Edelstein , S. J. Allosteric mechanisms of signal transduction . Science 308 , 1424 – 1428 ( 2005 ). OpenUrl Abstract / FREE Full Text 28. ↵ Marcus , R. A. On the theory of oxidation-reduction reactions involving electron transfer . I. J. Chem. Phys . 24 , 966 – 978 ( 1956 ). OpenUrl CrossRef 29. Onsager , L. Reciprocal relations in irreversible processes. I . Phys. Rev . 37 , 405 – 426 ( 1931 ). OpenUrl CrossRef Web of Science 30. ↵ Onsager , L. Reciprocal relations in irreversible processes. II . Phys. Rev . 38 , 2265 – 2279 ( 1931 ). OpenUrl CrossRef Web of Science 31. ↵ Replogle , J. M. et al. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq . Cell 185 , 2559 – 2575 .e28 ( 2022 ). OpenUrl CrossRef PubMed 32. Norman , T. M. et al. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes . Science 365 , 786 – 793 ( 2019 ). OpenUrl Abstract / FREE Full Text 33. Frangieh , C. J. et al. Multimodal pooled Perturb-CITE-seq screens in patient models define mechanisms of cancer immune evasion . Nat. Genet . 53 , 332 – 341 ( 2021 ). OpenUrl CrossRef PubMed 34. Tian , R. et al. CRISPR interference-based platform for multimodal genetic screens in human iPSC-derived neurons . Neuron 104 , 239 – 255 .e12 ( 2019 ). OpenUrl CrossRef PubMed 35. ↵ Tian , R. et al. Genome-wide CRISPRi/a screens in human neurons link lysosomal failure to ferroptosis . Nat. Neurosci . 24 , 1020 – 1034 ( 2021 ). OpenUrl CrossRef PubMed 36. ↵ Kim , D. et al. Gene regulatory network reconstruction: harnessing the power of single-cell multi-omic data . NPJ Syst. Biol. Appl . 9 , 51 ( 2023 ). OpenUrl CrossRef PubMed 37. Sanguinetti , G. & Huynh-Thu , V. A. Huynh-Thu , V. A. & Sanguinetti , G. Gene Regulatory Network Inference: An Introductory Survey. in Gene Regulatory Networks: Methods and Protocols (eds. Sanguinetti , G. & Huynh-Thu , V. A. ) 1 – 23 ( Springer New York , New York, NY , 2019 ). 38. ↵ Margolin , A. A. et al. ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context . BMC Bioinformatics 7 Suppl 1, S7 ( 2006 ). 39. ↵ Freedman , S. L. , Xu , B. , Goyal , S. & Mani , M. A dynamical systems treatment of transcriptomic trajectories in hematopoiesis . Development 150 , ( 2023 ). 40. ↵ Phillips , R. Napoleon is in equilibrium . Annu. Rev. Condens. Matter Phys . 6 , 85 – 111 ( 2015 ). OpenUrl CrossRef PubMed 41. ↵ Bintu , L. et al. Transcriptional regulation by the numbers: models . Curr. Opin. Genet. Dev . 15 , 116 – 124 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 42. ↵ May , R. M. Will a large complex system be stable? Nature 238 , 413 – 414 ( 1972 ). OpenUrl CrossRef PubMed Web of Science 43. ↵ Hari , K. , Ullanat , V. , Balasubramanian , A. , Gopalan , A. & Jolly , M. K. Landscape of epithelial-mesenchymal plasticity as an emergent property of coordinated teams in regulatory networks . Elife 11 , ( 2022 ). 44. ↵ Chauhan , L. , Ram , U. , Hari , K. & Jolly , M. K. Topological signatures in regulatory network enable phenotypic heterogeneity in small cell lung cancer . Elife 10 , ( 2021 ). 45. ↵ Waddington , C. H. The Strategy of the Genes . ( Routledge , 2014 ). 46. ↵ Chapman , A. R. et al. Correlated gene modules uncovered by high-precision single-cell transcriptomics . Proc. Natl. Acad. Sci. U. S. A . 119 , e2206938119 ( 2022 ). OpenUrl CrossRef PubMed 47. ↵ Russo , C. J. , Husain , K. & Murugan , A. Soft modes as a predictive framework for low-dimensional biological systems across scales . Annu. Rev. Biophys . 54 , 401 – 426 ( 2025 ). OpenUrl CrossRef PubMed 48. ↵ Ashburner , M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium . Nat. Genet . 25 , 25 – 29 ( 2000 ). OpenUrl CrossRef PubMed Web of Science 49. ↵ Yamagishi , J. F. & Hatakeyama , T. S. Linear response theory of evolved metabolic systems . Phys. Rev. Lett . 131 , 028401 ( 2023 ). OpenUrl CrossRef PubMed 50. Kaneko , K. , Furusawa , C. & Yomo , T. Universal Relationship in Gene-Expression Changes for Cells in Steady-Growth State . Phys. Rev. X 5 , 011014 ( 2015 ). OpenUrl 51. ↵ Schuh , L. et al. Gene Networks with Transcriptional Bursting Recapitulate Rare Transient Coordinated High Expression States in Cancer . Cell Syst 10 , 363 – 378 .e12 ( 2020 ). OpenUrl CrossRef PubMed 52. ↵ Mellis , I. A. , Melzer , M. E. , Bodkin , N. & Goyal , Y. Prevalence of and gene regulatory constraints on transcriptional adaptation in single cells . Genome Biol . 25 , 217 ( 2024 ). OpenUrl CrossRef PubMed 53. ↵ Lopez , R. , Regier , J. , Cole , M. B. , Jordan , M. I. & Yosef , N. Deep generative modeling for single-cell transcriptomics . Nat. Methods 15 , 1053 – 1058 ( 2018 ). OpenUrl CrossRef PubMed 54. Xu , C. et al. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models . Mol. Syst. Biol . 17 , e9620 ( 2021 ). OpenUrl CrossRef PubMed 55. Lotfollahi , M. et al. Predicting cellular responses to complex perturbations in high-throughput screens . Mol. Syst. Biol . 19 , e11517 ( 2023 ). OpenUrl CrossRef PubMed 56. Korsunsky , I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony . Nat. Methods 16 , 1289 – 1296 ( 2019 ). OpenUrl CrossRef PubMed 57. Singhal , V. et al. BANKSY unifies cell typing and tissue domain segmentation for scalable spatial omics data analysis . Nat. Genet . 56 , 431 – 441 ( 2024 ). OpenUrl CrossRef PubMed 58. Dumitrascu , B. , Villar , S. , Mixon , D. G. & Engelhardt , B. E. Optimal marker gene selection for cell type discrimination in single cell analyses . Nat. Commun . 12 , 1186 ( 2021 ). OpenUrl CrossRef PubMed 59. ↵ Klein , D. et al. Mapping cells through time and space with moscot . Nature 638 , 1065 – 1075 ( 2025 ). OpenUrl CrossRef PubMed 60. ↵ Dotan , E. , Jaschek , G. , Pupko , T. & Belinkov , Y. Effect of tokenization on transformers for biological sequences . Bioinformatics 40 , btae196 ( 2024 ). OpenUrl CrossRef PubMed 61. ↵ Medvedev , A. et al. BioToken and BioFM – biologically-informed tokenization enables accurate and efficient genomic foundation models . bioRxiv 2025.03.27.645711 ( 2025 ) doi: 10.1101/2025.03.27.645711 . OpenUrl Abstract / FREE Full Text 62. ↵ Goyal , Y. et al. Diverse clonal fates emerge upon drug treatment of homogeneous cancer cells . Nature 620 , 651 – 659 ( 2023 ). OpenUrl CrossRef PubMed 63. Zhang , Z. et al. Synthetic DNA barcodes identify singlets in scRNA-seq datasets and evaluate doublet algorithms . Cell Genom 100592 ( 2024 ). 64. ↵ Sankaran , V. G. , Weissman , J. S. & Zon , L. I. Cellular barcoding to decipher clonal dynamics in disease . Science 378 , eabm5874 ( 2022 ). OpenUrl CrossRef PubMed 65. ↵ Sun , H. , Kumari , N. , Melzer , M. E. , Backman , V. & Goyal , Y. Navigating nongenetic plasticity in cancer drug resistance . Annu. Rev. Cancer Biol . 9 , 245 – 266 ( 2025 ). OpenUrl CrossRef 66. ↵ Kuznets-Speck , B. , Ogonor , B. K. , Wytock , T. P. & Motter , A. E. Generative prediction of causal gene sets responsible for complex traits . Proc. Natl. Acad. Sci. U. S. A . 122 , e2415071122 ( 2025 ). OpenUrl CrossRef PubMed References [1]. ↵ Robert Zwanzig . Nonequilibrium statistical mechanics . Oxford university press , 2001 . [2]. ↵ Nicolaas Godfried Van Kampen . Stochastic processes in physics and chemistry , volume 1 . Elsevier , 1992 . [3]. ↵ Simon L Freedman , Bingxian Xu , Sidhartha Goyal , and Madhav Mani . A dynamical systems treatment of transcriptomic trajectories in hematopoiesis . Development , 150 ( 11 ): dev201280 , 2023 . OpenUrl CrossRef PubMed [4]. ↵ Jack Carr . Applications of centre manifold theory , volume 35 . Springer Science & Business Media , 2012 . [5]. ↵ Chandrakala Meena , Chittaranjan Hens , Suman Acharyya , Simcha Haber , Stefano Boccaletti , and Baruch Barzel . Emergent stability in complex network dynamics . Nature Physics , 19 ( 7 ): 1033 – 1042 , 2023 . OpenUrl CrossRef [6]. ↵ Andrew Gelman , John B Carlin , Hal S Stern , and Donald B Rubin . Bayesian data analysis . Chapman and Hall/CRC , 1995 . [7]. ↵ Carlos M Carvalho , Nicholas G Polson , and James G Scott . Handling sparsity via the horseshoe . In Artificial intelligence and statistics, pages 73 – 80 . PMLR , 2009 . [8]. ↵ Matthew D Hoffman , Andrew Gelman , et al. The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo . J. Mach. Learn. Res ., 15 ( 1 ): 1593 – 1623 , 2014 . OpenUrl CrossRef [9]. ↵ Orly Alter , Patrick O Brown , and David Botstein . Singular value decomposition for genome-wide expression data processing and modeling . Proceedings of the National Academy of Sciences , 97 ( 18 ): 10101 – 10106 , 2000 . OpenUrl Abstract / FREE Full Text [10]. ↵ Gene Ontology Consortium . The gene ontology (go) database and informatics resource . Nucleic acids research , 32 (suppl_1): D258 – D261 , 2004 . OpenUrl CrossRef PubMed Web of Science [11]. ↵ Thomas M Cover and Joy A Thomas . Information theory and statistics . Elements of information theory , 1 ( 1 ): 279 – 335 , 1991 . OpenUrl [12]. ↵ Stefan Peidli , Tessa D Green , Ciyue Shen , Torsten Gross , Joseph Min , Samuele Garda , Bo Yuan , Linus J Schumacher , Jake P Taylor-King , Debora S Marks , et al. scperturb: harmonized single-cell perturbation data . Nature Methods , 21 ( 3 ): 531 – 540 , 2024 . OpenUrl CrossRef PubMed [13]. ↵ F Alexander Wolf , Philipp Angerer , and Fabian J Theis . Scanpy: large-scale single-cell gene expression data analysis . Genome biology , 19 : 1 – 5 , 2018 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted July 01, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Fluctuation structure predicts genome-wide perturbation outcomes Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Fluctuation structure predicts genome-wide perturbation outcomes Benjamin Kuznets-Speck , Leon Schwartz , Hanxiao Sun , Madeline E. Melzer , Nitu Kumari , Benjamin Haley , Ekta Prashnani , Suriyanarayanan Vaikuntanathan , Yogesh Goyal bioRxiv 2025.06.27.661814; doi: https://doi.org/10.1101/2025.06.27.661814 Share This Article: Copy Citation Tools Fluctuation structure predicts genome-wide perturbation outcomes Benjamin Kuznets-Speck , Leon Schwartz , Hanxiao Sun , Madeline E. Melzer , Nitu Kumari , Benjamin Haley , Ekta Prashnani , Suriyanarayanan Vaikuntanathan , Yogesh Goyal bioRxiv 2025.06.27.661814; doi: https://doi.org/10.1101/2025.06.27.661814 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41970) Biophysics (21463) Cancer Biology (18605) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19911) Epidemiology (2067) Evolutionary Biology (24329) Genetics (15615) Genomics (22514) Immunology (17743) Microbiology (40424) Molecular Biology (17194) Neuroscience (88650) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4827) Physiology (7648) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)