Auxiliary Markov Chain Research Articles

Kinetic Network in Milestoning: Clustering, Reduction, and Transition Path Analysis.

We present a reduction of the Milestoning (ReM) algorithm to analyze the high-dimensional Milestoning kinetic network. The algorithm reduces the Milestoning network to low dimensions but preserves essential kinetic information, such as local residence time, exit time, and mean first passage time between any two states. This is achieved in three steps. First, nodes (milestones) in the high-dimensional Milestoning network are grouped into clusters based on the metastability identified by an auxiliary continuous-time Markov chain. Our clustering method is applicable not only to time-reversible networks but also to nonreversible networks generated from practical simulations with statistical fluctuations. Second, a reduced network is established via network transformation, containing only the core sets of clusters as nodes. Finally, transition pathways are analyzed in the reduced network based on the transition path theory. The algorithm is illustrated using a toy model and a solvated alanine dipeptide in two and four dihedral angles.

Stoechiometric and dynamical autocatalysis for diluted chemical reaction networks.

Autocatalysis underlies the ability of chemical and biochemical systems to replicate. Recently, Blokhuis et al. (PNAS 117(41):25230-25236, 2020) gave a stoechiometric definition of autocatalysis for reaction networks, stating the existence of a combination of reactions such that the balance for all autocatalytic species is strictly positive, and investigated minimal autocatalytic networks, called autocatalytic cores. By contrast, spontaneous autocatalysis-namely, exponential amplification of all species internal to a reaction network, starting from a diluted regime, i.e. low concentrations-is a dynamical property. We introduce here a topological condition (Top) for autocatalysis, namely: restricting the reaction network description to highly diluted species, we assume existence of a strongly connected component possessing at least one reaction with multiple products (including multiple copies of a single species). We find this condition to be necessary and sufficient for stoechiometric autocatalysis. When degradation reactions have small enough rates, the topological condition further ensures dynamical autocatalysis, characterized by a strictly positive Lyapunov exponent giving the instantaneous exponential growth rate of the system. The proof is generally based on the study of auxiliary Markov chains. We provide as examples general autocatalytic cores of Type I and Type III in the typology of Blokhuis et al. (PNAS 117(41):25230-25236, 2020) . In a companion article (Unterberger in Dynamical autocatalysis for autocatalytic cores, 2021), Lyapunov exponents and the behavior in the growth regime are studied quantitatively beyond the present diluted regime .

A Bayesian modified Ising model for identifying spatially variable genes from spatial transcriptomics data

A recent technology breakthrough in spatial molecular profiling (SMP) has enabled the comprehensive molecular characterizations of single cells while preserving spatial information. It provides new opportunities to delineate how cells from different origins form tissues with distinctive structures and functions. One immediate question in SMP data analysis is to identify genes whose expressions exhibit spatially correlated patterns, called spatially variable (SV) genes. Most current methods to identify SV genes are built upon the geostatistical model with Gaussian process to capture the spatial patterns. However, the Gaussian process models rely on ad hoc kernels that could limit the models' ability to identify complex spatial patterns. In order to overcome this challenge and capture more types of spatial patterns, we introduce a Bayesian approach to identify SV genes via a modified Ising model. The key idea is to use the energy interaction parameter of the Ising model to characterize spatial expression patterns. We use auxiliary variable Markov chain Monte Carlo algorithms to sample from the posterior distribution with an intractable normalizing constant in the model. Simulation studies using both simulated and synthetic data showed that the energy-based modeling approach led to higher accuracy in detecting SV genes than those kernel-based methods. When applied to two real spatial transcriptomics (ST) datasets, the proposed method discovered novel spatial patterns that shed light on the biological mechanisms. In summary, the proposed method presents a new perspective for analyzing ST data.

On generalized random environment INAR models of higher order: Estimation of random environment states

The behavior of a generalized random environment integer-valued autoregressive model of higher order with geometric marginal distribution and negative binomial thinning operator is dictated by a realization {zn}?,n=1 of an auxiliary Markov chain called random environment process. Element zn represents a state of the environment in moment n ? N and determines all parameters of the model in that moment. In order to apply the model, one first needs to estimate {zn}?,n=1, which was so far done by K-means data clustering. We argue that this approach ignores some information and performs poorly in certain situations. We propose a new method for estimating {zn}?,n=1, which includes the data transformation preceding the clustering, in order to reduce the information loss. To confirm its efficiency, we compare this new approach with the usual one when applied on the simulated and the real-life data, and notice all the benefits obtained from our method.

On the distribution of the hitting time for the [formula omitted]–urn Ehrenfest model

In this paper, we consider the N–urn Ehrenfest model. By utilizing an auxiliary continuous-time Markov chain, we obtain the explicit formula for the Laplace transform of the hitting time from a single state to a set A of states where A satisfies some symmetric properties. After obtaining the Laplace transform, we are able to compute the high-order moments (especially, variance) for the hitting time.

