Multicanonical sampling of rare events in random matrices

Probing the large deviations for the beta random walk in random medium

We consider a discrete-time random walk on a one-dimensional lattice with space- and time-dependent random jump probabilities, known as the beta random walk. We are interested in the probability that, for a given realization of the jump probabilities (a sample), a walker starting at the origin at time t=0 is at position beyond ξsqrt[T/2] at time T. This probability fluctuates from sample to sample and we study the large-deviation rate function, which characterizes the tails of its distribution at large time T≫1. It is argued that, up to a simple rescaling, this rate function is identical to the one recently obtained exactly by two of the authors for the continuum version of the model. That continuum model also appears in the macroscopic fluctuation theory of a class of lattice gases, e.g., in the so-called KMP model of heat transfer. An extensive numerical simulation of the beta random walk, based on an importance sampling algorithm, is found in good agreement with the detailed analytical predictions. A first-order transition in the tilted measure, predicted to occur in the continuum model, is also observed in the numerics.

Evolution of hierarchy and irreversibility in theoretical cell differentiation model

The process of cell differentiation in multicellular organisms is characterized by hierarchy and irreversibility in many cases. However, the conditions and selection pressures that give rise to these characteristics remain poorly understood. By using a mathematical model, here we show that the network of differentiation potency (differentiation diagram) becomes necessarily hierarchical and irreversible by increasing the number of terminally differentiated states under certain conditions. The mechanisms generating these characteristics are clarified using geometry in the cell state space. The results demonstrate that the hierarchical organization and irreversibility can manifest independently of direct selection pressures associated with these characteristics, instead they appear to evolve as byproducts of selective forces favoring a diversity of differentiated cell types. The study also provides a new perspective on the structure of gene regulatory networks that produce hierarchical and irreversible differentiation diagrams. These results indicate some constraints on cell differentiation, which are expected to provide a starting point for theoretical discussion of the implicit limits and directions of evolution in multicellular organisms.

Rare-event sampling analysis uncovers the fitness landscape of the genetic code

The genetic code refers to a rule that maps 64 codons to 20 amino acids. Nearly all organisms, with few exceptions, share the same genetic code, the standard genetic code (SGC). While it remains unclear why this universal code has arisen and been maintained during evolution, it may have been preserved under selection pressure. Theoretical studies comparing the SGC and numerically created hypothetical random genetic codes have suggested that the SGC has been subject to strong selection pressure for being robust against translation errors. However, these prior studies have searched for random genetic codes in only a small subspace of the possible code space due to limitations in computation time. Thus, how the genetic code has evolved, and the characteristics of the genetic code fitness landscape, remain unclear. By applying multicanonical Monte Carlo, an efficient rare-event sampling method, we efficiently sampled random codes from a much broader random ensemble of genetic codes than in previous studies, estimating that only one out of every 1020 random codes is more robust than the SGC. This estimate is significantly smaller than the previous estimate, one in a million. We also characterized the fitness landscape of the genetic code that has four major fitness peaks, one of which includes the SGC. Furthermore, genetic algorithm analysis revealed that evolution under such a multi-peaked fitness landscape could be strongly biased toward a narrow peak, in an evolutionary path-dependent manner.

Evolution enhances mutational robustness and suppresses the emergence of a new phenotype: A new computational approach for studying evolution

The aim of this paper is two-fold. First, we propose a new computational method to investigate the particularities of evolution. Second, we apply this method to a model of gene regulatory networks (GRNs) and explore the evolution of mutational robustness and bistability. Living systems have developed their functions through evolutionary processes. To understand the particularities of this process theoretically, evolutionary simulation (ES) alone is insufficient because the outcomes of ES depend on evolutionary pathways. We need a reference system for comparison. An appropriate reference system for this purpose is an ensemble of the randomly sampled genotypes. However, generating high-fitness genotypes by simple random sampling is difficult because such genotypes are rare. In this study, we used the multicanonical Monte Carlo method developed in statistical physics to construct a reference ensemble of GRNs and compared it with the outcomes of ES. We obtained the following results. First, mutational robustness was significantly higher in ES than in the reference ensemble at the same fitness level. Second, the emergence of a new phenotype, bistability, was delayed in evolution. Third, the bistable group of GRNs contains many mutationally fragile GRNs compared with those in the non-bistable group. This suggests that the delayed emergence of bistability is a consequence of the mutation-selection mechanism.

Living systems are products of evolution, and their present forms reflect their evolutionary history. Thus, to investigate the particularity of the evolutionary process by computer simulations, an appropriate reference system is needed for comparison with the outcomes of evolutionary simulations. In this study, we considered a model of gene regulatory networks (GRNs). Our idea was to construct a reference ensemble comprising randomly generated GRNs. To produce GRNs with high fitness values, which are rare, we employed a “rare event sampling” method developed in statistical physics. In particular, we focused on the evolution of mutational robustness. Living systems do not lose viability readily, even when some genes are mutated. This trait, called mutational robustness, has developed throughout evolution, along with functionality. Using the abovementioned method, we found that mutational robustness resulting from evolution exceeded that of the reference set. Therefore, mutational robustness is enhanced by evolution. We also found that the emergence of a new phenotype was significantly delayed in evolution. Our results suggest that this delay is a consequence of the fact that mutationally robust GRNs are favored by evolution.

Emergence of cooperative bistability and robustness of gene regulatory networks

