The stochastic modelling of kleptoparasitism using a Markov process

Kleptoparasitism, the stealing of food items from other animals, is a common behaviour observed across a huge variety of species, and has been subjected to significant modelling effort. Most such modelling has been deterministic, effectively assuming an infinite population, although recently some important stochastic models have been developed. In particular the model of Yates and Broom (Stochastic models of kleptoparasitism. J. Theor. Biol. 248 (2007), 480-489) introduced a stochastic version following the original model of Ruxton and Moody (The ideal free distribution with kleptoparasitism. J. Theor. Biol. 186 (1997), 449-458), and whilst they generated results of interest, they did not solve the model explicitly. In this paper, building on methods used already by van der Meer and Smallegange (A stochastic version of the Beddington-DeAngelis functional response: Modelling interference for a finite number of predators. J. Animal Ecol. 78 (2009) 134-142) we give an exact solution to the distribution of the population over the states for the Yates and Broom model and investigate the effects of some key biological parameters, especially for small populations where stochastic models can be expected to differ most from their deterministic equivalents.     

Mixture distributions in a stochastic gene expression model with delayed feedback: a WKB approximation approach

Journal of Mathematical Biology - Noise in gene expression can be substantively affected by the presence of production delay. Here we consider a mathematical model with bursty production of...     

Bayesian phylogenetic model selection using reversible jump Markov chain Monte Carlo

A common problem in molecular phylogenetics is choosing a model of DNA substitution that does a good job of explaining the DNA sequence alignment without introducing superfluous parameters. A number of methods have been used to choose among a small set of candidate substitution models, such as the likelihood ratio test, the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and Bayes factors. Current implementations of any of these criteria suffer from the limitation that only a small set of models are examined, or that the test does not allow easy comparison of non-nested models. In this article, we expand the pool of candidate substitution models to include all possible time-reversible models. This set includes seven models that have already been described. We show how Bayes factors can be calculated for these models using reversible jump Markov chain Monte Carlo, and apply the method to 16 DNA sequence alignments. For each data set, we compare the model with the best Bayes factor to the best models chosen using AIC and BIC. We find that the best model under any of these criteria is not necessarily the most complicated one; models with an intermediate number of substitution types typically do best. Moreover, almost all of the models that are chosen as best do not constrain a transition rate to be the same as a transversion rate, suggesting that it is the transition/transversion rate bias that plays the largest role in determining which models are selected. Importantly, the reversible jump Markov chain Monte Carlo algorithm described here allows estimation of phylogeny (and other phylogenetic model parameters) to be performed while accounting for uncertainty in the model of DNA substitution.     

Reversible jump Markov chain Monte Carlo computation and Bayesian model determination

Markov chain Monte Carlo methods for Bayesian computation haveuntil recently been restricted to problems where the joint distributionof all variables has a density with respect to some fixed standardunderlying measure. They have therefore not been available forapplication to Bayesian model determination, where the dimensionalityof the parameter vector is typically not fixed. This paper proposesa new framework for the construction of reversible Markov chainsamplers that jump between parameter subspaces of differingdimensionality, which is flexible and entirely constructive.It should therefore have wide applicability in model determinationproblems. The methodology is illustrated with applications tomultiple change-point analysis in one and two dimensions, andto a Bayesian comparison of binomial experiments.     

Hierarchical Dirichlet process model for gene expression clustering

Clustering is an important data processing tool for interpreting microarray data and genomic network inference. In this article, we propose a clustering algorithm based on the hierarchical Dirichlet processes (HDP). The HDP clustering introduces a hierarchical structure in the statistical model which captures the hierarchical features prevalent in biological data such as the gene express data. We develop a Gibbs sampling algorithm based on the Chinese restaurant metaphor for the HDP clustering. We apply the proposed HDP algorithm to both regulatory network segmentation and gene expression clustering. The HDP algorithm is shown to outperform several popular clustering algorithms by revealing the underlying hierarchical structure of the data. For the yeast cell cycle data, we compare the HDP result to the standard result and show that the HDP algorithm provides more information and reduces the unnecessary clustering fragments.     

A stochastic model of gene evolution with chaotic mutations

We develop here a new class of stochastic models of gene evolution in which the mutations are chaotic, i.e. a random subset of the 64 possible trinucleotides mutates at each evolutionary time t according to some substitution probabilities. Therefore, at each time t, the numbers and the types of mutable trinucleotides are unknown. Thus, the mutation matrix changes at each time t. The chaotic model developed generalizes the standard model in which all the trinucleotides mutate at each time t. It determines the occurrence probabilities at time t of trinucleotides which chaotically mutate according to three substitution parameters associated with the three trinucleotide sites. Two theorems prove that this chaotic model has a probability vector at each time t and that it converges to a uniform probability vector identical to that of the standard model. Furthermore, four applications of this chaotic model (with a uniform random strategy for the 64 trinucleotides and with a particular strategy for the three stop codons) allow an evolutionary study of the three circular codes identified in both eukaryotic and prokaryotic genes. A circular code is a particular set of trinucleotides whose main property is the retrieval of the frames in genes locally, i.e. anywhere in genes and particularly without start codons, and automatically with a window of a few nucleotides. After a certain evolutionary time and with particular values for the three substitution parameters, the chaotic models retrieve the main statistical properties of the three circular codes observed in genes. These applications also allow an evolutionary comparison between the standard and chaotic models.     

A stochastic gene evolution model with time dependent mutations

