Self-Organized Resonance during Search of a Diverse Chemical Space

Tal Kachman, Jeremy A. Owen, and Jeremy L. England Physics of Living Systems Group, Department of Physics, Massachusetts Institute of Technology, 400 Tech Square, Cambridge, Massachusetts 02139, USA Machine Learning for Healthcare and Life Sciences, IBM Research Laboratory, Haifa 3498825, Israel Department of Physics, Technion-Israel Institute of Technology, Haifa 320003, Israel (Received 25 February 2017; published 21 July 2017)

Matter is termed "active" when it experiences sustained inputs of work energy from chemical, mechanical, or other sources at the single-particle level.Even active matter systems of quite simple composition can exhibit a striking array of complex collective behaviors, some of which mimic dynamical patterns found in living organisms [1][2][3][4][5].Yet while it is known that being driven far from equilibrium is essential to many phenomena that emerge in active matter, it is challenging to make general thermodynamic statements that constrain or explain experimental observations.The need for such principles is felt especially acutely in examples where assembly outcomes are well-adapted to harvesting work energy from an available drive [5,6], suggesting a mechanism of structure formation that might readily be understood in thermodynamic terms.
New theoretical progress [7][8][9][10][11][12] in statistical mechanics has begun to point towards such a mechanism: given many interacting particles described by the configurational coordinates x and driven energy landscape H(x, t), held in contact with a heat bath at temperature T = 1/(βk B ), the relative probability π τ of evolving in time τ from a starting point x to configurations y or z is [7,10] e −βW x→z e −βW x→y (1) where W is the work absorbed from the drive, the * operator indicates time-symmetry-reversal of momenta and drive, and the average . . . is taken over repeated stochastic realizations of the transition from x to y or z (See Supplementary Material for derivation).
It is apparent in the above that several different factors can affect the relative likelihoods of dynamical outcomes in a given system of interest; like in the undriven case, states of lower energy are more likely, and so are states that are more likely to return to the starting state x in * tal.kachman@ibm.comthe same amount of time.Holding these energetic and kinetic factors fixed, however, reveals a new and uniquely nonequilibrium way that an outcome y can be made to have high probability: namely, by ensuring that as many of the likely microtrajectories leading to y from x as possible have exceptionally positive values of dissipated work W . Cases where this effect dominates have the potential to be quite important in systems where the ability to absorb work from an external drive is strongly dependent on system state (as, for example, in the case of mechanical resonance).Under such conditions, it has been proposed that likely dynamical outcomes should exhibit a dissipative adaptation phenomenon, in which their physical properties bear the signature of the highly work-absorbing states the system had to traverse during its history [10,11].
It must be emphasized at the outset that the character of this signature may differ depending on the detailed physical rules of a system under study.The fact that a system adopted highly work-absorbing configurations in the past only sometimes means that its likely outcome state will also be highly work-absorbing; indeed, alongside systems that tend to grow their rates of dissipation in "energy-seeking" fashion over time [5,9], there are also cases (one of which will be discussed below) where driving causes the opposite to happen (see Figure S8) .In either event, however, there is an opportunity to explain properties of the long-time-stable outcomes of driven dynamics once one recognizes that particular work-absorbing states have to be visited on the way into states of high kinetic stability.
We set out to establish and characterize clear examples of such dissipative adaptation in a system where both the dynamical and thermodynamic properties could be monitored and controlled with exquisite precision.Our first and main goal was to find a regime of physical rules in which "energy-seeking behavior" [5] could be observed, and then to investigate the thermodynamics underlying the emergence of this pattern of dynamics.
Accordingly, we devised a simulated "toy chemical" reaction in which a collection of idealized Newtonian par- (a) Twenty particles experience viscous drag as they move in one dimension.The only interaction between any two particles is a bonding reaction that follows two-state Arrhenius kinetics, as shown in (b), whereby a bond forms or breaks at a Poissonian rate determined by the distance ∆x between them.A bonding energy ε controls the tendency to form bonds at random.(c) In the absence of any drive, the particles tend to form a random graph at equilibrium (left).When a single particle (yellow) is subjected to an oscillatory external force, vibrations in the network change the rates of formation and breakage of certain bonds, and new structures are sampled in a nonequilibrium steady-state.A color gradient superimposed on the bonds indicates minimal connectivity distance from the driven particle (with reddest being closest).Particles are depicted with sizes proportional to their number of bonded neighbors in the graph.Produced using the ForceAtlas2 algorithm in Gephi.
ticles randomly form and break harmonic spring bonds with each other while subject to drag forces from a surrounding viscous medium held at temperature T .At each instant, a given pair of particles i and j separated by distance ∆x ij either share a single bond of internal energy 2 − ε, or do not (in which case, E ij = 0).In order to ensure that the model of the surrounding thermal bath is consistent with Newtonian mechanics, the formation and breakage of bonds is assumed to follow a simple two-state stochastic switching process whose bonding and unbonding rates are governed by an Arrhenius [13] is the height of the transition state barrier separating the bonded and unbonded states.The system may therefore be thought of loosely as a group of reacting "atoms" transiting randomly between differently bonded "molecular" states over time under the influence of thermal fluctuations.
The appeal of studying such a system is that we expect different patterns of bonding within a group of atoms to give rise to different mechanical response properties, because the spectrum of normal mode frequencies for the group is determined by the bond connectivity of the system as a whole.Moreover, the number of possible connec-tivities is large; even a group as small as 20 particles can be very roughly estimated to have 1 20! 2 ( 20 2 ) ∼ 10 39 bonding states available to it, comprising a vast space of different ways of hooking springs together to give rise to a distinct vibrational fingerprint.What remains is to determine how this space of fingerprints is explored along typical dynamical trajectories when the system is driven far from thermal equilibrium by a time-varying external field.
To answer this question, we carried out simulations [14] of 20 particles confined in a shallow, one-dimensional harmonic well under conditions of under-damping viscous drag, low temperature, and strong, sinusoidal forcing [15,16] of a single particle (labeled with index 0).In this setting, the diffusive motions of the particles due to fluctuating random forces from the thermal bath can be ignored (see Supplementary Material), and the i th particle obeys between stochastic events that form or break bonds.The Laplacian matrix L ij (t) depends on which particles are connected by springs at each moment in time, and changes according to the sequence of "reactions" between particles.The Arrhenius law governing such reactions contains an arbitrary rate constant r 0 , which we assume to be much less than the drive frequency so that the drive undergoes many full oscillations between each bonding or un-bonding event.In this regime, the driven long-time dynamics of the oscillator network for a given realization of L ij (t) has an exact analytical expression, from which the inter-particle distances ∆x ij (t) can be computed.Thus, the stochastic evolution of the spring connectivity is a Markov process whose transition rates to new bonding states are determined by the driven vibrational dynamics of the current graph connectivity.
To complete the dynamical rules for the system, it is necessary to specify how the transition state barrier heights B(∆x) vary with particle configuration.To start, we chose this function to have the simple, "catch bond" form B(∆x) = k∆x 2 to ensure that particles react only when they are close together in space.It should be noted that, while catch bonds can be realized on the mesoscale and above in various ways, and in general are a physically consistent way of describing interparticle interactions, they do not reproduce the nanoscale intuition that bonds under strain break more rapidly.Accordingly, the particle clusters in this case should only be thought of as "molecules" made of "atoms" to the extent that they exhibit diverse physical properties in novel combinations.
Using Gillespie's algorithm, we generated a large ensemble of trajectories which were initialized in the same connectivity state and then driven at frequency ω d for many iterations.At each iteration, we computed the normal modes of L ij and the associated rate of heat dis- sipation Q. See the Supplementary Material for details. Figure 2a plots the statistics of normal mode frequencies observed for driven and undriven (equilibrium) bonding connectivities.In the absence of drive, the thermal equilibrium ensemble of random graphs has a normal mode spectrum stereotypically shaped like a haystack [17], with a finite gap separating the lowest non-trivial vibrational frequency from zero.The striking effect when a sinusoidal driving force was applied to one of the particles in the network was that the steady-state shape of this normal mode spectrum changed, exhibiting an additional peak jutting out of the haystack at the drive frequency ω d .The emergence of this peak indicates a nonequilibrium enrichment in structures whose connectivities have natural frequencies of vibration that are better matched with the frequency of the drive.The peak generally appeared when the system was driven at frequencies in the haystack (Fig. 2b, S9), and in a parameter regime where damping was sufficiently low for most modes in the system to be under-damped.We further found that the reshaping of the normal mode spectrum brought about by these dynamics resulted in an altered rate of work absorption and heat dissipation by the system; in principle, this need not have been the case, since the normal modes with natural frequency matched to the drive might not have been strongly coupled to the motion of the driven particle, making resonance impossible.Fig. 3b establishes that, although the average rate of work absorption in the driven ensemble relaxes to steady-state more slowly than the shape of the normal mode spectrum, it is in any event always higher than for a random (thermal equilibrium) ensemble.Moreover, plotting the statistics of applied force amplitude in the normal mode basis as a function of frequency revealed that the spring networks that evolved in the presence of a drive were even biased towards connectivities that distributed the most drive amplitude to those modes with the greatest capacity to resonate, indicating a strong pressure on the networks to become organized in ways that increase work absorption (Fig. 2c).
Given the remarkably different physical properties of the driven nonequilibrium spring network ensemble relative to its equilibrium counterpart, we undertook to characterize the topologies of these networks to understand what structural differences enable them to be more "finely-tuned" to the drive than random graphs.In general, driven and equilibrium graph ensembles had different average numbers of bonds (Supplementary Material).However, the main effect of this coarse structural difference was simply to shift the center of the normal mode spectrum slightly in the direction of the drive frequency (Fig. 2a); thus, the change in bond number could not explain the appearance of a new peak in the spectrum precisely at the chosen value of drive frequency.
By eye, the topologies of these graphs were difficult to distinguish from random graphs obtained for an equilibrium ensemble of the same average bond number (Fig. 1a), except in the special case of ω d = 1.Although the driven ensemble generally had a different average number of bonds per particle than the equilibrium ensemble in the range of drive frequencies where adaptation was observed, the edges per node distribution of the driven ensemble of graphs was not qualitatively different than that of an equilibrium ensemble with the same average number of edges.However, we then computed the largest eigenvalue of the adjacency matrix for each network and compared across groups of fixed number of bonds in order to quantify the degree of interconnectedness of the whole network.By doing so, we discovered a noticeable decrease in interconnectedness of driven ensembles, reflecting system-wide topological correlations that are less likely to appear in the equilibrium ensemble (see Supplementary Material).In order to understand the dynamical mechanism for the observed self-organization, we characterized the thermodynamic properties of the driven steady-state (Figure 3).As one might expect, the increased resonance of the driven ensemble leads to a rate of work absorption far above that of the random graph ensemble that is sampled at thermal equilibrium (Fig. 3b).What causes this "energy-seeking" behavior [5]?By scatter-plotting the rates of work absorption for graphs at adjacent iteration steps in dynamical trajectories sam- pled from the driven steady-state, we uncovered a strong tendency in the spring networks that already had atypically high Q to transit to other states with similarly high values (Fig. 3a).Thus, the adaptive resonance observed with driving resulted from the extra kinetic stabilization of high Q structures achieved through driving.
The intuition behind this effect is that resonant structures stretch out more during a drive cycle, and thus the particle motions most associated with resonance tend to be associated with springs that spend more time in a stretched state.Since catch bond springs like the ones simulated here are less likely to break while stretched, resonant structures stay trapped in a high Q part of configuration space that is only sampled very rarely at equilibrium.
The emergence of resonance in this particular "toy chemistry" scenario is surprising, but can it be explained as an instance of dissipative adaptation?The signature of a process of dissipative adaptation is the changes in system configuration that are more irreversible (and thus give rise to the observed far-from-equilibrium organization) occur when the system passes through or into states that can absorb and dissipate especially large amounts of drive work [10,11].
To investigate this question, we coarse-grained the space of graph topologies into three macrostates (A, B, and C) according to whether their rate of work absorption Q at one particular drive frequency lay between certain thresholds.By following long trajectories under driving, we sampled many transitions between each pair of macrostates in the steady-state, and computed the stochastic entropy production [18] (∆S = ln r fwd rrev ) which in this case sets a lower bound on the heat dissipated during the transition.As Figure 3c shows, transitions that led to an increase in Q had a very similar distribution of ∆S, skewing significantly towards positive values.In contrast, transitions that brought about a decrease in Q were less likely to be highly irreversible, and thus had less positive values of ∆S.
Motivated by the dissipative adaptation hypothesis, we sought to test whether the difference in irreversibility between these two types of transitions might originate in their differing tendency to harvest work from the drive.In making this comparison, it is important to note that it is the extra work above the basal rate of Q absorbed during the transient relaxation from one graph connectivity's particular solution to another's that can contribute to statistical irreversibility in transitions between connectivities [11,19].As Figure 3d shows, transitions that are accompanied by a rise in the dissipated power Q do indeed tend to result in elevated levels of transient work being absorbed and dissipated during the transition.We surmise that the observed statistical tendency of the system to adopt structures with higher-than-equilibrium rates of work absorption does follow a dissipative adaptation mechanism, whereby the highly irreversible transitions that sustain the system's nonequilibrium bias towards resonant structures occur because the resonance helps them harvest more work from the external drive.
The self-organized energy-seeking behavior in this particular system is striking, but how general is it?Put another way, why does the fundamental relationship between statistical irreversibility and the absorption and dissipation of work not always imply a tendency towards elevated levels of work absorption through self-organized resonance?
It is best to answer this question by modifying the model to invert the observed effect.In this particular system, it turns out that the use of catch bonds that are kinetically stabilized by stretching is essential to the mechanism of increased resonance from driving.To show how the specific physical characteristics of a nonequilibrium steady-state can turn out differently under alternate conditions, we simulated the same 20 particle system with a different energy landscape of transition state barriers, so that the stretching of formed bonds gave rise to an increased (rather than decreased) rate of spontaneous snapping.Under these new conditions, the steady-state normal mode spectrum of the driven ensemble exhibits the opposite behavior: modes with frequencies close to resonance with the drive are strongly suppressed, while those that resonate less remain (see Figure S8).Thus, the steady-state dissipated power of the driven ensemble turns out to be lower than in the equilibrium state.One may think of the driven steady-state in this case as the "shattered" remnants of structures that resonated, absorbed work from the drive, and broke apart irreversibly.
Taken together, these findings illustrate that both spontaneous self-organization into states that absorb extra drive work, and ejection into shattered states that absorb less, are scenarios that can be made physically sensible.In a given nonequilibrium system of it may be the case that both such effects are at play to some extent [4]; here, however, the self-organized resonance we have characterized provides a clean, idealized in silico illustration of a potentially more general thermodynamic mechanism for emergent energy-seeking behavior in a many-body mixture with a diverse space of "chemical" combinations.We are thus encouraged to hope that future investigations of dissipative adaptation in more complex self-assembly reactions far from equilibrium may go significantly further towards a physical elucidation of life-like patterns of collective molecular behavior.

