Biosystems

Rule-based languages have proven effective for modeling systems of interacting structured entities as typically encountered in chemistry and molecular biology. We present PyKappa, a rule-based modeling package written in Python whose interpreted nature enables interactive simulation and analysis, including by agentic AI. The package seeks to broaden the base of developers by utilizing a widely known programming language and serves as an easy-to-deploy teaching tool. Using PyKappa, we conduct a case study of phase separation.

The circadian clock allows living systems to anticipate environmental daily variations by scheduling cellular mechanisms according to the time of the day. The reciprocal interplay between metabolism and circadian clock is essential for energy homeostasis, and the disruption of metabolic clock inputs contributes to dysfunction of clocks in peripheral tissues. However, the molecular mechanisms through which the feeding/fasting cycle entrains peripheral clocks remain incompletely understood. Here, we investigate experimentally and theoretically the role of AMP-dependent kinase (AMPK), a key fasting sensor, in metabolic regulation of the circadian clock. Using a luciferase reporter driven by the Bmal1 promoter to monitor molecular clock activity in human U2OS cells, and thanks to a signal processing method based on the Hilbert transform to retrieve the instantaneous phase and period of circadian signals, we show that pharmacological activation of AMPK by AICAR significantly elevates Bmal1 promoter activity and markedly lengthens the clock period in a dose-dependent manner. Conversely, the clock period is shortened by the pharmacological inhibition by SR18292 of Peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PGC-1-alpha). Importantly, the action of AICAR is abolished by FK866, a nicotinamide phosphoribosyltransferase (NAMPT) inhibitor, demonstrating that AMPK-dependent modulation of Bmal1 promoter activity requires NAD+ availability. These experimental results are recapitulated by numerical simulations of our previously published mathematical model, where PGC-1-alpha plays a key role in mediating AMPK activity to the clock. Moreover, the comparison of experimentally measured and theoretically predicted phase response curves, following AMPK activation at different circadian phases, provides an additional validation of the model. By integrating biological experiments with mathematical modeling, our results identify a key mechanism through which metabolic factors can entrain and regulate the circadian clock via NAD+-dependent AMPK-mediated regulation of Bmal1 promoter activity.

We show that the p53-Mdm2 gene interaction can be understood as a form of antithetic feedback, but with mixed feedback loops. In particular, we show that this viewpoint sheds light on the emergence and robustness properties of p53's oscillatory behaviour, as observed experimentally.

Dysfunctions in transcription cause severe human diseases. In eukaryotes, transcription initiation by RNA polymerase (Pol) II requires the formation of the pre-initiation complex (PIC) composed of general transcription factors (GTFs), Pol II, and an essential multiprotein coactivator, Mediator. In vitro experiments with GTFs and Pol II have depicted a linear sequence of PIC assembly. However, PIC assembly paths in vivo at the genomic scale and Mediator's role in these paths remain poorly understood. One strategy characterizes ChIP-seq occupancy alterations in Mediator mutants, but analysis of these data still lacks a formalized interpretation framework. In this study, we address this need by developing a quantitative framework that scores candidate assembly paths according to their ability to explain these alterations. This framework combines a model of PIC assembly based on a system of differential equations with parameter optimization based on a penalized criterion. We applied this new framework to datasets from two yeast Mediator mutants, in the Med17 and Med10 subunits. Our results show that the best-supported PIC assembly paths differ from those inferred in vitro. In particular, we show that it is important to distinguish the core TFIIH module from the TFIIK kinase module, and that similar paths starting with Mediator or TBP receive different levels of support. Furthermore, although multiple paths fit the data equally well, they differ markedly for specific gene groups. Taken together, the results show that the path can depend on the gene, and are also compatible with more than one path being operational in vivo for a single gene. This work improves our understanding of PIC assembly paths and transcription initiation, and opens perspectives for future modeling and experimental studies.