We develop here a new class of gene evolution models in which the nucleotide mutations are time dependent. These models allow to study nonlinear gene evolution by accelerating or decelerating the mutation rates at different evolutionary times. They generalize the previous ones which are based on constant mutation rates. The stochastic model developed in this class determines at some time t the occurrence probabilities of trinucleotides mutating according to 3 time dependent substitution parameters associated with the 3 trinucleotide sites. Therefore, it allows to simulate the evolution of the circular code recently observed in genes. By varying the class of function for the substitution parameters, 1 among 12 models retrieves after mutation the statistical properties of the observed circular code in the 3 frames of actual genes. In this model, the mutation rate in the 3rd trinucleotide site increases during gene evolution while the mutation rates in the 1st and 2nd sites decrease. This property agrees with the actual degeneracy of the genetic code. This approach can easily be generalized to study evolution of motifs of various lengths, e.g., dicodons, etc., with time dependent mutations.     

A simple stochastic gene substitution model

If the fitnesses of n haploid alleles in a finite population are assigned at random and if the alleles can mutate to one another, and if the population is initially fixed for the kth most fit allele, then the mean number of substitutions that will occur before the most fit allele is fixed is shown to be

$\text{1}\text{2}\text{+}\text{1}\text{k}\text{+}\text{∑}\text{i=2}\text{k?1}\text{(i+3)}\text{(2i(i+1))}$

when selection is strong and mutation is weak. This result is independent of the parameters that went into the model. The result is used to provide a partial explanation for the large variance observed in the rates of molecular evolution.     

A spatiotemporal stochastic process model for species spread

We use a spatiotemporal Markov process to model the spread of an ecological population through its environment over time. Available habitat is divided into sites, and a parametric function of spatial variables is used to model the probability that one site is colonized from another. This allows us both to make predictions about the future spread of a population, and to determine which are the important factors governing colonizations. The model evolves in discrete time, allowing the population distribution to change seasonally in accordance with breeding patterns. Discrete time formulations are natural for ecological populations, but are problematic due to difficulties of fitting and predicting over irregular time intervals. The model described here can accommodate years of missing data and can therefore fit and predict at irregular intervals. Two methods of approximating the likelihood are described and applied to ornithological survey data for the woodlark, Lullula arborea, from Thetford Forest in the U.K.     

Colored extrinsic fluctuations and stochastic gene expression

Stochasticity is both exploited and controlled by cells. Although the intrinsic stochasticity inherent in biochemistry is relatively well understood, cellular variation, or ‘noise’, is predominantly generated by interactions of the system of interest with other stochastic systems in the cell or its environment. Such extrinsic fluctuations are nonspecific, affecting many system components, and have a substantial lifetime, comparable to the cell cycle (they are ‘colored’). Here, we extend the standard stochastic simulation algorithm to include extrinsic fluctuations. We show that these fluctuations affect mean protein numbers and intrinsic noise, can speed up typical network response times, and can explain trends in high‐throughput measurements of variation. If extrinsic fluctuations in two components of the network are correlated, they may combine constructively (amplifying each other) or destructively (attenuating each other). Consequently, we predict that incoherent feedforward loops attenuate stochasticity, while coherent feedforwards amplify it. Our results demonstrate that both the timescales of extrinsic fluctuations and their nonspecificity substantially affect the function and performance of biochemical networks.     

Dimension reduction strategies for analyzing global gene expression data with a response

The analysis of global gene expression data from microarrays is breaking new ground in genetics research, while confronting modelers and statisticians with many critical issues. In this paper, we consider data sets in which a categorical or continuous response is recorded, along with gene expression, on a given number of experimental samples. Data of this type are usually employed to create a prediction mechanism for the response based on gene expression, and to identify a subset of relevant genes. This defines a regression setting characterized by a dramatic under-resolution with respect to the predictors (genes), whose number exceeds by orders of magnitude the number of available observations (samples). We present a dimension reduction strategy that, under appropriate assumptions, allows us to restrict attention to a few linear combinations of the original expression profiles, and thus to overcome under-resolution. These linear combinations can then be used to build and validate a regression model with standard techniques. Moreover, they can be used to rank original predictors, and ultimately to select a subset of them through comparison with a background 'chance scenario' based on a number of independent randomizations. We apply this strategy to publicly available data on leukemia classification.     

A stochastic Markov chain model to describe lung cancer growth and metastasis

A stochastic Markov chain model for metastatic progression is developed for primary lung cancer based on a network construction of metastatic sites with dynamics modeled as an ensemble of random walkers on the network. We calculate a transition matrix, with entries (transition probabilities) interpreted as random variables, and use it to construct a circular bi-directional network of primary and metastatic locations based on postmortem tissue analysis of 3827 autopsies on untreated patients documenting all primary tumor locations and metastatic sites from this population. The resulting 50 potential metastatic sites are connected by directed edges with distributed weightings, where the site connections and weightings are obtained by calculating the entries of an ensemble of transition matrices so that the steady-state distribution obtained from the long-time limit of the Markov chain dynamical system corresponds to the ensemble metastatic distribution obtained from the autopsy data set. We condition our search for a transition matrix on an initial distribution of metastatic tumors obtained from the data set. Through an iterative numerical search procedure, we adjust the entries of a sequence of approximations until a transition matrix with the correct steady-state is found (up to a numerical threshold). Since this constrained linear optimization problem is underdetermined, we characterize the statistical variance of the ensemble of transition matrices calculated using the means and variances of their singular value distributions as a diagnostic tool. We interpret the ensemble averaged transition probabilities as (approximately) normally distributed random variables. The model allows us to simulate and quantify disease progression pathways and timescales of progression from the lung position to other sites and we highlight several key findings based on the model.