Optimal transport reveals dynamic gene regulatory networks via gene velocity estimation.

Inferring gene regulatory networks from gene expression data is an important and challenging problem in the biology community. We propose OTVelo, a methodology that takes time-stamped single-cell gene expression data as input and predicts gene regulation across two time points. It is known that the rate of change of gene expression, which we will refer to as gene velocity, provides crucial information that enhances such inference; however, this information is not always available due to the limitations in sequencing depth. Our algorithm overcomes this limitation by estimating gene velocities using optimal transport. We then infer gene regulation using time-lagged correlation and Granger causality via regularized linear regression. Instead of providing an aggregated network across all time points, our method uncovers the underlying dynamical mechanism across time points. We validate our algorithm on 13 simulated datasets with both synthetic and curated networks and demonstrate its efficacy on 4 experimental data sets.

Parameter Identifiability in PDE Models of Fluorescence Recovery After Photobleaching.

Identifying unique parameters for mathematical models describing biological data can be challenging and often impossible. Parameter identifiability for partial differential equations models in cell biology is especially difficult given that many established in vivo measurements of protein dynamics average out the spatial dimensions. Here, we are motivated by recent experiments on the binding dynamics of the RNA-binding protein PTBP3 in RNP granules of frog oocytes based on fluorescence recovery after photobleaching (FRAP) measurements. FRAP is a widely-used experimental technique for probing protein dynamics in living cells, and is often modeled using simple reaction-diffusion models of the protein dynamics. We show that current methods of structural and practical parameter identifiability provide limited insights into identifiability of kinetic parameters for these PDE models and spatially-averaged FRAP data. We thus propose a pipeline for assessing parameter identifiability and for learning parameter combinations based on re-parametrization and profile likelihoods analysis. We show that this method is able to recover parameter combinations for synthetic FRAP datasets and investigate its application to real experimental data.

Topological data analysis of spatial patterning in heterogeneous cell populations: clustering and sorting with varying cell-cell adhesion.

Different cell types aggregate and sort into hierarchical architectures during the formation of animal tissues. The resulting spatial organization depends (in part) on the strength of adhesion of one cell type to itself relative to other cell types. However, automated and unsupervised classification of these multicellular spatial patterns remains challenging, particularly given their structural diversity and biological variability. Recent developments based on topological data analysis are intriguing to reveal similarities in tissue architecture, but these methods remain computationally expensive. In this article, we show that multicellular patterns organized from two interacting cell types can be efficiently represented through persistence images. Our optimized combination of dimensionality reduction via autoencoders, combined with hierarchical clustering, achieved high classification accuracy for simulations with constant cell numbers. We further demonstrate that persistence images can be normalized to improve classification for simulations with varying cell numbers due to proliferation. Finally, we systematically consider the importance of incorporating different topological features as well as information about each cell type to improve classification accuracy. We envision that topological machine learning based on persistence images will enable versatile and robust classification of complex tissue architectures that occur in development and disease.

SCOTv2: Single-Cell Multiomic Alignment with Disproportionate Cell-Type Representation.

Multiomic single-cell data allow us to perform integrated analysis to understand genomic regulation of biological processes. However, most single-cell sequencing assays are performed on separately sampled cell populations, as applying them to the same single-cell is challenging. Existing unsupervised single-cell alignment algorithms have been primarily benchmarked on coassay experiments. Our investigation revealed that these methods do not perform well for noncoassay single-cell experiments when there is disproportionate cell-type representation across measurement domains. Therefore, we extend our previous work-Single Cell alignment using Optimal Transport (SCOT)-by using unbalanced Gromov-Wasserstein optimal transport to handle disproportionate cell-type representation and differing sample sizes across single-cell measurements. Our method, SCOTv2, gives state-of-the-art alignment performance across five non-coassay data sets (simulated and real world). It can also integrate multiple (M≥2) single-cell measurements while preserving the self-tuning capabilities and computational tractability of its original version.

Single-Cell Multiomics Integration by SCOT.

