Modeling stochastic dynamics in biochemical systems with feedback using maximum caliber

Avoiding matrix exponentials for large transition rate matrices

Exact methods for the exponentiation of matrices of dimension N can be computationally expensive in terms of execution time (N3) and memory requirements (N2), not to mention numerical precision issues. A matrix often exponentiated in the natural sciences is the rate matrix. Here, we explore five methods to exponentiate rate matrices, some of which apply more broadly to other matrix types. Three of the methods leverage a mathematical analogy between computing matrix elements of a matrix exponential process and computing transition probabilities of a dynamical process (technically a Markov jump process, MJP, typically simulated using Gillespie). In doing so, we identify a novel MJP-based method relying on restricting the number of "trajectory" jumps that incurs improved computational scaling. We then discuss this method's downstream implications on mixing properties of Monte Carlo posterior samplers. We also benchmark two other methods of matrix exponentiation valid for any matrix (beyond rate matrices and, more generally, positive definite matrices) related to solving differential equations: Runge-Kutta integrators and Krylov subspace methods. Under conditions where both the largest matrix element and the number of non-vanishing elements scale linearly with N-reasonable conditions for rate matrices often exponentiated-computational time scaling with the most competitive methods (Krylov and one of the MJP-based methods) reduces to N2 with total memory requirements of N.

Critical Comparison of MaxCal and Other Stochastic Modeling Approaches in Analysis of Gene Networks

Learning the underlying details of a gene network with feedback is critical in designing new synthetic circuits. Yet, quantitative characterization of these circuits remains limited. This is due to the fact that experiments can only measure partial information from which the details of the circuit must be inferred. One potentially useful avenue is to harness hidden information from single-cell stochastic gene expression time trajectories measured for long periods of time—recorded at frequent intervals—over multiple cells. This raises the feasibility vs. accuracy dilemma while deciding between different models of mining these stochastic trajectories. We demonstrate that inference based on the Maximum Caliber (MaxCal) principle is the method of choice by critically evaluating its computational efficiency and accuracy against two other typical modeling approaches: (i) a detailed model (DM) with explicit consideration of multiple molecules including protein-promoter interaction, and (ii) a coarse-grain model (CGM) using Hill type functions to model feedback. MaxCal provides a reasonably accurate model while being significantly more computationally efficient than DM and CGM. Furthermore, MaxCal requires minimal assumptions since it is a top-down approach and allows systematic model improvement by including constraints of higher order, in contrast to traditional bottom-up approaches that require more parameters or ad hoc assumptions. Thus, based on efficiency, accuracy, and ability to build minimal models, we propose MaxCal as a superior alternative to traditional approaches (DM, CGM) when inferring underlying details of gene circuits with feedback from limited data.

Inferring a network from dynamical signals at its nodes

We give an approximate solution to the difficult inverse problem of inferring the topology of an unknown network from given time-dependent signals at the nodes. For example, we measure signals from individual neurons in the brain, and infer how they are inter-connected. We use Maximum Caliber as an inference principle. The combinatorial challenge of high-dimensional data is handled using two different approximations to the pairwise couplings. We show two proofs of principle: in a nonlinear genetic toggle switch circuit, and in a toy neural network.

Thermodynamics from first principles: Correlations and nonextensivity

