Neural dynamics as sampling: a model for stochastic computation in recurrent networks of spiking neurons

The organization of computations in networks of spiking neurons in the brain is still largely unknown, in particular in view of the inherently stochastic features of their firing activity and the experimentally observed trial-to-trial variability of neural systems in the brain. In principle there exists a powerful computational framework for stochastic computations, probabilistic inference by sampling, which can explain a large number of macroscopic experimental data in neuroscience and cognitive science. But it has turned out to be surprisingly difficult to create a link between these abstract models for stochastic computations and more detailed models of the dynamics of networks of spiking neurons. Here we create such a link and show that under some conditions the stochastic firing activity of networks of spiking neurons can be interpreted as probabilistic inference via Markov chain Monte Carlo (MCMC) sampling. Since common methods for MCMC sampling in distributed systems, such as Gibbs sampling, are inconsistent with the dynamics of spiking neurons, we introduce a different approach based on non-reversible Markov chains that is able to reflect inherent temporal processes of spiking neuronal activity through a suitable choice of random variables. We propose a neural network model and show by a rigorous theoretical analysis that its neural activity implements MCMC sampling of a given distribution, both for the case of discrete and continuous time. This provides a step towards closing the gap between abstract functional models of cortical computation and more detailed models of networks of spiking neurons.     

Smoothing spline ANOVA frailty model for recurrent event data

Gap time hazard estimation is of particular interest in recurrent event data. This article proposes a fully nonparametric approach for estimating the gap time hazard. Smoothing spline analysis of variance (ANOVA) decompositions are used to model the log gap time hazard as a joint function of gap time and covariates, and general frailty is introduced to account for between-subject heterogeneity and within-subject correlation. We estimate the nonparametric gap time hazard function and parameters in the frailty distribution using a combination of the Newton-Raphson procedure, the stochastic approximation algorithm (SAA), and the Markov chain Monte Carlo (MCMC) method. The convergence of the algorithm is guaranteed by decreasing the step size of parameter update and/or increasing the MCMC sample size along iterations. Model selection procedure is also developed to identify negligible components in a functional ANOVA decomposition of the log gap time hazard. We evaluate the proposed methods with simulation studies and illustrate its use through the analysis of bladder tumor data.     

MCMC estimation of Markov models for ion channels

Ion channels are characterized by inherently stochastic behavior which can be represented by continuous-time Markov models (CTMM). Although methods for collecting data from single ion channels are available, translating a time series of open and closed channels to a CTMM remains a challenge. Bayesian statistics combined with Markov chain Monte Carlo (MCMC) sampling provide means for estimating the rate constants of a CTMM directly from single channel data. In this article, different approaches for the MCMC sampling of Markov models are combined. This method, new to our knowledge, detects overparameterizations and gives more accurate results than existing MCMC methods. It shows similar performance as QuB-MIL, which indicates that it also compares well with maximum likelihood estimators. Data collected from an inositol trisphosphate receptor is used to demonstrate how the best model for a given data set can be found in practice.     

Inference for nonlinear epidemiological models using genealogies and time series

Phylodynamics - the field aiming to quantitatively integrate the ecological and evolutionary dynamics of rapidly evolving populations like those of RNA viruses - increasingly relies upon coalescent approaches to infer past population dynamics from reconstructed genealogies. As sequence data have become more abundant, these approaches are beginning to be used on populations undergoing rapid and rather complex dynamics. In such cases, the simple demographic models that current phylodynamic methods employ can be limiting. First, these models are not ideal for yielding biological insight into the processes that drive the dynamics of the populations of interest. Second, these models differ in form from mechanistic and often stochastic population dynamic models that are currently widely used when fitting models to time series data. As such, their use does not allow for both genealogical data and time series data to be considered in tandem when conducting inference. Here, we present a flexible statistical framework for phylodynamic inference that goes beyond these current limitations. The framework we present employs a recently developed method known as particle MCMC to fit stochastic, nonlinear mechanistic models for complex population dynamics to gene genealogies and time series data in a Bayesian framework. We demonstrate our approach using a nonlinear Susceptible-Infected-Recovered (SIR) model for the transmission dynamics of an infectious disease and show through simulations that it provides accurate estimates of past disease dynamics and key epidemiological parameters from genealogies with or without accompanying time series data.     

Bayesian sequential inference for stochastic kinetic biochemical network models.

As postgenomic biology becomes more predictive, the ability to infer rate parameters of genetic and biochemical networks will become increasingly important. In this paper, we explore the Bayesian estimation of stochastic kinetic rate constants governing dynamic models of intracellular processes. The underlying model is replaced by a diffusion approximation where a noise term represents intrinsic stochastic behavior and the model is identified using discrete-time (and often incomplete) data that is subject to measurement error. Sequential MCMC methods are then used to sample the model parameters on-line in several data-poor contexts. The methodology is illustrated by applying it to the estimation of parameters in a simple prokaryotic auto-regulatory gene network.     

Accelerated maximum likelihood parameter estimation for stochastic biochemical systems