Although the availability of various sequencing technologies allows us to capture different genome properties at single-cell resolution, with the exception of a few co-assaying technologies, applying different sequencing assays on the same single cell is impossible. Single-cell alignment using optimal transport (SCOT) is an unsupervised algorithm that addresses this limitation by using optimal transport to align single-cell multiomics data. First, it preserves the local geometry by constructing a k-nearest neighbor (k-NN) graph for each data set (or domain) to capture the intra-domain distances. SCOT then finds a probabilistic coupling matrix that minimizes the discrepancy between the intra-domain distance matrices. Finally, it uses the coupling matrix to project one single-cell data set onto another through barycentric projection, thus aligning them. SCOT requires tuning only two hyperparameters and is robust to the choice of one. Furthermore, the Gromov-Wasserstein distance in the algorithm can guide SCOT's hyperparameter tuning in a fully unsupervised setting when no orthogonal alignment information is available. Thus, SCOT is a fast and accurate alignment method that provides a heuristic for hyperparameter selection in a real-world unsupervised single-cell data alignment scenario. We provide a tutorial for SCOT and make its source code publicly available on GitHub.

SCOT: Single-Cell Multi-Omics Alignment with Optimal Transport.

Recent advances in sequencing technologies have allowed us to capture various aspects of the genome at single-cell resolution. However, with the exception of a few of co-assaying technologies, it is not possible to simultaneously apply different sequencing assays on the same single cell. In this scenario, computational integration of multi-omic measurements is crucial to enable joint analyses. This integration task is particularly challenging due to the lack of sample-wise or feature-wise correspondences. We present single-cell alignment with optimal transport (SCOT), an unsupervised algorithm that uses the Gromov-Wasserstein optimal transport to align single-cell multi-omics data sets. SCOT performs on par with the current state-of-the-art unsupervised alignment methods, is faster, and requires tuning of fewer hyperparameters. More importantly, SCOT uses a self-tuning heuristic to guide hyperparameter selection based on the Gromov-Wasserstein distance. Thus, in the fully unsupervised setting, SCOT aligns single-cell data sets better than the existing methods without requiring any orthogonal correspondence information.

Influenza spread on context-specific networks lifted from interaction-based diary data.

Studying the spread of infections is an important tool in limiting or preventing future outbreaks. A first step in understanding disease dynamics is constructing networks that reproduce features of real-world interactions. In this paper, we generate networks that maintain some features of the partial interaction networks that were recorded in an existing diary-based survey at the University of Warwick. To preserve realistic structure in our artificial networks, we use a context-specific approach. In particular, we propose different algorithms for producing larger home, work and social networks. Our networks are able to maintain much of the interaction structure in the original diary-based survey and provide a means of accounting for the interactions of survey participants with non-participants. Simulating a discrete susceptible-infected-recovered model on the full network produces epidemic behaviour which shares characteristics with previous influenza seasons. Our approach allows us to explore how disease transmission and dynamic responses to infection differ depending on interaction context. We find that, while social interactions may be the first to be reduced after influenza infection, limiting work and school encounters may be significantly more effective in controlling the overall severity of the epidemic.

Detecting Shared Genetic Architecture Among Multiple Phenotypes by Hierarchical Clustering of Gene-Level Association Statistics.

Emerging large-scale biobanks pairing genotype data with phenotype data present new opportunities to prioritize shared genetic associations across multiple phenotypes for molecular validation. Past research, by our group and others, has shown gene-level tests of association produce biologically interpretable characterization of the genetic architecture of a given phenotype. Here, we present a new method, Ward clustering to identify Internal Node branch length outliers using Gene Scores (WINGS), for identifying shared genetic architecture among multiple phenotypes. The objective of WINGS is to identify groups of phenotypes, or "clusters," sharing a core set of genes enriched for mutations in cases. We validate WINGS using extensive simulation studies and then combine gene-level association tests with WINGS to identify shared genetic architecture among 81 case-control and seven quantitative phenotypes in 349,468 European-ancestry individuals from the UK Biobank. We identify eight prioritized phenotype clusters and recover multiple published gene-level associations within prioritized clusters.

Topological data analysis of zebrafish patterns.