Introducing a Novel Multi-Level Method for Simulating the pH Dependence of Charge State Fluctuations and Conformational Ensembles of Intrinsically Disordered Proteins

Simulations of intrinsically disordered proteins (IDPs) use fixed charge models where the charge states of ionizable residues are set using model compound pKa values. Recent experiments have questioned the validity of fixed charge models. For a 100-residue IDP without sequence biases there can be roughly 103 or more unique possible charge states. Accurate and fully ergodic simulations require converged conformational sampling of each of the charge states. To facilitate this we have developed a multi-level Monte Carlo and free energy based method that leverages the ABSINTH implicit solvation model to enable efficient simulations of the conformational ensembles of all possible charge states for a given sequence. This is achieved by grouping charge states based on their net charge and adding auxiliary Markov chains that enable protonation/deprotonation of ionizable residues. At a given pH, the populations of thermodynamically accessible charge states and conformations are calculated using the Bennett Acceptance Ratio method and data from simulations obtained for all possible charge states. Key features of the ABSINTH model allow reductions in the number of replicas needed for free energy calculations. We also estimate shifts in free energies upon protonation / deprotonation relative to model compounds. Using these estimates, we quantify charge state fluctuations - known as charge regulation - to estimate pKa shifts within IDPs and extract the pH dependencies of the populations of charge states and conformational ensembles. In peptides excised from IDP sequences we find that local sequence contexts cause shifts in pKa values that can be larger than one pH unit. The method, with the incorporation of heuristics to enable increased efficiencies, is being deployed for simulations of charge states and conformational ensembles of a spectrum of IDPs.

Layered nested Markov chain Monte Carlo.

A configurational sampling algorithm based on nested layerings of Markov chains (Layered Nested Markov Chain Monte Carlo or L-NMCMC) is presented for simulations of systems characterized by rugged free energy landscapes. The layerings are generated using a set of auxiliary potential energy surfaces. The implementation of the method is demonstrated in the context of a rugged, two-dimensional potential energy surface. The versatility of the algorithm is next demonstrated on a simple, many-body system, namely, a canonical Lennard-Jones fluid in the liquid state. In that example, different layering schemes and auxiliary potentials are used, including variable cutoff distances and excluded-volume tempering. In addition to calculating a variety of properties of the system, it is also shown that L-NMCMC, when combined with a free-energy perturbation formalism, provides a straightforward means to construct approximate free-energy surfaces at no additional computational cost using the sampling distributions of each auxiliary Markov chain. The proposed L-NMCMC scheme is general in that it could be complementary to any number of methods that rely on sampling from a target distribution or methods that exploit a hierarchy of time scales and/or length scales through decomposition of the potential energy.

Ancestral Lineages and Limit Theorems for Branching Markov Chains in Varying Environment

We consider branching processes in discrete time for structured population in varying environment. Each individual has a trait which belongs to some general state space and both its reproduction law and the trait inherited by its offsprings may depend on its trait and the environment. We study the long-time behavior of the population and the ancestral lineage of typical individuals under general assumptions. We focus on the mean growth rate and the trait distribution among the population. The approach relies on many-to-one formulae and the analysis of the genealogy, and a key role is played by well-chosen (possibly non-homogeneous) Markov chains. The applications use large deviations principles and the Harris ergodicity for these auxiliary Markov chains.

Minimal auxiliary Markov chains through sequential elimination of states

ABSTRACTWhen using an auxiliary Markov chain to compute the distribution of a pattern statistic, the computational complexity is directly related to the number of Markov chain states. Theory related to minimal deterministic finite automata have been applied to large state spaces to reduce the number of Markov chain states so that only a minimal set remains. In this paper, a characterization of equivalent states is given so that extraneous states are deleted during the process of forming the state space, improving computational efficiency. The theory extends the applicability of Markov chain based methods for computing the distribution of pattern statistics.

Auxiliary Parameter MCMC for Exponential Random Graph Models

Exponential random graph models (ERGMs) are a well-established family of statistical models for analyzing social networks. Computational complexity has so far limited the appeal of ERGMs for the analysis of large social networks. Efficient computational methods are highly desirable in order to extend the empirical scope of ERGMs. In this paper we report results of a research project on the development of snowball sampling methods for ERGMs. We propose an auxiliary parameter Markov chain Monte Carlo (MCMC) algorithm for sampling from the relevant probability distributions. The method is designed to decrease the number of allowed network states without worsening the mixing of the Markov chains, and suggests a new approach for the developments of MCMC samplers for ERGMs. We demonstrate the method on both simulated and actual (empirical) network data and show that it reduces CPU time for parameter estimation by an order of magnitude compared to current MCMC methods.

