Olfactory neurogenesis occurs throughout the lives of vertebrates, including in humans, and relies on the continuous differentiation and integration of neurons into a complex network. How progenitor cells convert fluctuations in cell-cell signaling into streamlined fate decisions over both space and time is poorly understood. Here, we track multicellular dynamics in the zebrafish olfactory epithelium, undertake targeted perturbations, and find that neurogenesis is driven by mutual antagonism between Notch signaling and insulinoma-associated 1a (Insm1a) that is responsive to inter-organ retinoic acid signaling. Single-cell analysis reveals that olfactory neurons emerge from transient groups of cells termed cellular neighborhoods. Stochastic modeling shows that neighborhood self-assembly is maintained by a tightly regulated bistable toggle switch. Differentiating cells migrate apically in response to brain-derived neurotrophic factor (BDNF) to take up residence as mature sensory neurons. Cumulatively, these findings reveal how stochastic signaling networks spatiotemporally regulate a balance between progenitors and derivatives, driving sustained neurogenesis in an intricate organ system.

Biochemical reaction networks are intrinsically stochastic. Such stochastic effects are pronounced when molecules have low copy numbers. In an enzymatic cascade reaction network, the formation of product molecules is dictated by the enzymatic signal. Here the noise level in the signaling molecules can have counterintuitive effects. Specifically, increased noise level can reduce the uncertainty in amplified responses of product formation. This phenomenon has been termed Stochastic Focusing. Of particular interest is the sensitivity of product formation to changes in the distribution of signaling enzyme molecules.

SUMMARY Olfactory neurogenesis occurs continuously throughout the lives of vertebrates, including in humans, and relies on the rapid, unceasing differentiation and integration of neurons into a complex multicellular network. The system-wide regulation of this intricate choreography is poorly understood; in particular, it is unclear how progenitor cells convert stochastic fluctuations in cell-cell signaling, over both space and time, into streamlined fate decisions. Here, we track single-cell level multicellular dynamics in the developing zebrafish olfactory epithelium, perturb signaling pathways with temporal specificity, and find that the continuous generation of neurons is driven by the spatially-restricted self-assembly of transient groups of progenitor cells, i.e. cellular neighborhoods. Stochastic modeling and validation of the underlying genetic circuit reveals that neighborhood self-assembly is driven by a tightly regulated bistable toggle switch between Notch signaling and the transcription factor Insulinoma-associated 1a that is responsive to inter-organ retinoic acid signaling. Newly differentiating neurons emerge from neighborhoods and, in response to brain-derived neurotrophic factor signaling, migrate across the olfactory epithelium to take up residence as apically-located, mature sensory neurons. After developmental olfactory neurogenesis is complete, inducing injury results in a robust expansion of neighborhoods, followed by neuroregeneration. Taken together, these findings provide new insights into how stochastic signaling networks spatially pattern and regulate a delicate balance between progenitors and their neuronal derivatives to drive sustained neurogenesis during both development and regeneration.

Biomolecular species such as proteins, DNAs, and RNA interact in cell to form reaction networks that regulate cellular functions. At small copy numbers, stochasticity arises from thermal fluctuations and plays important roles on cellular phenotypes. The underlying stochastic chemical kinetics (SCK) is governed by the discrete chemical master equation (dCME). Recent developments of the ACME algorithm enables the exact solution to dCME. With the constructions of exactly computed time-evolving and steady state probability landscapes, it is now possible to investigate the landscape properties of many stochastic networks, including where the probability basins are located, their significance, how they are connected, and where cycles of various dimensions appear.

