Systemic sclerosis (SSc) is an autoimmune disease marked by excessive extracellular matrix (ECM) deposited by myofibroblasts. The disease carries the risk of pathologically progressing to internal organs, particularly the lung. Since abnormally large densities of myofibroblasts are associated with SSc, clinical studies consider trials that reduce the density of myofibroblasts: imatinib treatment ( N ), which promotes apoptosis in myofibroblasts, and SSc by fresolimumab ( S ), which inhibits TGF- β , a key growth factor of myofibroblasts. In this paper, we develop a mathematical model of SSc by a system of partial differential equations, and use it to explore a range of treatment strategies with N and S . For example, we considered administering S in fractions three weeks apart. In this case, we determined the smallest amount of S such that ρ ( t ) , the density of the ECM at time t , will continuously and oscillatingly decrease from a disease level 2 ρ 0 , where ρ 0 is the density ρ in health. Since SSc has no cure, ρ ( t ) cannot decrease below ρ 0 . We found that with the smallest amount of S , ρ ( t ) decreases over a few months to 1.19 ρ 0 and remains nearly stable thereafter. The results of the paper could be useful in the design of future clinical trials aimed to decrease the excessive extracellular matrix in SSc patients.
Extrusomes are extrusive organelles found in a variety of organisms, including cnidarians, dinoflagellates, and protists. These organelles typically contain a fluid-filled capsule that houses a structure, such as a coiled tubule or barb, which is rapidly ejected toward a target. Although the ejection mechanisms and morphology vary widely, a common feature is the rapid acceleration of the ejected structure through a fluid. In this paper, we develop an idealized model to simulate the collapse of an extrusome capsule, which enables the ejection of a barb or other internal structure. Specifically, we used the immersed boundary method to numerically simulate the collapse of a simplified capsule, modeled as the two sides of an elliptical shell with a flat plate along the bottom. As the capsule collapses, the elliptical sides straighten, ejecting both the internal fluid and the enclosed structure towards either free fluid or a flexible target. We investigate the effects of key model parameters, such as the size of the capsule opening. We also explore the role of the Reynolds number (Re), to consider the fluid dynamics across a range of regimes, from inertial-dominated flows relevant to some extrusome firings to viscous-dominated flows characteristic of cellular-scale processes. Our results demonstrate that decreasing the capsule opening gap size leads to increased firing velocity and shorter ejection times. Similarly, increasing the capsule's minor axis reduces the time it takes the barb to reach the target as a larger volume of fluid is moved. The relationship between Re and the time to contact is nonmonotonic, but higher Re values generally result in faster target contact, even at longer initial distances. Furthermore, we observe that higher Re values enhance the robustness of the target contact in different configurations. Finally, we quantify how the stiffness of the barb affects its ability to reach the target. We find that the large deformations of the flexible barbs slow their trajectories and that the stiffer barbs hit their targets sooner. These findings provide a foundational understanding of the biomechanics and fluid dynamics of extrusome ejection through the collapse of a capsule. The insights gained may contribute to the development of microinjectors in drug delivery, where precise and rapid mechanical movements are critical.
This paper provides a rigorous mathematical resolution of the open global stability problem for a "shock-and-kill" model of HIV-1/SIV infection in brain reservoirs recently formulated by Roda et al. (2021). The model explicitly incorporates the effects of latency-reversing agents and enhanced immune clearance of reactivated cells. We derive an explicit formula for the basic reproduction number R 0 , which serves as the sole threshold parameter governing viral eradication versus persistence and integrates infection pathways from both productive and latent compartments. By combining the next-generation matrix approach with an extended graph-theoretic Lyapunov method for multigraphs with parallel arcs, we rigorously establish that the disease-free equilibrium is globally asymptotically stable when R 0 ≤ 1 , whereas a unique productive equilibrium exists and is globally asymptotically stable when R 0 > 1 . To resolve the sign-indefinite quadratic perturbations induced by structurally distinct parallel transmission arcs-a fundamental bottleneck of classical graph-theoretic Lyapunov schemes-we develop a refined composite Lyapunov framework equipped with hierarchically calibrated parameters. Systematic asymptotic scaling and multi-parameter tuning eliminate indefinite cyclic quadratic interactions, securing strict negative definiteness of the Lyapunov derivative and overcoming key limitations of conventional graph-based methods. These global stability results provide a definitive mathematical answer to whether therapeutic interventions guarantee viral eradication or lead to persistent brain-reservoir infection. Furthermore, they furnish a rigorous theoretical foundation for the "shock-and-kill" strategy and establish mathematically precise conditions to guide the design of safe, effective interventions for eliminating HIV-1/SIV from CNS reservoirs.
Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naïve model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.
Cancer is one of the major global health challenges and predicting the response of malignant tumors to chemotherapy is particularly difficult due to the complex interaction between malignant cell growth and drug effects across time and space. Mathematical modeling plays an important role in understanding these dynamics, as it provides a systematic way to simulate malignant tumor behavior under treatment and to evaluate potential outcomes of different therapeutic strategies. In this work, a reaction-diffusion model is considered to analyze the dynamics of malignant tumors during chemotherapy treatment. The temporal derivative of the model equation is discretized using the Crank-Nicolson scheme, while the spatial derivatives are approximated using an improvised cubic B-spline collocation method. Nonlinearities in the governing model are treated with the Rubin-Graves method, which transforms the nonlinear terms into a tractable linear form. This approach maintains diagonal dominance of the system matrix while ensuring consistent enforcement of boundary conditions. To evaluate stability, we apply the Fourier spectral analysis to the proposed scheme. Four numerical experiments are conducted to test the effectiveness of the scheme under different treatment parameters. The computed numerical solutions and graphical representations are accurate and provide good agreement with previously reported results. Our approach is both easy to implement and computationally efficient, making it a reliable approach for biomedical simulations of tumor dynamics.
The dynamics of tumor-immune interactions within a complex tumor microenvironment are typically modeled using a system of ordinary differential equations or partial differential equations. These models introduce some unknown parameters that need to be estimated accurately and efficiently from the limited, noisy experimental data. Moreover, due to the intricate biological complexity and limitations in experimental measurements, tumor-immune dynamics are not fully understood, and therefore, only partial knowledge of the underlying physics may be available, resulting in unknown or missing terms within the system of equations. Thus, there are twofold challenges in modeling tumor dynamics: (i) accurate estimation of model parameters and (ii) discovery of the mathematical equations governing the physical and biological systems. These types of problems are referred to as gray-box identification areas, where both experimental data and partial system knowledge are used to recover unknown parameters and missing components. In this study, we develop a cancer biology-informed neural network model (CBINN) to infer the unknown parameters in the system of equations as well as to discover the missing mechanisms from sparse and noisy measurements. We test the performance of the CBINN model on three distinct nonlinear compartmental tumour-immune models and evaluate its robustness across multiple synthetic noise levels. By harnessing these highly nonlinear dynamics, our CBINN framework effectively estimates the unknown model parameters and uncovers the underlying physical laws or mathematical structures that govern these biological systems, from scattered and noisy measurements. The models chosen here represent the dynamic patterns commonly observed in compartmental models of tumor-immune interactions, thereby validating the generalizability and efficacy of our methodology. Structural and practical identifiablility of the model parameters are also discussed using computational and Fisher information matrix based analysis. This work provides valuable guidance for researchers addressing inverse problems and gray-box identification challenges in complex dynamical systems.
The establishment of polyploid populations is constrained by minority cytotype exclusion, whereby newly formed polyploids are lost when rare in predominantly diploid populations. Here, we consider this problem through a continuous-time approximation of a discrete model of tetraploid establishment. The spatio-temporal dynamics of sexually reproducing mixed-ploidy populations is then formally investigated using a reaction-diffusion framework, which allows us to determine the conditions for spatial invasion of tetraploids. We first study the local dynamics of gamete frequencies in populations composed of diploid, triploid, and tetraploid cytotypes, where υ denotes the per-generation proportion of unreduced (diploid) gametes and ϕ denotes the relative viable-gamete contribution of triploid cytotypes. The conjugacy between models of cytotype and gamete frequency dynamics is formally established. Then, we show that the system admits a bistable structure and characterize its equilibria. In one spatial dimension, the continuous-time approximation yields a closed-form traveling-wave solution and a wave speed that scales with the standard deviation of distances between mother and offspring birth locations, σ . Extending the analysis to radially symmetric geometry, we show that successful establishment from a localized tetraploid patch requires a critical nucleus of radius R c , for which we derive an asymptotic approximation R c ∼ υ - 1 σ ( 1 - ϕ ) / 2 for small υ . We confirm these analytical predictions through numerical simulations of the continuous-time model. Our analyses reveal that sufficiently large founding patches can nucleate expanding bistable waves without intrinsic tetraploid fitness advantages, or reproductive strategies that circumvent frequency-dependent selection in sexually reproducing populations.
Understanding how stable developmental patterns emerge from gene regulatory networks remains a central problem in developmental biology. Here, we study how classical homeotic mutations reshape the epigenetic landscape of the floral gene regulatory network of Arabidopsis thaliana. We represent this landscape as an Epigenetic Forest: a collection of rooted in-arborescences induced by the state transition graph of a Boolean gene regulatory network, where each tree is the basin of attraction of a stable gene expression pattern associated with a floral or meristematic identity. We apply this framework to the wild-type network, three single homeotic mutants (ap1, pi, and ag), and three double mutants (ap3-pi, ag-pi, and ap1-ag). For each genotype, we quantify landscape organization using complementary descriptors of basin structure, convergence depth, fate diversity, dominance, inequality, and Jensen-Shannon divergence from wild type. The resulting landscapes reveal distinct modes of mutant-induced deformation. The ap1 mutant restricts fate accessibility and concentrates trajectories into dominant basins, whereas pi eliminates B-function-dependent identities while largely preserving global basin organization. In contrast, ag increases effective fate diversity despite the loss of reproductive identity, reflecting defective meristem termination and convergence to WUS-associated states. Double mutants exhibit non-additive deformation: ap3-pi is indistinguishable from pi across all reported descriptors, consistent with logical saturation of the AND-like B-function module, whereas ag-pi and ap1-ag produce distinct redistributions of fate accessibility. The framework recovers canonical floral identities and experimentally observed mutant phenotypes while treating the epigenetic landscape as a finite, computable object determined by regulatory logic. The framework thus provides a topology-based description of developmental robustness, epistasis, and mutant-induced landscape deformation.
Bacteria often develop distinct phenotypes to adapt to environmental stress. In particular, they can produce biofilms, dense communities of bacteria that live in a complex extracellular matrix. While previous studies have investigated how bacterial biofilms are regulated under laboratory conditions, they have not considered (1) the data requirements necessary to estimate model parameters and (2) how bacteria respond to recurring stressors in their natural habitats. To address (1), we adapted a mechanistic population model to explore the dynamics of biofilm formation in the presence of predator stress, using synthetic data. We used a Maximum Likelihood Estimation framework to measure crucial parameters underpinning the biofilm formation dynamics. We used genetic algorithms to propose an optimal data collection schedule that minimised parameter identifiability confidence interval widths. Our sensitivity analysis revealed that, within the explored regimes, we could simplify the binding dynamics and eliminate biofilm detachment. To address (2), we proposed a structured version of our model to capture the long-term behaviour and evolutionary selection. In our extended model, the subpopulations feature different maximal rates of biofilm formation. We compared the selection under different predator types and amounts and identified key parameters that affected the speed of selection via sensitivity analysis.
Lonafarnib (LNF) is an investigational drug targeting hepatitis delta virus (HDV) but not hepatitis B virus (HBV), providing a unique opportunity to model HDV kinetics and how changes in HDV affect HBV. We performed a detailed kinetic analysis and developed a mathematical model to explain serum HBV DNA, HDV RNA and hepatitis B surface antigen (HBsAg) kinetics in 15 HBV/HDV coinfected patients receiving LNF-based treatment. After a delay of 0-2 days, patients experienced a rapid 1st-phase HDV-decline followed by either a viral plateau, 2nd slower-decline phase, or viral breakthrough (VB). LNF monotherapy led to a flat-partial-response (often followed by VB), while LNF combination therapy with ritonavir or pegylated interferon-α (PEG-IFN α ) was associated with a biphasic HDV decline (without VB). All treatments except LNF + PEG-IFN α had at least one patient experiencing an increase in HBV on-treatment. Our model successfully reproduced the observed HDV and HBV kinetics. We estimated an HDV RNA half-life of 1.26 days [95% confidence interval, CI 1.05-1.47] in serum and treatment efficacy of 94% in inhibiting HDV RNA production across all treatments [95% CI 89-97%], as reflected by the 1st phase HDV decline. The 2nd phase of HDV decline was explained by a time-dependent increase in efficacy, reaching a maximum of 98.9%. The model explained the increase in serum HBV DNA by a median four-fold [interquartile range, IQR: 1-28] increase in HBV DNA production rate when HDV declined below an inhibitory threshold. The stability of serum HBsAg was explained by a constant number of HBsAg-producing cells.
Infectious disease superspreading caused by heterogeneity in contact behavior has been observed to be an important determinant of epidemic dynamics and size in both empirical and theoretical settings. However, it has also been observed that the importance of this type of superspreading changes throughout an epidemic, generally in a decreasing manner as infections cascade from individuals with many contacts to those with fewer contacts. We provide an exact mathematical formulation of this phenomenon in strongly-immunizing (SIR) epidemics on static contact networks. Building on the edge-based modeling framework, we construct three metrics to track how superspreading changes through the course of an epidemic, respectively measuring infected nodes' contacts, exposures, and transmissions: (1) the mean degree of infected nodes, (2) the mean number of susceptible neighbors of infected nodes, and (3) the mean number of secondary cases that will be caused by newly infected nodes. We prove results about the behaviors of these metrics, highlighting the fact that their peak times all occur at less than half the time it takes for population-level infection prevalence to peak. This suggests that the importance of superspreading will be low when an epidemic is already near its peak, so contact-based control strategies are best employed as early in an outbreak as possible. We discuss implications for accurately measuring epidemiological parameters from incidence, mobility, contact tracing, and transmission data.
Haar-like wavelets sparsify the phylogenetic covariance matrices of large, uniformly random k-regular trees with overwhelmingly high probability. This motivates the Haar-like distance, a β -diversity metric that implicitly ranks the splits of a reference phylogeny by their relevance in differentiating two microbial environments, offering an interpretation as to why the environments differ compositionally. Nevertheless, uniform binary trees exhibit statistical features distinct from those of the trees used by practitioners, leaving the extent of sparsification and the practical validity of the implied Haar-like distance speculative. To address this, our manuscript examines the sparsification of phylogenetic covariance matrices of large critical beta-splitting random trees, a model introduced to better reflect real-world phylogenies. By obtaining sharp asymptotic estimates of the first and second moments of the external path length in this ensemble, we demonstrate that the Haar-like basis also pseudo-diagonalizes the phylogenetic covariance matrix of most large trees in this more realistic framework. Additionally, we devise a test to assess the statistical significance of splits in the reference phylogeny identified by the Haar-like distance. We apply the test to a well-studied microbial mat to further substantiate the presumption that the identified splits represent genuine biological signals differentiating the top and bottom layers of the mat.
Hair follicles, the organs that produce hair, go through a constant cycle composed of phases of growth, regression, and rest. During this cycle, matrix keratinocytes (MKs), the cells responsible for hair fiber synthesis, proliferate for several years and then undergo spontaneous apoptosis. Damage to MKs and perturbations in their normal dynamics result in a shortened growth phase of the hair cycle, leading to hair loss. The most common factors causing such disruption are hormonal imbalance and attacks by the immune system. Androgenetic alopecia (AGA) is a form of hair loss caused by high sensitivity to androgens, and alopecia areata (AA) is a condition where hair loss is caused by an autoimmune reaction against MKs. In this study, we inform a mathematical model for the human hair cycle with experimental data for the lengths of hair cycle phases available from male control subjects and subjects with AGA. We also connect a mathematical model for AA with estimates for the duration of hair cycle phases obtained from the literature. Subsequently, with each model we perform parameter screening, uncertainty quantification, and global sensitivity analysis, and we compare the results across control, AGA, and AA conditions. The findings reveal that, in AGA subjects, there is greater uncertainty associated with the duration of hair growth than in control subjects. Additionally, compared to control and AGA conditions, in AA it is more certain that longer hair growth phase could not be expected. The global sensitivity analysis results show that, in AGA conditions, synthesis of regulatory molecules in the dermal papilla and stem cell input to the MK population have high impact on hair growth duration, which agrees with physiological understanding for the effect of androgens on hair follicles in AGA.
Collective systems that self-organise to maximise the group's ability to collect and distribute information can be successful in environments with high spatial and temporal variation. Such organisations are abundant in nature, as sharing information is a key benefit of many biological collective systems, and have been influential in the design of many artificial collectives such as swarm robotics. Understanding how these systems may be spatially distributed to optimise their collective potential is therefore of importance in both ecology and in collective systems design. Here, we develop a mathematical model which uses an optimisation framework to determine the higher-order spatial structure of a collective that optimises group-level knowledge transfer. The domain of the objective function is a set of weighted hypergraphs, which can fully represent the spatial structure from a topological perspective. By varying the parameters within the objective function and the constraints, we determine how the optimal spatial structure may vary when individuals differ in their information gathering ability and how this variation differs in the context of resource constraints. Our key findings are that the amount of resources in the environment can lead to specific subgroup sizes being optimal for the group as a whole when individuals are homogeneous in their information gathering abilities. Further, when there is variation in information gathering abilities, our model implies that the sharing of space between smaller subgroups of the population, rather than the whole population, is optimal for collective knowledge sharing. Our results have applications across diverse contexts from behavioural ecology to bio-inspired collective systems design.
Building upon the ODE model describing the dynamics of healthy and leukemic cells introduced in Kumar et al. (2024); Stiehl and Marciniak-Czochra (2012), we propose an extended framework that incorporates a control variable representing the effects of chemotherapy. This extension aims to provide a more refined mathematical basis for investigating anti-cancer strategies. First, we perform a stability analysis of the equilibria associated with healthy and leukemic states, partly estimated from clinical data. This analysis reveals a complex structure, including the emergence of a continuum of coexistence states and bifurcation thresholds that play a key role in the subsequent optimization stage. Based on these findings, we investigate an optimal control problem to minimize leukemia stem cells while limiting drug toxicity. Pontryagin's Maximum Principle provides necessary conditions for optimality, and direct numerical optimization confirms the predicted structures, motivating the study of the static problem. This static formulation reveals an unconventional feature: the cost functional becomes set-valued due to the continuum of equilibria, placing the problem outside the scope of standard methods. Simulations reveal a turnpike phenomenon, where over long time horizons the dynamic trajectories closely approximate the ideal static structure. Finally, a sensitivity analysis of the performance criterion with respect to key parameters complements the study, providing preliminary insights into which biological mechanisms may influence the optimal therapeutic outcomes. We conclude with a discussion of these findings.
The mitochondrial dicarboxylate carrier SLC25A10 mediates reversible exchange among succinate, malate, and phosphate, contributing to mitochondrial metabolic regulation. Structural studies establish a ping-pong mechanism, but most mathematical models still assume sequential binding, lacking mechanistic justification and overlooking the alternation of a single binding site. Here, we present the first mechanistically derived and thermodynamically consistent model of SLC25A10 based on a ping-pong framework. The model incorporates competitive binding of succinate, malate, and phosphate, heteroexchange, reversibility, and electroneutrality, and is calibrated using experimental datasets from intact mitochondria and reconstituted proteoliposomes. To estimate kinetic parameters and quantify their uncertainty, we employed Bayesian inference, enabling statistically rigorous calibration to uptake and competition assays. The model introduces new terms that quantify which substrate and from which side of the membrane is most likely to start the transport cycle. Beyond reproducing experimentally observed exchange kinetics, the model resolves non-equilibrium transport dynamics that are difficult to access directly in classical uptake assays. In particular, the simulations reveal a two-phase response in which an initial phosphate-driven high-flux uptake regime for malate and succinate is followed by a slower redistribution phase in which the two dicarboxylates continue to readjust primarily against each other. The model also predicts that mitochondrial morphology modulates early transport behaviour, with matrix swelling increasing and matrix condensation decreasing the initial SLC25A10 flux magnitude. More broadly, the framework provides a quantitative basis for studying how substrate competition, thermodynamic driving forces, and compartment geometry shape SLC25A10-mediated exchange, and it offers a transferable modelling strategy for other carriers in the SLC25 family.
Homeostasis is widely observed in biological systems and refers to their ability to maintain an output quantity approximately constant despite variations in external disturbances. Mathematically, homeostasis can be formulated through an input-output function mapping an external parameter to an output variable. Infinitesimal homeostasis occurs at isolated points where the derivative of this input-output function vanishes, allowing tools from singularity theory and combinatorial matrix theory to characterize and classify homeostatic mechanisms in terms of network topology. Although the theoretical framework allows homeostasis subnetworks to be identified directly from combinatorial structures of the input-output network without numerical simulation, the required combinatorial enumeration becomes increasingly intractable as network size grows. Moreover, the reliance on advanced graph-theoretic concepts limits its broader accessibility and practical use across disciplines, particularly in biological applications. To overcome these limitations, we develop a Python-based algorithm that automates the identification of homeostasis subnetworks and their associated homeostasis conditions directly from network topology. Given an input-output network specified solely by its connectivity structure and the designation of input and output nodes, the algorithm automatically identifies the relevant graph-theoretical structures and enumerates all homeostatic mechanisms. We demonstrate the applicability of the algorithm across a range of biological examples, including small and large networks, networks with a single input parameter (with single or multiple input nodes), multiple input parameters, and cases where input and output coincide. This wide applicability stems from our extension of the theoretical framework from single-input-single-output networks to networks with multiple input nodes through an augmented single-input-node representation. The resulting computational framework provides a scalable and systematic approach to classifying homeostatic mechanisms in complex biological networks, facilitating the application of advanced mathematical theory to a broad range of biological systems.
The risk of vector-borne disease is highly dependent on the local community composition of hosts and vectors, as well as the means by which disease is introduced into a susceptible population. The mosquito species Aedes aegypti and Aedes albopictus are both vectors of dengue virus, but differ in their biting preferences and ability to transmit the disease. The two species compete for habitat at the larval stage and their spatial distributions are highly heterogeneous, due in part to variability in factors such as temperature and resource quality affecting the outcome of competition and the resulting abundance of each species. In addition to affecting vector population dynamics, temperature also strongly affects multiple aspects of the disease transmission process. We present the basic reproduction number R 0 for a deterministic temperature-dependent transmission model between humans, Ae. aegypti, and Ae. albopictus, then develop a stochastic continuous-time Markov chain transmission model and determine the probability of disease extinction P 0 for introduction by exposed or infectious humans, Ae. aegypti, or Ae. albopictus. We explore how both R 0 and P 0 depend on a number of variables, including temperature, vector species composition, vector-host ratio, and mosquito biting behavior. We discuss our results in the context of changes in climate and neighborhood-level spread of mosquito populations and dengue.
Polyploidy occurs in plants and animals, and is an important force in speciation and genome evolution. The main focus of this paper is the following fundamental question that was recently posed by Huber and Maher: Given the ploidy numbers of a collection of extant species, or their ploidy profile, what is the smallest number of hybridizations needed in any evolutionary history for these species to completely represent these numbers? In this paper, we shall show that this question can be rephrased in terms of addition chains and the closely related addition sequences, which have been studied for over a century in mathematics and computer science. These are sequences of natural numbers that start with 1, so that each number in the sequence larger than 1 is the sum of two other numbers arising earlier in the sequence. In our first main result, we show that finding the smallest number of hybridization events to explain a ploidy profile, or the hybrid number, is equivalent to solving the so-called addition sequence problem. This immediately implies that computing the hybridization number is computationally intractable. Even so, it also leads to new connections to representing polyploid evolution using networks. More specifically, in our second main result we show that ploidy profiles representable by tree-child networks are exactly the addition chains, implying a polynomial-time algorithm for identifying these profiles. We then consider beaded tree-child networks, which permit the representation of autopolyploidy events, and in our third main result we provide a greedy polynomial-time algorithm to decide whether a given profile can be realized by such a network. We expect that our results can be leveraged in future work through, for example, making use of known algorithms for computing short addition sequences to give bounds for the hybrid number, and in guiding network reconstruction for polyploid species.
Species sharing a habitat will co-evolve to make use of the available resources, as consumption is modulated by competition and negative feedback loops between consumers and resources. The dietary range of a given species determines the resources it has access to and thus the other species with which it competes. A narrow dietary range avoids competition at the cost of over-reliance on a small selection of resources; conversely a wide dietary range provides more alternatives but also more chance of competition with other species. Here, we investigate the evolution of dietary range within a mathematical model of niche formation. We find highly path dependent co-evolution dynamics characterised by long-lived quasi-stable states. Ultimately, stochastic effects drive the evolution of generalist diets, as we uncover in our analysis and simulations.