In many phylogenetic applications, such as cancer and virus evolution, time trees, evolutionary histories where speciation events are timed, are inferred. Of particular interest are clock-like trees, where all leaves are sampled at the same time and have equal distance to the root. One popular approach to model clock-like trees is coalescent theory, which is used in various tree inference software packages. Methodologically, phylogenetic inference methods require a tree space over which the inference is performed, and the geometry of this space plays an important role in statistical and computational aspects of tree inference algorithms. It has recently been shown that coalescent tree spaces possess a unique geometry, different from that of classical phylogenetic tree spaces. Here we introduce and study a space of discrete coalescent trees. They assume that time is discrete, which is natural in many computational applications. This tree space is a generalisation of the previously studied ranked nearest neighbour interchange space, and is built upon tree-rearrangement operations. We generalise existing results about ranked trees, including an algorithm for computing distances in polynomial time, and in particular provide new results for both the space of discrete coalescent trees and the space of ranked trees. We establish several geometrical properties of these spaces and show how these properties impact various algorithms used in phylogenetic analyses. Our tree space is a discretisation of a previously introduced time tree space, called t-space, and hence our results can be used to approximate solutions to various open problems in t-space.

Differentiation is the process whereby a cell acquires a specific phenotype, by differential gene expression as a function of time. This is thought to result from the dynamical functioning of an underlying Gene Regulatory Network (GRN). The precise path from the stochastic GRN behavior to the resulting cell state is still an open question. In this work we propose to reduce a stochastic model of gene expression, where a cell is represented by a vector in a continuous space of gene expression, to a discrete coarse-grained model on a limited number of cell types. We develop analytical results and numerical tools to perform this reduction for a specific model characterizing the evolution of a cell by a system of piecewise deterministic Markov processes (PDMP). Solving a spectral problem, we find the explicit variational form of the rate function associated to a large deviations principle, for any number of genes. The resulting Lagrangian dynamics allows us to define a deterministic limit of which the basins of attraction can be identified to cellular types. In this context the quasipotential, describing the transitions between these basins in the weak noise limit, can be defined as the unique solution of an Hamilton-Jacobi equation under a particular constraint. We develop a numerical method for approximating the coarse-grained model parameters, and show its accuracy for a symmetric toggle-switch network. We deduce from the reduced model an approximation of the stationary distribution of the PDMP system, which appears as a Beta mixture. Altogether those results establish a rigorous frame for connecting GRN behavior to the resulting cellular behavior, including the calculation of the probability of jumps between cell types.

Inherent stochasticity in gene expression leads to distributions of mRNA copy numbers in a population of identical cells. These distributions are determined primarily by the multitude of states of a gene promoter, each driving transcription at a different rate. In an era where single-cell mRNA copy number data are more and more available, there is an increasing need for fast computations of mRNA distributions. In this paper, we present a method for computing separate distributions for each species of mRNA molecules, i.e. mRNAs that have been either partially or fully processed post-transcription. The method involves the integration over all possible realizations of promoter states, which we cast into a set of linear ordinary differential equations of dimension [Formula: see text], where M is the number of available promoter states and [Formula: see text] is the mRNA copy number of species j up to which one wishes to compute the probability distribution. This approach is superior to solving the Master equation (ME) directly in two ways: (a) the number of coupled differential equations in the ME approach is [Formula: see text], where [Formula: see text] is the cutoff for the probability of the jth species of mRNA; and (b) the ME must be solved up to the cutoffs [Formula: see text], which must be selected a priori. In our approach, the equation for the probability to observe n mRNAs of any species depends only on the the probability of observing [Formula: see text] mRNAs of that species, thus yielding a correct probability distribution up to an arbitrary n. To demonstrate the validity of our derivations, we compare our results with Gillespie simulations for ten randomly selected system parameters.

Jellyfish have been called one of the most energy-efficient animals in the world due to the ease in which they move through their fluid environment, by product of their bell kinematics coupled with their morphological, muscular, material properties. We investigated jellyfish locomotion by conducting in silico comparative studies and explored swimming performance across different fluid scales (i.e., Reynolds Number), bell contraction frequencies, and contraction phase kinematics (duty cycle) for a jellyfish with a fineness ratio of 1 (ratio of bell height to bell diameter). To study these relationships, an open source implementation of the immersed boundary method was used (IB2d) to solve the fully coupled fluid-structure interaction problem of a flexible jellyfish bell in a viscous fluid. Thorough 2D parameter subspace explorations illustrated optimal parameter combinations in which give rise to enhanced swimming performance. All performance metrics indicated a higher sensitivity to bell actuation frequency than fluid scale or duty cycle, via Sobol sensitivity analysis, on a higher performance parameter subspace. Moreover, Pareto-like fronts were identified in the overall performance space involving the cost of transport and forward swimming speed. Patterns emerged within these performance spaces when highlighting different parameter regions, which complemented the global sensitivity results. Lastly, an open source computational model for jellyfish locomotion is offered to the science community that can be used as a starting place for future numerical experimentation.

We rigorously prove the passage from a Lotka-Volterra reaction-diffusion system towards a cross-diffusion system at the fast reaction limit. The system models a competition of two species, where one species has a more diverse diet than the other. The resulting limit gives a cross-diffusion system of a starvation driven type. We investigate the linear stability of homogeneous equilibria of those systems and rule out the possibility of cross-diffusion induced instability (Turing instability). Numerical simulations are included which are compatible with the theoretical results.