The dynamics of reaction coordinates during barrier-crossing are key to understanding activated processes in complex systems such as proteins. The default assumption from Kramers' physical intuition is that of a diffusion process. However, the dynamics of barrier-crossing in natural complex molecules are largely unexplored. Here we investigate the transition dynamics of alanine dipeptide isomerization, the simplest complex system with a large number of non-reaction coordinates that can serve as an adequate thermal bath feeding energy into the reaction coordinates. We separate conformations along the time axis and construct the dynamic probability surface of reaction. We quantify its topological structure and rotational flux using persistent homology and differential form. Our results uncovered a region with a strong reactive vortex in the configuration-time space, where the highest probability peak and the transition state ensemble are located. This reactive region contains strong rotational fluxes: Most reactive trajectories swirl multiple times around this region in the subspace of the two most important reaction coordinates. Furthermore, the rotational fluxes result from cooperative movement along the isocommitter surfaces and orthogonal barrier-crossing. Overall, our findings offer a first glimpse into the reactive vortex regions that characterize the non-diffusive dynamics of barrier-crossing of a naturally occurring activation process.

Single-cell RNA sequencing is a powerful method that helps delineate the regulatory mechanisms shaping the diverse cellular populations. Heterogeneous cell populations consist of individual cells, and the expression of distinct sets of genes can differentiate one sub-population of cells from another, as they are responsible for the emergence of distinct cellular phenotypes. Of particular importance are cells at transition states that bridge these different cellular phenotypes. In this study, we develop a method to identify the cells at transition states bridging different cellular phenotypes. Our approach is based on persistent homology, which enabled us to identify the group of cells located on the boundaries between different sub-populations of cells. We applied this method to study the reprogramming of human fibroblasts toward induced pluripotent stem cells using single-cell time-course data. Even though only the data that is representative of the early stages of the reprogramming process are analyzed, we are able to uncover transient cells bridging different cell sub-populations. The most prominent group of transient cells are found to be enriched for NANOG, which is a known stem cell transcription factor that takes part in the maintenance of pluripotency and other stem cell marker genes. Overall, our method can identify cells in transient states bridging major cellular phenotypes, even though they are only a small fraction of the overall cell population. We also discuss how this approach can link the topology of the surface of cellular transcripts and bring order to the transition between cellular states and how it automatically uncovers the underlying time process.

Dynamics of reaction coordinates during barrier-crossing are key to understand activated processes in complex systems such as proteins. The default assumption from Kramers physical intuition is that of a diffusion process. However, the dynamics of barrier-crossing in natural complex molecules are largely unexplored. Here we investigate the transition dynamics of alanine-dipeptide isomerization, the simplest complex system with a large number of non-reaction coordinates that can serve as an adequate thermal bath feeding energy into the reaction coordinates. We separate conformations along the time axis and construct the dynamic probability surface of reaction. We quantify its topological structure and rotational flux using persistent homology and differential form. Our results uncovered a region with strong reactive vortex in the configuration-time space, where the highest probability peak and the transition state ensemble are located. This reactive region contains strong rotational fluxes: Most reactive trajectories swirl multiple times around this region in the subspace of the two most-important reaction coordinates. Furthermore, the rotational fluxes result from cooperative movement along the isocommitter surfaces and orthogonal barrier-crossing. Overall, our findings offer a first glimpse into the reactive vortex regions that characterize the non-diffusive dynamics of barrier-crossing of a naturally occurring activation process.

To gain insight into the reaction mechanism of activated processes, we introduce an exact approach for quantifying the topology of high-dimensional probability surfaces of the underlying dynamic processes. Instead of Morse indexes, we study the homology groups of a sequence of superlevel sets of the probability surface over high-dimensional configuration spaces using persistent homology. For alanine-dipeptide isomerization, a prototype of activated processes, we identify locations of probability peaks and connecting ridges, along with measures of their global prominence. Instead of a saddle point, the transition state ensemble (TSE) of conformations is at the most prominent probability peak after reactants/products, when proper reaction coordinates are included. Intuition-based models, even those exhibiting a double-well, fail to capture the dynamics of the activated process. Peak occurrence, prominence, and locations can be distorted upon subspace projection. While principal component analysis accounts for conformational variance, it inflates the complexity of the surface topology and destroys the dynamic properties of the topological features. In contrast, TSE emerges naturally as the most prominent peak beyond the reactant/product basins, when projected to a subspace of minimum dimension containing the reaction coordinates. Our approach is general and can be applied to investigate the topology of high-dimensional probability surfaces of other activated processes.

