Multiscale topology classifies cells in subcellular spatial transcriptomics.

Spatial transcriptomics measures in situ gene expression at millions of locations within a tissue1, hitherto with some trade-off between transcriptome depth, spatial resolution and sample size2. Although integration of image-based segmentation has enabled impactful work in this context, it is limited by imaging quality and tissue heterogeneity. By contrast, recent array-based technologies offer the ability to measure the entire transcriptome at subcellular resolution across large samples3-6. Presently, there exist no approaches for cell type identification that directly leverage this information to annotate individual cells. Here we propose a multiscale approach to automatically classify cell types at this subcellular level, using both transcriptomic information and spatial context. We showcase this on both targeted and whole-transcriptome spatial platforms, improving cell classification and morphology for human kidney tissue and pinpointing individual sparsely distributed renal mouse immune cells without reliance on image data. By integrating these predictions into a topological pipeline based on multiparameter persistent homology7-9, we identify cell spatial relationships characteristic of a mouse model of lupus nephritis, which we validate experimentally by immunofluorescence. The proposed framework readily generalizes to new platforms, providing a comprehensive pipeline bridging different levels of biological organization from genes through to tissues.

Multiparameter persistent homology landscapes identify immune cell spatial patterns in tumors.

Highly resolved spatial data of complex systems encode rich and nonlinear information. Quantification of heterogeneous and noisy data-often with outliers, artifacts, and mislabeled points-such as those from tissues, remains a challenge. The mathematical field that extracts information from the shape of data, topological data analysis (TDA), has expanded its capability for analyzing real-world datasets in recent years by extending theory, statistics, and computation. An extension to the standard theory to handle heterogeneous data is multiparameter persistent homology (MPH). Here we provide an application of MPH landscapes, a statistical tool with theoretical underpinnings. MPH landscapes, computed for (noisy) data from agent-based model simulations of immune cells infiltrating into a spheroid, are shown to surpass existing spatial statistics and one-parameter persistent homology. We then apply MPH landscapes to study immune cell location in digital histology images from head and neck cancer. We quantify intratumoral immune cells and find that infiltrating regulatory T cells have more prominent voids in their spatial patterns than macrophages. Finally, we consider how TDA can integrate and interrogate data of different types and scales, e.g., immune cell locations and regions with differing levels of oxygenation. This work highlights the power of MPH landscapes for quantifying, characterizing, and comparing features within the tumor microenvironment in synthetic and real datasets.

Topological approximate Bayesian computation for parameter inference of an angiogenesis model.

MOTIVATION: Inferring the parameters of models describing biological systems is an important problem in the reverse engineering of the mechanisms underlying these systems. Much work has focused on parameter inference of stochastic and ordinary differential equation models using Approximate Bayesian Computation (ABC). While there is some recent work on inference in spatial models, this remains an open problem. Simultaneously, advances in topological data analysis (TDA), a field of computational mathematics, have enabled spatial patterns in data to be characterized. RESULTS: Here, we focus on recent work using TDA to study different regimes of parameter space for a well-studied model of angiogenesis. We propose a method for combining TDA with ABC to infer parameters in the Anderson-Chaplain model of angiogenesis. We demonstrate that this topological approach outperforms ABC approaches that use simpler statistics based on spatial features of the data. This is a first step toward a general framework of spatial parameter inference for biological systems, for which there may be a variety of filtrations, vectorizations and summary statistics to be considered. AVAILABILITY AND IMPLEMENTATION: All code used to produce our results is available as a Snakemake workflow from github.com/tt104/tabc_angio.

Geometric anomaly detection in data.

The quest for low-dimensional models which approximate high-dimensional data is pervasive across the physical, natural, and social sciences. The dominant paradigm underlying most standard modeling techniques assumes that the data are concentrated near a single unknown manifold of relatively small intrinsic dimension. Here, we present a systematic framework for detecting interfaces and related anomalies in data which may fail to satisfy the manifold hypothesis. By computing the local topology of small regions around each data point, we are able to partition a given dataset into disjoint classes, each of which can be individually approximated by a single manifold. Since these manifolds may have different intrinsic dimensions, local topology discovers singular regions in data even when none of the points have been sampled precisely from the singularities. We showcase this method by identifying the intersection of two surfaces in the 24-dimensional space of cyclo-octane conformations and by locating all of the self-intersections of a Henneberg minimal surface immersed in 3-dimensional space. Due to the local nature of the topological computations, the algorithmic burden of performing such data stratification is readily distributable across several processors.