Cell-fate transition can be modeled by ordinary differential equations (ODEs) which describe the behavior of several molecules in interaction, and for which each stable equilibrium corresponds to a possible phenotype (or 'biological trait'). In this paper, we focus on simple ODE systems modeling two molecules which each negatively (or positively) regulate the other. It is well-known that such models may lead to monostability or multistability, depending on the selected parameters. However, extensive numerical simulations have led systems biologists to conjecture that in the vast majority of cases, there cannot be more than two stable points. Our main result is a proof of this conjecture. More specifically, we provide a criterion ensuring at most bistability, which is indeed satisfied by most commonly used functions. This includes Hill functions, but also a wide family of convex and sigmoid functions. We also determine which parameters lead to monostability, and which lead to bistability, by developing a more general framework encompassing all our results.

In this paper we consider mathematical modeling of the dynamics of self-organized patterning of spatially confined human embryonic stem cells (hESCs) treated with BMP4 (gastruloids) described in recent experimental works (Warmflash in Nat Methods 11:847-854, 2014; Chhabra in PloS Biol 17: 3000498, 2019). In the first part of the paper we use the activator-inhibitor equations of Gierer and Meinhardt to identify 3 reaction-diffusion regimes for each of the three morphogenic proteins, BMP4, Wnt and Nodal, based on the characteristic features of the dynamic patterning. We identify appropriate boundary conditions which correspond to the experimental setup and perform numerical simulations of the reaction-diffusion (RD) systems, using the finite element approximation, to confirm that the RD systems in these regimes produce realistic dynamics of the protein concentrations. In the second part of the paper we use analytic tools to address the questions of the existence and stability of non-homogeneous steady states for the reaction-diffusion systems of the type considered in the first part of the paper.

Hybrid models of genetic regulatory networks allow for a simpler analysis with respect to fully detailed quantitative models, still maintaining the main dynamical features of interest. In this paper we consider a piecewise affine model of a genetic regulatory network, in which the parameters describing the production function are affected by polytopic uncertainties. In the first part of the paper, after recalling how the problem of finding a Lyapunov function is solved in the nominal case, we present the considered polytopic uncertain system and then, after describing how to deal with sliding mode solutions, we prove a result of existence of a parameter dependent Lyapunov function subject to the solution of a feasibility linear matrix inequalities problem. In the second part of the paper, based on the previously described Lyapunov function, we are able to determine a set of domains where the system is guaranteed to converge, with the exception of a zero measure set of times, independently from the uncertainty realization. Finally a three nodes network example shows the validity of the results.

Linear functions of the site frequency spectrum (SFS) play a major role for understanding and investigating genetic diversity. Estimators of the mutation rate (e.g. based on the total number of segregating sites or average of the pairwise differences) and tests for neutrality (e.g. Tajima's D) are perhaps the most well-known examples. The distribution of linear functions of the SFS is important for constructing confidence intervals for the estimators, and to determine significance thresholds for neutrality tests. These distributions are often approximated using simulation procedures. In this paper we use multivariate phase-type theory to specify, characterize and calculate the distribution of linear functions of the site frequency spectrum. In particular, we show that many of the classical estimators of the mutation rate are distributed according to a discrete phase-type distribution. Neutrality tests, however, are generally not discrete phase-type distributed. For neutrality tests we derive the probability generating function using continuous multivariate phase-type theory, and numerically invert the function to obtain the distribution. A main result is an analytically tractable formula for the probability generating function of the SFS. Software implementation of the phase-type methodology is available in the R package PhaseTypeR, and R code for the reproduction of our results is available as an accompanying vignette.

In this paper, a stochastic eco-epidemiological system with patchy structure and transport-related infection is proposed and the stochastic dynamical behaviors are investigated. Firstly, by constructing suitable Lyapunov functions, it is revealed that there is a unique globally positive solution starting from the positive initial value. Secondly, it is proved that the presented system is stochastically ultimately bounded and the average in time of the second moment of solution is bounded. Thirdly, we prove that the large enough stochastic perturbations may lead the predator population and the diseases in the predator to be extinct while it is persistent in the deterministic system. Finally, some numerical simulations are given to test our theoretical results.

When modeling infectious diseases, it is common to assume that infection-derived immunity is either (1) non-existent or (2) perfect and lifelong. However there are many diseases in which infection-derived immunity is known to be present but imperfect. There are various ways in which infection-derived immunity can fail, which can ultimately impact the probability that an individual be reinfected by the same pathogen, as well as the long-run population-level prevalence of the pathogen. Here we discuss seven different models of imperfect infection-derived immunity, including waning, leaky and all-or-nothing immunity. For each model we derive the probability that an infected individual becomes reinfected during their lifetime, given that the system is at endemic equilibrium. This can be thought of as the impact that each of these infection-derived immunity failures have on reinfection. This measure is useful because it provides us with a way to compare different modes of failure of infection-derived immunity.

Motivated by modelling epidemics like COVID-19, this paper proposes a generalized chain binomial process which integrates two types of infectives, those with symptoms and those without. Testing of infectives and vaccination of susceptibles are then incorporated as preventive protective measures. Our interest relates to the distribution of the state of the population at the end of infection and to the reproduction number [Formula: see text] with the associated extinction condition. The method uses the construction of a family of martingales and a branching approximation for large populations, respectively. A more general branching process for epidemics is also constructed and studied. Finally, some results obtained are illustrated by numerical examples.