In network biology, analyzing pathways is critical to study mechanisms or sequences of biological events. The pathway co-membership approach aims to link proteins belonging to the same pathway inside a network that is a clique. Although this method is great for projecting pathway knowledge onto proteins, it remains unsuitable for some downstream applications relying on graph distances like random walks with restart. To address this limitation, we introduce a weighted pathway co-membership method that incorporates biological distances between proteins. We achieve this by developing an abstraction framework for Reactome pathways, which captures both their hierarchical organization and sequential relationships. By assigning weights to these abstractions, we define a scoring system for proteins within the pathway co-membership clique. We illustrate our approach using a Reactome pathway ("Signaling by TGF-beta receptor complex", R-HSA-170834), and explain how we generate the weighted co-membership cliques of the 29 Reactome top-level pathways. Next, we construct the entire weighted pathway co membership network of the Reactome knowledge base. Validation and evaluation of our method is performed using cancer hallmark gene sets and a random walk with restart experiment, comparing our approach to state-of-the-art methods. We show that our approach outperforms state-of-the-art methods for gene prioritization on cancer hallmarks.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Artificial intelligence (AI) assistants can improve individual performance but may weaken unaided performance and limit interest diversity, raising questions about their long-term effects on collective intelligence (CI). Here, we consider an evolutionary model in which individuals predict, consult AI, and imitate better-performing peers. Classic reward structures that promote CI without AI can fail under Chatbot AI, whose predictions depend on users' interests, by reducing interest diversity or propagating bias. Directly rewarding or penalizing AI use does not reliably mitigate these effects. Instead, a balanced incentive that preserves minority expertise and rewards collective contributions can promote CI, but an additional incentive is needed to restore unaided collective accuracy. Our framework offers an explanation for emerging evidence of AI's societal impact and clarifies its long-term implications.

In our past study on human cilia from primary ciliary dyskinesia (PCD) patients using cryo-electron tomography and mass spectrometry, defects in the DNAH5 gene were proven to cause the loss of the entire outer dynein arm along with other proteins. This phenomenon was not observed in the unicellular green alga Chlamydomonas. In this study, we examined the loss of ciliary components caused by defects in various genes, based on proteomic and structural analyses of motile cilia from multiple PCD patients. The absence of single outer arm dynein genes causes a loss or significant decrease in other inner and outer arm dyneins. Interestingly, some patients with outer dynein defects show decreased intraflagellar transport (IFT) proteins, suggesting an influence of cargo components or assembly on transport. These phenomena were not observed in Chlamydomonas, suggesting a more complex mechanism of ciliogenesis in humans despite the high similarity in the final 3D architecture of the 9+2 axoneme. Defects in certain central pair proteins unexpectedly reduced components on the doublets as well. Our results demonstrate protein-protein interactions during stages of human motile ciliogenesis that are more complex than those in unicellular organisms.

In this study, we develop a mathematical model of the hypothalamic pituitary adrenal (HPA) axis and systematically investigate the model complexity required to reproduce observed hormone dynamics. The model is fitted to ACTH and cortisol measurements from five healthy volunteers reported in the literature. We show that simultaneous reproduction of both hormone profiles requires a delay differential equation model with at least five state variables among the model architectures considered. Following fitting the volunteer data, we investigate the response of the model to sustained stress using bifurcation analysis. The coupled system exhibits bifurcations over a range of stress intensities. This analysis reveals bifurcations whose critical stress thresholds depend strongly on the modulatory state of the system. These results provide a mathematical framework for investigating transitions between physiological and dysregulated HPA axis dynamics, which may be related to states of health and disease.