Faster exact distributions of pattern statistics through sequential elimination of states

When using an auxiliary Markov chain (AMC) to compute sampling distributions, the computational complexity is directly related to the number of Markov chain states. For certain complex pattern statistics, minimal deterministic finite automata (DFA) have been used to facilitate efficient computation by reducing the number of AMC states. For example, when statistics of overlapping pattern occurrences are counted differently than non-overlapping occurrences, a DFA consisting of prefixes of patterns extended to overlapping occurrences has been generated and then minimized to form an AMC. However, there are situations where forming such a DFA is computationally expensive, e.g., with computing the distribution of spaced seed coverage. In dealing with this problem, we develop a method to obtain a small set of states during the state generation process without forming a DFA, and show that a huge reduction in the size of the AMC can be attained.

A Decomposition Based Proof for Fast Mixing of a Markov Chain over Balanced Realizations of a Joint Degree Matrix

A joint degree matrix (JDM) specifies the number of connections between nodes of given degrees in a graph, for all degree pairs, and uniquely determines the degree sequence of the graph. We consider the space of all balanced realizations of an arbitrary JDM, realizations in which the links between any two fixed-degree groups of nodes are placed as uniformly as possible. We prove that a swap Markov chain Monte Carlo algorithm in the space of all balanced realizations of an arbitrary graphical JDM mixes rapidly, i.e., the relaxation time of the chain is bounded from above by a polynomial in the number of nodes $n$. To prove fast mixing, we first prove a general factorization theorem similar to the Martin--Randall method for disjoint decompositions (partitions). This theorem can be used to bound from below the spectral gap with the help of fast mixing subchains within every partition and a bound on an auxiliary Markov chain between the partitions. Our proof of the general factorization theorem is direct and uses conductance based methods (Cheeger inequality).

Hamiltonian Switch Metropolis Monte Carlo Simulations for Improved Conformational Sampling of Intrinsically Disordered Regions Tethered to Ordered Domains of Proteins

There is growing interest in the topic of intrinsically disordered proteins (IDPs). Atomistic Metropolis Monte Carlo (MMC) simulations based on novel implicit solvation models have yielded useful insights regarding sequence-ensemble relationships for IDPs modeled as autonomous units. However, a majority of naturally occurring IDPs are tethered to ordered domains. Tethering introduces additional energy scales and this creates the challenge of broken ergodicity for standard MMC sampling or molecular dynamics that cannot be readily alleviated by using generalized tempering methods. We have designed, deployed, and tested our adaptation of the Nested Markov Chain Monte Carlo sampling algorithm. We refer to our adaptation as Hamiltonian Switch Metropolis Monte Carlo (HS-MMC) sampling. In this method, transitions out of energetic traps are enabled by the introduction of an auxiliary Markov chain that draws conformations for the disordered region from a Boltzmann distribution that is governed by an alternative potential function that only includes short-range steric repulsions and conformational restraints on the ordered domain. We show using multiple, independent runs that the HS-MMC method yields conformational distributions that have similar and reproducible statistical properties, which is in direct contrast to standard MMC for equivalent amounts of sampling. The method is efficient and can be deployed for simulations of a range of biologically relevant disordered regions that are tethered to ordered domains.

A Monte Carlo Metropolis-Hastings Algorithm for Sampling from Distributions with Intractable Normalizing Constants

Simulating from distributions with intractable normalizing constants has been a long-standing problem in machine learning. In this letter, we propose a new algorithm, the Monte Carlo Metropolis-Hastings (MCMH) algorithm, for tackling this problem. The MCMH algorithm is a Monte Carlo version of the Metropolis-Hastings algorithm. It replaces the unknown normalizing constant ratio by a Monte Carlo estimate in simulations, while still converges, as shown in the letter, to the desired target distribution under mild conditions. The MCMH algorithm is illustrated with spatial autologistic models and exponential random graph models. Unlike other auxiliary variable Markov chain Monte Carlo (MCMC) algorithms, such as the Møller and exchange algorithms, the MCMH algorithm avoids the requirement for perfect sampling, and thus can be applied to many statistical models for which perfect sampling is not available or very expensive. The MCMH algorithm can also be applied to Bayesian inference for random effect models and missing data problems that involve simulations from a distribution with intractable integrals.

Almost sure stability and instability for switching-jump-diffusion systems with state-dependent switching

This work develops almost sure stability criteria for switching-jump-diffusion systems with state-dependent switching. By means of introducing certain auxiliary Markov chains and constructing certain order-preserving couplings, upper and lower “stability envelops” are constructed, which lead to systems with “upper and lower” approximating Markov chains. Using these approximations, sufficient conditions that are relatively easily verifiable for the almost sure stability and instability are obtained. When the jump process is missing, it is demonstrated that the techniques work equally well and provide a way to analyze the corresponding switching diffusion systems with x-dependent switching. In addition, stochastic stabilization and destabilization are examined. Moreover, illustrative examples are provided for demonstration.