Self-organized pattern behavior is ubiquitous throughout nature, from fish schooling to collective cell dynamics during organism development. Qualitatively these patterns display impressive consistency, yet variability inevitably exists within pattern-forming systems on both microscopic and macroscopic scales. Quantifying variability and measuring pattern features can inform the underlying agent interactions and allow for predictive analyses. Nevertheless, current methods for analyzing patterns that arise from collective behavior capture only macroscopic features or rely on either manual inspection or smoothing algorithms that lose the underlying agent-based nature of the data. Here we introduce methods based on topological data analysis and interpretable machine learning for quantifying both agent-level features and global pattern attributes on a large scale. Because the zebrafish is a model organism for skin pattern formation, we focus specifically on analyzing its skin patterns as a means of illustrating our approach. Using a recent agent-based model, we simulate thousands of wild-type and mutant zebrafish patterns and apply our methodology to better understand pattern variability in zebrafish. Our methodology is able to quantify the differential impact of stochasticity in cell interactions on wild-type and mutant patterns, and we use our methods to predict stripe and spot statistics as a function of varying cellular communication. Our work provides an approach to automatically quantifying biological patterns and analyzing agent-based dynamics so that we can now answer critical questions in pattern formation at a much larger scale.

A Functional Data Method for Causal Dynamic Network Modeling of Task-Related fMRI.

Functional MRI (fMRI) is a popular approach to investigate brain connections and activations when human subjects perform tasks. Because fMRI measures the indirect and convoluted signals of brain activities at a lower temporal resolution, complex differential equation modeling methods (e.g., Dynamic Causal Modeling) are usually employed to infer the neuronal processes and to fit the resulting fMRI signals. However, this modeling strategy is computationally expensive and remains to be mostly a confirmatory or hypothesis-driven approach. One major statistical challenge here is to infer, in a data-driven fashion, the underlying differential equation models from fMRI data. In this paper, we propose a causal dynamic network (CDN) method to estimate brain activations and connections simultaneously. Our method links the observed fMRI data with the latent neuronal states modeled by an ordinary differential equation (ODE) model. Using the basis function expansion approach in functional data analysis, we develop an optimization-based criterion that combines data-fitting errors and ODE fitting errors. We also develop and implement a block coordinate-descent algorithm to compute the ODE parameters efficiently. We illustrate the numerical advantages of our approach using data from realistic simulations and two task-related fMRI experiments. Compared with various effective connectivity methods, our method achieves higher estimation accuracy while improving the computational speed by from tens to thousands of times. Though our method is developed for task-related fMRI, we also demonstrate the potential applicability of our method (with a simple modification) to resting-state fMRI, by analyzing both simulated and real data from medium-sized networks.

Geometry-Dependent Arrhythmias in Electrically Excitable Tissues.

Little is known about how individual cells sense the macroscopic geometry of their tissue environment. Here, we explore whether long-range electrical signaling can convey information on tissue geometry to individual cells. First, we studied an engineered electrically excitable cell line. Cells grown in patterned islands of different shapes showed remarkably diverse firing patterns under otherwise identical conditions, including regular spiking, period-doubling alternans, and arrhythmic firing. A Hodgkin-Huxley numerical model quantitatively reproduced these effects, showing how the macroscopic geometry affected the single-cell electrophysiology via the influence of gap junction-mediated electrical coupling. Qualitatively similar geometry-dependent dynamics were observed in human induced pluripotent stem cell (iPSC)-derived cardiomyocytes. The cardiac results urge caution in translating observations of arrhythmia in vitro to predictions in vivo, where the tissue geometry is very different. We study how to extrapolate electrophysiological measurements between tissues with different geometries and different gap junction couplings.

Iridophores as a source of robustness in zebrafish stripes and variability in Danio patterns.

Zebrafish (Danio rerio) feature black and yellow stripes, while related Danios display different patterns. All these patterns form due to the interactions of pigment cells, which self-organize on the fish skin. Until recently, research focused on two cell types (melanophores and xanthophores), but newer work has uncovered the leading role of a third type, iridophores: by carefully orchestrated transitions in form, iridophores instruct the other cells, but little is known about what drives their form changes. Here we address this question from a mathematical perspective: we develop a model (based on known interactions between the original two cell types) that allows us to assess potential iridophore behavior. We identify a set of mechanisms governing iridophore form that is consistent across a range of empirical data. Our model also suggests that the complex cues iridophores receive may act as a key source of redundancy, enabling both robust patterning and variability within Danio.

Modeling microtubule-based transport and anchoring of mRNA.

Localization of messenger RNA (mRNA) at the vegetal cortex plays an important role in the early development of Xenopus laevis oocytes. While it is known that molecular motors are responsible for the transport of mRNA cargo along microtubules to the cortex, the mechanisms of localization remain unclear. We model cargo transport along microtubules using partial differential equations with spatially-dependent rates. A theoretical analysis of reduced versions of our model predicts effective velocity and diffusion rates for the cargo and shows that randomness of microtubule networks enhances effective transport. A more complex model using parameters estimated from fluorescence microscopy data reproduces the spatial and timescales of mRNA localization observed in Xenopus oocytes, corroborates experimental hypotheses that anchoring may be necessary to achieve complete localization, and shows that anchoring of mRNA complexes actively transported to the cortex is most effective in achieving robust accumulation at the cortex.