Coinfections are common in a nature and often involve competition among coinfecting parasites to extract resources from the host. As a consequence of this competition, parasites can face trade-offs between investing in competition with other parasites and exploiting the host for their own growth and transmission. Such within-host interactions among parasites can either increase or decrease virulence depending on how they alter parasite density, resource use, and host damage. We tested the inhibitory effect of different species of the entomopathogenic fungus Metarhizium on one another in vitro and the impact of coinfection on virulence towards Tribolium castaneum. We found that Metarhizium pemphigi is able to inhibit all other isolates of Metarhizium tested, whereas Metarhizium pingshaense was relatively resistant to inhibition by other Metarhizium strains. Virulence towards the host, Tribolium castaneum, differed very little among isolates tested. We did, however, observe differences between virulence in single and coinfections in two cases. Where there were virulence differences, these were both consistent with in vitro inhibition patterns and interestingly resulted in reduced virulence in coinfection. This is inconsistent with broader patterns of virulence in coinfection, where virulence seems to be driven by the more virulent parasite. It is, however, consistent with data from other systems where there is evidence of interference competition. Our data thus support the idea that, in contrast to other forms of between parasite interaction, interference competition results in reduced virulence in coinfection.

Game theory in biology started as an attempt to model animal contests by making assumptions about the fitness costs individuals pay in aggressive interactions. Among the different approaches, a notable one is to construct mechanistic models, using assumptions about behaviours and cognitive processes. Reinforcement learning is an important cognitive process, and we use it here to model pairwise and coalitional contests. We study situations where territorial neighbours become so-called dear enemies. A defending individual can get help from a neighbour through a defender-neighbour coalition against a challenger attempting territory takeover. In our model, coalition members have an advantage in contests in terms of costs of aggression, with the challenger being exposed to aggression from both members and each member only receiving part of the aggression from the challenger. We find that notable costs of border conflicts between neighbours, together with substantial advantages for coalitions, favour intervention. Neighbours intervene when they can help weaker defenders, in this way avoiding costs of border renegotiation with a new and potentially stronger territory owner. We also introduce different forms of perceived costs into the model, which we show correspond to the range from pure mutual assessment to partial self-assessment.

Trigger waves are positive feedback-mediated traveling waves that propagate without slowing down or diminishing in amplitude, enabling biochemical coordination over distances where diffusion alone would be insufficient. Apoptosis can spread as trigger waves of caspase activity, as demonstrated in Xenopus eggs and egg extracts. We measured apoptotic wave speeds in extracts at two cell cycle states and found that waves propagated more slowly in M-phase arrested cytostatic factor (CSF) extracts than in interphase extracts (21 versus 27 m min-1). Luther's equation, v=2{surd}D{kappa}, predicts wave speed from the diffusion coefficient D and the positive-feedback doubling rate {kappa}. However, D and {kappa} were comparable between the two extracts and could not explain the speed difference. This discrepancy highlighted a key limitation of the classical equation: it does not account for an activation threshold in trigger wave initiation, and in fact the threshold required to trigger apoptosis was found to differ between cell cycle states. We therefore extended Luther's equation to include such a threshold, obtaining v{propto}{surd}(AD/Uc ), where A is the reaction rate and Uc is the activation threshold. Perturbing apoptotic sensitivity confirmed the predicted threshold-speed relationship. These findings establish the activation threshold as a key determinant of trigger wave speed and provide a broadly applicable scaling law for predicting the spread of switch-like biochemical signals.

Recent biological research has shown significant interest and progress towards adoption of agentic AI systems. Domain specific agentic platforms with access to large collections of biological tools with code generation and/or querying capabilities, such as Biomni (illustrative case in point), are representative of the capability to orchestrate automated research workflows. Such platforms enable collaborative AI workspace for drug discovery, genomics, bioinformatics, literature research, and experimental planning for biologists and researchers. When such platforms are used for the aforementioned workflows, individually validating every agentic decision and the research process by which it was reached becomes challenging due to the scale of decisions involved and oftentimes, insufficient disclosure of assumptions and methodological choices, limited traceability of evidence, and interpretive judgments embedded throughout the workflow, giving rise to verification debt. Verification debt in agentic environments could result in unexpected failures and consequences in downstream experiments. This study focuses on empirical qualitative characterization of verification debt in agentic (prompt to outcomes) workflows. Biomni was chosen the testing AI agentic platform in this study mainly due to the descriptive traces it provides to track every agentic step. Verification debt is examined through two complementary angles: the individual agentic decision steps that carry claims and evidence, and the underlying research process by which the agent arrives at and reports them. This study attempts to characterize verification debt with 10 test cases (covering proteomics, genomics and drug discovery) with varying degree of biological nuance across three identical and independent runs. Verification debt is represented through 9 decision-step metrics and 13 research-process metrics, and also propose the VERIFY framework for evaluating scientific verifiability and validity in agentic workflows.