Use of SAMC for Bayesian analysis of statistical models with intractable normalizing constants

Statistical inference for the models with intractable normalizing constants has attracted much attention. During the past two decades, various approximation- or simulation-based methods have been proposed for the problem, such as the Monte Carlo maximum likelihood method and the auxiliary variable Markov chain Monte Carlo methods. The Bayesian stochastic approximation Monte Carlo algorithm specifically addresses this problem: It works by sampling from a sequence of approximate distributions with their average converging to the target posterior distribution, where the approximate distributions can be achieved using the stochastic approximation Monte Carlo algorithm. A strong law of large numbers is established for the Bayesian stochastic approximation Monte Carlo estimator under mild conditions. Compared to the Monte Carlo maximum likelihood method, the Bayesian stochastic approximation Monte Carlo algorithm is more robust to the initial guess of model parameters. Compared to the auxiliary variable MCMC methods, the Bayesian stochastic approximation Monte Carlo algorithm avoids the requirement for perfect samples, and thus can be applied to many models for which perfect sampling is not available or very expensive. The Bayesian stochastic approximation Monte Carlo algorithm also provides a general framework for approximate Bayesian analysis.

Stability of a Markov-modulated Markov Chain, with Application to a Wireless Network Governed by two Protocols

We consider a discrete-time Markov chain (Xt,Yt), t = 0, 1, 2,…, where the X-component forms a Markov chain itself. Assume that (Xt) is Harris-ergodic and consider an auxiliary Markov chain [Formula: see text] whose transition probabilities are the averages of transition probabilities of the Y-component of the (X, Y)-chain, where the averaging is weighted by the stationary distribution of the X-component. We first provide natural conditions in terms of test functions ensuring that the [Formula: see text]-chain is positive recurrent and then prove that these conditions are also sufficient for positive recurrence of the original chain (Xt, Yt). The we prove a “multi-dimensional” extension of the result obtained. In the second part of the paper, we apply our results to two versions of a multi-access wireless model governed by two randomised protocols.

Stability of a Markov-modulated Markov chain, with application to a wireless network governed by two protocols

We consider a discrete-time Markov chain (Xt,Yt), t = 0, 1, 2,…, where the X-component forms a Markov chain itself. Assume that (Xt) is Harris-ergodic and consider an auxiliary Markov chain {Yt} whose transition probabilities are the averages of transition probabilities of the Y-component of the (X, Y)-chain, where the averaging is weighted by the stationary distribution of the X-component. We first provide natural conditions in terms of test functions ensuring that the Y-chain is positive recurrent and then prove that these conditions are also sufficient for positive recurrence of the original chain (Xt, Yt). The we prove a “multi-dimensional” extension of the result obtained. In the second part of the paper, we apply our results to two versions of a multi-access wireless model governed by two randomised protocols.

Excited Brownian motions as limits of excited random walks

We obtain the convergence in law of a sequence of excited (also called cookies) random walks toward an excited Brownian motion. This last process is a continuous semi-martingale whose drift is a function, say φ, of its local time. It was introduced by Norris, Rogers and Williams as a simplified version of Brownian polymers, and then recently further studied by the authors. To get our results we need to renormalize together the sequence of cookies, the time and the space in a convenient way. The proof follows a general approach already taken by Toth and his coauthors in multiple occasions, which goes through Ray-Knight type results. Namely we first prove, when φ is bounded and Lipschitz, that the convergence holds at the level of the local time processes. This is done via a careful study of the transition kernel of an auxiliary Markov chain which describes the local time at a given level. Then we prove a tightness result and deduce the convergence at the level of the full processes.

Waiting time distribution of generalized later patterns

In this paper the concept of later waiting time distributions for patterns in multi-state trials is generalized to cover a collection of compound patterns that must all be counted pattern-specific numbers of times, and a practical method is given to compute the generalized distribution. The solution given applies to overlapping counting and two types of non-overlapping counting, and the underlying sequences are assumed to be Markovian of a general order. Patterns are allowed to be weighted so that an occurrence is counted multiple times, and patterns may be completely included in longer patterns. Probabilities are computed through an auxiliary Markov chain. As the state space associated with the auxiliary chain can be quite large if its setup is handled in a naïve fashion, an algorithm is given for generating a “minimal” state space that leaves out states that can never be reached. For the case of overlapping counting, a formula that relates probabilities for intersections of events to probabilities for unions of subsets of the events is also used, so that the distribution is also computed in terms of probabilities for competing patterns. A detailed example is given to illustrate the methodology.