Topological data analysis distinguishes parameter regimes in the Anderson-Chaplain model of angiogenesis.

Angiogenesis is the process by which blood vessels form from pre-existing vessels. It plays a key role in many biological processes, including embryonic development and wound healing, and contributes to many diseases including cancer and rheumatoid arthritis. The structure of the resulting vessel networks determines their ability to deliver nutrients and remove waste products from biological tissues. Here we simulate the Anderson-Chaplain model of angiogenesis at different parameter values and quantify the vessel architectures of the resulting synthetic data. Specifically, we propose a topological data analysis (TDA) pipeline for systematic analysis of the model. TDA is a vibrant and relatively new field of computational mathematics for studying the shape of data. We compute topological and standard descriptors of model simulations generated by different parameter values. We show that TDA of model simulation data stratifies parameter space into regions with similar vessel morphology. The methodologies proposed here are widely applicable to other synthetic and experimental data including wound healing, development, and plant biology.

Algebra, Geometry and Topology of ERK Kinetics.

The MEK/ERK signalling pathway is involved in cell division, cell specialisation, survival and cell death (Shaul and Seger in Biochim Biophys Acta (BBA)-Mol Cell Res 1773(8):1213-1226, 2007). Here we study a polynomial dynamical system describing the dynamics of MEK/ERK proposed by Yeung et al. (Curr Biol 2019, https://doi.org/10.1016/j.cub.2019.12.052 ) with their experimental setup, data and known biological information. The experimental dataset is a time-course of ERK measurements in different phosphorylation states following activation of either wild-type MEK or MEK mutations associated with cancer or developmental defects. We demonstrate how methods from computational algebraic geometry, differential algebra, Bayesian statistics and computational algebraic topology can inform the model reduction, identification and parameter inference of MEK variants, respectively. Throughout, we show how this algebraic viewpoint offers a rigorous and systematic analysis of such models.

Multiscale Methods for Signal Selection in Single-Cell Data.

Analysis of single-cell transcriptomics often relies on clustering cells and then performing differential gene expression (DGE) to identify genes that vary between these clusters. These discrete analyses successfully determine cell types and markers; however, continuous variation within and between cell types may not be detected. We propose three topologically motivated mathematical methods for unsupervised feature selection that consider discrete and continuous transcriptional patterns on an equal footing across multiple scales simultaneously. Eigenscores (eigi) rank signals or genes based on their correspondence to low-frequency intrinsic patterning in the data using the spectral decomposition of the Laplacian graph. The multiscale Laplacian score (MLS) is an unsupervised method for locating relevant scales in data and selecting the genes that are coherently expressed at these respective scales. The persistent Rayleigh quotient (PRQ) takes data equipped with a filtration, allowing the separation of genes with different roles in a bifurcation process (e.g., pseudo-time). We demonstrate the utility of these techniques by applying them to published single-cell transcriptomics data sets. The methods validate previously identified genes and detect additional biologically meaningful genes with coherent expression patterns. By studying the interaction between gene signals and the geometry of the underlying space, the three methods give multidimensional rankings of the genes and visualisation of relationships between them.

Coloured Noise from Stochastic Inflows in Reaction-Diffusion Systems.

In this paper, we present a framework for investigating coloured noise in reaction-diffusion systems. We start by considering a deterministic reaction-diffusion equation and show how external forcing can cause temporally correlated or coloured noise. Here, the main source of external noise is considered to be fluctuations in the parameter values representing the inflow of particles to the system. First, we determine which reaction systems, driven by extrinsic noise, can admit only one steady state, so that effects, such as stochastic switching, are precluded from our analysis. To analyse the steady-state behaviour of reaction systems, even if the parameter values are changing, necessitates a parameter-free approach, which has been central to algebraic analysis in chemical reaction network theory. To identify suitable models, we use tools from real algebraic geometry that link the network structure to its dynamical properties. We then make a connection to internal noise models and show how power spectral methods can be used to predict stochastically driven patterns in systems with coloured noise. In simple cases, we show that the power spectrum of the coloured noise process and the power spectrum of the reaction-diffusion system modelled with white noise multiply to give the power spectrum of the coloured noise reaction-diffusion system.

Joining and decomposing reaction networks.

In systems and synthetic biology, much research has focused on the behavior and design of single pathways, while, more recently, experimental efforts have focused on how cross-talk (coupling two or more pathways) or inhibiting molecular function (isolating one part of the pathway) affects systems-level behavior. However, the theory for tackling these larger systems in general has lagged behind. Here, we analyze how joining networks (e.g., cross-talk) or decomposing networks (e.g., inhibition or knock-outs) affects three properties that reaction networks may possess-identifiability (recoverability of parameter values from data), steady-state invariants (relationships among species concentrations at steady state, used in model selection), and multistationarity (capacity for multiple steady states, which correspond to multiple cell decisions). Specifically, we prove results that clarify, for a network obtained by joining two smaller networks, how properties of the smaller networks can be inferred from or can imply similar properties of the original network. Our proofs use techniques from computational algebraic geometry, including elimination theory and differential algebra.

The role of the Hes1 crosstalk hub in Notch-Wnt interactions of the intestinal crypt.

The Notch pathway plays a vital role in determining whether cells in the intestinal epithelium adopt a secretory or an absorptive phenotype. Cell fate specification is coordinated via Notch's interaction with the canonical Wnt pathway. Here, we propose a new mathematical model of the Notch and Wnt pathways, in which the Hes1 promoter acts as a hub for pathway crosstalk. Computational simulations of the model can assist in understanding how healthy intestinal tissue is maintained, and predict the likely consequences of biochemical knockouts upon cell fate selection processes. Chemical reaction network theory (CRNT) is a powerful, generalised framework which assesses the capacity of our model for monostability or multistability, by analysing properties of the underlying network structure without recourse to specific parameter values or functional forms for reaction rates. CRNT highlights the role of ß-catenin in stabilising the Notch pathway and damping oscillations, demonstrating that Wnt-mediated actions on the Hes1 promoter can induce dynamic transitions in the Notch system, from multistability to monostability. Time-dependent model simulations of cell pairs reveal the stabilising influence of Wnt upon the Notch pathway, in which ß-catenin- and Dsh-mediated action on the Hes1 promoter are key in shaping the subcellular dynamics. Where Notch-mediated transcription of Hes1 dominates, there is Notch oscillation and maintenance of fate flexibility; Wnt-mediated transcription of Hes1 favours bistability akin to cell fate selection. Cells could therefore regulate the proportion of Wnt- and Notch-mediated control of the Hes1 promoter to coordinate the timing of cell fate selection as they migrate through the intestinal epithelium and are subject to reduced Wnt stimuli. Furthermore, mutant cells characterised by hyperstimulation of the Wnt pathway may, through coupling with Notch, invert cell fate in neighbouring healthy cells, enabling an aberrant cell to maintain its neighbours in mitotically active states.

Parameter-free methods distinguish Wnt pathway models and guide design of experiments.

The canonical Wnt signaling pathway, mediated by ß-catenin, is crucially involved in development, adult stem cell tissue maintenance, and a host of diseases including cancer. We analyze existing mathematical models of Wnt and compare them to a new Wnt signaling model that targets spatial localization; our aim is to distinguish between the models and distill biological insight from them. Using Bayesian methods we infer parameters for each model from mammalian Wnt signaling data and find that all models can fit this time course. We appeal to algebraic methods (concepts from chemical reaction network theory and matroid theory) to analyze the models without recourse to specific parameter values. These approaches provide insight into aspects of Wnt regulation: the new model, via control of shuttling and degradation parameters, permits multiple stable steady states corresponding to stem-like vs. committed cell states in the differentiation hierarchy. Our analysis also identifies groups of variables that should be measured to fully characterize and discriminate between competing models, and thus serves as a guide for performing minimal experiments for model comparison.

Nanog Fluctuations in Embryonic Stem Cells Highlight the Problem of Measurement in Cell Biology.

A number of important pluripotency regulators, including the transcription factor Nanog, are observed to fluctuate stochastically in individual embryonic stem cells. By transiently priming cells for commitment to different lineages, these fluctuations are thought to be important to the maintenance of, and exit from, pluripotency. However, because temporal changes in intracellular protein abundances cannot be measured directly in live cells, fluctuations are typically assessed using genetically engineered reporter cell lines that produce a fluorescent signal as a proxy for protein expression. Here, using a combination of mathematical modeling and experiment, we show that there are unforeseen ways in which widely used reporter strategies can systematically disturb the dynamics they are intended to monitor, sometimes giving profoundly misleading results. In the case of Nanog, we show how genetic reporters can compromise the behavior of important pluripotency-sustaining positive feedback loops, and induce a bifurcation in the underlying dynamics that gives rise to heterogeneous Nanog expression patterns in reporter cell lines that are not representative of the wild-type. These findings help explain the range of published observations of Nanog variability and highlight the problem of measurement in live cells.

Nuclear to cytoplasmic shuttling of ERK promotes differentiation of muscle stem/progenitor cells.

The transition between the proliferation and differentiation of progenitor cells is a key step in organogenesis, and alterations in this process can lead to developmental disorders. The extracellular signal-regulated kinase 1/2 (ERK) signaling pathway is one of the most intensively studied signaling mechanisms that regulates both proliferation and differentiation. How a single molecule (e.g. ERK) can regulate two opposing cellular outcomes is still a mystery. Using both chick and mouse models, we shed light on the mechanism responsible for the switch from proliferation to differentiation of head muscle progenitors and implicate ERK subcellular localization. Manipulation of the fibroblast growth factor (FGF)-ERK signaling pathway in chick embryos in vitro and in vivo demonstrated that blockage of this pathway accelerated myogenic differentiation, whereas its activation diminished it. We next examined whether the spatial subcellular localization of ERK could act as a switch between proliferation (nuclear ERK) and differentiation (cytoplasmic ERK) of muscle progenitors. A myristoylated peptide that blocks importin 7-mediated ERK nuclear translocation induced robust myogenic differentiation of muscle progenitor/stem cells in both head and trunk. In the mouse, analysis of Sprouty mutant embryos revealed that increased ERK signaling suppressed both head and trunk myogenesis. Our findings, corroborated by mathematical modeling, suggest that ERK shuttling between the nucleus and the cytoplasm provides a switch-like transition between proliferation and differentiation of muscle progenitors.

Persistent homology of time-dependent functional networks constructed from coupled time series.

We use topological data analysis to study "functional networks" that we construct from time-series data from both experimental and synthetic sources. We use persistent homology with a weight rank clique filtration to gain insights into these functional networks, and we use persistence landscapes to interpret our results. Our first example uses time-series output from networks of coupled Kuramoto oscillators. Our second example consists of biological data in the form of functional magnetic resonance imaging data that were acquired from human subjects during a simple motor-learning task in which subjects were monitored for three days during a five-day period. With these examples, we demonstrate that (1) using persistent homology to study functional networks provides fascinating insights into their properties and (2) the position of the features in a filtration can sometimes play a more vital role than persistence in the interpretation of topological features, even though conventionally the latter is used to distinguish between signal and noise. We find that persistent homology can detect differences in synchronization patterns in our data sets over time, giving insight both on changes in community structure in the networks and on increased synchronization between brain regions that form loops in a functional network during motor learning. For the motor-learning data, persistence landscapes also reveal that on average the majority of changes in the network loops take place on the second of the three days of the learning process.

Reduction of dimension for nonlinear dynamical systems.

We consider reduction of dimension for nonlinear dynamical systems. We demonstrate that in some cases, one can reduce a nonlinear system of equations into a single equation for one of the state variables, and this can be useful for computing the solution when using a variety of analytical approaches. In the case where this reduction is possible, we employ differential elimination to obtain the reduced system. While analytical, the approach is algorithmic and is implemented in symbolic software such as MAPLE or SageMath. In other cases, the reduction cannot be performed strictly in terms of differential operators, and one obtains integro-differential operators, which may still be useful. In either case, one can use the reduced equation to both approximate solutions for the state variables and perform chaos diagnostics more efficiently than could be done for the original higher-dimensional system, as well as to construct Lyapunov functions which help in the large-time study of the state variables. A number of chaotic and hyperchaotic dynamical systems are used as examples in order to motivate the approach.

Evasion of the Innate Immune Type I Interferon System by Monkeypox Virus.

UNLABELLED: The vaccinia virus (VACV) E3 protein has been shown to be important for blocking activation of the cellular innate immune system and allowing viral replication to occur unhindered. Mutation or deletion of E3L severely affects viral host range and pathogenesis. While the monkeypox virus (MPXV) genome encodes a homologue of the VACV E3 protein, encoded by the F3L gene, the MPXV gene is predicted to encode a protein with a truncation of 37 N-terminal amino acids. VACV with a genome encoding a similarly truncated E3L protein (VACV-E3LΔ37N) has been shown to be attenuated in mouse models, and infection with VACV-E3LΔ37N has been shown to lead to activation of the host antiviral protein kinase R pathway. In this report, we present data demonstrating that, despite containing a truncated E3 homologue, MPXV phenotypically resembles a wild-type (wt) VACV rather than VACV-E3LΔ37N. Thus, MPXV appears to contain a gene or genes that can suppress the phenotypes associated with an N-terminal truncation in E3. The suppression maps to sequences outside F3L, suggesting that the suppression is extragenic in nature. Thus, MPXV appears to have evolved mechanisms to minimize the effects of partial inactivation of its E3 homologue. IMPORTANCE: Poxviruses have evolved to have many mechanisms to evade host antiviral innate immunity; these mechanisms may allow these viruses to cause disease. Within the family of poxviruses, variola virus (which causes smallpox) is the most pathogenic, while monkeypox virus is intermediate in pathogenicity between vaccinia virus and variola virus. Understanding the mechanisms of monkeypox virus innate immune evasion will help us to understand the evolution of poxvirus innate immune evasion capabilities, providing a better understanding of how poxviruses cause disease.

Role of seasonality on predator-prey-subsidy population dynamics.

The role of seasonality on predator-prey interactions in the presence of a resource subsidy is examined using a system of non-autonomous ordinary differential equations (ODEs). The problem is motivated by the Arctic, inhabited by the ecological system of arctic foxes (predator), lemmings (prey), and seal carrion (subsidy). We construct two nonlinear, nonautonomous systems of ODEs named the Primary Model, and the n-Patch Model. The Primary Model considers spatial factors implicitly, and the n-Patch Model considers space explicitly as a "Stepping Stone" system. We establish the boundedness of the dynamics, as well as the necessity of sufficiently nutritional food for the survival of the predator. We investigate the importance of including the resource subsidy explicitly in the model, and the importance of accounting for predator mortality during migration. We find a variety of non-equilibrium dynamics for both systems, obtaining both limit cycles and chaotic oscillations. We were then able to discuss relevant implications for biologically interesting predator-prey systems including subsidy under seasonal effects. Notably, we can observe the extinction or persistence of a species when the corresponding autonomous system might predict the opposite.

Algebraic Systems Biology: A Case Study for the Wnt Pathway.

Steady-state analysis of dynamical systems for biological networks gives rise to algebraic varieties in high-dimensional spaces whose study is of interest in their own right. We demonstrate this for the shuttle model of the Wnt signaling pathway. Here, the variety is described by a polynomial system in 19 unknowns and 36 parameters. It has degree 9 over the parameter space. This case study explores multistationarity, model comparison, dynamics within regions of the state space, identifiability, and parameter estimation, from a geometric point of view. We employ current methods from computational algebraic geometry, polyhedral geometry, and combinatorics.

Clinical drug resistance linked to interconvertible phenotypic and functional states of tumor-propagating cells in multiple myeloma.

The phenotype and function of cells enriched in tumor-propagating activity and their relationship to the phenotypic architecture in multiple myeloma (MM) are controversial. Here, in a cohort of 30 patients, we show that MM composes 4 hierarchically organized, clonally related subpopulations, which, although phenotypically distinct, share the same oncogenic chromosomal abnormalities as well as immunoglobulin heavy chain complementarity region 3 area sequence. Assessed in xenograft assays, myeloma-propagating activity is the exclusive property of a population characterized by its ability for bidirectional transition between the dominant CD19(-)CD138(+) plasma cell (PC) and a low frequency CD19(-)CD138(-) subpopulation (termed Pre-PC); in addition, Pre-PCs are more quiescent and unlike PCs, are primarily localized at extramedullary sites. As shown by gene expression profiling, compared with PCs, Pre-PCs are enriched in epigenetic regulators, suggesting that epigenetic plasticity underpins the phenotypic diversification of myeloma-propagating cells. Prospective assessment in paired, pretreatment, and posttreatment bone marrow samples shows that Pre-PCs are up to 300-fold more drug-resistant than PCs. Thus, clinical drug resistance in MM is linked to reversible, bidirectional phenotypic transition of myeloma-propagating cells. These novel biologic insights have important clinical implications in relation to assessment of minimal residual disease and development of alternative therapeutic strategies in MM.

Parameter-free model discrimination criterion based on steady-state coplanarity.

We introduce a procedure for deciding when a mass-action model is incompatible with observed steady-state data that does not require any parameter estimation. Thus, we avoid the difficulties of nonlinear optimization typically associated with methods based on parameter fitting. Instead, we borrow ideas from algebraic geometry to construct a transformation of the model variables such that any set of steady states of the model under that transformation lies on a common plane, irrespective of the values of the model parameters. Model rejection can then be performed by assessing the degree to which the transformed data deviate from coplanarity. We demonstrate our method by applying it to models of multisite phosphorylation and cell death signaling. Our framework offers a parameter-free perspective on the statistical model selection problem, which can complement conventional statistical methods in certain classes of problems where inference has to be based on steady-state data and the model structures allow for suitable algebraic relationships among the steady-state solutions.