The most studied symbiosis between land vertebrates and algae, is the one between Ambystoma maculatum and Oophila amblystomatis, where the latter colonizes the amphibian''s eggs. Although clear benefit given by algal growth for the salamander embryo has been observed, mutualism is not widely accepted. Mutualistic hypotheses posit benefits from the host to the algae which include protection, N- and C-sources, but no experimental evidence has been obtained. In fact, recent studies suggest that the algae are stressed during this relationship and propose fermentation rather than photosynthesis as main driver for algal growth. In this work, we analyzed embryonic host developmental rate, algal and embryo growth and photosynthetic function of eggs housed at different light intensities. We provide evidence that photosynthetic growth is the main driver for Oophila and that higher algal growth yields faster developmental rate and embryo size. Also, we provide strong evidence that rates of photosynthetic electron transport (and the concomitant O2 evolution) must be higher than those of consumption by the embryo for this symbiosis to be mutualistic. Since photosynthetic rate is dependent on light, this work illustrates the importance of light intensity in this symbiotic relationship.

Oncogenic mutations are known to occur at non-random rates and specific hotspots, yet the structural and thermodynamic factors underlying these hotspots remain poorly understood. This study presents molecular dynamics simulations of 22 cancer hotspots across 10 oncogene families, simulated as histone-DNA complexes. Across 18 of 22 mutations, thermodynamic effects were predominantly localized to within 6 [A] of the DNA-histone interface (mean capture 104%), with van der Waals interactions driving the effect at the thermodynamic extremes, indicating that oncogenic mutations alter nucleosome binding through precise local contact changes rather than global structural rearrangements. Further, analysis revealed a gene-family correlated pattern in thermodynamic stability. RAS family mutations showed a consistent trend toward nucleosome stabilization (mean {Delta}{Delta}G = -4.82 kcal mol-1), while kinase domain mutations trended toward destabilization (mean {Delta}{Delta}G = +40.12 kcal mol-1), a difference reaching statistical significance in this exploratory analysis (Mann-Whitney U, p = 0.008). This interface-specific mechanism, combined with the gene-family-correlated thermodynamic pattern, provides a biophysical framework for understanding nucleosome-level contributions to cancer hotspot mutation biology.

Large shear deformation in cartilage can lead to tissue failure and injury, which are often precursors to diseases such as post-traumatic osteoarthritis. Under shear, cartilage strain-stiffens prior to failure, which makes stiffening an important indicator of changes to the loading of the extracellular matrix that gives rise to this tissue's properties. While a variety of stiffening mechanisms ranging from uncrimping to fiber-alignment have been explored, whether such mechanisms apply to the depth dependent shear strain stiffening in articular cartilage remains poorly understood. Here, we present experiments and theory that support the idea that the strain-stiffening behavior is caused by a mechanical phase transition from a random isotropic network whose mechanics arise primarily from bending-dominated modes to one where the mechanics arise from stretching-dominated modes. Indeed, consistent with network model predictions, we find that stress-strain data for all depths can be collapsed onto a single universal curve that diverges at a depth-dependent critical strain. These results show that increasing the local concentration of collagen, which constitutes the primary load-bearing network in cartilage, decreases the critical strain at which the strain-stiffening transition occurs, while increasing the low-strain elastic modulus. Furthermore, we demonstrate that bulk compression leads to a higher critical strain and decreased low-strain elastic modulus. Network simulations capture this rich behavior, and indicate that under compression the network collapses, which allows for larger strains before the transition from bending-dominated to stretching-dominated mechanics occurs. Collectively, these findings shed light on the roles that tissue composition and axial strain play in the mechanisms that organize how articular cartilage stiffens during the lead-up to injury.