Gene regulatory networks (GRNs) are complex systems in which many genes regulate mutually to adapt the cell state to environmental conditions. In addition to function, the GRNs possess several kinds of robustness. This robustness means that systems do not lose their functionality when exposed to disturbances such as mutations or noise, and is widely observed at many levels in living systems. Both function and robustness have been acquired through evolution. In this respect, GRNs utilized in living systems are rare among all possible GRNs. In this study, we explored the fitness landscape of GRNs and investigated how robustness emerged in highly-fit GRNs. We considered a toy model of GRNs with one input gene and one output gene. The difference in the expression level of the output gene between two input states, "on" and "off", was considered as fitness. Thus, the determination of the fitness of a GRN was based on how sensitively it responded to the input. We employed the multicanonical Monte Carlo method, which can sample GRNs randomly in a wide range of fitness levels, and classified the GRNs according to their fitness. As a result, the following properties were found: (1) Highly-fit GRNs exhibited bistability for intermediate input between "on" and "off". This means that such GRNs responded to two input states by using different fixed points of dynamics. This bistability emerges necessarily as fitness increases. (2) These highly-fit GRNs were robust against noise because of their bistability. In other words, noise robustness is a byproduct of high fitness. (3) GRNs that were robust against mutations were not extremely rare among the highly-fit GRNs. This implies that mutational robustness is readily acquired through the evolutionary process. These properties are universal irrespective of the evolutionary pathway, because the results do not rely on evolutionary simulation.

Probing large deviations of the Kardar-Parisi-Zhang equation at short times with an importance sampling of directed polymers in random media

The one-point distribution of the height for the continuum Kardar-Parisi-Zhang equation is determined numerically using the mapping to the directed polymer in a random potential at high temperature. Using an importance sampling approach, the distribution is obtained over a large range of values, down to a probability density as small as 10^{-1000} in the tails. The short-time behavior is investigated and compared with recent analytical predictions for the large-deviation forms of the probability of rare fluctuations, showing a spectacular agreement with the analytical expressions. The flat and stationary initial conditions are studied in the full space, together with the droplet initial condition in the half-space.

Theory for the conditioned spectral density of noninvariant random matrices

We develop a theoretical approach to compute the conditioned spectral density of N×N noninvariant random matrices in the limit N→∞. This large deviation observable, defined as the eigenvalue distribution conditioned to have a fixed fraction k of eigenvalues smaller than x∈R, provides the spectrum of random matrix samples that deviate atypically from the average behavior. We apply our theory to sparse random matrices and unveil strikingly different and generic properties, namely, (i) their conditioned spectral density has compact support, (ii) it does not experience any abrupt transition for k around its typical value, and (iii) its eigenvalues do not accumulate at x. Moreover, our work points towards other types of transitions in the conditioned spectral density for values of k away from its typical value. These properties follow from the weak or absent eigenvalue repulsion in sparse ensembles and they are in sharp contrast to those displayed by classic or rotationally invariant random matrices. The exactness of our theoretical findings are confirmed through numerical diagonalization of finite random matrices.

A single-walker approach for studying quasi-nonergodic systems

The jump-walking Monte-Carlo algorithm is revisited and updated to study the equilibrium properties of systems exhibiting quasi-nonergodicity. It is designed for a single processing thread as opposed to currently predominant algorithms for large parallel processing systems. The updated algorithm is tested on the Ising model and applied to the lattice-gas model for sorption in aerogel at low temperatures, when dynamics of the system is critically slowed down. It is demonstrated that the updated jump-walking simulations are able to produce equilibrium isotherms which are typically hidden by the hysteresis effect characteristic of the standard single-flip simulations.

Phase transitions in the condition-number distribution of Gaussian random matrices

We study the statistics of the condition number κ=λ_{max}/λ_{min} (the ratio between largest and smallest squared singular values) of N×M Gaussian random matrices. Using a Coulomb fluid technique, we derive analytically and for large N the cumulative P(κ<x) and tail-cumulative P(κ>x) distributions of κ. We find that these distributions decay as P(κ<x)≈exp[-βN^{2}Φ_{-}(x)] and P(κ>x)≈exp[-βNΦ_{+}(x)], where β is the Dyson index of the ensemble. The left and right rate functions Φ_{±}(x) are independent of β and calculated exactly for any choice of the rectangularity parameter α=M/N-1>0. Interestingly, they show a weak nonanalytic behavior at their minimum 〈κ〉 (corresponding to the average condition number), a direct consequence of a phase transition in the associated Coulomb fluid problem. Matching the behavior of the rate functions around 〈κ〉, we determine exactly the scale of typical fluctuations ∼O(N^{-2/3}) and the tails of the limiting distribution of κ. The analytical results are in excellent agreement with numerical simulations.

How many eigenvalues of a Gaussian random matrix are positive?

We study the probability distribution of the index N(+), i.e., the number of positive eigenvalues of an N×N Gaussian random matrix. We show analytically that, for large N and large N(+) with the fraction 0≤c=N(+)/N≤1 of positive eigenvalues fixed, the index distribution P(N(+)=cN,N)~exp[-βN(2)Φ(c)] where β is the Dyson index characterizing the Gaussian ensemble. The associated large deviation rate function Φ(c) is computed explicitly for all 0≤c≤1. It is independent of β and displays a quadratic form modulated by a logarithmic singularity around c=1/2. As a consequence, the distribution of the index has a Gaussian form near the peak, but with a variance Δ(N) of index fluctuations growing as Δ(N)~lnN/βπ(2) for large N. For β=2, this result is independently confirmed against an exact finite-N formula, yielding Δ(N)=lnN/2π(2)+C+O(N(-1)) for large N, where the constant C for even N has the nontrivial value C=(γ+1+3ln2)/2π(2)≃0.185 248… and γ=0.5772… is the Euler constant. We also determine for large N the probability that the interval [ζ(1),ζ(2)] is free of eigenvalues. Some of these results have been announced in a recent letter [Phys. Rev. Lett. 103, 220603 (2009)].