The standard formulation of thermostatistics, being based on the Boltzmann-Gibbs distribution and logarithmic Shannon entropy, describes idealized uncorrelated systems with extensive energies and short-range interactions. In this Rapid Communication, we use the fundamental principles of ergodicity (via Liouville's theorem), the self-similarity of correlations, and the existence of the thermodynamic limit to derive generalized forms of the equilibrium distribution for long-range-interacting systems. Significantly, our formalism provides a justification for the well-studied nonextensive thermostatistics characterized by the Tsallis distribution, which it includes as a special case. We also give the complementary maximum entropy derivation of the same distributions by constrained maximization of the Gibbs-Shannon entropy. The consistency between the ergodic and maximum entropy approaches clarifies the use of the latter in the study of correlations and nonextensive thermodynamics.

The Maximum Caliber Variational Principle for Nonequilibria

Ever since Clausius in 1865 and Boltzmann in 1877, the concepts of entropy and of its maximization have been the foundations for predicting how material equilibria derive from microscopic properties. But, despite much work, there has been no equally satisfactory general variational principle for nonequilibrium situations. However, in 1980, a new avenue was opened by E.T. Jaynes and by Shore and Johnson. We review here maximum caliber, which is a maximum-entropy-like principle that can infer distributions of flows over pathways, given dynamical constraints. This approach is providing new insights, particularly into few-particle complex systems, such as gene circuits, protein conformational reaction coordinates, network traffic, bird flocking, cell motility, and neuronal firing.

Quantitative Kinetic Models from Intravital Microscopy: A Case Study Using Hepatic Transport

The liver performs critical physiological functions, including metabolizing and removing substances, such as toxins and drugs, from the bloodstream. Hepatotoxicity itself is intimately linked to abnormal hepatic transport, and hepatotoxicity remains the primary reason drugs in development fail and approved drugs are withdrawn from the market. For this reason, we propose to analyze, across liver compartments, the transport kinetics of fluorescein-a fluorescent marker used as a proxy for drug molecules-using intravital microscopy data. To resolve the transport kinetics quantitatively from fluorescence data, we account for the effect that different liver compartments (with different chemical properties) have on fluorescein's emission rate. To do so, we develop ordinary differential equation transport models from the data where the kinetics is related to the observable fluorescence levels by "measurement parameters" that vary across different liver compartments. On account of the steep non-linearities in the kinetics and stochasticity inherent to the model, we infer kinetic and measurement parameters by generalizing the method of parameter cascades. For this application, the method of parameter cascades ensures fast and precise parameter estimates from noisy time traces.

Maximum Caliber Can Build and Infer Models of Oscillation in a Three-Gene Feedback Network

Single-cell protein expression time trajectories provide rich temporal data quantifying cellular variability and its role in dictating fitness. However, theoretical models to analyze and fully extract information from these measurements remain limited for three reasons: (i) gene expression profiles are noisy, rendering models of averages inapplicable, (ii) experiments typically measure only a few protein species while leaving other molecular actors-necessary to build traditional bottom-up models-unnoticed, and (iii) measured data are in fluorescence, not particle number. We recently addressed these challenges in an alternate top-down approach using the principle of Maximum Caliber (MaxCal) to model genetic switches with one and two protein species. In the present work we address scalability and broader applicability of MaxCal by extending to a three-gene (A, B, C) feedback network that exhibits oscillation, commonly known as the repressilator. We test MaxCal's inferential power by using synthetic data of noisy protein number time traces-serving as a proxy for experimental data-generated from a known underlying model. We notice that the minimal MaxCal model-accounting for production, degradation, and only one type of symmetric coupling between all three species-reasonably infers several underlying features of the circuit such as the effective production rate, degradation rate, frequency of oscillation, and protein number distribution. Next, we build models of higher complexity including different levels of coupling between A, B, and C and rigorously assess their relative performance. While the minimal model (with four parameters) performs remarkably well, we note that the most complex model (with six parameters) allowing all possible forms of crosstalk between A, B, and C slightly improves prediction of rates, but avoids ad hoc assumption of all the other models. It is also the model of choice based on Bayesian information criteria. We further analyzed time trajectories in arbitrary fluorescence (using synthetic trajectories) to mimic realistic data. We conclude that even with a three-protein system including both fluorescence noise and intrinsic gene expression fluctuations, MaxCal can faithfully infer underlying details of the network, opening future directions to model other network motifs with many species.

Perspective: Maximum caliber is a general variational principle for dynamical systems

We review here Maximum Caliber (Max Cal), a general variational principle for inferring distributions of paths in dynamical processes and networks. Max Cal is to dynamical trajectories what the principle of maximum entropy is to equilibrium states or stationary populations. In Max Cal, you maximize a path entropy over all possible pathways, subject to dynamical constraints, in order to predict relative path weights. Many well-known relationships of non-equilibrium statistical physics-such as the Green-Kubo fluctuation-dissipation relations, Onsager's reciprocal relations, and Prigogine's minimum entropy production-are limited to near-equilibrium processes. Max Cal is more general. While it can readily derive these results under those limits, Max Cal is also applicable far from equilibrium. We give examples of Max Cal as a method of inference about trajectory distributions from limited data, finding reaction coordinates in bio-molecular simulations, and modeling the complex dynamics of non-thermal systems such as gene regulatory networks or the collective firing of neurons. We also survey its basis in principle and some limitations.

An implementation of the maximum-caliber principle by replica-averaged time-resolved restrained simulations

Inferential methods can be used to integrate experimental informations and molecular simulations. The maximum entropy principle provides a framework for using equilibrium experimental data, and it has been shown that replica-averaged simulations, restrained using a static potential, are a practical and powerful implementation of such a principle. Here we show that replica-averaged simulations restrained using a time-dependent potential are equivalent to the principle of maximum caliber, the dynamic version of the principle of maximum entropy, and thus may allow us to integrate time-resolved data in molecular dynamics simulations. We provide an analytical proof of the equivalence as well as a computational validation making use of simple models and synthetic data. Some limitations and possible solutions are also discussed.

An introduction to the maximum entropy approach and its application to inference problems in biology

A cornerstone of statistical inference, the maximum entropy framework is being increasingly applied to construct descriptive and predictive models of biological systems, especially complex biological networks, from large experimental data sets. Both its broad applicability and the success it obtained in different contexts hinge upon its conceptual simplicity and mathematical soundness. Here we try to concisely review the basic elements of the maximum entropy principle, starting from the notion of 'entropy', and describe its usefulness for the analysis of biological systems. As examples, we focus specifically on the problem of reconstructing gene interaction networks from expression data and on recent work attempting to expand our system-level understanding of bacterial metabolism. Finally, we highlight some extensions and potential limitations of the maximum entropy approach, and point to more recent developments that are likely to play a key role in the upcoming challenges of extracting structures and information from increasingly rich, high-throughput biological data.

Maximum Caliber Can Characterize Genetic Switches with Multiple Hidden Species

Gene networks with feedback often involve interactions between multiple species of biomolecules, much more than experiments can actually monitor. Coupled with this is the challenge that experiments often measure gene expression in noisy fluorescence instead of protein numbers. How do we infer biophysical information and characterize the underlying circuits from this limited and convoluted data? We address this by building stochastic models using the principle of Maximum Caliber (MaxCal). MaxCal uses the basic information on synthesis, degradation, and feedback-without invoking any other auxiliary species and ad hoc reactions-to generate stochastic trajectories similar to those typically measured in experiments. MaxCal in conjunction with Maximum Likelihood (ML) can infer parameters of the model using fluctuating trajectories of protein expression over time. We demonstrate the success of the MaxCal + ML methodology using synthetic data generated from known circuits of different genetic switches: (i) a single-gene autoactivating circuit involving five species (including mRNA), (ii) a mutually repressing two-gene circuit (toggle switch) with seven species (including mRNA) considering stochastic time traces of two proteins, and (iii) the same toggle switch circuit considering stochastic time traces of only one of the two proteins. To further challenge the MaxCal + ML inference scheme, we repeat our analysis for the second and third scenario with traces expressed in noisy fluorescence instead of protein number to closely mimic typical experiments. We show that, for all of these models with increasing complexity and obfuscation, the minimal model of MaxCal is still able to capture the fluctuations of the trajectory and infer basic underlying rate parameters when benchmarked against the known values used to generate the synthetic data. Importantly, the model also yields an effective feedback parameter that can be used to quantify interactions within these circuits. These applications show the promise of MaxCal's ability to characterize circuits with limited data, and its utility to better understand evolution and advance design strategies for specific functions.

Building Predictive Models of Genetic Circuits Using the Principle of Maximum Caliber

Learning the underlying details of a gene network is a major challenge in cellular and synthetic biology. We address this challenge by building a chemical kinetic model that utilizes information encoded in the stochastic protein expression trajectories typically measured in experiments. The applicability of the proposed method is demonstrated in an auto-activating genetic circuit, a common motif in natural and synthetic gene networks. Our approach is based on the principle of maximum caliber (MaxCal)-a dynamical analog of the principle of maximum entropy-and builds a minimal model using only three constraints: 1) protein synthesis, 2) protein degradation, and 3) positive feedback. The MaxCal-generated model (described with four parameters) was benchmarked against synthetic data generated using a Gillespie algorithm on a known reaction network (with seven parameters). MaxCal accurately predicts underlying rate parameters of protein synthesis and degradation as well as experimental observables such as protein number and dwell-time distributions. Furthermore, MaxCal yields an effective feedback parameter that can be useful for circuit design. We also extend our methodology and demonstrate how to analyze trajectories that are not in protein numbers but in arbitrary fluorescence units, a more typical condition in experiments. This "top-down" methodology based on minimal information-in contrast to traditional "bottom-up" approaches that require ad hoc knowledge of circuit details-provides a powerful tool to accurately infer underlying details of feedback circuits that are not otherwise visible in experiments and to help guide circuit design.