Biochemical reaction networks often display different behaviors and exhibit different phenotypes. In non-equilibrium networks, switching between different phenotypes can occur spontaneously under certain conditions. While there have been significant efforts towards understanding the mechanism of dynamic phenotype switching, a full understanding in many cases remains lacking. Calculation of the rotational probability flux can help elucidating the mechanism of dynamic switching. However, when the copy numbers of molecular species are small, the assumption of a system with continuous state is invalid, and the calculation of rotational probability flux becomes difficult. To the best of our knowledge, rotational probability flux has not been formulated in discrete state space. Here, we first develop a theoretical framework of discrete differential forms for stochastic reaction kinetics. We then introduce the concept of discrete rotational probability flux in high dimensional reaction networks and show how it can be computed. Results of discrete rotational probability flux in the toggle switch network are presented and compared with the amount of entropy production to facilitate understanding of dynamic phenotype switching.

Feed-forward loops (FFLs) are among the most ubiquitously found motifs of reaction networks in nature. However, little is known about their stochastic behavior and the variety of network phenotypes they can exhibit. In this study, we provide full characterizations of the properties of stochastic multimodality of FFLs, and how switching between different network phenotypes are controlled. We have computed the exact steady-state probability landscapes of all eight types of coherent and incoherent FFLs using the finite-butter Accurate Chemical Master Equation (ACME) algorithm, and quantified the exact topological features of their high-dimensional probability landscapes using persistent homology. Through analysis of the degree of multimodality for each of a set of 10,812 probability landscapes, where each landscape resides over 105-106 microstates, we have constructed comprehensive phase diagrams of all relevant behavior of FFL multimodality over broad ranges of input and regulation intensities, as well as different regimes of promoter binding dynamics. In addition, we have quantified the topological sensitivity of the multimodality of the landscapes to regulation intensities. Our results show that with slow binding and unbinding dynamics of transcription factor to promoter, FFLs exhibit strong stochastic behavior that is very different from what would be inferred from deterministic models. In addition, input intensity play major roles in the phenotypes of FFLs: At weak input intensity, FFL exhibit monomodality, but strong input intensity may result in up to 6 stable phenotypes. Furthermore, we found that gene duplication can enlarge stable regions of specific multimodalities and enrich the phenotypic diversity of FFL networks, providing means for cells toward better adaptation to changing environment. Our results are directly applicable to analysis of behavior of FFLs in biological processes such as stem cell differentiation and for design of synthetic networks when certain phenotypic behavior is desired.

The lysogeny-lysis switch in bacteriophage lambda serves as a model for understanding cell fate decisions. The molecular network controlling this switch has been explored through extensive experimental and computational studies. Yet, the specific role of protein-DNA interactions, like the binding of CI2 and Cro2 proteins to operator sites, in regulating lysogeny stability during prophage induction remains less understood. This study employs a minimalistic model and the Accurate Chemical Master Equation (ACME) method to construct detailed probability landscapes of the network's behavior under varying conditions, such as different dissociation and CI2 degradation rates which simulate UV irradiation effects. Our findings indicate that Cro2 binding at OR3 and CI2 at OR1 significantly influence lysogeny stability, with the former destabilizing and the latter stabilizing it. Conversely, interactions at OR1 by Cro2 and OR3 by CI show minimal impact on this stability. Through the ACME approach, we could examine the network's global behavior under conditions unapproachable by conventional stochastic simulations. This study highlights the critical roles of specific protein-DNA interactions in maintaining lysogeny and provides insight into the broader dynamics of the lysogeny-lysis switch under various physiological stresses. ### Competing Interest Statement The authors have declared no competing interest.