ABSTRACT: BACKGROUND: A prerequisite for the mechanistic simulation of a biochemical system is detailed knowledge of its kinetic parameters. Despite recent experimental advances, the estimation of unknown parameter values from observed data is still a bottleneck for obtaining accurate simulation results. Many methods exist for parameter estimation in deterministic biochemical systems; methods for discrete stochastic systems are less well developed. Given the probabilistic nature of stochastic biochemical models, a natural approach is to choose parameter values that maximize the probability of the observed data with respect to the unknown parameters, a.k.a. the maximum likelihood parameter estimates (MLEs). MLE computation for all but the simplest models requires the simulation of many system trajectories that are consistent with experimental data. For models with unknown parameters, this presents a computational challenge, as the generation of consistent trajectories can be an extremely rare occurrence. RESULTS: We have developed Monte Carlo Expectation-Maximization with Modified Cross-Entropy Method (MCEM2): an accelerated method for calculating MLEs that combines advances in rare event simulation with a computationally efficient version of the Monte Carlo expectation-maximization (MCEM) algorithm. Our method requires no prior knowledge regarding parameter values, and it automatically provides a multivariate parameter uncertainty estimate. We applied the method to five stochastic systems of increasing complexity, progressing from an analytically tractable pure-birth model to a computationally demanding model of yeast-polarization. Our results demonstrate that MCEM2 substantially accelerates MLE computation on all tested models when compared to a stand-alone version of MCEM. Additionally, we show how our method identifies parameter values for certain classes of models more accurately than two recently proposed computationally efficient methods. CONCLUSIONS: This work provides a novel, accelerated version of a likelihood-based parameter estimation method that can be readily applied to stochastic biochemical systems. In addition, our results suggest opportunities for added efficiency improvements that will further enhance our ability to mechanistically simulate biological processes.     

Towards linkage analysis with markers in linkage disequilibrium by graphical modelling

We review recent developments of MCMC integration methods for computations on graphical models for two applications in statistical genetics: modelling allelic association and pedigree based linkage analysis. We discuss and illustrate estimation of graphical models from haploid and diploid genotypes, and the importance of MCMC updating schemes beyond what is strictly necessary for irreducibility. We then outline an approach combining these methods to compute linkage statistics when alleles at the marker loci are in linkage disequilibrium. Other extensions suitable for analysis of SNP genotype data in pedigrees are also discussed and programs that implement these methods, and which are available from the author's web site, are described. We conclude with a discussion of how this still experimental approach might be further developed.     

Variable selection for discriminant analysis with Markov random field priors for the analysis of microarray data

MOTIVATION: Discriminant analysis is an effective tool for the classification of experimental units into groups. Here, we consider the typical problem of classifying subjects according to phenotypes via gene expression data and propose a method that incorporates variable selection into the inferential procedure, for the identification of the important biomarkers. To achieve this goal, we build upon a conjugate normal discriminant model, both linear and quadratic, and include a stochastic search variable selection procedure via an MCMC algorithm. Furthermore, we incorporate into the model prior information on the relationships among the genes as described by a gene-gene network. We use a Markov random field (MRF) prior to map the network connections among genes. Our prior model assumes that neighboring genes in the network are more likely to have a joint effect on the relevant biological processes. RESULTS: We use simulated data to assess performances of our method. In particular, we compare the MRF prior to a situation where independent Bernoulli priors are chosen for the individual predictors. We also illustrate the method on benchmark datasets for gene expression. Our simulation studies show that employing the MRF prior improves on selection accuracy. In real data applications, in addition to identifying markers and improving prediction accuracy, we show how the integration of existing biological knowledge into the prior model results in an increased ability to identify genes with strong discriminatory power and also aids the interpretation of the results.     

Exploratory Bayesian model selection for serial genetics data

Characterizing the process by which molecular and cellular level changes occur over time will have broad implications for clinical decision making and help further our knowledge of disease etiology across many complex diseases. However, this presents an analytic challenge due to the large number of potentially relevant biomarkers and the complex, uncharacterized relationships among them. We propose an exploratory Bayesian model selection procedure that searches for model simplicity through independence testing of multiple discrete biomarkers measured over time. Bayes factor calculations are used to identify and compare models that are best supported by the data. For large model spaces, i.e., a large number of multi-leveled biomarkers, we propose a Markov chain Monte Carlo (MCMC) stochastic search algorithm for finding promising models. We apply our procedure to explore the extent to which HIV-1 genetic changes occur independently over time.     

Performance of Markov chain-Monte Carlo approaches for mapping genes in oligogenic models with an unknown number of loci

