Different definitions of keystone taxa and functions agree that they have an outsized role in maintaining community composition, and thus, the keystone concept continues to attract attention even 50 years after its introduction. In this Perspective, we base our definition of microbial keystones on the original one to explore its implications and limitations. We review different mechanisms behind keystoneness, cover the strengths and weaknesses of current keystone prediction methods and present findings on keystones discovered in recent experiments. In addition, we suggest a new prediction method for keystones based on metabolic modelling. Finally, we discuss the role of keystones in community control strategies. Overall, the development of new prediction methods and the insights from recent experiments illustrate the continued relevance of the keystone concept for microbial communities. In this Perspective, Garza et al. revisit the keystone concept and discuss its relevance for microbial keystone taxa and functions. For this, they explore the different mechanisms behind keystoneness, cover the strengths and weaknesses of current keystone prediction methods, present recent findings on keystone taxa and functions and explore the role of keystones in community control.

Secondary cell wall thickening is essential for plant structural development, providing the mechanical strength and rigidity required for upright growth. However, direct observation of this process in its endogenous developmental context within living plants has remained limited. Cellulose, the predominant component of secondary walls and the most abundant biopolymer on Earth, is synthesized at the plasma membrane by complexes containing CELLULOSE SYNTHASE A (CESA) proteins. Despite its central role, the precise timing and regulation of cellulose deposition during plant development remain unclear. To address this gap, we developed a real-time bioluminescence imaging system in the model grass Brachypodium distachyon using a luciferase transcriptional reporter driven by the CESA8 cis-regulatory region. Bioluminescence imaging revealed a consistent spatial pattern of CESA8 expression within elongating internodes, coinciding with regions undergoing secondary wall deposition and progressive increases in cellulose crystallinity. Time-lapse imaging showed that expression follows a robust daily rhythm driven by temperature cycles, independent of light or endogenous circadian signals. Temperature-pulse experiments uncovered rapid, transient inverse responses that were accurately predicted by a mathematical model based on an incoherent feedforward loop. CESA8 expression correlated strongly with stem elongation, linking structural reinforcement with temperature-driven shoot growth in grasses.

The generation of robust circadian oscillations in the suprachiasmatic nucleus relies on intercellular coupling and synchronization of cell-autonomous clocks. When uncoupled, individual cells display considerable heterogeneity in amplitude and period, and a large fraction of cells undergo damped oscillations. We review here modeling approaches and summarize theoretical findings regarding synchronization mechanisms and entrainment properties of the multi-cellular circadian network. We also discuss recent experimental efforts aiming to characterize the synchronization dynamics.

During cell differentiation, identical pluripotent cells undergo a specification process marked by changes in the expression of key genes, regulated by transcription factors that can inhibit the transcription of a competing gene or activate their own transcription. This specification is orchestrated by gene regulatory networks (GRNs), encompassing transcription factors, biochemical reactions, and signalling cascades. Mathematical models for these GRNs have been proposed in various contexts, to replicate observed robustness in differentiation properties. This includes reproducible proportions of differentiated cells with respect to parametric or stochastic noise and the avoidance of transitions between differentiated states. Understanding the GRN components controlling these features is crucial. Our study thoroughly explored an extended version of the Toggle Switch model with auto-activation loops. This model represents cells evolving from common progenitors in one out of two fates (A or B, bistable regime) or, additionally, remaining in their progenitor state (C, tristable regime). Such a differentiation into populations with three distinct cell fates is observed during blastocyst formation in mammals, where inner cell mass cells can remain in that state or differentiate into epiblast cells or primitive endoderm. Systematic analysis revealed that the existence of a stable non-differentiated state significantly impacts the GRN's robustness against parametric variations and stochastic noise. This state reduces the sensitivity of cell populations to parameters controlling key gene expression asymmetry and prevents cells from making transitions after acquiring a new identity. Stochastic noise enhances robustness by decreasing sensitivity to initial expression levels and helping the system escape from the non-differentiated state to differentiated cell fates, making the differentiation more efficient.

Several human-associated microbial communities exist in multiple configurations and can change their composition in response to perturbations, remaining in an altered state even after the perturbation ends. Multistability has been previously proposed to explain this behavior for gut microbiota in particular, but has not been clearly demonstrated experimentally. Here, we first investigate the life history strategies of three common human gut bacteria to identify mechanisms driving alternative states. We then use this data to build and parameterize a kinetic model, which predicts that alternative states emerge due to phenotype switching between subpopulations of the same species. Perturbation experiments support these predictions, and confirm the existence of alternative states. Finally, simulations show that phenotype switching can also explain alternative states in larger communities. Thus, a transient perturbation combined with metabolic flexibility is sufficient for alternative communities to emerge.