FIG. 1 .
FIG. 1.(a) Twenty particles experience viscous drag as they move in one dimension.The only interaction between any two particles is a bonding reaction that follows two-state Arrhenius kinetics, as shown in (b), whereby a bond forms or breaks at a Poissonian rate determined by the distance ∆x between them.A bonding energy ε controls the tendency to form bonds at random.(c) In the absence of any drive, the particles tend to form a random graph at equilibrium (left).When a single particle (yellow) is subjected to an oscillatory external force, vibrations in the network change the rates of formation and breakage of certain bonds, and new structures are sampled in a nonequilibrium steady-state.A color gradient superimposed on the bonds indicates minimal connectivity distance from the driven particle (with reddest being closest).Particles are depicted with sizes proportional to their number of bonded neighbors in the graph.Produced using the ForceAtlas2 algorithm in Gephi.

FIG. 2 .
FIG. 2. (a)The normal mode spectrum P(ω) is sampled for an equilibrium ensemble of spring networks (blue) and exhibits a stereotypical haystack shape.When driven at frequency ω d = 1.2, the normal mode spectrum rearranges (red) in two noticeable respects: First, the average number of springs decreases; indeed, an equilibrium normal mode spectrum computed for a value of ε that leads to the same average number of bonds (orange) has a very similar shape.However, in the nonequilibrium spectrum, there is an additional peak at the chosen drive frequency.(b) The drive induces a peak in the nonequilibrium normal mode spectrum at any chosen drive frequency within the haystack region of the spectrum.(c) The average forcing magnitude distributed to normal modes of a given frequency is plotted for different drive frequencies.The diagonal line indicates that modes with greater resonant response to the drive frequency tend to be more strongly forced.

FIG. 3 .
FIG. 3. (a) Scatter plot of rates of dissipation Q for spring networks at adjacent iteration steps along a driven Markov microtrajectory of the network graph.States with atypically high values of Q tend to jump to states that retain this property, indicating that this subset of spring networks is kinetically stabilized by driving.(b) Mean dissipation rate Q over time of the ensemble of random graphs, initialized in the equilibrium and driven with an oscillatory drive.Across a wide range of drive frequencies ω d , the dissipation rate Q increases monotonically from its value in the initial equilibrium distribution (c) Distributions of stochastic entropy ∆S = ln r fwd rrev for transitions in the driven steady-state.Markov transitions in spring network connectivity that lead to an increase in Q (red) show a significantly greater tendency toward positive ∆S than those that lead to a decrease, and are thus more irreversible.(d) Distributions of transient work ∆W for transitions in the driven steady-state.This quantity measures the extra work performed by the drive during dissipative relaxation after formation or breakage of a bond.Transitions that lead to an increase in Q typically absorb and dissipate more work from the drive relative to those that lead to a decrease in Q.