# Spontaneous fine-tuning to environment in many-species chemical reaction networks

See allHide authors and affiliations

Edited by Stanislas Leibler, The Rockefeller University, New York, NY, and approved June 12, 2017 (received for review January 12, 2017)

## Significance

A qualitatively more diverse range of possible behaviors emerge in many-particle systems once external drives are allowed to push the system far from equilibrium; nonetheless, general thermodynamic principles governing nonequilibrium pattern formation and self-assembly have remained elusive, despite intense interest from researchers across disciplines. Here, we use the example of a randomly wired driven chemical reaction network to identify a key thermodynamic feature of a complex, driven system that characterizes the “specialness” of its dynamical attractor behavior. We show that the network’s fixed points are biased toward the extremization of external forcing, causing them to become kinetically stabilized in rare corners of chemical space that are either atypically weakly or strongly coupled to external environmental drives.

## Abstract

A chemical mixture that continually absorbs work from its environment may exhibit steady-state chemical concentrations that deviate from their equilibrium values. Such behavior is particularly interesting in a scenario where the environmental work sources are relatively difficult to access, so that only the proper orchestration of many distinct catalytic actors can power the dissipative flux required to maintain a stable, far-from-equilibrium steady state. In this article, we study the dynamics of an in silico chemical network with random connectivity in an environment that makes strong thermodynamic forcing available only to rare combinations of chemical concentrations. We find that the long-time dynamics of such systems are biased toward states that exhibit a fine-tuned extremization of environmental forcing.

An arrangement of matter may be said to be finely tuned to its environment if it is configured to interact with that environment in a way that is highly atypical for random rearrangements of its components. Whereas some of the most striking examples of such matching between system and environment are found in the architecture of living things (1, 2), the category is in principle far broader, and one might therefore hope to provide an account of emergent fine-tuning in general physical terms.

Recent progress in theoretical nonequilibrium statistical mechanics (3⇓⇓⇓–7) has helped to clarify a general relationship between the likelihood of a system adopting a particular microscopic configuration and the amount of energy absorbed and dissipated during the system’s dynamical history. Essentially, the irreversibility of being driven into a kinetically stable state must be compensated by the release of heat in the surroundings, which frequently occurs when the system absorbs and dissipates energy from an environmental source of work. The implications of this thermodynamic frame are particularly clear in systems where the specific configuration adopted has a large influence on the accessibility of energy from external drives, such that strongly driven states are rare. In such a scenario, having a highly dissipative history becomes tantamount to having a configurational history that includes the rare, finely tuned states that make dissipation possible. This observation suggests that there may be a broad class of nonequilibrium systems for which state-dependent drive strength is the relevant thermodynamic parameter for describing the bias in their exploration of configuration space.

Motivated by these considerations, we carried out simulations of a complexly driven, randomly wired, many-species chemical reaction network, with the aim of identifying a thermodynamic principle governing self-organization far from equilibrium. We found that such networks have a pronounced bias toward settling into rare states in chemical space with extremal forcing, suggesting that emergent fine-tuning is a general tendency of their dynamical behavior.

## Two-Species Model

In this study, we are ultimately interested in the thermodynamic properties of the nonequilibrium fixed points that emerge when a many-species chemical system is driven in a way that varies complexly as a function of the system state. However, to illustrate the relationship between dynamics and thermodynamics in a simple example, we first consider a chemical network with only two possible species *A*). The introduction of this new reaction branch may seem surprising at first, but it is quite natural. It simply reflects that there was always another reaction

The fact that there is a force

The rate of dissipation in a fixed point stabilized by such driving would be positive, and the amount of dissipation per reaction would be determined by the chemical potential drop powering the driven reaction. However, the detailed question of whether such a nonequilibrium fixed point might lie at high or low chemical forcing would entirely depend on the assumption of how

## Many-Species Model

We set out to study the dynamics of a randomly wired chemical network of many reacting species in a dilute, well-stirred mixture, driven out of equilibrium by complexly structured thermodynamic forces induced by the surroundings, schematically illustrated in Fig. 1*B*.

As detailed in *Materials and Methods*, species were allowed to participate in one-body, two-body, and catalyzed three-body reaction channels (

Thus described, such a chemical reaction network has a single dynamical fixed point (known as its chemical equilibrium point) where the concentrations of all different species are equal. To drive the system away from equilibrium, we introduce a collection of additional reaction channels biased by various thermodynamic forces

To this end, we chose forcing functions *C*) (*Materials and Methods*). Note that whereas the forcing function itself contains no explicit time dependence, it does depend strongly on the system state. To achieve this manner of forcing in an explicit model would in general require a complex collection of hidden degrees of freedom that more rapidly relax to a nonequilibrium steady state determined both by the system state and by the fixed chemical potentials of various external particle baths. However, because of the fundamental relationships that hold between kinetic and thermodynamic quantities, a consistent picture of the thermodynamic fluxes flowing through these hidden processes can be inferred from the system kinetics alone.

