Biological applications of the theory of birth-and-death processes

GeTeSEPdb: A comprehensive database and online tool for the identification and analysis of gene profiles with temporal-specific expression patterns

Gene expression is dynamic and varies at different stages of processes. The identification of gene profiles with temporal-specific expression patterns can provide valuable insights into ongoing biological processes, such as the cell cycle, cell development, circadian rhythms, or responses to external stimuli such as drug treatments or viral infections. However, currently, no database defines, identifies or archives gene profiles with temporal-specific expression patterns. Here, using a high-throughput regression analysis approach, eight linear and nonlinear parametric models were fitted to gene expression profiles from time-series experiments to identify eight types of gene profiles with temporal-specific expression patterns. We curated 2684 time-series transcriptome datasets and identified 2644,370 gene profiles exhibiting temporal-specific expression patterns. The results were stored in the database GeTeSEPdb (gene profiles with temporal-specific expression patterns database, http://www.inbirg.com/GeTeSEPdb/). Moreover, we implemented an online tool to identify gene profiles with temporal-specific expression patterns from user-submitted data. In summary, GeTeSEPdb is a comprehensive web service that can be used to identify and analyse gene profiles with temporal-specific expression patterns. This approach facilitates the exploration of transcriptional changes and temporal patterns of responses. We firmly believe that GeTeSEPdb will become a valuable resource for biologists and bioinformaticians.

A birth-death model to understand bacterial antimicrobial heteroresistance from time-kill curves

Antimicrobial heteroresistance refers to the presence of different subpopulations with heterogeneous antimicrobial responses within the same bacterial isolate, so they show reduced susceptibility compared with the main population. Though it is widely accepted that heteroresistance can play a crucial role in the outcome of antimicrobial treatments, predictive Antimicrobial Resistance (AMR) models accounting for bacterial heteroresistance are still scarce and need to be refined as the techniques to measure heteroresistance become standardised and consistent conclusions are drawn from data. In this work, we propose a multivariate Birth-Death (BD) model of bacterial heteroresistance and analyse its properties in detail. Stochasticity in the population dynamics is considered since heteroresistance is often characterised by low initial frequencies of the less susceptible subpopulations, those mediating AMR transmission and potentially leading to treatment failure. We also discuss the utility of the heteroresistance model for practical applications and calibration under realistic conditions, demonstrating that it is possible to infer the model parameters and heteroresistance distribution from time-kill data, i.e., by measuring total cell counts alone and without performing any heteroresistance test.

Resilience of weighted networks with dynamical behavior against multi-node removal

In many real-world networks, interactions between nodes are weighted to reflect their strength, such as predator-prey interactions in the ecological network and passenger numbers in airline networks. These weighted networks are prone to cascading effects caused by minor perturbations, which can lead to catastrophic outcomes. This vulnerability highlights the importance of studying weighted network resilience to prevent system collapses. However, due to many variables and weight parameters coupled together, predicting the behavior of such a system governed by a multi-dimensional rate equation is challenging. To address this, we propose a dimension reduction technique that simplifies a multi-dimensional system into a one-dimensional state space. We applied this methodology to explore the impact of weights on the resilience of four dynamics whose weights are assigned by three weight assignment methods. The four dynamical systems are the biochemical dynamical system (B), the epidemic dynamical system (E), the regulatory dynamical system (R), and the birth-death dynamical system (BD). The results show that regardless of the weight distribution, for B, the weights are negatively correlated with the activities of the network, while for E, R, and BD, there is a positive correlation between the weights and the activities of the network. Interestingly, for B, R, and BD, the change in the weights of the system has little impact on the resilience of the system. However, for the E system, the greater the weights the more resilient the system. This study not only simplifies the complexity inherent in weighted networks but also enhances our understanding of their resilience and response to perturbations.

Fundamental patterns of signal propagation in complex networks

Various disasters stem from minor perturbations, such as the spread of infectious diseases and cascading failure in power grids. Analyzing perturbations is crucial for both theoretical and application fields. Previous researchers have proposed basic propagation patterns for perturbation and explored the impact of basic network motifs on the collective response to these perturbations. However, the current framework is limited in its ability to decouple interactions and, therefore, cannot analyze more complex structures. In this article, we establish an effective, robust, and powerful propagation framework under a general dynamic model. This framework reveals classical and dense network motifs that exert critical acceleration on signal propagation, often reducing orders of magnitude compared with conclusions generated by previous work. Moreover, our framework provides a new approach to understand the fundamental principles of complex systems and the negative feedback mechanism, which is of great significance for researching system controlling and network resilience.

Clonal selection of hematopoietic stem cells after gene therapy for sickle cell disease

Gene therapy (GT) provides a potentially curative treatment option for patients with sickle cell disease (SCD); however, the occurrence of myeloid malignancies in GT clinical trials has prompted concern, with several postulated mechanisms. Here, we used whole-genome sequencing to track hematopoietic stem cells (HSCs) from six patients with SCD at pre- and post-GT time points to map the somatic mutation and clonal landscape of gene-modified and unmodified HSCs. Pre-GT, phylogenetic trees were highly polyclonal and mutation burdens per cell were elevated in some, but not all, patients. Post-GT, no clonal expansions were identified among gene-modified or unmodified cells; however, an increased frequency of potential driver mutations associated with myeloid neoplasms or clonal hematopoiesis (DNMT3A- and EZH2-mutated clones in particular) was observed in both genetically modified and unmodified cells, suggesting positive selection of mutant clones during GT. This work sheds light on HSC clonal dynamics and the mutational landscape after GT in SCD, highlighting the enhanced fitness of some HSCs harboring pre-existing driver mutations. Future studies should define the long-term fate of mutant clones, including any contribution to expansions associated with myeloid neoplasms.

A zero-agnostic model for copy number evolution in cancer

Motivation

New low-coverage single-cell DNA sequencing technologies enable the measurement of copy number profiles from thousands of individual cells within tumors. From this data, one can infer the evolutionary history of the tumor by modeling transformations of the genome via copy number aberrations. A widely used model to infer such copy number phylogenies is the copy number transformation (CNT) model in which a genome is represented by an integer vector and a copy number aberration is an event that either increases or decreases the number of copies of a contiguous segment of the genome. The CNT distance between a pair of copy number profiles is the minimum number of events required to transform one profile to another. While this distance can be computed efficiently, no efficient algorithm has been developed to find the most parsimonious phylogeny under the CNT model.

Results

We introduce the zero-agnostic copy number transformation (ZCNT) model, a simplification of the CNT model that allows the amplification or deletion of regions with zero copies. We derive a closed form expression for the ZCNT distance between two copy number profiles and show that, unlike the CNT distance, the ZCNT distance forms a metric. We leverage the closed-form expression for the ZCNT distance and an alternative characterization of copy number profiles to derive polynomial time algorithms for two natural relaxations of the small parsimony problem on copy number profiles. While the alteration of zero copy number regions allowed under the ZCNT model is not biologically realistic, we show on both simulated and real datasets that the ZCNT distance is a close approximation to the CNT distance. Extending our polynomial time algorithm for the ZCNT small parsimony problem, we develop an algorithm, Lazac, for solving the large parsimony problem on copy number profiles. We demonstrate that Lazac outperforms existing methods for inferring copy number phylogenies on both simulated and real data.

Gibbs process distinguishes survival and reveals contact-inhibition genes in Glioblastoma multiforme

Tumor growth is a spatiotemporal birth-and-death process with loss of heterotypic contact-inhibition of locomotion (CIL) of tumor cells promoting invasion and metastasis. Therefore, representing tumor cells as two-dimensional points, we can expect the tumor tissues in histology slides to reflect realizations of spatial birth-and-death process which can be mathematically modeled to reveal molecular mechanisms of CIL, provided the mathematics models the inhibitory interactions. Gibbs process as an inhibitory point process is a natural choice since it is an equilibrium process of the spatial birth-and-death process. That is if the tumor cells maintain homotypic contact inhibition, the spatial distributions of tumor cells will result in Gibbs hard core process over long time scales. In order to verify if this is the case, we applied the Gibbs process to 411 TCGA Glioblastoma multiforme patient images. Our imaging dataset included all cases for which diagnostic slide images were available. The model revealed two groups of patients, one of which - the "Gibbs group," showed the convergence of the Gibbs process with significant survival difference. Further smoothing the discretized (and noisy) inhibition metric, for both increasing and randomized survival time, we found a significant association of the patients in the Gibbs group with increasing survival time. The mean inhibition metric also revealed the point at which the homotypic CIL establishes in tumor cells. Besides, RNAseq analysis between patients with loss of heterotypic CIL and intact homotypic CIL in the Gibbs group unveiled cell movement gene signatures and differences in Actin cytoskeleton and RhoA signaling pathways as key molecular alterations. These genes and pathways have established roles in CIL. Taken together, our integrated analysis of patient images and RNAseq data provides for the first time a mathematical basis for CIL in tumors, explains survival as well as uncovers the underlying molecular landscape for this key tumor invasion and metastatic phenomenon.

Nonlinear Markov chains with finite state space: Invariant distributions and long-term behaviour

Nonlinear Markov chains with finite state space were introduced by Kolokoltsov (Nonlinear Markov Processes and Kinetic Equations, 2010). The characteristic property of these processes is that the transition probabilities depend not only on the state, but also on the distribution of the process. Here we provide first results regarding their invariant distributions and long-term behaviour: we show that under a continuity assumption an invariant distribution exists and provide a sufficient criterion for the uniqueness of the invariant distribution. Moreover, we present examples of peculiar limit behaviour that cannot occur for classical linear Markov chains. Finally, we present for the case of small state spaces sufficient (and easy-to-verify) criteria for the ergodicity of the process.

Epidemiological inference from pathogen genomes: A review of phylodynamic models and applications

Phylodynamics requires an interdisciplinary understanding of phylogenetics, epidemiology, and statistical inference. It has also experienced more intense application than ever before amid the SARS-CoV-2 pandemic. In light of this, we present a review of phylodynamic models beginning with foundational models and assumptions. Our target audience is public health researchers, epidemiologists, and biologists seeking a working knowledge of the links between epidemiology, evolutionary models, and resulting epidemiological inference. We discuss the assumptions linking evolutionary models of pathogen population size to epidemiological models of the infected population size. We then describe statistical inference for phylodynamic models and list how output parameters can be rearranged for epidemiological interpretation. We go on to cover more sophisticated models and finish by highlighting future directions.

Generalized fluctuation-dissipation theorem for non-Markovian reaction networks

Intracellular biochemical networks often display large fluctuations in the molecule numbers or the concentrations of reactive species, making molecular approaches necessary for system descriptions. For Markovian reaction networks, the fluctuation-dissipation theorem (FDT) has been well established and extensively used in fast evaluation of fluctuations in reactive species. For non-Markovian reaction networks, however, the similar FDT has not been established so far. Here, we present a generalized FDT (gFDT) for a large class of non-Markovian reaction networks where general intrinsic-event waiting-time distributions account for the effect of intrinsic noise and general stochastic reaction delays represent the impact of extrinsic noise from environmental perturbations. The starting point is a generalized chemical master equation (gCME), which describes the probabilistic behavior of an equivalent Markovian reaction network and identifies the structure of the original non-Markovian reaction network in terms of stoichiometries and effective transition rates (extensions of common reaction propensity functions). From this formulation follows directly the solution of the linear noise approximation of the stationary gCME for all the components in the non-Markovian reaction network. While the gFDT can quickly trace noisy sources in non-Markovian reaction networks, example analysis verifies its effectiveness.

First detection of threshold crossing events under intermittent sensing

The time taken by a random variable to cross a threshold for the first time, known as the first passage time, is of interest in many areas of sciences and engineering. Conventionally, there is an implicit assumption that the notional "sensor" monitoring the threshold crossing event is always active. In many realistic scenarios, the sensor monitoring the stochastic process works intermittently. Then, the relevant quantity of interest is the first detection time, which denotes the time when the sensor detects the random variable to be above the threshold for the first time. In this Letter, a birth-death process monitored by a random intermittent sensor is studied for which the first detection time distribution is obtained. In general, it is shown that the first detection time is related to and is obtainable from the first passage time distribution. Our analytical results display an excellent agreement with simulations. Furthermore, this framework is demonstrated in several applications-the susceptible infected susceptible compartmental and logistic models and birth-death processes with resetting. Finally, we solve the practically relevant problem of inferring the first passage time distribution from the first detection time.

Nearly reducible finite Markov chains: Theory and algorithms

Finite Markov chains, memoryless random walks on complex networks, appear commonly as models for stochastic dynamics in condensed matter physics, biophysics, ecology, epidemiology, economics, and elsewhere. Here, we review exact numerical methods for the analysis of arbitrary discrete- and continuous-time Markovian networks. We focus on numerically stable methods that are required to treat nearly reducible Markov chains, which exhibit a separation of characteristic timescales and are therefore ill-conditioned. In this metastable regime, dense linear algebra methods are afflicted by propagation of error in the finite precision arithmetic, and the kinetic Monte Carlo algorithm to simulate paths is unfeasibly inefficient. Furthermore, iterative eigendecomposition methods fail to converge without the use of nontrivial and system-specific preconditioning techniques. An alternative approach is provided by state reduction procedures, which do not require additional a priori knowledge of the Markov chain. Macroscopic dynamical quantities, such as moments of the first passage time distribution for a transition to an absorbing state, and microscopic properties, such as the stationary, committor, and visitation probabilities for nodes, can be computed robustly using state reduction algorithms. The related kinetic path sampling algorithm allows for efficient sampling of trajectories on a nearly reducible Markov chain. Thus, all of the information required to determine the kinetically relevant transition mechanisms, and to identify the states that have a dominant effect on the global dynamics, can be computed reliably even for computationally challenging models. Rare events are a ubiquitous feature of realistic dynamical systems, and so the methods described herein are valuable in many practical applications.

Exclusion of the fittest predicts microbial community diversity in fluctuating environments

Microorganisms live in environments that inevitably fluctuate between mild and harsh conditions. As harsh conditions may cause extinctions, the rate at which fluctuations occur can shape microbial communities and their diversity, but we still lack an intuition on how. Here, we build a mathematical model describing two microbial species living in an environment where substrate supplies randomly switch between abundant and scarce. We then vary the rate of switching as well as different properties of the interacting species, and measure the probability of the weaker species driving the stronger one extinct. We find that this probability increases with the strength of demographic noise under harsh conditions and peaks at either low, high, or intermediate switching rates depending on both species' ability to withstand the harsh environment. This complex relationship shows why finding patterns between environmental fluctuations and diversity has historically been difficult. In parameter ranges where the fittest species was most likely to be excluded, however, the beta diversity in larger communities also peaked. In sum, how environmental fluctuations affect interactions between a few species pairs predicts their effect on the beta diversity of the whole community.

Stochastic response of bacterial cells to antibiotics: its mechanisms and implications for population and evolutionary dynamics

The effectiveness of antibiotics against bacterial infections has been declining due to the emergence of resistance. Precisely understanding the response of bacteria to antibiotics is critical to maximizing antibiotic-induced bacterial eradication while minimizing the emergence of antibiotic resistance. Cell-to-cell heterogeneity in antibiotic susceptibility is observed across various bacterial species for a wide range of antibiotics. Heterogeneity in antibiotic susceptibility is not always due to the genetic differences. Rather, it can be caused by non-genetic mechanisms such as stochastic gene expression and biased partitioning upon cell division. Heterogeneous susceptibility leads to the stochastic growth and death of individual cells and stochastic fluctuations in population size. These fluctuations have important implications for the eradication of bacterial populations and the emergence of genotypic resistance.

Percolation on feature-enriched interconnected systems

Percolation is an emblematic model to assess the robustness of interconnected systems when some of their components are corrupted. It is usually investigated in simple scenarios, such as the removal of the system's units in random order, or sequentially ordered by specific topological descriptors. However, in the vast majority of empirical applications, it is required to dismantle the network following more sophisticated protocols, for instance, by combining topological properties and non-topological node metadata. We propose a novel mathematical framework to fill this gap: networks are enriched with features and their nodes are removed according to the importance in the feature space. We consider features of different nature, from ones related to the network construction to ones related to dynamical processes such as epidemic spreading. Our framework not only provides a natural generalization of percolation but, more importantly, offers an accurate way to test the robustness of networks in realistic scenarios.

MEDALT: single-cell copy number lineage tracing enabling gene discovery

We present a Minimal Event Distance Aneuploidy Lineage Tree (MEDALT) algorithm that infers the evolution history of a cell population based on single-cell copy number (SCCN) profiles, and a statistical routine named lineage speciation analysis (LSA), whichty facilitates discovery of fitness-associated alterations and genes from SCCN lineage trees. MEDALT appears more accurate than phylogenetics approaches in reconstructing copy number lineage. From data from 20 triple-negative breast cancer patients, our approaches effectively prioritize genes that are essential for breast cancer cell fitness and predict patient survival, including those implicating convergent evolution.The source code of our study is available at https://github.com/KChen-lab/MEDALT .

Predicting the Future Course of Opioid Overdose Mortality: An Example From Two US States

BACKGROUND

The rapid growth of opioid abuse and the related mortality across the United States has spurred the development of predictive models for the allocation of public health resources. These models should characterize heterogeneous growth across states using a drug epidemic framework that enables assessments of epidemic onset, rates of growth, and limited capacities for epidemic growth.

METHODS

We used opioid overdose mortality data for 146 North and South Carolina counties from 2001 through 2014 to compare the retrodictive and predictive performance of a logistic growth model that parameterizes onsets, growth, and carrying capacity within a traditional Bayesian Poisson space-time model.

RESULTS

In fitting the models to past data, the performance of the logistic growth model was superior to the standard Bayesian Poisson space-time model (deviance information criterion: 8,088 vs. 8,256), with reduced spatial and independent errors. Predictively, the logistic model more accurately estimated fatality rates 1, 2, and 3 years in the future (root mean squared error medians were lower for 95.7% of counties from 2012 to 2014). Capacity limits were higher in counties with greater population size, percent population age 45-64, and percent white population. Epidemic onset was associated with greater same-year and past-year incidence of overdose hospitalizations.

CONCLUSION

Growth in annual rates of opioid fatalities was capacity limited, heterogeneous across counties, and spatially correlated, requiring spatial epidemic models for the accurate and reliable prediction of future outcomes related to opioid abuse. Indicators of risk are identifiable and can be used to predict future mortality outcomes.

Model-Based Detection of Whole-Genome Duplications in a Phylogeny

Ancient whole-genome duplications (WGDs) leave signatures in comparative genomic data sets that can be harnessed to detect these events of presumed evolutionary importance. Current statistical approaches for the detection of ancient WGDs in a phylogenetic context have two main drawbacks. The first is that unwarranted restrictive assumptions on the "background" gene duplication and loss rates make inferences unreliable in the face of model violations. The second is that most methods can only be used to examine a limited set of a priori selected WGD hypotheses and cannot be used to discover WGDs in a phylogeny. In this study, we develop an approach for WGD inference using gene count data that seeks to overcome both issues. We employ a phylogenetic birth-death model that includes WGD in a flexible hierarchical Bayesian approach and use reversible-jump Markov chain Monte Carlo to perform Bayesian inference of branch-specific duplication, loss, and WGD retention rates across the space of WGD configurations. We evaluate the proposed method using simulations, apply it to data sets from flowering plants, and discuss the statistical intricacies of model-based WGD inference.

Stochastic dynamics of Francisella tularensis infection and replication

We study the pathogenesis of Francisella tularensis infection with an experimental mouse model, agent-based computation and mathematical analysis. Following inhalational exposure to Francisella tularensis SCHU S4, a small initial number of bacteria enter lung host cells and proliferate inside them, eventually destroying the host cell and releasing numerous copies that infect other cells. Our analysis of disease progression is based on a stochastic model of a population of infectious agents inside one host cell, extending the birth-and-death process by the occurrence of catastrophes: cell rupture events that affect all bacteria in a cell simultaneously. Closed expressions are obtained for the survival function of an infected cell, the number of bacteria released as a function of time after infection, and the total bacterial load. We compare our mathematical analysis with the results of agent-based computation and, making use of approximate Bayesian statistical inference, with experimental measurements carried out after murine aerosol infection with the virulent SCHU S4 strain of the bacterium Francisella tularensis, that infects alveolar macrophages. The posterior distribution of the rate of replication of intracellular bacteria is consistent with the estimate that the time between rounds of bacterial division is less than 6 hours in vivo.

Infecting a host cell is required for the replication of many types of bacteria and viruses. In some cases, infected cells release new infectious agents continuously over their lifetime. In others, such as the Francisella tularensis bacterium studied here, they are released in a single burst that coincides with the cell’s death. We show how a stochastic model, the birth-and-death process with catastrophe, can be used to characterise infection in a single cell, thereby allowing us to account for burst events and quantify the kinetics of pathogenesis in the lung, the initial site of infection, as well as in other organs that the infection spreads to. We learn about the parameters of the mathematical model of Francisella tularensis infection making use of the experimental measurements of bacterial loads, together with approximate Bayesian statistical inference methods. The most important parameter describing the pathogenesis is the rate of replication of intracellular bacteria.

Parameter estimation for discretely observed linear birth-and-death processes

Birth-and-death processes are widely used to model the development of biological populations. Although they are relatively simple models, their parameters can be challenging to estimate, as the likelihood can become numerically unstable when data arise from the most common sampling schemes, such as annual population censuses. A further difficulty arises when the discrete observations are not equi-spaced, for example, when census data are unavailable for some years. We present two approaches to estimating the birth, death, and growth rates of a discretely observed linear birth-and-death process: via an embedded Galton-Watson process and by maximizing a saddlepoint approximation to the likelihood. We study asymptotic properties of the estimators, compare them on numerical examples, and apply the methodology to data on monitored populations.

Social Cooperativity of Bacteria during Reversible Surface Attachment in Young Biofilms: a Quantitative Comparison of Pseudomonas aeruginosa PA14 and PAO1

The initial pivotal phase of bacterial biofilm formation known as reversible attachment, where cells undergo a period of transient surface attachment, is at once universal and poorly understood. What is more, although we know that reversible attachment culminates ultimately in irreversible attachment, it is not clear how reversible attachment progresses phenotypically, as bacterial surface-sensing circuits fundamentally alter cellular behavior. We analyze diverse observed bacterial behavior one family at a time (defined as a full lineage of cells related to one another by division) using a unifying stochastic model and show that our findings lead to insights on the time evolution of reversible attachment and the social cooperative dimension of surface attachment in PAO1 and PA14 strains.

What are bacteria doing during “reversible attachment,” the period of transient surface attachment when they initially engage a surface, besides attaching themselves to the surface? Can an attaching cell help any other cell attach? If so, does it help all cells or employ a more selective strategy to help either nearby cells (spatial neighbors) or its progeny (temporal neighbors)? Using community tracking methods at the single-cell resolution, we suggest answers to these questions based on how reversible attachment progresses during surface sensing for Pseudomonas aeruginosa strains PAO1 and PA14. Although PAO1 and PA14 exhibit similar trends of surface cell population increase, they show unanticipated differences when cells are considered at the lineage level and interpreted using the quantitative framework of an exactly solvable stochastic model. Reversible attachment comprises two regimes of behavior, processive and nonprocessive, corresponding to whether cells of the lineage stay on the surface long enough to divide, or not, before detaching. Stark differences between PAO1 and PA14 in the processive regime of reversible attachment suggest the existence of two surface colonization strategies. PAO1 lineages commit quickly to a surface compared to PA14 lineages, with early c-di-GMP-mediated exopolysaccharide (EPS) production that can facilitate the attachment of neighbors. PA14 lineages modulate their motility via cyclic AMP (cAMP) and retain memory of the surface so that their progeny are primed for improved subsequent surface attachment. Based on the findings of previous studies, we propose that the differences between PAO1 and PA14 are potentially rooted in downstream differences between Wsp-based and Pil-Chp-based surface-sensing systems, respectively.

Angiosperm-Wide and Family-Level Analyses of AP2/ERF Genes Reveal Differential Retention and Sequence Divergence After Whole-Genome Duplication

Plants are immobile and often face stressful environmental conditions, prompting the evolution of genes regulating environmental responses. Such evolution is achieved largely through gene duplication and subsequent divergence. One of the most important gene families involved in regulating plant environmental responses and development is the AP2/ERF superfamily; however, the evolutionary history of these genes is unclear across angiosperms and in major angiosperm families adapted to various ecological niches. Specifically, the impact on gene copy number of whole-genome duplication events occurring around the time of the origins of several plant families is unknown. Here, we present the first angiosperm-wide comparative study of AP2/ERF genes, identifying 75 Angiosperm OrthoGroups (AOGs), each derived from an ancestral angiosperm gene copy. Among these AOGs, 21 retain duplicates with increased copy number in many angiosperm lineages, while the remaining 54 AOGs tend to maintain low copy number. Further analyses of multiple species in the Brassicaceae family indicated that family-specific duplicates experienced differential selective pressures in coding regions, with some paralogs showing signs of positive selection. Further, cis regulatory elements also exhibit extensive divergence between duplicates in Arabidopsis. Moreover, comparison of expression levels suggested that AP2/ERF genes with frequently retained duplicates are enriched for broad expression patterns, offering increased opportunities for functional diversification via changes in expression patterns, and providing a mechanism for repeated duplicate retention in some AOGs. Our results represent the most comprehensive evolutionary history of the AP2/ERF gene family, and support the hypothesis that AP2/ERF genes with broader expression patterns are more likely to be retained as duplicates than those with narrower expression profiles, which could lead to a higher chance of duplicate gene subfunctionalization. The greater tendency of some AOGs to retain duplicates, allowing expression and functional divergence, may facilitate the evolution of complex signaling networks in response to new environmental conditions.

Grammar of protein domain architectures

Genomes appear similar to natural language texts, and protein domains can be treated as analogs of words. To investigate the linguistic properties of genomes further, we calculated the complexity of the “protein languages” in all major branches of life and identified a nearly universal value of information gain associated with the transition from a random domain arrangement to the current protein domain architecture. An exploration of the evolutionary relationship of the protein languages identified the domain combinations that discriminate between the major branches of cellular life. We conclude that there exists a “quasi-universal grammar” of protein domains and that the nearly constant information gain we identified corresponds to the minimal complexity required to maintain a functional cell.

From an abstract, informational perspective, protein domains appear analogous to words in natural languages in which the rules of word association are dictated by linguistic rules, or grammar. Such rules exist for protein domains as well, because only a small fraction of all possible domain combinations is viable in evolution. We employ a popular linguistic technique, n-gram analysis, to probe the “proteome grammar”—that is, the rules of association of domains that generate various domain architectures of proteins. Comparison of the complexity measures of “protein languages” in major branches of life shows that the relative entropy difference (information gain) between the observed domain architectures and random domain combinations is highly conserved in evolution and is close to being a universal constant, at ∼1.2 bits. Substantial deviations from this constant are observed in only two major groups of organisms: a subset of Archaea that appears to be cells simplified to the limit, and animals that display extreme complexity. We also identify the n-grams that represent signatures of the major branches of cellular life. The results of this analysis bolster the analogy between genomes and natural language and show that a “quasi-universal grammar” underlies the evolution of domain architectures in all divisions of cellular life. The nearly universal value of information gain by the domain architectures could reflect the minimum complexity of signal processing that is required to maintain a functioning cell.

Environmental fitness heterogeneity in the Moran process

Many mathematical models of evolution assume that all individuals experience the same environment. Here, we study the Moran process in heterogeneous environments. The population is of finite size with two competing types, which are exposed to a fixed number of environmental conditions. Reproductive rate is determined by both the type and the environment. We first calculate the condition for selection to favour the mutant relative to the resident wild-type. In large populations, the mutant is favoured if and only if the mutant's spatial average reproductive rate exceeds that of the resident. But environmental heterogeneity elucidates an interesting asymmetry between the mutant and the resident. Specifically, mutant heterogeneity suppresses its fixation probability; if this heterogeneity is strong enough, it can even completely offset the effects of selection (including in large populations). By contrast, resident heterogeneity has no effect on a mutant's fixation probability in large populations and can amplify it in small populations.

Hidden heterogeneity and circadian-controlled cell fate inferred from single cell lineages

The origin of lineage correlations among single cells and the extent of heterogeneity in their intermitotic times (IMT) and apoptosis times (AT) remain incompletely understood. Here we developed single cell lineage-tracking experiments and computational algorithms to uncover correlations and heterogeneity in the IMT and AT of a colon cancer cell line before and during cisplatin treatment. These correlations could not be explained using simple protein production/degradation models. Sister cell fates were similar regardless of whether they divided before or after cisplatin administration and did not arise from proximity-related factors, suggesting fate determination early in a cell's lifetime. Based on these findings, we developed a theoretical model explaining how the observed correlation structure can arise from oscillatory mechanisms underlying cell fate control. Our model recapitulated the data only with very specific oscillation periods that fit measured circadian rhythms, thereby suggesting an important role of the circadian clock in controlling cellular fates.

Incorporating extrinsic noise into the stochastic simulation of biochemical reactions: A comparison of approaches

The stochastic simulation algorithm (SSA) has been widely used for simulating biochemical reaction networks. SSA is able to capture the inherently intrinsic noise of the biological system, which is due to the discreteness of species population and to the randomness of their reciprocal interactions. However, SSA does not consider other sources of heterogeneity in biochemical reaction systems, which are referred to as extrinsic noise. Here, we extend two simulation approaches, namely, the integration-based method and the rejection-based method, to take extrinsic noise into account by allowing the reaction propensities to vary in time and state dependent manner. For both methods, new efficient implementations are introduced and their efficiency and applicability to biological models are investigated. Our numerical results suggest that the rejection-based method performs better than the integration-based method when the extrinsic noise is considered.

Duration of a minor epidemic

Disease outbreaks in stochastic SIR epidemic models are characterized as either minor or major. When ℛ0<1, all epidemics are minor, whereas if ℛ0>1, they can be minor or major. In 1955, Whittle derived formulas for the probability of a minor or a major epidemic. A minor epidemic is distinguished from a major one in that a minor epidemic is generally of shorter duration and has substantially fewer cases than a major epidemic. In this investigation, analytical formulas are derived that approximate the probability density, the mean, and the higher-order moments for the duration of a minor epidemic. These analytical results are applicable to minor epidemics in stochastic SIR, SIS, and SIRS models with a single infected class. The probability density for minor epidemics in more complex epidemic models can be computed numerically applying multitype branching processes and the backward Kolmogorov differential equations. When ℛ0 is close to one, minor epidemics are more common than major epidemics and their duration is significantly longer than when ℛ0≪1 or ℛ0≫1.

Antibiotic-induced population fluctuations and stochastic clearance of bacteria

Effective antibiotic use that minimizes treatment failures remains a challenge. A better understanding of how bacterial populations respond to antibiotics is necessary. Previous studies of large bacterial populations established the deterministic framework of pharmacodynamics. Here, characterizing the dynamics of population extinction, we demonstrated the stochastic nature of eradicating bacteria with antibiotics. Antibiotics known to kill bacteria (bactericidal) induced population fluctuations. Thus, at high antibiotic concentrations, the dynamics of bacterial clearance were heterogeneous. At low concentrations, clearance still occurred with a non-zero probability. These striking outcomes of population fluctuations were well captured by our probabilistic model. Our model further suggested a strategy to facilitate eradication by increasing extinction probability. We experimentally tested this prediction for antibiotic-susceptible and clinically-isolated resistant bacteria. This new knowledge exposes fundamental limits in our ability to predict bacterial eradication. Additionally, it demonstrates the potential of using antibiotic concentrations that were previously deemed inefficacious to eradicate bacteria.

Computational methods for birth-death processes

Many important stochastic counting models can be written as general birth-death processes (BDPs). BDPs are continuous-time Markov chains on the non-negative integers in which only jumps to adjacent states are allowed. BDPs can be used to easily parameterize a rich variety of probability distributions on the non-negative integers, and straightforward conditions guarantee that these distributions are proper. BDPs also provide a mechanistic interpretation - birth and death of actual particles or organisms - that has proven useful in evolution, ecology, physics, and chemistry. Although the theoretical properties of general BDPs are well understood, traditionally statistical work on BDPs has been limited to the simple linear (Kendall) process. Aside from a few simple cases, it remains impossible to find analytic expressions for the likelihood of a discretely-observed BDP, and computational difficulties have hindered development of tools for statistical inference. But the gap between BDP theory and practical methods for estimation has narrowed in recent years. There are now robust methods for evaluating likelihoods for realizations of BDPs: finite-time transition, first passage, equilibrium probabilities, and distributions of summary statistics that arise commonly in applications. Recent work has also exploited the connection between continuously- and discretely-observed BDPs to derive EM algorithms for maximum likelihood estimation. Likelihood-based inference for previously intractable BDPs is much easier than previously thought and regression approaches analogous to Poisson regression are straightforward to derive. In this review, we outline the basic mathematical theory for BDPs and demonstrate new tools for statistical inference using data from BDPs.

Birth/birth-death processes and their computable transition probabilities with biological applications

Birth-death processes track the size of a univariate population, but many biological systems involve interaction between populations, necessitating models for two or more populations simultaneously. A lack of efficient methods for evaluating finite-time transition probabilities of bivariate processes, however, has restricted statistical inference in these models. Researchers rely on computationally expensive methods such as matrix exponentiation or Monte Carlo approximation, restricting likelihood-based inference to small systems, or indirect methods such as approximate Bayesian computation. In this paper, we introduce the birth/birth-death process, a tractable bivariate extension of the birth-death process, where rates are allowed to be nonlinear. We develop an efficient algorithm to calculate its transition probabilities using a continued fraction representation of their Laplace transforms. Next, we identify several exemplary models arising in molecular epidemiology, macro-parasite evolution, and infectious disease modeling that fall within this class, and demonstrate advantages of our proposed method over existing approaches to inference in these models. Notably, the ubiquitous stochastic susceptible-infectious-removed (SIR) model falls within this class, and we emphasize that computable transition probabilities newly enable direct inference of parameters in the SIR model. We also propose a very fast method for approximating the transition probabilities under the SIR model via a novel branching process simplification, and compare it to the continued fraction representation method with application to the 17th century plague in Eyam. Although the two methods produce similar maximum a posteriori estimates, the branching process approximation fails to capture the correlation structure in the joint posterior distribution.

Stochastic Model of Maturation and Vesicular Exchange in Cellular Organelles

The dynamical organization of membrane-bound organelles along intracellular transport pathways relies on vesicular exchange between organelles and on the maturation of the organelle's composition by enzymatic reactions or exchange with the cytoplasm. The relative importance of each mechanism in controlling organelle dynamics remains controversial, in particular for transport through the Golgi apparatus. Using a stochastic model, we identify two classes of dynamical behavior that can lead to full maturation of membrane-bound compartments. In the first class, maturation corresponds to the stochastic escape from a steady state in which export is dominated by vesicular exchange, and is very unlikely for large compartments. In the second class, it occurs in a quasi-deterministic fashion and is almost size independent. Whether a system belongs to the first or second class is largely controlled by homotypic fusion.

Dynamic patterns of information flow in complex networks

Although networks are extensively used to visualize information flow in biological, social and technological systems, translating topology into dynamic flow continues to challenge us, as similar networks exhibit fundamentally different flow patterns, driven by different interaction mechanisms. To uncover a network’s actual flow patterns, here we use a perturbative formalism, analytically tracking the contribution of all nodes/paths to the flow of information, exposing the rules that link structure and dynamic information flow for a broad range of nonlinear systems. We find that the diversity of flow patterns can be mapped into a single universal function, characterizing the interplay between the system’s topology and its dynamics, ultimately allowing us to identify the network’s main arteries of information flow. Counter-intuitively, our formalism predicts a family of frequently encountered dynamics where the flow of information avoids the hubs, favoring the network’s peripheral pathways, a striking disparity between structure and dynamics.

Complex networks are a useful tool to investigate spreading processes but topology alone is insufficient to predict information flow. Here the authors propose a measure of information flow and predict its behavior from the interplay between structure and dynamics.

Modelling the evolution of transcription factor binding preferences in complex eukaryotes

Transcription factors (TFs) exert their regulatory action by binding to DNA with specific sequence preferences. However, different TFs can partially share their binding sequences due to their common evolutionary origin. This "redundancy" of binding defines a way of organizing TFs in "motif families" by grouping TFs with similar binding preferences. Since these ultimately define the TF target genes, the motif family organization entails information about the structure of transcriptional regulation as it has been shaped by evolution. Focusing on the human TF repertoire, we show that a one-parameter evolutionary model of the Birth-Death-Innovation type can explain the TF empirical repartition in motif families, and allows to highlight the relevant evolutionary forces at the origin of this organization. Moreover, the model allows to pinpoint few deviations from the neutral scenario it assumes: three over-expanded families (including HOX and FOX genes), a set of "singleton" TFs for which duplication seems to be selected against, and a higher-than-average rate of diversification of the binding preferences of TFs with a Zinc Finger DNA binding domain. Finally, a comparison of the TF motif family organization in different eukaryotic species suggests an increase of redundancy of binding with organism complexity.

Structure-based control of complex networks with nonlinear dynamics

What can we learn about controlling a system solely from its underlying network structure? Here we adapt a recently developed framework for control of networks governed by a broad class of nonlinear dynamics that includes the major dynamic models of biological, technological, and social processes. This feedback-based framework provides realizable node overrides that steer a system toward any of its natural long-term dynamic behaviors, regardless of the specific functional forms and system parameters. We use this framework on several real networks, identify the topological characteristics that underlie the predicted node overrides, and compare its predictions to those of structural controllability in control theory. Finally, we demonstrate this framework's applicability in dynamic models of gene regulatory networks and identify nodes whose override is necessary for control in the general case but not in specific model instances.

Disentangling the effects of selection and loss bias on gene dynamics

We combine mathematical modeling of genome evolution with comparative analysis of prokaryotic genomes to estimate the relative contributions of selection and intrinsic loss bias to the evolution of different functional classes of genes and mobile genetic elements (MGE). An exact solution for the dynamics of gene family size was obtained under a linear duplication-transfer-loss model with selection. With the exception of genes involved in information processing, particularly translation, which are maintained by strong selection, the average selection coefficient for most nonparasitic genes is low albeit positive, compatible with observed positive correlation between genome size and effective population size. Free-living microbes evolve under stronger selection for gene retention than parasites. Different classes of MGE show a broad range of fitness effects, from the nearly neutral transposons to prophages, which are actively eliminated by selection. Genes involved in antiparasite defense, on average, incur a fitness cost to the host that is at least as high as the cost of plasmids. This cost is probably due to the adverse effects of autoimmunity and curtailment of horizontal gene transfer caused by the defense systems and selfish behavior of some of these systems, such as toxin-antitoxin and restriction modification modules. Transposons follow a biphasic dynamics, with bursts of gene proliferation followed by decay in the copy number that is quantitatively captured by the model. The horizontal gene transfer to loss ratio, but not duplication to loss ratio, correlates with genome size, potentially explaining increased abundance of neutral and costly elements in larger genomes.

Social interactions among grazing reef fish drive material flux in a coral reef ecosystem

In human financial and social systems, exchanges of information among individuals cause speculative bubbles, behavioral cascades, and other correlated actions that profoundly influence system-level function. Exchanges of information are also widespread in ecological systems, but their effects on ecosystem-level processes are largely unknown. Herbivory is a critical ecological process in coral reefs, where diverse assemblages of fish maintain reef health by controlling the abundance of algae. Here, we show that social interactions have a major effect on fish grazing rates in a reef ecosystem. We combined a system for observing and manipulating large foraging areas in a coral reef with a class of dynamical decision-making models to reveal that reef fish use information about the density and actions of nearby fish to decide when to feed on algae and when to flee foraging areas. This "behavioral coupling" causes bursts of feeding activity that account for up to 68% of the fish community's consumption of algae. Moreover, correlations in fish behavior induce a feedback, whereby each fish spends less time feeding when fewer fish are present, suggesting that reducing fish stocks may not only reduce total algal consumption but could decrease the amount of algae each remaining fish consumes. Our results demonstrate that social interactions among consumers can have a dominant effect on the flux of energy and materials through ecosystems, and our methodology paves the way for rigorous in situ measurements of the behavioral rules that underlie ecological rates in other natural systems.

Mathematical Modeling of Avidity Distribution and Estimating General Binding Properties of Transcription Factors from Genome-Wide Binding Profiles

The shape of the experimental frequency distributions (EFD) of diverse molecular interaction events quantifying genome-wide binding is often skewed to the rare but abundant quantities. Such distributions are systematically deviated from standard power-law functions proposed by scale-free network models suggesting that more explanatory and predictive probabilistic model(s) are needed. Identification of the mechanism-based data-driven statistical distributions that provide an estimation and prediction of binding properties of transcription factors from genome-wide binding profiles is the goal of this analytical survey. Here, we review and develop an analytical framework for modeling, analysis, and prediction of transcription factor (TF) DNA binding properties detected at the genome scale. We introduce a mixture probabilistic model of binding avidity function that includes nonspecific and specific binding events. A method for decomposition of specific and nonspecific TF-DNA binding events is proposed. We show that the Kolmogorov-Waring (KW) probability function (PF), modeling the steady state TF binding-dissociation stochastic process, fits well with the EFD for diverse TF-DNA binding datasets. Furthermore, this distribution predicts total number of TF-DNA binding sites (BSs), estimating specificity and sensitivity as well as other basic statistical features of DNA-TF binding when the experimental datasets are noise-rich and essentially incomplete. The KW distribution fits equally well to TF-DNA binding activity for different TFs including ERE, CREB, STAT1, Nanog, and Oct4. Our analysis reveals that the KW distribution and its generalized form provides the family of power-law-like distributions given in terms of hypergeometric series functions, including standard and generalized Pareto and Waring distributions, providing flexible and common skewed forms of the transcription factor binding site (TFBS) avidity distribution function. We suggest that the skewed binding events may be due to a wide range of evolutionary processes of creating weak avidity TFBS associated with random mutations, while the rare high-avidity binding sites (i.e., high-avidity evolutionarily conserved canonical e-boxes) rarely occurred. These, however, may be positively selected in microevolution.

Evolution by gene loss

The recent increase in genomic data is revealing an unexpected perspective of gene loss as a pervasive source of genetic variation that can cause adaptive phenotypic diversity. This novel perspective of gene loss is raising new fundamental questions. How relevant has gene loss been in the divergence of phyla? How do genes change from being essential to dispensable and finally to being lost? Is gene loss mostly neutral, or can it be an effective way of adaptation? These questions are addressed, and insights are discussed from genomic studies of gene loss in populations and their relevance in evolutionary biology and biomedicine.

Coupling bounds for approximating birth-death processes by truncation

Birth-death processes are continuous-time Markov counting processes. Approximate moments can be computed by truncating the transition rate matrix. Using a coupling argument, we derive bounds for the total variation distance between the process and its finite approximation.

Closed-form stochastic solutions for non-equilibrium dynamics and inheritance of cellular components over many cell divisions

Stochastic dynamics govern many important processes in cellular biology, and an underlying theoretical approach describing these dynamics is desirable to address a wealth of questions in biology and medicine. Mathematical tools exist for treating several important examples of these stochastic processes, most notably gene expression and random partitioning at single-cell divisions or after a steady state has been reached. Comparatively little work exists exploring different and specific ways that repeated cell divisions can lead to stochastic inheritance of unequilibrated cellular populations. Here we introduce a mathematical formalism to describe cellular agents that are subject to random creation, replication and/or degradation, and are inherited according to a range of random dynamics at cell divisions. We obtain closed-form generating functions describing systems at any time after any number of cell divisions for binomial partitioning and divisions provoking a deterministic or random, subtractive or additive change in copy number, and show that these solutions agree exactly with stochastic simulation. We apply this general formalism to several example problems involving the dynamics of mitochondrial DNA during development and organismal lifetimes.

On Binomial Observations of Continuous-Time Markovian Population Models

In this paper we consider a class of stochastic processes based on binomial observations of continuous-time, Markovian population models. We derive the conditional probability mass function of the next binomial observation given a set of binomial observations. For this purpose, we first find the conditional probability mass function of the underlying continuous-time Markovian population model, given a set of binomial observations, by exploiting a conditional Bayes' theorem from filtering, and then use the law of total probability to find the former. This result paves the way for further study of the stochastic process introduced by the binomial observations. We utilize our results to show that binomial observations of the simple birth process are non-Markovian.

On Binomial Observations of Continuous-Time Markovian Population Models

In this paper we consider a class of stochastic processes based on binomial observations of continuous-time, Markovian population models. We derive the conditional probability mass function of the next binomial observation given a set of binomial observations. For this purpose, we first find the conditional probability mass function of the underlying continuous-time Markovian population model, given a set of binomial observations, by exploiting a conditional Bayes' theorem from filtering, and then use the law of total probability to find the former. This result paves the way for further study of the stochastic process introduced by the binomial observations. We utilize our results to show that binomial observations of the simple birth process are non-Markovian.

SEX, LIES AND SELF-REPORTED COUNTS: BAYESIAN MIXTURE MODELS FOR HEAPING IN LONGITUDINAL COUNT DATA VIA BIRTH-DEATH PROCESSES

Surveys often ask respondents to report non-negative counts, but respondents may misremember or round to a nearby multiple of 5 or 10. This phenomenon is called heaping, and the error inherent in heaped self-reported numbers can bias estimation. Heaped data may be collected cross-sectionally or longitudinally and there may be covariates that complicate the inferential task. Heaping is a well-known issue in many survey settings, and inference for heaped data is an important statistical problem. We propose a novel reporting distribution whose underlying parameters are readily interpretable as rates of misremembering and rounding. The process accommodates a variety of heaping grids and allows for quasi-heaping to values nearly but not equal to heaping multiples. We present a Bayesian hierarchical model for longitudinal samples with covariates to infer both the unobserved true distribution of counts and the parameters that control the heaping process. Finally, we apply our methods to longitudinal self-reported counts of sex partners in a study of high-risk behavior in HIV-positive youth.

Sustainable deployment of QTLs conferring quantitative resistance to crops: first lessons from a stochastic model

Quantitative plant disease resistance is believed to be more durable than qualitative resistance, since it exerts less selective pressure on the pathogens. However, the process of progressive pathogen adaptation to quantitative resistance is poorly understood, which makes it difficult to predict its durability or to derive principles for its sustainable deployment. Here, we study the dynamics of pathogen adaptation in response to quantitative plant resistance affecting pathogen reproduction rate and its colonizing capacity. We developed a stochastic model for the continuous evolution of a pathogen population within a quantitatively resistant host. We assumed that pathogen can adapt to a host by the progressive restoration of reproduction rate or of colonizing capacity, or of both. Our model suggests that a combination of quantitative trait loci (QTLs) affecting distinct pathogen traits was more durable if the evolution of repressed traits was antagonistic. Otherwise, quantitative resistance that depressed only pathogen reproduction was more durable. In order to decelerate the progressive pathogen adaptation, QTLs that decrease the pathogen's maximum capacity to colonize must be combined with QTLs that decrease the spore production per lesion or the infection efficiency or that increase the latent period. Our theoretical framework can help breeders to develop principles for sustainable deployment of QTLs.

Statistical properties of pairwise distances between leaves on a random Yule tree

A Yule tree is the result of a branching process with constant birth and death rates. Such a process serves as an instructive null model of many empirical systems, for instance, the evolution of species leading to a phylogenetic tree. However, often in phylogeny the only available information is the pairwise distances between a small fraction of extant species representing the leaves of the tree. In this article we study statistical properties of the pairwise distances in a Yule tree. Using a method based on a recursion, we derive an exact, analytic and compact formula for the expected number of pairs separated by a certain time distance. This number turns out to follow a increasing exponential function. This property of a Yule tree can serve as a simple test for empirical data to be well described by a Yule process. We further use this recursive method to calculate the expected number of the n-most closely related pairs of leaves and the number of cherries separated by a certain time distance. To make our results more useful for realistic scenarios, we explicitly take into account that the leaves of a tree may be incompletely sampled and derive a criterion for poorly sampled phylogenies. We show that our result can account for empirical data, using two families of birds species.

Inference of epidemiological dynamics based on simulated phylogenies using birth-death and coalescent models

Quantifying epidemiological dynamics is crucial for understanding and forecasting the spread of an epidemic. The coalescent and the birth-death model are used interchangeably to infer epidemiological parameters from the genealogical relationships of the pathogen population under study, which in turn are inferred from the pathogen genetic sequencing data. To compare the performance of these widely applied models, we performed a simulation study. We simulated phylogenetic trees under the constant rate birth-death model and the coalescent model with a deterministic exponentially growing infected population. For each tree, we re-estimated the epidemiological parameters using both a birth-death and a coalescent based method, implemented as an MCMC procedure in BEAST v2.0. In our analyses that estimate the growth rate of an epidemic based on simulated birth-death trees, the point estimates such as the maximum a posteriori/maximum likelihood estimates are not very different. However, the estimates of uncertainty are very different. The birth-death model had a higher coverage than the coalescent model, i.e. contained the true value in the highest posterior density (HPD) interval more often (2-13% vs. 31-75% error). The coverage of the coalescent decreases with decreasing basic reproductive ratio and increasing sampling probability of infecteds. We hypothesize that the biases in the coalescent are due to the assumption of deterministic rather than stochastic population size changes. Both methods performed reasonably well when analyzing trees simulated under the coalescent. The methods can also identify other key epidemiological parameters as long as one of the parameters is fixed to its true value. In summary, when using genetic data to estimate epidemic dynamics, our results suggest that the birth-death method will be less sensitive to population fluctuations of early outbreaks than the coalescent method that assumes a deterministic exponentially growing infected population.