High-Dimensional Mutant and Modular Thermodynamic Cycles, Molecular Switching, and Free Energy Transduction

Understanding how distinct parts of proteins produce coordinated behavior has driven and continues to drive advances in protein science and enzymology. However, despite consensus about the conceptual basis for allostery, the idiosyncratic nature of allosteric mechanisms resists general approaches. Computational methods can identify conformational transition states from structural changes, revealing common switching mechanisms that impose multistate behavior. Thermodynamic cycles use factorial perturbations to measure coupling energies between side chains in molecular switches that mediate shear during domain motion. Such cycles have now been complemented by modular cycles that measure energetic coupling between separable domains. For one model system, energetic coupling between domains has been shown to be quantitatively equivalent to that between dynamic side chains. Linkages between domain motion, switching residues, and catalysis make nucleoside triphosphate hydrolysis conditional on domain movement, confirming an essential yet neglected aspect of free energy transduction and suggesting the potential generality of these studies.

Competition enhances stochasticity in biochemical reactions

We study stochastic dynamics of two competing complexation reactions (i) A + B↔AB and (ii) A + C↔AC. Such reactions are common in biology where different reactants compete for common resources--examples range from binding enzyme kinetics to gene expression. On the other hand, stochasticity is inherent in biological systems due to small copy numbers. We investigate the complex interplay between competition and stochasticity, using coupled complexation reactions as the model system. Within the master equation formalism, we compute the exact distribution of the number of complexes to analyze equilibrium fluctuations of several observables. Our study reveals that the presence of competition offered by one reaction (say A + C↔AC) can significantly enhance the fluctuation in the other (A + B↔AB). We provide detailed quantitative estimates of this enhanced fluctuation for different combinations of rate constants and numbers of reactant molecules that are typical in biology. We notice that fluctuations can be significant even when two of the reactant molecules (say B and C) are infinite in number, maintaining a fixed stoichiometry, while the other reactant (A) is finite. This is purely due to the coupling mediated via resource sharing and is in stark contrast to the single reaction scenario, where large numbers of one of the components ensure zero fluctuation. Our detailed analysis further highlights regions where numerical estimates of mass action solutions can differ from the actual averages. These observations indicate that averages can be a poor representation of the system, hence analysis that is purely based on averages such as mass action laws can be potentially misleading in such noisy biological systems. We believe that the exhaustive study presented here will provide qualitative and quantitative insights into the role of noise and its enhancement in the presence of competition that will be relevant in many biological settings.