We implemented this forcing scheme on a many-species chemical reaction network by randomly selecting a subset of reaction channels for duplication and asymmetrization. Thus, for a given channel

Let us summarize that the motivation for our choice of the

## Results

With the reactions fixed, the dynamical evolution of the **Z** is constructed using mass-action kinetics. We numerically solved the reaction rate equations for many randomly generated chemical reaction networks, choosing initial conditions uniformly from the simplex of species concentrations with total *Materials and Methods*. We did not carry out a systematic search of parameter space, which is far too vast, but we found in general that the effects reported below required that nonequilibrium forcing be sufficiently strong, the number of catalyzed reactions be sufficiently large, and reaction graph connectivity be sufficiently sparse. Other than these qualitative constraints, the phenomena observed arose generically across the range of parameter values we examined.

In general, the dynamics evolved to one of many possible fixed points [i.e., to one of the solutions

Fig. 2*A* shows the concentration dynamics of the

We designed the chemical space of this system to contain rare states of extremal thermodynamic forcing, which might be recognized as examples of apparent fine-tuning. Thus, we undertook to analyze the total force magnitude on the network over the course of its dynamical evolution. Labeling the force on the

In the particular realization in Fig. 2*B*, we see the forcing tended to exhibit one of two behaviors. In one scenario, forcing stayed low or decreased; these outcomes corresponded to final concentration distributions that were relatively close to uniform, indicating a near-equilibrium behavior. Alternatively, it was also often the case that the forcing would rise to an extremely high value as the system approached a fixed point.

This observation is made precise in Fig. 2*C* by comparing a histogram of total forces obtained from a uniform sampling of chemical space with the histogram of the final total force at the end of the trajectories, where the dynamics have reached a fixed point. The motivation for this comparison is to detect whether the dynamical fixed points appear to be exceptional combinations of their constituents compared with the ensemble of random rearrangements achieved through uniform sampling of chemical space.

The final force distribution for the representative network studied in Fig. 2 displays a bimodal structure, with the final total forces being either exceptionally high or low, relative to the uniformly sampled set. Whereas other network realizations did not necessarily show this bimodal structure, their final total force distributions were peaked at either unusually high or low levels, essentially filling in one or the other of the humps of the bimodal distribution (Fig. S1). Thus, the dynamical evolution appeared to be biased toward landing on fixed points that have recognizably special and atypical thermodynamic relationships to their external environment (16).

To test for this bias more rigorously, we randomly realized 200 different reaction networks with distinct connectivity, unforced reaction timescales, and forcing parameters and then evolved their dynamics from five random initial conditions. To assess the atypicality of the final total force *C* for each randomly realized network and then assigning a percentile rank in this distribution for each fixed point reached. For example, a trajectory ending with final total force

As a further control, we uniformly sampled configuration space a second time and determined the percentile rank of the observed total force according to the same procedure by which the fixed points were ranked; Fig. 3*A*, *Inset* is a histogram of these percentile ranks for fixed points and for the control set. By computing the ratio between these two histograms we quantified the bias toward each percentile rank, which is displayed in Fig. 3*A*.

Strikingly, the final force percentile distribution separates into two types of behavior, paralleling the outcomes observed in the single reaction network studied above. Frequently, the fixed point of a random reaction network resides at low force, corresponding to a near-equilibrium fine-tuned outcome: The system gets stuck in a near-equilibrium state that is atypically decoupled from the external drive. However, there is also a heavy tail in the percentile distribution, and a significant fraction of all trajectories terminate at fixed points with percentile ranks above

We also computed the work rate per maximum current, *B*, despite the greater frequency of low-force fixed points. In aggregate, the high-force fixed points thus lead to such a pronounced increase in

Our observations lead us to ask how the dynamics of these random networks become biased toward having fixed points that are thermodynamically special. What is it about these rare configurations—which make up a minuscule fraction of the whole chemical space—that causes them so frequently to turn out to be the inevitable dynamical attractor?

To explain this, we note that in our model, all driven reactions were duplicated from an existing undriven version. Accordingly, we can divide the equations of motion into undriven and driven pieces,

To assess the impact of high forcing, we investigated how it shifts the Pearson correlation (or normalized inner product)

Low-force configurations have a correlation peaked near *B* demonstrates that the magnitudes of these two reaction vectors remain largely unchanged, and correlated, even at high force. With these two different vectors sampled quasi-independently, there is an increased likelihood that they will turn out to be equal and opposite, canceling out to a high-force fixed point. The key to this mechanism is the high dimensionality of the chemical space and the random strength and direction of forcing on reactions.

## Discussion

In the ensemble of random reaction networks studied here, a challenging driving environment made it possible to define fine-tuned order in terms of an atypically strong matching to available sources of work. As anticipated in previous theoretical works (7, 17), our central finding was that kinetically stable behaviors of such a system are indeed biased toward appearing to be finely tuned to the external drive: In other words, the long-time behavior of the system is enriched in outcomes that would be observed only with small likelihood in a random and uniform sampling of the whole space of possibilities.

The mechanism for emergent fine-tuning at high force appears to require that randomly interacting components lead to a high-dimensional space of possible configurations with distinct properties. To maintain a high-force fixed point, the system must experience enough forcing from the environment to give rise to driven reactive fluxes that can exactly counterbalance the undriven reactions occurring alongside them. Whereas such counterbalancing is by no means guaranteed to occur solely as the result of high forcing, the high-force regime is composed of different locations in chemical space that may be thought of as each being a quasi-independent attempt at achieving this counterbalancing by random accident. In a setting where the effective number of such attempts is large enough, the existence of at least one high-force fixed point can be surprisingly likely.

We might have suspected that our initial assumption that forced and unforced reaction branches were governed by the same reaction timescales (

We found in general that, once such an attractive fixed point exists, the system dynamics can find it on the timescale of the slowest unforced reaction. At first, this rate of relaxation to the steady state may seem unexpectedly rapid: A structure of correlations in concentration extending across the entire system turns out not to take significantly longer to develop than the time for a single slow reaction to occur at random. The reason for this effect is that we have modeled this chemical network using mass-action kinetics, which assume a limit of large particle numbers. Thus, implicitly, the system has the opportunity to “experiment” with many different chemical combinations in parallel at the level of individual stochastic molecular events, until a successfully cooperative one is discovered (18, 19).

It is instructive to be able to demonstrate that even a chemical mixture devoid of any obvious selective pressure via self-replication could discover and stay in specially structured combinations that appear to be good at meeting a complex challenge presented by the environment—they can not only do this, but also do so on timescales rapid enough to be accessible to direct observation. With further study, such knowledge might prove to be useful in explaining or promoting the emergence of life-like behaviors in initially lifeless material settings.

## Materials and Methods

Our model is a well-stirred, dilute mixture of

The reaction rates are given by mass-action kinetics (11),

We break detailed balance by introducing additional biased reactions driven by a collection of nonlinear, generalized thermodynamic forces. To this end, we randomly select a fraction

With the rates, the species-concentration dynamics are dictated by the macroscopic reaction rate equations*i*) *ii*) *iii*) *iv*) *v*) *vi*) *vii*) *viii*) *ix*)