In the preimplantation mammalian embryo, stochastic cell-to-cell expression heterogeneity is followed by signal reinforcement to initiate the specification of Inner Cell Mass (ICM) cells into Epiblast (Epi). The expression of NANOG, the key transcription factor for the Epi fate, is necessary but not sufficient: coincident expression of other factors is required. To identify possible Nanog-helper genes, we analyzed gene expression variability in five time-stamped single-cell transcriptomic datasets using differential entropy, a quantitative measure of cell-to-cell heterogeneity. The entropy of Nanog displays a peak-shaped temporal pattern from the 16-cell to the 64-cell stage, consistent with its key role in Epi specification. By estimating the entropy profiles of the 21 genes common to all five datasets, we identified three genes - Pecam1, Sox2, and Hnf4a - whose variability in expression patterns mirrors that of Nanog. We further performed gene regulatory network inference using CARDAMOM, an algorithm that exploits temporal dynamics and transcriptional bursting. The results revealed that these three genes exhibit reciprocal activation with Nanog at the 32-cell stage. This regulatory motif reinforces fate-switching decisions and co-expression states. Our innovative analysis of single-cell transcriptomic data thus uncovers a likely role for Pecam1, Sox2, and Hnf4a as key genes that, when coincidentally expressed with Nanog, initiate ICM differentiation.

Several human-associated microbial communities exist in multiple configurations and can change their composition in response to perturbations, remaining in an altered state even after the perturbation ends. Multistability has been previously proposed to explain this behavior for gut microbiota in particular, but has not been clearly demonstrated experimentally. Here, we first investigated the life history strategies of three common human gut bacteria to identify mechanisms driving alternative states. We then used this data to build and parameterize a kinetic model, which predicted that alternative states emerge due to phenotype switching between subpopulations of the same species. Perturbation experiments supported these predictions, and confirmed the existence of alternative states. Finally, simulations showed that phenotype switching can also explain alternative states in larger communities. Thus, a transient perturbation combined with metabolic flexibility is sufficient for alternative communities to emerge, implying that they are not necessarily explained by differences between individuals. ### Competing Interest Statement The authors have declared no competing interest.

Chronotherapy aims at optimising the time of day and dosing of drugs administration. This is a promising perspective because the toxicity and efficacy of many drugs show a dependence on the time of the day at which they are administrated. Efficient cancer chronotherapy requires a good understanding of the interplay between the cell cycle and the circadian clock. Computational models offer a way to study the dynamics resulting from the coupling between these two biological oscillators and to predict successful therapeutic protocols. We review here recent advances and highlight key challenges for further developments of predictive mathematical models.

Physiological processes are governed by intricate networks of transcriptional and post-translational regulations. Inter-cellular interactions and signaling pathways further modulate the response of the cells to environmental conditions. Understanding the dynamics of these systems in healthy conditions and their alterations in pathologic situations requires a “systems” approach. Computational models allow to formalize and to simulate the dynamics of complex networks. Here, we briefly illustrate, through a few selected examples, how modeling helps to answer non-trivial questions regarding rhythmic phenomena, signaling and decision-making in cellular systems. These examples relate to cell differentiation, metabolic regulation, chronopharmacology and calcium dynamics.

Background: Heterogeneity in gene expression among single cells is critical for cell differentiation. In the preimplantation mammalian embryo, stochastic cell-to-cell expression heterogeneity is followed by signal reinforcement to initiate the differentiation of the Inner Cell Mass (ICM) cells into Epiblast (Epi) or Primitive Endoderm. A key question is to identify genes that initiate and drive coordination of pluripotency factors and specification of Epi cells alongside the key transcription factor NANOG. Results: To address this question, we fit single-cell expression data to gamma distributions and compute their inter-cellular differential entropy profile. Using five single-cell transcriptomic datasets (RT-qPCR and RNA-seq), we recover genes that are well-known to play a key role in Epi specification and identify candidate genes, whose coordinated heterogeneity in expression, alongside Nanog drives ICM cells differentiation. We show that other typical measures of variability, such as the Fano factor or the coefficient of variation, show a less robust temporal profile than differential inter-cellular entropy. Moreover, candidate genes would not have been identified with principal component and correlation analyses. Conclusions: The identified candidate genes, with a differential entropy time profile similar to that of Nanog, could play a crucial role in the differentiation of ICM cells. The estimation of cell-to-cell variability in gene expression based on inter-cellular differential entropy is robust among datasets and reveals insights that other standard methods such as principal component analysis and Fano factors do not. Differential entropy thus provides a useful quantity to assess gene regulation during early embryonic development.