Analysis of Active Transport by Fluorescence Recovery after Photobleaching.

Fluorescence recovery after photobleaching (FRAP) is a well-established experimental technique to study binding and diffusion of molecules in cells. Although a large number of analytical and numerical models have been developed to extract binding and diffusion rates from FRAP recovery curves, active transport of molecules is typically not included in the existing models that are used to estimate these rates. Here we present a validated numerical method for estimating diffusion, binding/unbinding rates, and active transport velocities using FRAP data that captures intracellular dynamics through partial differential equation models. We apply these methods to transport and localization of mRNA molecules in Xenopus laevis oocytes, where active transport processes are essential to generate developmental polarity. By providing estimates of the effective velocities and diffusion, as well as expected run times and lengths, this approach can help quantify dynamical properties of localizing and nonlocalizing RNA. Our results confirm the distinct transport dynamics in different regions of the cytoplasm, and suggest that RNA movement in both the animal and vegetal directions may influence the timescale of RNA localization in Xenopus oocytes. We also show that model initial conditions extracted from FRAP postbleach intensities prevent underestimation of diffusion, which can arise from the instantaneous bleaching assumption. The numerical and modeling approach presented here to estimate parameters using FRAP recovery data is a broadly applicable tool for systems where intracellular transport is a key molecular mechanism.

Using in vivo imaging to measure RNA mobility in Xenopus laevis oocytes.

RNA localization in the Xenopus oocyte is responsible for the establishment of polarity during oogenesis as well as the specification of germ layers during embryogenesis. However, the inability to monitor mRNA localization in live vertebrate oocytes has posed a major barrier to understanding the mechanisms driving directional transport. Here we describe a method for imaging MS2 tagged RNA in live Xenopus oocytes to study the dynamics of RNA localization. We also focus on methods for implementing and analyzing FRAP data. This protocol is optimized for imaging of the RNAs in stage II oocytes but it can be adapted to study dynamics of other molecules during oogenesis. Using this approach, mobility can be measured in different regions of the oocyte, enabling the direct observation of molecular dynamics throughout the oocyte.

Modelling stripe formation in zebrafish: an agent-based approach.

Zebrafish have distinctive black stripes and yellow interstripes that form owing to the interaction of different pigment cells. We present a two-population agent-based model for the development and regeneration of these stripes and interstripes informed by recent experimental results. Our model describes stripe pattern formation, laser ablation and mutations. We find that fish growth shortens the necessary scale for long-range interactions and that iridophores, a third type of pigment cell, help align stripes and interstripes.

Invariant manifolds and global bifurcations.

Invariant manifolds are key objects in describing how trajectories partition the phase spaces of a dynamical system. Examples include stable, unstable, and center manifolds of equilibria and periodic orbits, quasiperiodic invariant tori, and slow manifolds of systems with multiple timescales. Changes in these objects and their intersections with variation of system parameters give rise to global bifurcations. Bifurcation manifolds in the parameter spaces of multi-parameter families of dynamical systems also play a prominent role in dynamical systems theory. Much progress has been made in developing theory and computational methods for invariant manifolds during the past 25 years. This article highlights some of these achievements and remaining open problems.

Theory of passive harmonic mode-locking using waveguide arrays.

A comprehensive theoretical treatment is given of the phenomenon of harmonic mode-locking in a laser cavity mode-locked by the nonlinear mode-coupling behavior in a waveguide array. The theoretical model completely characterizes oscillatory instabilities and the transition from M to M+1 pulses as a function of increased gain.

Curvature effects on spiral spectra: generation of point eigenvalues near branch points.

It is shown that spiral waves may possess many isolated point eigenvalues that appear near branch points of the linear dispersion relation. These eigenvalues are created by the same mechanism that leads to infinitely many bound states for selfadjoint Schrödinger operators with sufficiently weakly decaying long-range potentials. For spirals, the weak decay of the potential is due to the curvature effects on the profile of the spiral in an intermediate spatial range that separates the spiral core from the far field.