Numerical solutions were obtained using the NDSolve function in Wolfram Mathematica version 10.1.0.0, using the “StiffnessSwitching” method with accuracy goal set to infinity. Parameters were as follows: For Fig. 2, *i*) *ii*) *iii*) *iv*)

This section includes additional analyses of the forcing statistics in driven chemical reaction network realizations. In Fig. S1, the dynamics and forcing statistics are presented for different randomly realized networks that exhibit distinct properties in their fixed-point ensembles. Figs. S2 and S3 reproduce the analysis of Figs. 3 and 4 of the main text for an ensemble of chemical reaction networks in which the driven reaction time constants are independently sampled. Figs. S2 and S3 demonstrate that the effects are robust to additional disorder in the reaction rates.

## Acknowledgments

The authors are grateful to Dan Goldman, Pankaj Mehta, and Robert Marsland for thoughtful comments on this article. J.M.H. and J.L.E. are supported by the Gordon and Betty Moore Foundation through Grant GBMF4343. J.L.E. further acknowledges the Cabot family for their generous support of Massachusetts Institute of Technology.

## Footnotes

↵

^{1}J.M.H. and J.L.E. contributed equally to this work.- ↵
^{2}To whom correspondence should be addressed. Email: jengland{at}mit.edu.

Author contributions: J.M.H. and J.L.E. designed research, performed research, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1700617114/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵.
- Milo R, et al.

- ↵
- ↵.
- Seifert U

- ↵
- ↵
- ↵.
- Perunov P,
- Marsland RA,
- England JL

- ↵
- ↵
- ↵.
- Poletinni M,
- Esposito M

- ↵.
- Van Kampen NG

- ↵.
- Nowak MA

- ↵
- ↵
- ↵.
- Esposito M

- ↵
- ↵.
- England JL

- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics