MOTIVATION: Stochastic gene expression and cell-to-cell heterogeneity have attracted increased interest in recent years, enabled by advances in single-cell measurement technologies. These studies are also increasingly complemented by quantitative biophysical modeling, often using the framework of stochastic biochemical kinetic models. However, inferring parameters for such models (i.e., the kinetic rates of biochemical reactions) remains a technical and computational challenge, particularly doing so in a manner that can leverage high-throughput single-cell sequencing data. RESULTS: In this work, we develop a chemical master equation model reference library-based computational pipeline to infer kinetic parameters describing noisy mRNA distributions from single-cell RNA sequencing data, using the commonly applied stochastic telegraph model. The approach fits kinetic parameters via steady-state distributions, as measured across a population of cells in snapshot data. Our pipeline also serves as a tool for comprehensive analysis of parameter identifiability, in both a priori (studying model properties in the absence of data) and a posteriori (in the context of a particular dataset) use-cases. The pipeline can perform both of these tasks, i.e. inference and identifiability analysis, in an efficient and scalable manner, and also serves to disentangle contributions to uncertainty in inferred parameters from experimental noise versus structural properties of the model. We found that for the telegraph model, the majority of the parameter space is not practically identifiable from single-cell RNA sequencing data, and low experimental capture rates worsen the identifiability. Our methodological framework could be extended to other data types in the fitting of small biochemical network models. AVAILABILITY AND IMPLEMENTATION: All code relevant to this work is available at https://github.com/Read-Lab-UCI/TelegraphLikelihoodInfer, archival DOI: https://doi.org/10.5281/zenodo.16915450.

DNA methylation is a heritable epigenetic modification that plays an essential role in mammalian development. Genomic methylation patterns are dynamically maintained, with DNA methyltransferases mediating inheritance of methyl marks onto nascent DNA over cycles of replication. A recently developed experimental technique employing immunoprecipitation of bromodeoxyuridine labeled nascent DNA followed by bisulfite sequencing (Repli-BS) measures post-replication temporal evolution of cytosine methylation, thus enabling genome-wide monitoring of methylation maintenance. In this work, we combine statistical analysis and stochastic mathematical modeling to analyze Repli-BS data from human embryonic stem cells. We estimate site-specific kinetic rate constants for the restoration of methyl marks on >10 million uniquely mapped cytosines within the CpG (cytosine-phosphate-guanine) dinucleotide context across the genome using Maximum Likelihood Estimation. We find that post-replication remethylation rate constants span approximately two orders of magnitude, with half-lives of per-site recovery of steady-state methylation levels ranging from shorter than ten minutes to five hours and longer. Furthermore, we find that kinetic constants of maintenance methylation are correlated among neighboring CpG sites. Stochastic mathematical modeling provides insight to the biological mechanisms underlying the inference results, suggesting that enzyme processivity and/or collaboration can produce the observed kinetic correlations. Our combined statistical/mathematical modeling approach expands the utility of genomic datasets and disentangles heterogeneity in methylation patterns arising from replication-associated temporal dynamics versus stable cell-to-cell differences.

In many biological settings, two or more cells come into physical contact to form a cell-cell interface. In some cases, the cell-cell contact must be transient, forming on timescales of seconds. One example is offered by the T cell, an immune cell which must attach to the surface of other cells in order to decipher information about disease. The aspect ratio of these interfaces (tens of nanometers thick and tens of micrometers in diameter) puts them into the thin-layer limit, or "lubrication limit", of fluid dynamics. A key question is how the receptors and ligands on opposing cells come into contact. What are the relative roles of thermal undulations of the plasma membrane and deterministic forces from active filopodia? We use a computational fluid dynamics algorithm capable of simulating 10-nanometer-scale fluid-structure interactions with thermal fluctuations up to seconds- and microns-scales. We use this to simulate two opposing membranes, variously including thermal fluctuations, active forces, and membrane permeability. In some regimes dominated by thermal fluctuations, proximity is a rare event, which we capture by computing mean first-passage times using a Weighted Ensemble rare-event computational method. Our results demonstrate a parameter regime in which the time it takes for an active force to drive local contact actually increases if the cells are being held closer together (e.g., by nonspecific adhesion), a phenomenon we attribute to the thin-layer effect. This leads to an optimal initial cell-cell separation for fastest receptor-ligand binding, which could have relevance for the role of cellular protrusions like microvilli. We reproduce a previous experimental observation that fluctuation spatial scales are largely unaffected, but timescales are dramatically slowed, by the thin-layer effect. We also find that membrane permeability would need to be above physiological levels to abrogate the thin-layer effect.

Stochastic simulation has been a powerful tool for studying the dynamics of gene regulatory networks, particularly in terms of understanding how cell-phenotype stability and fate-transitions are impacted by noisy gene expression. However, gene networks often have dynamics characterized by multiple attractors. Stochastic simulation is often inefficient for such systems, because most of the simulation time is spent waiting for rare, barrier-crossing events to occur. We present a rare-event simulation-based method for computing epigenetic landscapes and phenotype-transitions in metastable gene networks. Our computational pipeline was inspired by studies of metastability and barrier-crossing in protein folding, and provides an automated means of computing and visualizing essential stationary and dynamic information that is generally inaccessible to conventional simulation. Applied to a network model of pluripotency in Embryonic Stem Cells, our simulations revealed rare phenotypes and approximately Markovian transitions among phenotype-states, occurring with a broad range of timescales. The relative probabilities of phenotypes and the transition paths linking pluripotency and differentiation are sensitive to global kinetic parameters governing transcription factor-DNA binding kinetics. Our approach significantly expands the capability of stochastic simulation to investigate gene regulatory network dynamics, which may help guide rational cell reprogramming strategies. Our approach is also generalizable to other types of molecular networks and stochastic dynamics frameworks.