The discrete Chemical Master Equation (dCME) provides a general framework to study stochastic biochemical networks. As some molecular species may have large copy numbers of molecules while others are in small copy numbers, this divergence in copy numbers makes the explicit construction of the state space necessary for exact computation of probability distribution challenging, as the large copy number of species makes the enumeration of microstates intractable. One solution to this problem is to truncate the state space where there is little probability mass. However, to ensure the accuracy of the results, it is important to obtain an upper bound for the truncation error and minimize it a priori. Here, we expand the theoretical framework previously developed for state space truncation and provide the error bounds for truncating the state space below a minimum copy number and over a maximum copy number of each molecular specie in the network. In addition, we provide error bound for truncating inner-specie regions so microstates located in these regions can also be truncated when the probability mass is negligible. We compare the probability distribution computed on the truncated state space with that computed on the full state space. We show that a priori estimated errors can bound actual errors, when the state space is truncated from the appropriate upper, lower, and inner-specie regions. We give examples of results on the isomerization network and the feedback network, where truncations based on a priori estimated errors result in the reduction of more than 70 percent of the state space with negligible error. Our approach is applicable to more complex reaction networks.

During the process of tissue formation and regeneration, cells migrate collectively while remaining connected through intercellular adhesions. However, the roles of cell–substrate and cell–cell mechanical interactions in regulating collective cell migration are still unclear. In this study, we employ a newly developed finite element cellular model to study collective cell migration by exploring the effects of mechanical feedback between cell and substrate and mechanical signal transmission between adjacent cells. Our viscoelastic model of cells consists many triangular elements and is of high resolution. Cadherin adhesion between cells is modeled explicitly as linear springs at subcellular level. In addition, we incorporate a mechano-chemical feedback loop between cell–substrate mechanics and Rac-mediated cell protrusion. Our model can reproduce a number of experimentally observed patterns of collective cell migration during wound healing, including cell migration persistence, separation distance between cell pairs and migration direction. Moreover, we demonstrate that cell protrusion determined by the cell–substrate mechanics plays an important role in guiding persistent and oriented collective cell migration. Furthermore, this guidance cue can be maintained and transmitted to submarginal cells of long distance through intercellular adhesions. Our study illustrates that our finite element cellular model can be employed to study broad problems of complex tissue in dynamic changes at subcellular level.

Coagulation and fragmentation (CF) is a fundamental process by which particles attach to each other to form clusters while existing clusters break up into smaller ones. It is a ubiquitous process that plays a key role in many physical and biological phenomena. CF is typically a stochastic process that often occurs in confined spaces with a limited number of available particles. In this study, we use the discrete Chemical Master Equation (dCME) to describe the CF process. Using the newly developed Accurate Chemical Master Equation (ACME) method, we calculate the time-dependent behavior of the CF system. We investigate the effects of a number important factors that influence the overall behavior of the system, including the dimensionality, the ratio of attachment to detachment rates among clusters, and the initial conditions. By comparing CF in one and three dimensions we conclude that systems in higher dimensions are more likely to form large clusters. We also demonstrate how the ratio of the attachment to detachment rates affect the dynamics and the steady-state of the system. Finally, we demonstrate the relationship between the formation of large clusters and the initial condition.

Coagulation and fragmentation (CF) is a process in which particles aggregate into clusters and clusters break down into smaller clusters or particles. This ubiquitous process plays important roles in life-threatening diseases such as Alzheimer's disease and brain shrinkage. This process often happens in confined space with limited number of particles, and thus the behavior of the system is highly stochastic. A fundamental approach to study CF and its stochasticity is through solving the corresponding Discrete Chemical Master Equation (dCME), which provides exact descriptions of the time-evolving and the steady states of the system of biochemical reaction networks. However, current models of CF have limitations. For instance, stochastic simulation algorithm has difficult in sampling rare events and its convergence is difficult to determine. The error in the reconstructed steady-state probability distribution is also often unknown. Recent theoretical models consider the process as one dimensional and do not account for attachment, detachment, synthesis and degradation together. Here we study CF by solving the dCME exactly using the newly developed ACME method [1][2], and investigate systematically how synthesis, degradation, attachment and detachment of molecular particles affect the behavior of a CF system. In addition, we examine the effects of dimensionality of clusters. We demonstrate how these factors may have profound impacts on system behavior. Systematic analysis of the CF process as described in this study can help in further understanding their roles in biological problems such as amyloid-beta aggregation in neurodegenerative disease and actin cable formation. [1] Youfang Cao, Anna Terebus, and Jie Liang, SIAM Multiscale Modeling and Simulation, 2016. 14(2):923-963. [2] Youfang Cao, Anna Terebus, and Jie Liang, Bulletin of Mathematical Biology, 2016. 78(4):617-661.