Prediabetes, a subclinical state of high glucose, carries a risk of transition to diabetes. One cause of prediabetes is insulin resistance, which impairs the ability of insulin to control blood glucose. However, many individuals with high insulin resistance retain normal glucose due to compensation by enhanced insulin secretion by beta cells. Individuals seem to differ in their maximum compensation level, termed beta cell carrying capacity, such that low carrying capacity is associated with a higher risk of prediabetes and diabetes. Carrying capacity has not been quantified using a mathematical model and cannot be estimated directly from measured glucose and insulin levels in patients, unlike insulin resistance and beta cell function which can be estimated using HOMA-IR and HOMA-B formula.Here we present a mathematical model of beta cell compensation and carrying capacity, and develop a new formula called HOMA-C to estimate it from glucose and insulin measurements. HOMA-C estimates the maximal potential beta cell function of an individual, rather than the current beta cell function. We test this approach using longitudinal cohorts of prediabetic people, finding 10-fold variation in carrying capacity. Low carrying capacity is associated with higher risk of transitioning to diabetes. We estimate the timescales of beta cell compensation and insulin resistance using large datasets, showing that, unlike previous mathematical models, the new model can explain the slow rise in glucose over decades. Our mathematical understanding of beta cell carrying capacity may help to assess the risk of prediabetes in each individual.### Competing Interest StatementThe authors have declared no competing interest.

We present a comprehensive dynamical analysis of a chaotic chemical model referred to as the autocatalator, when subject to a periodic administration of one substrate. Our investigation encompasses the dynamical characterization of both unforced and forced systems utilizing isospikes and largest Lyapunov exponents-based parameter planes, bifurcation diagrams, and analysis of complex oscillations. Additionally, we present a phase diagram showing the effect of the period and amplitude of the forcing signal on the system’s behavior. Furthermore, we show how the landscapes of parameter planes are altered in response to forcing application. This analysis contributes to a deeper understanding of the intricate dynamics induced by the periodic forcing of a chaotic system.

Ecological memory refers to the influence of past events on an ecosystem’s response to external or endogenous changes. Whereas memory’s contribution to ecosystem dynamics is well recognized, its role in microbial community dynamics has gained limited attention. Quantitative models of microbial communities have recently shown how resilience and multistability can emerge from interactions between community members under environmental perturbations. Here, we propose a general approach to incorporate memory into dynamical ecological models, and apply it to a topical microbial community model to demonstrate how memory: (1) slows down the overall dynamics, thus enhancing resistance while reducing resilience to perturbations; (2) can prevent state shifts and mitigate hysteresis, thus increasing resilience; (3) favours long transient dynamics, such as long-standing oscillations and delayed regime shifts; and (4) contributes to the emergence and persistence of alternative stable states. These results highlight the fundamental role of memory in key aspects of microbial community dynamics.

We present a simple nonlinear model describing the population dynamics of cancerous, healthy, and effector cells submitted to ionizing radiation. We examine the situations in which intermittently high doses of radiation affect the cells as it occurs in radiotherapy. We characterize the effect of various irradiation parameters, such as the total delivered dose for a given standard dose rate and the number of treatment sessions seeking situations that optimize normal cells survival. In other words, we are interested in enhancing the effectiveness of a radiotherapy treatment by adjusting the temporal pattern of the dose deliverance through an exhaustive dynamical characterization of the system. This includes the study of the parameter planes, the populations' dynamical behavior going from fixed points to chaotic oscillations, and the basins of attraction. This work leads to counter-intuitive predictions and constitutes a step towards a model of chronomodulated radiotherapy.

The dynamics of a community of four planktonic bacterial strains isolated from river water was followed in R2 broth for 72 h in batch experiments. These strains were identified as Janthinobacterium sp., Brevundimonas sp., Flavobacterium sp. and Variovorax sp. 16S rRNA gene sequencing and flow cytometry analyses were combined to monitor the change in abundance of each individual strain in bi-cultures and quadri-culture. Two interaction networks were constructed that summarize the impact of the strains on each other’s growth rate in exponential phase and carrying capacity in stationary phase. The networks agree on the absence of positive interactions but also show differences, implying that ecological interactions can be specific to particular growth phases. Janthinobacterium sp. was the fastest growing strain and dominated the co-cultures. However, its growth rate was negatively affected by the presence of other strains 10 to 100 times less abundant than Janthinobacterium sp. In general, we saw a positive correlation between growth rate and carrying capacity in this system. In addition, growth rate in monoculture was predictive of carrying capacity in co-culture. Taken together, our results highlight the necessity to take growth phases into account when measuring interactions within a microbial community. In addition, evidence that a minor strain can greatly influence the dynamics of a dominant one underlines the necessity to choose population models that do not assume a linear dependency of interaction strength to abundance of other species for accurate parameterization from such empirical data.