Data-driven quantification of the robustness and sensitivity of cell signaling networks

Robustness and sensitivity of responses generated by cell signaling networks has been associated with survival and evolvability of organisms. However, existing methods analyzing robustness and sensitivity of signaling networks ignore the experimentally observed cell-to-cell variations of protein abundances and cell functions or contain ad hoc assumptions. We propose and apply a data-driven maximum entropy based method to quantify robustness and sensitivity of Escherichia coli (E. coli) chemotaxis signaling network. Our analysis correctly rank orders different models of E. coli chemotaxis based on their robustness and suggests that parameters regulating cell signaling are evolutionary selected to vary in individual cells according to their abilities to perturb cell functions. Furthermore, predictions from our approach regarding distribution of protein abundances and properties of chemotactic responses in individual cells based on cell population averaged data are in excellent agreement with their experimental counterparts. Our approach is general and can be used to evaluate robustness as well as generate predictions of single cell properties based on population averaged experimental data in a wide range of cell signaling systems.

Microcanonical origin of the maximum entropy principle for open systems

There are two distinct approaches for deriving the canonical ensemble. The canonical ensemble either follows as a special limit of the microcanonical ensemble or alternatively follows from the maximum entropy principle. We show the equivalence of these two approaches by applying the maximum entropy formulation to a closed universe consisting of an open system plus bath. We show that the target function for deriving the canonical distribution emerges as a natural consequence of partial maximization of the entropy over the bath degrees of freedom alone. By extending this mathematical formalism to dynamical paths rather than equilibrium ensembles, the result provides an alternative justification for the principle of path entropy maximization as well.

Markov processes follow from the principle of maximum caliber

Markov models are widely used to describe stochastic dynamics. Here, we show that Markov models follow directly from the dynamical principle of maximum caliber (Max Cal). Max Cal is a method of deriving dynamical models based on maximizing the path entropy subject to dynamical constraints. We give three different cases. First, we show that if constraints (or data) are given in the form of singlet statistics (average occupation probabilities), then maximizing the caliber predicts a time-independent process that is modeled by identical, independently distributed random variables. Second, we show that if constraints are given in the form of sequential pairwise statistics, then maximizing the caliber dictates that the kinetic process will be Markovian with a uniform initial distribution. Third, if the initial distribution is known and is not uniform we show that the only process that maximizes the path entropy is still the Markov process. We give an example of how Max Cal can be used to discriminate between different dynamical models given data.