Markov chain-Monte Carlo (MCMC) techniques for multipoint mapping of quantitative trait loci have been developed on nuclear-family and extended-pedigree data. These methods are based on repeated sampling-peeling and gene dropping of genotype vectors and random sampling of each of the model parameters from their full conditional distributions, given phenotypes, markers, and other model parameters. We further refine such approaches by improving the efficiency of the marker haplotype-updating algorithm and by adopting a new proposal for adding loci. Incorporating these refinements, we have performed an extensive simulation study on simulated nuclear-family data, varying the number of trait loci, family size, displacement, and other segregation parameters. Our simulation studies show that our MCMC algorithm identifies the locations of the true trait loci and estimates their segregation parameters well-provided that the total number of sibship pairs in the pedigree data is reasonably large, heritability of each individual trait locus is not too low, and the loci are not too close together. Our MCMC algorithm was shown to be significantly more efficient than LOKI (Heath 1997) in our simulation study using nuclear-family data.     

Neural fate decisions mediated by combinatorial regulation of Hes1 and miR-9

In the nervous system, Hes1 shows an oscillatory manner in neural progenitors but a persistent one in neurons. Many models involving Hes1 have been provided for the study of neural differentiation but few of them take the role of microRNA into account. It is known that a microRNA, miR-9, plays crucial roles in modulating Hes1 oscillations. However, the roles of miR-9 in controlling Hes1 oscillations and inducing transition between different cell fates still need to be further explored. Here we provide a mathematical model to show the interaction between miR-9 and Hes1, with the aim of understanding how the Hes1 oscillations are produced, how they are controlled, and further, how they are terminated. Based on the experimental findings, the model demonstrates the essential roles of Hes1 and miR-9 in regulating the dynamics of the system. In particular, the model suggests that the balance between miR-9 and Hes1 plays important roles in the choice between progenitor maintenance and neural differentiation. In addition, the synergistic (or antagonistic) effects of several important regulations are investigated so as to elucidate the effects of combinatorial regulation in neural decision-making. Our model provides a qualitative mechanism for understanding the process in neural fate decisions regulated by Hes1 and miR-9.     

Linkage analysis of quantitative trait loci in multiple line crosses

Simple line crosses, for example, backcross and F₂, are commonly used in mapping quantitative trait loci (QTL). However, these simple crosses are rarely used alone in commercial plant breeding; rather, crosses involving multiple inbred lines or several simple crosses but connected by shared inbred lines may be common in plant breeding. Mapping QTL using crosses of multiple lines is more relevant to plant breeding. Unfortunately, current statistical methods and computer programs of QTL mapping are all designed for simple line crosses or multiple line crosses but under a regular mating system. It is not straightforward to extend the existing methods to handle multiple line crosses under irregular and complicated mating designs. The major hurdle comes from irregular inbreeding, multiple generations, and multiple alleles. In this study, we develop a Bayesian method implemented via the Markov chain Monte Carlo (MCMC) algorithm for mapping QTL using complicated multiple line crosses. With the MCMC algorithm, we are able to draw a complete path of the gene flow from founder alleles to their descendents via a recursive process. This has greatly simplified the problem caused by irregular mating and inbreeding in the mapping population. Adopting the reversible jump MCMC algorithm, we are able to simultaneously search for multiple QTL along the genome. We can even infer the posterior distribution of the number of QTL, one of the most important parameters in QTL study. Application of the new MCMC based QTL mapping procedure is demonstrated using two different mating designs. Design I involves two inbred lines and their derived F₁, F₂, and BC populations. Design II is a half-diallel cross involving three inbred lines. The two designs appear different, but can be handled with the same robust computer program.     

A comparison of strategies for Markov chain Monte Carlo computation in quantitative genetics

In quantitative genetics, Markov chain Monte Carlo (MCMC) methods are indispensable for statistical inference in non-standard models like generalized linear models with genetic random effects or models with genetically structured variance heterogeneity. A particular challenge for MCMC applications in quantitative genetics is to obtain efficient updates of the high-dimensional vectors of genetic random effects and the associated covariance parameters. We discuss various strategies to approach this problem including reparameterization, Langevin-Hastings updates, and updates based on normal approximations. The methods are compared in applications to Bayesian inference for three data sets using a model with genetically structured variance heterogeneity.     

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.     

A new method of identifying a systems based on the minimal quadratic discrepancy criterion for biophysics problems

A wide range of biophysical systems are described by nonlinear dynamic models mathematically presented as a set of ordinary differential equations in the Cauchy explicit form: [formula: see text] Fij(X1(t),..,XN(t),t), (i = 1,...,N, j = 1,...,M), where Fij (X1(t), ..., XN(t), t) is a set of basis functions satisfying the Lipschitz condition. We investigate the problem of evaluation of model constants aij (the system identification) using experimental data about the time dependence of the dynamic parameters of the system Xi(t). A new method of system identification for the class of similar nonlinear dynamic models is proposed. It is shown that the problem of identifying an initial nonlinear model can be reduced to the solution of a system of linear equations for the matrix of the dynamic model constants [aj]i. It is proposed to determine the set of dynamic model constants aij using the criterion of minimal quadratic discrepancy for the time dependence of the set of dynamic parameters Xi(t). An important special case of the nonlinear model, the quadratic model, is considered. Test problems of identification using this method are presented for two nonlinear systems: the Van der Pol type multiparametric nonlinear oscillator and the strange attractor of Ressler, a widely known example of dynamic systems showing the stochastic behavior.