The genetic code is the rule that (apart from minor modifications) is used by all contemporary life on Earth to translate genes into proteins. It is summarized by a 4x4x4 table which maps the 64 possible triplet RNA codons to the 20 canonical amino acids ("AAs"). Understanding the origin of this code is a fundamental open problem in biology. We review the idea of a comma-free code, and put forward a general reason to suspect the primordial code was comma-free (even though the modern code is not). We explain why this primordial code should then persist as a comma-free subcode within the modern code, and identify the five optimal candidates for such a subcode. We then introduce the idea of "compatible comma-free codes" and explain how they allow multiple codes to evolve in parallel, circumventing a standard argument suggesting that code evolution should effectively freeze when the associated genome becomes sufficiently complex. In particular, we highlight the possibility that the original code was a comma-free doublet code, which was then overtaken by a compatible comma-free triplet code, which finally evolved (after the evolution of the start codon and ribosome) to the modern saturated triplet code. To assess this idea, we classify all inequivalent compatible comma-free doublet-triplet codes and, among these, identify the three most promising predecessors of the modern code. Consistency between these three best codes and the "wobble rules" in gene translation provide non-trivial evidence for this possibility.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Evolution has given rise to the diversity of life on Earth, but it can also cause species extinctions. While major extinction events are usually associated with external drivers, in this study, we show mechanistically how evolution can trigger abrupt extinction avalanches in ecological networks without any environmental change. In particular, we investigate an eco-evolutionary simulation model of a tripartite mutualistic-antagonistic network, such as a plant-pollinator-herbivore network, with additional intraguild competition. The prevalence of self-organised extinction avalanches depends on a subtle balance of the different interaction types, where antagonism and similarity-based competition are needed for initial diversification while mutualism and trait-independent competition foster evolutionarily driven extinction events. Extinction avalanches are triggered by a mutualist evolving to interact with plants that are not subject to sufficient antagonistic control. Subsequently, the mutualist, its plant partners and their main antagonist potentially outcompete the rest of the community. Our findings point out that evolutionary murder may claim many victims simultaneously and highlight the need to take different interaction types into account to gain a comprehensive understanding of the links between evolution and biodiversity in ecological communities.

It is important to determine whether native state dynamics drive the misfolding and oligomerization of the prion protein, which are important events in prion disease, and how they are modulated by conformational conversion. Native (N) mouse prion protein (moPrP) is known to form small (OS) and large (OL) oligomers rich in {beta}-sheet, and in this study, photoinduced electron transfer-fluorescence correlation spectroscopy (PET-FCS) has been used to characterize intramolecular dynamics within individual monomeric units in both isolated OS and OL, as well as the diffusion properties of the oligomers. It is estimated that OS and OL comprise of about 15 and 55 monomeric units, respectively. Microsecond dynamics at each of the two regions that are the 1-3 and 2-3 interfaces of native protein are distinct in N, OS and OL, although they occur on very similar timescales. Analysis of the evolution of the distribution of diffusion times, determined using the maximum entropy method, indicates heterogeneity in the oligomerization reaction. Analysis of the change in the fluctuations which occur in two different timescales in the native state ensemble shows that they are damped more at the erstwhile 1-3 interface than the erstwhile 2-3 interface. The difference in the extent of damping at the erstwhile 1-3 and 2-3 interfaces can be explained on the basis of the structural changes known to occur across each region. The changes in dynamics occur concurrently in both regions, indicating that the structural changes accompanying conformational conversion also occur simultaneously during the oligomerization of moPrP.