The mammalian cell cycle is governed by a network of cyclin/Cdk complexes which signal the progression into the successive phases of the cell division cycle. Once coupled to the circadian clock, this network produces oscillations with a 24 h period such that the progression into each phase of the cell cycle is synchronized to the day-night cycle. Here, we use a computational model for the circadian clock control of the cell cycle to investigate the entrainment in a population of cells characterized by some variability in the kinetic parameters. Our numerical simulations showed that successful entrainment and synchronization are only possible with a sufficient circadian amplitude and an autonomous period close to 24 h. Cellular heterogeneity, however, introduces some variability in the entrainment phase of the cells. Many cancer cells have a disrupted clock or compromised clock control. In these conditions, the cell cycle runs independently of the circadian clock, leading to a lack of synchronization of cancer cells. When the coupling is weak, entrainment is largely impacted, but cells maintain a tendency to divide at specific times of day. These differential entrainment features between healthy and cancer cells can be exploited to optimize the timing of anti-cancer drug administration in order to minimize their toxicity and to maximize their efficacy. We then used our model to simulate such chronotherapeutic treatments and to predict the optimal timing for anti-cancer drugs targeting specific phases of the cell cycle. Although qualitative, the model highlights the need to better characterize cellular heterogeneity and synchronization in cell populations as well as their consequences for circadian entrainment in order to design successful chronopharmacological protocols.

Bacterial growth often alters the environment, which in turn can impact interspecies interactions among bacteria. Here, we used an in vitro batch system containing mucin beads to emulate the dynamic host environment and to study its impact on the interactions between two abundant and prevalent human gut bacteria, the primary fermenter Bacteroides thetaiotaomicron and the butyrate producer Roseburia intestinalis. By combining machine learning and flow cytometry, we found that the number of viable B. thetaiotaomicron cells decreases with glucose consumption due to acid production, while R. intestinalis survives post-glucose depletion by entering a slow growth mode. Both species attach to mucin beads, but only viable cell counts of B. thetaiotaomicron increase significantly. The number of viable co-culture cells varies significantly over time compared to those of monocultures. A combination of targeted metabolomics and RNA-seq showed that the slow growth mode of R. intestinalis represents a diauxic shift towards acetate and lactate consumption, whereas B. thetaiotaomicron survives glucose depletion and low pH by foraging on mucin sugars. In addition, most of the mucin monosaccharides we tested inhibited the growth of R. intestinalis but not B. thetaiotaomicron. We encoded these causal relationships in a kinetic model, which reproduced the observed dynamics. In summary, we explored how R. intestinalis and B. thetaiotaomicron respond to nutrient scarcity and how this affects their dynamics. We highlight the importance of understanding bacterial metabolic strategies to effectively modulate microbial dynamics in changing conditions.

The human gut is a complex ecosystem consisting of hundreds of microbial species interacting with each other and with the human host. Mathematical models of the gut microbiome integrate our knowledge of this system and help to formulate hypotheses to explain observations. The generalized Lotka-Volterra model has been widely used for this purpose, but it does not describe interaction mechanisms and thus does not account for metabolic flexibility. Recently, models that explicitly describe gut microbial metabolite production and consumption have become popular. These models have been used to investigate the factors that shape gut microbial composition and to link specific gut microorganisms to changes in metabolite concentrations found in diseases. Here, we review how such models are built and what we have learned so far from their application to human gut microbiome data. In addition, we discuss current challenges of these models and how these can be addressed in the future.

To survive and reproduce, a cell must process information from its environment and its own internal state and respond accordingly, in terms of metabolic activity, gene expression, movement, growth, division and differentiation. These signal–response decisions are made by complex networks of interacting genes and proteins, which function as biochemical switches and clocks, and other recognizable information-processing circuitry. This theme issue of Interface Focus (in two parts) brings together articles on time-keeping and decision-making in living cells—work that uses precise mathematical modelling of underlying molecular regulatory networks to understand important features of cell physiology. Part I focuses on time-keeping: mechanisms and dynamics of biological oscillators and modes of synchronization and entrainment of oscillators, with special attention to circadian clocks.

During development, cells from a population of common progenitors evolve towards different fates characterized by distinct levels of specific transcription factors, a process known as cell differentiation. This evolution is governed by gene regulatory networks modulated by intercellular signalling. In order to evolve towards distinct fates, cells forming the population of common progenitors must display some heterogeneity. We applied a modelling approach to obtain insights into the possible sources of cell-to-cell variability initiating the specification of cells of the inner cell mass into epiblast or primitive endoderm cells in early mammalian embryo. At the single-cell level, these cell fates correspond to three possible steady states of the model. A combination of numerical simulations and bifurcation analyses predicts that the behaviour of the model is preserved with respect to the source of variability and that cell–cell coupling induces the emergence of multiple steady states associated with various cell fate configurations, and to a distribution of the levels of expression of key transcription factors. Statistical analysis of these time-dependent distributions reveals differences in the evolutions of the variance-to-mean ratios of key variables of the system, depending on the simulated source of variability, and, by comparison with experimental data, points to the rate of synthesis of the key transcription factor NANOG as a likely initial source of heterogeneity.