Axonal microtubules are dynamically instable bundles in the interior part of the axon. The dynamics of these bundles are of vital importance in the behavior of axon such as their degeneration. Each axon typically contains 10~100 microtubule bundles with average length of 4μm. These bundles are coated with cytoplasm and are cross linked with random number of tau proteins. In some circumstances such as acceleration or deceleration of head in space or during the strike, they are placed in tension which may cause rupture of these bundles or disconnection of tau protein cross links. Mechanical behavior and rupture modality of microtubule bundles are becoming more and more important recently. In our model, viscoelastic microtubule bundles constituted from several discrete masses connected to the neighboring mass with a standard linear solid (SLS), a spring damper model. In addition we take into account the effect of cytoplasm by Dissipative Particle Dynamic (DPD) to investigate the rupture nature and mechanical behavior of these bundles and the effect of cytoplasm on their mechanical behavior. We obtain these results for various amounts of suddenly applied end forces to the group of axonal microtubule bundles.

Axon is an important part of the neuronal cells and axonal microtubules are bundles in axons. In axons, microtubules are coated with microtubule-associated protein tau, a natively unfolded filamentous protein in the central nervous system. These proteins are responsible for cross-linking axonal microtubule bundles. Through complimentary dimerization with other tau proteins, bridges are formed between nearby microtubules creating bundles. Formation of bundles of microtubules causes their transverse reinforcement and has been shown to enhance their ability to bear compressive loads. Though microtubules are conventionally regarded as bearing compressive loads, in certain circumstances during traumatic brain injuries, they are placed in tension. In our model, microtubule bundles were formed from a large number of discrete masses. We employed Standard Linear Solid model (SLS), a viscoelastic model, to computationally simulate microtubules. In this study, we investigated the dynamic responses of two dimensional axonal microtubules under suddenly applied end forces by implementing discrete masses connected to their neighboring masses with a Standard Linear Solid unit. We also investigated the effect of the applied force rate and magnitude on the deformation of bundles. Under tension, a microtubule fiber may rupture as a result of a sudden force. Using the developed model, we could predict the critical regions of the axonal microtubule bundles in the presence of varying end forces. We finally analyzed the nature of microtubular failure under varying mechanical stresses.

Axon is a filament in neuronal system and axonal microtubules are bundles in axons. In axons, microtubules are coated with microtubule-associated protein tau, a natively unfolded profuse filamentous protein in the central nervous system. These proteins are responsible for the cross-linked structure of the axonal microtubule bundles. Through complimentary dimerization with other tau proteins, bridges are formed to nearby microtubules to create bundles. The transverse reinforcement of microtubules by cross-linking to the cytoskeleton has been shown to enhance their ability to bear compressive loads. Though microtubules are conventionally regarded as bearing compressive loads, in certain circumstances such as in traumatic stretch injury, they are placed in tension. We employ Standard Linear Solid, a viscoelastic model, to computationally simulate microtubules. This study investigates the dynamic response of two dimensional axonal microtubules under suddenly applied end forces. We obtain the results for steady state behavior of axonal microtubule for different forces.

Axon is a filament in neuronal system and axonal microtubules are bundles in axons. In axons, microtubules are coated with microtubule-associated protein tau, a natively unfolded profuse filamentous protein in the central nervous system. These proteins are responsible for the cross-linked structure of the axonal microtubule bundles. Through complimentary dimerization with other tau proteins, bridges are formed to nearby microtubules to create bundles. The transverse reinforcement of microtubules by cross-linking to the cytoskeleton has been shown to enhance their ability to bear compressive loads. Though microtubules are conventionally regarded as bearing compressive loads, in certain circumstances such as in traumatic stretch injury, they are placed in tension. We employ Standard Linear Solid, a viscoelastic model, to computationally simulate microtubules. This study investigates the dynamic response of two dimensional axonal microtubules under suddenly applied end forces. We obtain the results for steady state behavior of axonal microtubule for different forces.