Typical Stochastic Paths in the Transient ... - American Chemical Society

May 7, 2019 - (e.g., fiber length, number, connectivity) from those of the ...... accurately identify the time of extremal values in the number and le...
1 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B 2019, 123, 4792−4802

Typical Stochastic Paths in the Transient Assembly of Fibrous Materials Published as part of The Journal of Physical Chemistry virtual special issue “Deciphering Molecular Complexity in Dynamics and Kinetics from the Single Molecule to the Single Cell Level”. Schuyler B. Nicholson and Rebecca A. Bone

Downloaded via NOTTINGHAM TRENT UNIV on August 15, 2019 at 00:38:11 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Department of Chemistry, University of Massachusetts Boston, Boston, Massachusetts 02125, United States

Jason R. Green* Department of Chemistry, University of Massachusetts Boston, Boston, Massachusetts 02125, United States Department of Physics, University of Massachusetts Boston, Boston, Massachusetts 02125, United States Center for Quantum and Nonequilibrium Systems, University of Massachusetts Boston, Boston, Massachusetts 02125, United States ABSTRACT: When chemically fueled, molecular self-assembly can sustain dynamic aggregates of polymeric fibers hydrogelswith tunable properties. If the fuel supply is finite, the hydrogel is transient, as competing reactions switch molecular subunits between active and inactive states, drive fiber growth and collapse, and dissipate energy. Because the process is away from equilibrium, the structure and mechanical properties can reflect the history of preparation. As a result, the formation of these active materials is not readily susceptible to a statistical treatment in which the configuration and properties of the molecular building blocks specify the resulting material structure. Here, we illustrate a stochastic−thermodynamic and information−theoretic framework for this purpose and apply it to these self-annihilating materials. Among the possible paths, the framework variationally identifies those that are typicalloosely, the minimum number with the majority of the probability. We derive these paths from computer simulations of experimentally-informed stochastic chemical kinetics and a physical kinetics model for the growth of an active hydrogel. The model reproduces features observed by confocal microscopy, including the fiber length, lifetime, and abundance as well as the observation of fast fiber growth and stochastic fiber collapse. The typical mesoscopic paths we extract are less than 0.23% of those possible, but they accurately reproduce material properties such as mean fiber length. select,22 self-heal,9 self-erase,23 self-organize,11 and self-replicate.24 Underlying these behaviors of nonequilibrium materials is that they depend not only on the structure and morphology but also on the history and formation mechanism.25 Out of equilibrium, it also remains a challenge to self-assemble arbitrary structures with prespecified features in high yield.26−29 Recent experiments have demonstrated syntheses where the structures formed are sensitive to the pathways that building blocks traverse.30,31 For instance, changing the solvent environment can influence assembly pathways and produce peptide amphiphiles with different morphologies32 and imposing a multistep assembly process can reduce the number of allowed assembly pathways, giving control of multicomponent micelle formation.33 Collectively, these results suggest the material

I. INTRODUCTION Dissipative self-assembly1 can create, sustain, regulate, and destroy the structure of supramolecular materials.2−4 Because of this ability to dynamically control structure, there have been efforts to synthesize molecular architectures that are responsive to chemistry,5 light,6 and pH.7−10 For example, the chemically fueled assembly of hydrogels endows these materials with tunable properties that make them potential drug delivery systems.11 However, these fibrous gels only exist while they consume fuel and dissipate energy, which raises fundamental questions about how to predict the material structure from the nature and interactions of the basic building blocks. This fact, coupled with a growing awareness of mechanisms far from equilibrium,1,12−16 highlights the challenge of choosing building blocks for supramolecular assembly with a well-defined set of pathways that lead to a target structure.17−21 Away from thermodynamic equilibrium, self-assembly can produce structural phases that are able to self-annihilate,9 self© 2019 American Chemical Society

Received: March 25, 2019 Revised: May 6, 2019 Published: May 7, 2019 4792

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

produces an active material but only while dissipating energy. Chemical kinetics and pH control the formation and annihilation of these fibrous materials and their mechanical properties.11 Boekhoven et al.9 use inactive precursors, such as N,N′dibenzoyl-L-cystine, which are alkylated at a carbonyl site. Alkylation by dimethyl sulfate [DMS, (CH3)2SO4], the chemical fuel, forms the active material; it removes the repulsive interactions and transforms inactive monomers into active monomers that can assemble into fibers. Under basic conditions, however, this activation reaction competes with deactivation by ester hydrolysis. Deactivation of active monomers can occur in solution and in fibers, which can destabilize and dismantle as a result. During the whole assembly process, the system passes through a succession of nonequilibrium states where fibers dominate the observed structure. Individual fibers, though, are dynamic, with some growing and others stochastically collapsing. Without fuel, fibers cannot persist. Here, we model this dissipative self-assembly process,9 incorporating both experimentally informed chemical kinetics and stochastic fiber dynamics. Two mechanisms are necessary to analyze the sequence of configurations traversed during the dissipative assembly process: the chemical kinetics, which evolves the concentrations of molecules in time, and the physical kinetics of the assembly, which controls the growth and decay of fibers. The chemical kinetic component is the set of chemical reactions and their associated rate coefficients:

structure may be controlled by designing molecular building blocks that follow particular pathways.32−37 However, away from equilibrium, it is not entirely clear how the thermodynamic and kinetic forces under experimental control, if they are known,19 imprint on the dynamic histories of self-assembly processes.18 In dissipative materials, including self-annihilating hydrogels, there have been efforts to understand how external energy sources control the form and the dynamic function of emergent structure.38 However, the basic principles of statistical dynamics needed to describe the dissipative nature of these syntheses and the transient existence of structure are still under development. Most lacking is a theory to predict yields and material properties (e.g., fiber length, number, connectivity) from those of the constituent building blocks. Experimental and theoretical work has led to some design principles and the assembly of intricate three-dimensional structures.39,40 However, most advances are grounded in an equilibrium thermodynamic framework,41 relying on the formation of structure through the paths carved out by thermodynamic gradients and that lead to a free energy minimum. Far less is understood far from equilibrium. The challenge in theoretically treating the kinetic pathways that lead to assembled structures is the vastness of the space in which they reside. Put simply, it is necessary to know which pathways of assembling systems are statistically significant and which are critical to controlling the products that are accessible in high yield. Our goal here is to determine which sequences of configuration states or paths are typical and which are rare, using stochastic chemical kinetics that model experiments.

k0

DMS + H 2O → MMS− + MeOH + H+

II. MODEL FOR THE DISSIPATIVE SELF-ASSEMBLY OF HYDROGELS As a proof-of-principle, we consider the dissipative assembly of hydrogels introduced by Boekhoven, Hendriksen, Koper, Eelkema, and van Esch.9 The assembly process is chemically fueled, operating through the action of multiple reactions on a switchable molecular building block. A generic scheme is shown in Figure 1. When a chemical fuel is available, monomeric building blocks can interconvert between an inactive precursor (Mi) and an active state (Ma). Molecules in the active state can self-assemble into a supramolecular material. Their deactivation by a separate reaction produces waste. The overall reaction cycle

k1

DMS + M i → +MMS− + M a k2

Ma + OH− → + M i + MeOH

(1a) (1b) (1c)

We use rate coefficients for the molecular gelator determined through HPLC and 1H NMR in ref 9 (Table S3, pH 11). The background reaction of the DMS fuel with H2O, reaction 1a, dominates early in the process. As the fuel is consumed, the activation reaction, reaction 1b, becomes more probable, and inactive molecules Mi transition to active molecules Ma. Active molecules can also react with OH− to return to the inactive state through the deactivation reaction, reaction 1c. Together, these reactions and their rates describe the concentrations of all molecules in time but not the growth and decay of fibers. To model the assembly of activated monomers into fibers, we add “chemical species” to this mechanism for active monomers that are bound to fibers of varying length. We denote these species Ma ⊗ Fj to distinguish active monomers in fibers from those in solution, Ma. These species, and the population of fibers, are necessary to model the chemical reaction network and its influence on the material properties. Our simulations track individual molecules undergoing reactions and their location in solution or fibers. Each fiber we track is assigned to a bin spanning a range of fiber lengths. The set of N − 2 (here N = 20) bins Fj = jΔF are equally spaced between 2 and 15,000 molecules. The range of fiber bins is fixed in time but depends on the volume of simulations. These coarse-grained fiber states make the total set of states x ∈ X, with

Figure 1. A schematic chemical reaction network capable of driving the assembly of a supramolecular material from a switchable molecular building block, M. An essential design element is the two irreversible reactions. One reaction consumes a chemical fuel to activate molecular precursors Mi (e.g., by methylation of a molecular gelator). The other reaction reverts active building blocks Ma into the precursor (e.g., ester hydrolysis). Only activated building blocks can assemble into polymeric fibers. The net result of the process is the conversion of fuel into waste and the transient existence of structure.

X = {M i, M a, Ma ⊗ Fj}

(2)

and j = 1, 2, ..., N − 2. These extra states contribute to the chemical kinetics through two additional reactions: 4793

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B k2

Ma ⊗ Fj + OH− → M i + M a ⊗ Fl + MeOH k3

2Ma → Ma ⊗ F1

(3a) (3b)

The first reaction is the deactivation of molecules that belong to fibers, akin to the deactivation of free active monomers in reaction 1c. It accounts for active molecules in a fiber Ma ⊗ Fj becoming inactive Mi and the fiber in bin Fj (potentially) transitioning to a different fiber bin Fl. The second reaction accounts for the birth of new fibers. New fibers consist of at least two uncharged (active) molecules bound together. The rate constant k3 is an unknown parameter that in practice can be estimated from collision theory.42 We set k3 = 3 × 10−3 M−1 s−1 to be consistent with the experimentally determined rate coefficients. The experimental rate coefficients, along with k3, are held fixed throughout. Together, reactions 1 and 3 are a minimal mechanism for the stochastic chemical kinetics of the fuel-driven assembly of molecular gelators. Fiber growth, decay, and lifetime are pH-dependent. To account for variable rates of fiber growth and decay, we expanded the model with exponential growth rates of the fiber length F. From the time of fiber birth tg0, each fiber grows as

Figure 2. Trajectory-averaged concentration (mM) of inactive building blocks (Mi, blue), free active building blocks (Ma, orange), and all N−2 building blocks bound in fibers (∑ j = 1 M a ⊗ Fj , gray) as a function of time (h) from stochastic simulations of the chemical kinetics. Data are an average over 6 × 104 trajectories.

activation and deactivation are rare, and the inactive pool approaches its original size.

III. RESULTS AND DISCUSSION A. Material Properties. Three observables of the fibers are of immediate interest: the concentration of monomers (inactive, free active, and fiber-bound active), the number of fibers, and the length of fibers. Monomers can interconvert between inactive Mi and active Ma forms, but the total number of monomers is conserved. First, the average concentration of the species m is

g

F(t ) = 2e(t − t0 )/ τg

(4)

The coefficient of 2 represents the initial number of molecules in a fiber. If the deactivation reaction, reaction 3a, occurs at a time t0d, a fiber is randomly selected and subsequently decays according to

⟨Cmt⟩ =



∑ χ(mlt , m)

(6)

l=1

We define the number of molecules in all simulations N̅ = NsNm from the number of molecules Nm and the number of simulations N s . Here, N A is Avogadro’s number and

d

F(t ) = F0d e−(t − t0 )/ τd

1 NsNAV

(5)

where Fd0 is the length of the selected fiber at td0. We set τd = τg/ 300 to agree with ref 9. The time scale τg modulates the growth and decay of fibers, while the chemical kinetic model governs the pools of molecules that fuel fiber formation. We chose the rate coefficient for fiber creation k3 and the fiber growth τg and decay τd parameters to match experimental time scales. To model the effect of pH on the material properties and lifetime, we vary τg = {35, 63, 70, 92, 150}. Observations of the dissipative assembly process using confocal microscopy show three regimes of behavior: an initial growth in the number of fibers under fuel-rich conditions, a period of simultaneous growth and stochastic collapse, and, finally, decay of the material upon depletion of the fuel DMS. To reproduce the stochastic and ephemeral nature of the fibers, we simulate the model mechanism with kinetic Monte Carlo using the Gillespie algorithm.42 In the infinite sample limit, the probability of a mixture composition is governed by the chemical master equation.43 Each simulation has an initial 50 mM concentration of inactive monomers, Nm = 301,150 molecules, and a simulation volume of V = 1 × 10−17 L. Simulations of this mechanism reproduce chemical features of the transient fibrous materials observed in experiments.9 Figure 2 shows concentration profiles (averaged over an ensemble of 6 × 104 trajectories) for inactive monomers, free active monomers, and active monomers in a fiber. Data shown are for τg = 70. The conversion of inactive monomers to active monomers leads to a growth in the pool of active molecules, both those that are free and those that are bound. The fuel concentration (not shown) is consumed steadily, and its concentration decays exponentially. Around 3 h after the start of the reaction, the fuel is nearly depleted, chemical cycles of

mlt ∈

{M , M , ∑ M i

a

j

a

}

⊗ Fj refers to the lth inactive or active

monomer. The indicator function χ is one if mtl is in state m and zero otherwise. Second, given the number of fibers in the jth bin at time t, |F jt |, the average ⟨NFt ⟩ over an ensemble of Ns simulations is ⟨NFt ⟩ =

1 Ns

Ns N − 2

∑∑

F jt (l) (7)

l=1 j=1

The average is over all simulations of the assembly process. Third, we define the average fiber length from a molecular perspective. The average is the typical length of fiber to which a molecule belongs: N−2

⟨Fmt ⟩ =

∑ Pr(m ∈ F jt )F jt j=1

=

1 N̅

N − 2 N̅

∑ ∑ χ(mlt , F jt )F jt (8)

j=1 l=1

Now, χ(mtl, Ftj) indicates whether molecule mtl belongs to in bin Ftj. These three observables are readily calculated

a fiber in our

simulations. The parameter τg directly controls the fiber growth rate and the decay rate through our definition of τd. Though indirect, τg also helps to control the length and lifetime of fibers by influencing the size of the active monomer pool. Regardless of the value of τg, the average number of fibers grows to a global maximum, undergoes weak oscillations as fibers grow/collapse, 4794

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

Figure 3. (a) The average number of fibers ⟨NtF⟩ as a function of time (in h) oscillates between periods of net growth and net decay. The growth rate parameter τg controls the maximum number of molecules in fibers assembled. (b) Up to the maximum, there is apparent power law growth ⟨NFt ⟩ ∝ t 2.8, as shown by ln⟨NFt ⟩ against ln t. The initial rate of fiber birth (and the mean number of fibers) is independent of τg, which instead determines the time and magnitude of the maximum number of average fibers. Comparing τg = 35 to τg = 150, faster fiber growth suppresses the maximum number of fibers formed. (c) By contrast, the first peak in the average fiber length does not grow monotonically with τg. In all cases, averages are over at least 6 × 104 stochastic trajectories.

Figure 4. (a) State diagram of a time-inhomogeneous Markov model for the fuel-driven, dissipative assembly of hydrogels. Molecules can occupy states that are inactive or active. Active states are further partitioned between a free state and k = 1, 2, ..., N − 2 bound states. Each bound state, Ma ⊗ Fk, represents active molecules in fibers with a range of lengths. Larger k values correspond to fiber states with longer fibers. Empirical transition probabilities for the Markov model are estimated from counts in the ensemble of stochastic trajectories. (b) The assembly process is irreversible: the instantaneous, global current shows clear deviations from zero as a function of time (in h) for the time-inhomogeneous Markov state model. At long times (>2 h), the system approaches equilibrium, where there are few fibers and the fuel level is low. By t ≈ 3 h, there is zero current, and detailed balance is satisfied. At intermediate times, however, fibers simultaneously grow and collapse. There are periods of positive (negative) current that correspond to net fiber growth (collapse). The magnitude of the nonzero current measures the degree to which the kinetics break detailed balance. Snapshots of the system are a projection of the (0 + 1) kinetics onto a 2D lattice to illustrate the fibrous structures that form (τg = 70).

The average fiber length shown in Figure 3c exhibits a more complicated behavior than the average number of fibers. For τg = 35, the rate of fiber growth significantly outpaces the rate at which the active pool is replenished. The disparity between these rates tends to result in only a few fibers being born. These fibers grow until the active pool is depleted. Only once the active pool is replenished will more fibers be spawned. When τg is large, such as τg = 150, the delayed exhaustion of the active pool gives time for generations of new fibers to be born and grow, which creates a broad distribution of fiber lengths. For τg = 63, 70, and 92, the fiber growth rate does not exhaust the active pool as quickly as τg = 35, and there are fewer fibers on average than τg = 150. These rates of fiber growth produce longer fibers (at their peak length) than both τg = 35 and τg = 150. The transient activation of the cycles stimulating fiber growth and collapse present a theoretical challenge: making predictions of material properties, such as the longest fibers that will form or the peak number of fibers, without asymptotic formulas or assumptions of steady state. The steady state for this system is reached once the fuel is consumed and the gel has dissipated, leaving an isotropic solution.

and, ultimately, decays to zero, Figure 3a. Reactions 3a and 3b both draw from the same pool of active monomers. However, the log−log plot in Figure 3b shows that the rate of fiber birth is independent of τg up to the maximum of the mean number of molecules in fibers. Comparing to the experimental results in ref 9, we see more pronounced oscillations in fiber length. From our data, the global maximum occurs when the fuel level is still high and at a time when there is both a depleted active pool and a large number of fiber-bound monomers. With the pool of active monomers exhausted, the number of fibers undergoes a period of decay as active monomers in fibers are deactivated, reaction 3a. Deactivation generates inactive precursors, which can be activated again by DMS fuel and continue the reaction cycle. The oscillations in the mean number of fibers are the result of these periods of depleted active molecules and fiber decay, followed by a buildup of free active monomers and fiber growth. When the fuel is depleted, fibers and reactions become rarer, and by around 3 h, the fuel is mostly consumed, signaling the end of the reaction. Overall then, these pronounced oscillations are due to the finite number of molecules and finite simulation volume. 4795

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

C. Measuring Irreversibility of Dissipative Assembly. The fibrous hydrogel forms by expending chemical fuel and producing waste. The energy dissipated drives the system through nonequilibrium states.45 The irreversibility of the assembly process is readily apparent in the current from stochastic thermodynamics.46 Using the time-dependent transition matrices, Rtxy, the current Jtxy is the difference between the single-step joint-probability of transitions forward and backward between x and y

B. Assembly Pathways and Their Probabilities. With the experimentally informed model of the kinetics and quantitative measures of the material properties, we sought to identify the sequences of configuration states or paths that are typical (and, consequently, which are rare). At a high level, we first introduce a coarse-grained symbolization to analyze paths. Paths, here, are generated by each molecule in the assembling system. Each molecule traverses a sequence of symbols during the assembly process as it transitions between coarse-grained states (inactive, free active, and the N − 2 fiber-bound active molecule states). From these state-to-state transitions and ensembles of kinetic Monte Carlo simulations, we populate a discrete-space, discrete-time, time-inhomogeneous Markovstate model.44 From the transition matrices, we can compute and analyze the coarse-grained dynamics with our informationtheoretic framework. More specifically, molecular-level paths are generated by tracking molecules as they transition between coarse-grained states xt ∈ X, eq 2, illustrated in Figure 4a. Over the course of an assembly process, each molecule will pass through a sequence of states xt̂ ≡ xt , xt − δt , ..., xt0+ δt , xt0

t t Jxyt = R xy py (t ) − R yx px (t )

The current for each state is an element of the matrix Jt. Because this matrix is antisymmetric, summing over the lower diagonal, J t = ∑x < y Jxyt , gives the total current. The process is away from equilibrium when the current is nonzero.47 The total current Jt is a global measure of the degree to which detailed balance is broken. Despite Jt measuring irreversibly, a larger value of Jt does not imply a greater distance from equilibrium because Jt is not a metric.48 Figure 4b shows the current as a function of time for our model system. To visualize the fibers, we map the molecule counts from stochastic kinetics simulations to a 2D lattice with periodic boundary conditions. Initial fiber positions are randomly selected and grow according to eq 4. We do not enforce area exclusion, so fibers can overlap. Figure 4 shows configurations from one reaction where the number and length of fibers change substantially during the nonequilibrium portion of the reaction. These configurations illustrate that fibers grow and decay when Jt ≠ 0. Once the system has consumed most of the fuel, t > 3 h, the current goes to zero. While the size of the current is not a quantitative measure of distance from equilibrium, the sign of Jt does give additional information. In our convention for the stateto-state transitions, fiber growth is dominant when the global current is positive Jt > 0, though individual fibers may be growing or collapsing. When Jt < 0, fibers are collapsing on average. The oscillations in Figure 4 between positive and negative current mark the cycles from net growth to net decay. Overall, the current demonstrates that the assembly and disassembly of fibers driven by the chemical reaction network is a timedependent nonequilibrium process. D. Variationally Typical Assembly Paths. Relative fluctuations of macroscopic quantities at equilibrium are exceedingly small,49 and the ensemble average coincides with the most probable or “typical” value. For transient or small systems, fluctuations need not be negligible. The question is then how to predict material properties that are average observables from the myriad paths in the dissipative assembly of the hydrogel? Despite significant fluctuations, however, there is often still a concentration of probability in a typical set of paths. This concentration of probability makes it possible to calculate macroscopic quantities by averaging over only a subset of paths. In other words, the problem of efficiently computing macroscopic observables comes down to finding the set of typical paths. The asymptotic equipartition property (AEP) states that infinite sequences of random variables often converge to a typical set of high probability.50,51 While the asymptotic results are elegant and powerful,52−56 for the transient systems here, we will use the less emphasized fact that concentration to a typical set can occur for all length sequences. The asymptotic typical sequences all converge to the same probability, e−nhμ, meaning they are a uniform subset of all sequences with the asymptotic

(9)

with x̂t ∈ X̂ t. The joint probability μ(x̂t) that a molecule follows the path x̂t ∈ X̂ t out of the set of possible paths up to time t, X̂ t, is normalized so that ∑x ̂ μ(xt̂ ) = 1. Knowledge of μ(x̂t) enables t direct calculation of instantaneous or path averages of observables and all marginal distributions for all times up to t. We employ a Markovian approximation of the joint probabilities for paths: ρ(xt̂ ) = R xt −txtδ−tδtR xt −t−2δtδxtt−2δt ... R xt0t +δtxt p(xt0) 0

(10)

0

(11)

t

The time-dependent transition matrix R is populated from tracking individual molecules in an ensemble of st o c has t ic a s se mb ly pr oc e s s e s . E a c h e l e m e n t i s R xt −ty δt ≡ Pr[X t = xt |X t − δt = yt − δt ] and normalized with the t − δt

convention ∑x R xt −ty δt = 1. The marginal probability px(t) is t − δt t normalized ∑x px(t) = 1 and found through px (t ) = ∑y R xt −ty δt py (t − δt ). We monitor the assembly process t − δt

in each simulation over 3.6 × 104 s, and we partition this time into 1 × 104 time windows of δt = 3.6 s. If all entries on a column of Rt are zero at any time, we set the diagonal elements to be one, which ensures probability conservation. Even with a Markov state model, it is challenging to estimate the probability of a path that both spans the entire assembly process and resolves the transient features of the fiber dynamics. To overcome this challenge, we estimate joint probabilities over coarse-grained time windows, τ = t − t0. We apply a sliding window to extract paths of length τ and calculate the joint probability for all paths ρτ(x̂t) = Pr(xt, xt−δt, ..., xt0). The subscript on ρτ(x̂t) indicates that the joint probabilities are for paths over a window τ. For times t > τ, we slide the joint probability forward one time step, t′ > t, ρτ(x̂t′) = Pr(xt′, xt′−δt, ..., xt′−τ). A sufficiently long block length will show the effect of the chemistry on the paths while still being short enough for the joint distribution to be computationally tractable. Here, we fix blocks to be n = 10 time steps, or τ = 10δt = 36 s of the 36,000 s assembly process. With these details, we analyze the relationship between pathways at this coarse-grained level and the inherent irreversibly of the assembly process.

4796

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

36 Figure 5. (a) The size of the variational typical set |Aϵ* | for τg = 35, 70, and 150 undergoes periods of expansion and contraction over the assembly process (color legend in part d). Less than 0.23% of all nonzero paths are typical, at most a couple thousand, but they hold >99% of the joint probability. (b) The joint entropy rate hτ measures the mean uncertainty about the next state, given the previous τ = 36 s history. The greater the diversity in length and number of fibers, the greater the uncertainty in the average sequence of states a molecule will follow. (inset) For times after the maximum, the entropy rate decays at a near exponential rate as the reaction slows and the number of paths with nonzero probability decreases. (c) Average concentration of inactive precursors and active building blocks in free and fiber bound states from the stochastic simulations (solid black) and from the Markov model using only typical paths (dashed color). (d) Average fiber length from stochastic simulation (solid black) and computed using only typical paths from the Markov model (dashed color).

entropy rate hμ.57 Joint probability distributions of finite-length typical sequences, however, are often not uniformly distributed. They instead lie within the exponential neighborhood of a finitetime entropy rate hτ(t). What is needed to approximate averages is this neighborhood and the sequences within itthe finite typical setthat hold a majority of the probability. To make the typical set precise for finite-length sequences, the joint entropy rate at t over a block length τ is defined as hτ (t ) = −

1 τ

We use the normalized size of the typical set W ϵt = |Aϵt | /|Ωt | and the number of all paths with nonzero probability at t, Ωt. The optimal ϵ = ϵ* in eq 13 yields the variational typical set Aϵτ*. The size of the typical set | Aϵτ*| is the number of paths that are necessary for an accurate (high probability) description. The size of the typical set for a block length of τ = 36 s is shown for three values of τg = 35, 70, and 150 in Figure 5a. There are at 36 |, corresponding to ≤0.23% of most a few thousand paths in | Aϵ* all paths with nonzero probability. Despite their small number, these typical paths account for >99% of the joint probability. If probability is concentrated onto a subset of states, only a portion of the distribution is necessary to estimate averages. Using the paths in the typical set, the average is

∑ ρτ (X̂t ) ln ρτ (X̂t ) X̂ t

(12)

Loosely, the entropy rate is the mean uncertainty about the next coarse-grained state that will be visited along a sequence of states. The time dependence of the entropy rate is shown in Figure 5b. Beyond the largest peak, h 36 (t) decreases approximately exponentially. This decay coincides with the slowing down of the reaction when fewer paths are possible and there is less uncertainty about future states; inspecting the paths in the typical set gives further support. At late times, the only τlength path that is typical is a sequence of states with all inactive monomers, x = Mi. The entropy rate hτ(t) fixes the location of the neighborhood that contains typical paths. The typical sequences are those that lie within the ϵ-neighborhood around e−τhτ(t):51 Aϵt ≡ {xt̂ : e−τ(hτ (t )+ϵ) ≤ ρτ (xt̂ ) ≤ e−τ(hτ (t )−ϵ)}

⟨O(X̂τ (t ))⟩Aϵ* =

px (t ) =

(13)

∑ {xt̂ ∈ A ϵ36* | xt = x}

In the limit, lim t , τ →∞ Pr(xt̂ ∈ 1 in probability. For finite length symbol sequences, there is no guarantee of such strong convergence. In the cases we have considered, there is still a subset that holds the majority of the probability because of the dynamics. The question then becomes how to separate the typical from the atypical paths.58 The variational typical set59,60 characterizes macroscopic behavior through a collection of micro- or mesoscopic finitetime typical sequences of states. A natural variational principle is to maximize the probability of typical paths while also minimizing the number of paths at a given τ. We implement this principle by making the ϵ-neighborhood around hτ(t) timedependent and finding the ϵ*(t) that satisfies ϵ

ρτ (X̂τ (t ))O(X̂τ (t )) (15)

For the system considered here, the most probable path always corresponds to the sequence of states for a monomeric unit that remains inactive throughout the time interval. Thus, average observables that correspond to material properties cannot be calculated solely using the most probable path. Instead, the typical average is necessary to accurately calculate the concentration and fiber length/number. We can also contract the joint probability

p Aϵt ) →

ϵ*(t ) ≡ arg max[Pr(Aϵt ) − W ϵt ]

∑ X̂ τ (t ) ∈ A ϵt*

ρτ (xt xt − δt ... xt − 36) (16)

to obtain the marginal probability that each typical path ends in state x. Using only the typical paths here gives an approximation of the marginal distribution at time t. Converting this probability to mM units, Figure 5c shows that the typical paths closely approximate the concentration from the stochastic kinetics, as illustrated with τg = 70. The typical marginal distribution can also be used to compute averages of fiber observables. Figure 5d shows the average fiber length calculated exactly and the typical averages for τg = 35, 70, and 150. Despite only using a small percentage of all paths, we find good agreement between the actual and typical averages. The current can be decomposed to give deeper insight into the assembly process and used to assign dynamical processes to features of the total current. Figure 6 shows the current Jt decomposed into different (time-dependent) components. Figure 6b shows the irreversibility due to inactive/active

(14) 4797

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

monomers. The peaks in Jt are a combination of fiber birth, t J32 = J Ma ⊗ F1 in Figure 6c, and fiber growth and decay, J - in Figure 6d, - = {Ma ⊗ Fj , j = 2, 3, ..., N − 2}. When there is overall fiber growth, J - is positive. Likewise, a negative J corresponds to net fiber decay. The small peak in Jt just before 0.5 h is a combination of a shoulder in the average birth rate of fibers, Figure 6c, and a decrease in the rate of fiber collapse, Figure 6d. What is clear from decomposing the current is that local maxima correspond to fiber birth/growth and local minima correspond to fiber collapse. By distinguishing each contribution to Jt, we can link the information measures of the typical set h36 36 and | Aϵ* | to the chemistry and the fiber dynamics. Figure 7 36 t |. All of shows J as a function of time together with h36 and | Aϵ* these observables are normalized so that their maximum value is one. Each vertical line designates a peak in one measure, showing that peak locations of the information-theoretic measures coincide with those in the current. At later times, peaks in Jt can align to shoulders in information measures and vice versa. From the decomposition of the current, we can see that increases in the number and length of fibers correspond to an increase in the uncertainty about what state will be visited next (on average) along the path. Likewise, Jt < 0 tends to correspond to a dip in h36, as overall fiber collapse leads to a reduction in the uncertainty. The increase in uncertainty correlates with the size of the typical set, second row of Figure 7. Again, we see that an increase in the number of typical paths corresponds to an increase in fiber growth and number. The fact 36 | correspond to peaks in Jt that the maxima in both h36 and | Aϵ* means we can correlate these peaks with fiber birth and growth. So far, we have seen that gross behavior of fiber growth and decay can be linked to information measures. This behavior is captured by the fiber observables, such as the number and length of individual fibers. Since the typical set can be interpreted as a mechanistic representation of the assembly process, we would expect the information measures, | Aϵτ*| and hτ, to correlate with macroscopic features of the material. If we normalize the fiber measures, ⟨NtF⟩ and ⟨Ftm⟩, along with the information measures

Figure 6. Decomposition of the current into contributions from the chemical kinetics and fiber growth. (a) The current J t = ∑x < y Jxyt measures the degree to which all transitions in the Markov model break t t + J32 + J - . It has local components from detailed balance: J t = J21 the (b) activation−deactivation reaction Jt21, (c) the birth of fibers Jt32 = JMa⊗F1, and (d) the set of transitions for fiber growth and collapse J - , - = {M a ⊗ Fj , j = 2, 3, ..., N − 2}. The current associated with fiber birth decays exponentially as the finite amount of fuel is consumed.

transitions Jt21, which steadily decreases with the concomitant decline in DMS fuel that drives the activation/deactivation of

Figure 7. Total current Jt for τg = 35, 70, 92, and 150 normalized so that the maximum is one (color). Temporal location of peaks in J (solid vertical 36 lines) aligned with those of (a) h36 and (b) |Aϵ* | (dashed lines) also normalized to one. Maxima (minima) in the current correspond to high rates of growth (decay) in fiber number and length. The corresponding maxima (minima) in these information measures have the same physical origins. 4798

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B

Figure 8. Each column of rows a−d corresponds to a different value of τg = 35, 70, 92, and 150. Each row compares ⟨NF⟩ or ⟨Ftm⟩ (colors from left to right) to |Aϵτ*| or h36 (black), after normalizing such that the maximum of each function is one. Dashed lines mark the peak of each informationtheoretic quantity. Colored lines mark peaks in the fiber measures. For dominant peaks, there is strong correlation between these locations. Row d shows the times of peaks in fiber measures and information measures for all places where a dashed and solid line are near each other, i.e., neglecting shoulders. Peak locations are well-correlated with a slope of one. Note, while τg = 63 is not shown in parts a−d, it is included in the correlation plots (part e).

new fibers is enhanced and more fibers are produced, Figure 3b. The entropy rate, as a global measure of the dynamics, shifts, reflecting its dependence on both the number and length of fibers. The size of the typical set is equally sensitive to peaks in both the length and number of fibers, Figure 8c,d. The peak locations for both information measures correlate strongly with peaks in fiber observables. Figure 8e shows that the times of corresponding peaks are strongly correlated, with a slope of one. The correlation in peaks suggests that the information measures accurately identify the time of extremal values in the number and length of fibers.

such that the maximum of each is again one, we can directly compare how well the information measures capture the fiber properties. Rows a and b of Figure 8 compare the dynamical entropy rate to ⟨NtF⟩ and ⟨Ftm⟩ for τg = 35, 70, 92, and 150. Again, we find agreement between peaks in uncertainty, h36, and peaks in ⟨Ftm⟩ and ⟨NtF⟩. The agreement is especially good for ⟨NtF⟩ and h36 with τg = 35. The entropy rate is dependent on both the number and size of fibers, but for τg = 35, fibers tend to grow and then truncate at around the same length, meaning that there is less uncertainty in the length of fibers and h36 depends more strongly on the number of fibers. Early in the assembly process, fibers can be born and grow simultaneously. As τg grows, the birth rate of 4799

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

The Journal of Physical Chemistry B

jij1N × 1 jj jj N × 1 j0 7 = jjjj jj∂ jj jj j0N ×1 k

IV. CONCLUSIONS In the dissipative self-assembly of hydrogels, myriad kinetic pathways lead to a multitude of assembled states. The breadth of pathways makes this nonequilibrium process both a challenge to developing predictive theory and a promising avenue for functional materials. The massive space of possible pathways makes understanding and controlling those that are most likely an open challenge. Here, we have introduced an experimentally informed model that reproduces key features of recent molecular gelators, including the observation via confocal microscopy of fast fiber growth and stochastic fiber collapse. This model could be adapted to the design of other stochastic kinetic systems for chemically fueled dissipative assembly. To calculate average observables away from equilibrium, it may be tempting to use the most probable path as a representative of the overall system behavior. For the dissipative self-assembling hydrogel here, however, the most probable path does not describe the formation of the fibrous structure. Another approach to calculating observables is to use all paths with positive weight. Here, the exponential growth in the number of paths leads to considerable computational expense. The variational typical set of pathways in a nonequilibrium selfassembling system, however, balances accuracy and efficiency in computing time-dependent averages. We showed that the local maximum length and number of fibers is correlated with maximums in the entropy rate and size of the typical set. These results demonstrate the typical set of paths is a mechanistic description of the chemistry and structure formation from which we can calculate the average fiber length and concentrations. The set of typical paths are 99% of the path probability, making them a high-fidelity and theoretically efficient representation of the dissipative assembly mechanism.



ρηn− 1 = ρηn 7 ρηn + 1 = [|R1n − 1⟩⟨ρ1,nη − 1| , ..., |RNn− 1⟩⟨ρNn , η − 1|]



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jason R. Green: 0000-0003-2572-2838 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Acknowledgment is made to the donors of The American Chemical Society Petroleum Research Fund (ACS PRF No. 55195-DNI6) for support of this research. We also acknowledge support through the Joseph P. Healey Research Grant Program and the use of the supercomputing facilities managed by the Research Computing Group at the University of Massachusetts Boston as well as the University of Massachusetts Green HighPerformance Computing Cluster.



REFERENCES

(1) Whitesides, G. M.; Grzybowski, B. A. Self-assembly at all scales. Science 2002, 295, 2418−2421. (2) van Rossum, S. A. P.; Tena-Solsona, M.; van Esch, J. H.; Eelkema, R.; Boekhoven, J. Dissipative out-of-equilibrium assembly of man-made supramolecular materials. Chem. Soc. Rev. 2017, 46, 5519−5535. (3) Stupp, S. I.; LeBonheur, V.; Walker, K.; Li, L. S.; Huggins, K. E.; Keser, M.; Amstutz, A. Supramolecular Materials: Self-Organized Nanostructures. Science 1997, 276, 384−389. (4) Amabilino, D. B.; Smith, D. K.; Steed, J. W. Supramolecular materials. Chem. Soc. Rev. 2017, 46, 2404−2420. (5) Boekhoven, J.; Poolman, J. M.; Maity, C.; Li, F.; van der Mee, L.; Minkenberg, C. B.; Mendes, E.; van Esch, J. H.; Eelkema, R. Catalytic control over supramolecular gel formation. Nat. Chem. 2013, 5, 433− 437. (6) van Leeuwen, T.; Pijper, T. C.; Areephong, J.; Feringa, B. L.; Browne, W. R.; Katsonis, N. Reversible photochemical control of cholesteric liquid crystals with a diamine-based diarylethene chiroptical switch. J. Mater. Chem. 2011, 21, 3142. (7) Aggeli, A.; Bell, M.; Boden, N.; Keen, J. N.; Knowles, P. F.; McLeish, T. C. B.; Pitkeathly, M.; Radford, S. E. Responsive gels formed by the spontaneous self-assembly of peptides into polymeric β-sheet tapes. Nature 1997, 386, 259−262.

(17)

For n ≤ η, the joint is evolved as (18)

The joint probability at n over a window length η is ρnη. To evolve n from ρnη to ρn+1 η , we must project down from ρη(Xη, Xη−1, ..., X1) n to ρη−1(Xη, Xη−1, ..., X2) through the projection operator 7(N

(19)

Repeating these steps, the joint probability of block length η can be evolved up to any arbitrary n. To further increase computational efficiency, we only calculate the nonzero values of ρnη. In order to track the associated path, we also store the set of sequence labels for each joint probability. The set of joint probabilities and sequences {ρnη, X̂ η} with nonzero probabilities here is vastly smaller than an ordered list of all sequences and their corresponding probabilities.

In this section, we outline an efficient calculation of the joint probability distribution for time windows of length τ. For convenience, we will use time steps, n = t/δt, instead of time t. The window length is then η = τ/δt. For N states and paths of length n, the entire joint distribution ρ can be represented by a matrix of at most Nn values. Here, we will represent ρn at n as an N × Nn−1 matrix. The joint probability at n = 2, ρ2 = ρ(X2, X1) is written element by element as ρ(x 2 , x1) = R x12x1px (1). The outer product of the jth column of Rn−1 and the jth row of ρn−1 generates an N × Nn−2 block of entries in ρn. Concatenating each of the N blocks gives ρn. For example, ρ3 is

ρn = [|R1n − 1⟩⟨ρ1n − 1| , |R 2n − 1⟩⟨ρ2n − 1| , ..., |RNn− 1⟩⟨ρNn− 1|]

0 N × 1 μ 0 N × 1zyz zz zz zz 1N × 1 ∂ zz z ∂ ∏ 0 N × 1zzzz zz 0 N × 1 μ 1N × 1 z{

where 1N×1 is an N × 1 vector of ones and 0N×1 is an N × 1 vector of zeros. The general evolution from n to n + 1 over a block length η can then be seen as a computation over two steps

APPENDIX

ρ3 = [|R12⟩⟨ρ12 | , |R 2 2⟩⟨ρ2 2 | , ..., |RN 2⟩⟨ρN 2 |]

Article

n−1

) × (N n − 2)

4800

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B (8) Tagliazucchi, M.; Weiss, E. A.; Szleifer, I. Dissipative self-assembly of particles interacting through time-oscillatory potentials. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 9751−9756. (9) Boekhoven, J.; Hendriksen, W. E.; Koper, G. J. M.; Eelkema, R.; van Esch, J. H. Transient assembly of active materials fueled by a chemical reaction. Science 2015, 349, 1075−1079. (10) Longo, G. S.; Pérez-Chávez, N. A.; Szleifer, I. How protonation modulates the interaction between proteins and pH-responsive hydrogel films. Curr. Opin. Colloid Interface Sci. 2019, 41, 27−39. (11) Rieß, B.; Boekhoven, J. Applications of Dissipative Supramolecular Materials with a Tunable Lifetime. ChemNanoMat 2018, 4, 710−719. (12) Grzelczak, M.; Vermant, J.; Furst, E. M.; Liz-Marzán, L. M. Directed Self-Assembly of Nanoparticles. ACS Nano 2010, 4, 3591− 3605. (13) Grant, J.; Jack, R. L.; Whitelam, S. Analyzing mechanisms and microscopic reversibility of self-assembly. J. Chem. Phys. 2011, 135, 214505. (14) Whitelam, S.; Schulman, R.; Hedges, L. Self-Assembly of Multicomponent Structures In and Out of Equilibrium. Phys. Rev. Lett. 2012, 109, 265506. (15) Green, J. R.; Costa, A. B.; Grzybowski, B. A.; Szleifer, I. Relationship between dynamical entropy and energy dissipation far from thermodynamic equilibrium. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 16339−16343. (16) Whitelam, S.; Hedges, L. O.; Schmit, J. D. Self-Assembly at a Nonequilibrium Critical Point. Phys. Rev. Lett. 2014, 112, 155504. (17) Jankowski, E.; Glotzer, S. C. Screening and designing patchy particles for optimized self-assembly propensity through assembly pathway engineering. Soft Matter 2012, 8, 2852−2859. (18) Haxton, T. K.; Whitelam, S. Do hierarchical structures assemble best via hierarchical pathways? Soft Matter 2013, 9, 6851−6861. (19) Zenk, J.; Schulman, R. An Assembly Funnel Makes Biomolecular Complex Assembly Efficient. PLoS One 2014, 9, e111233. (20) Whitelam, S.; Jack, R. L. The Statistical Mechanics of Dynamic Pathways to Self-Assembly. Annu. Rev. Phys. Chem. 2015, 66, 143−163. (21) Jacobs, W. M.; Reinhardt, A.; Frenkel, D. Rational design of selfassembly pathways for complex multicomponent structures. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 6313−6318. (22) Tena-Solsona, M.; Wanzke, C.; Riess, B.; Bausch, A. R.; Boekhoven, J. Self-selection of dissipative assemblies driven by primitive chemical reaction networks. Nat. Commun. 2018, 9, 2044. (23) Klajn, R.; Wesson, P. J.; Bishop, K. J. M.; Grzybowski, B. A. Writing self-erasing images using metastable nanoparticle “inks. Angew. Chem., Int. Ed. 2009, 48, 7035−7039. (24) Colomer, I.; Morrow, S. M.; Fletcher, S. P. A transient selfassembling self-replicator. Nat. Commun. 2018, 9, 2239. (25) Amabilino, D. B.; Smith, D. K.; Steed, J. W. Supramolecular materials. Chem. Soc. Rev. 2017, 46, 2404−2420. (26) Samanta, D.; Klajn, R. Aqueous light-controlled self-assembly of nanoparticles. Adv. Opt. Mater. 2016, 4, 1373−1377. (27) Kim, Y.; Tamaoki, N. A photoresponsive planar chiral azobenzene dopant with high helical twisting power. J. Mater. Chem. C 2014, 2, 9258. (28) Kundu, P. K.; Samanta, D.; Leizrowice, R.; Margulis, B.; Zhao, H.; Börner, M.; Udayabhaskararao, T.; Manna, D.; Klajn, R. Lightcontrolled self-assembly of non-photoresponsive nanoparticles. Nat. Chem. 2015, 7, 646−652. (29) Cademartiri, L.; Bishop, K. J. M. Programmable Self-Assembly. Nat. Mater. 2015, 14, 2−9. (30) Dambenieks, A. K.; Vu, P. H. Q.; Fyles, T. M. Dissipative assembly of a membrane transport system. Chem. Sci. 2014, 5, 3396. (31) Korevaar, P. A.; George, S. J.; Markvoort, A. J.; Smulders, M. M. J.; Hilbers, P. A. J.; Schenning, A. P. H. J.; Greef, T. F. A. D.; Meijer, E. W. Pathways complexity in supramolecular polymerization. Nature 2012, 481, 492−496. (32) Korevaar, P. A.; Newcomb, C. J.; Meijer, E. W.; Stupp, S. I. Pathway selection in Peptide Amphiphile Assembly. J. Am. Chem. Soc. 2014, 136, 8540−8543.

(33) Gröschel, A. H.; Schacher, F. H.; Schmalz, H.; Borisov, O. V.; Zhulina, E. B.; Walther, A.; Müller, A. H. E. Precise hierarchical selfassembly of multicompartment micelles. Nat. Commun. 2012, 3, 710. (34) Hirst, A. R.; Roy, S.; Arora, M.; Das, A. K.; Hodson, N.; Murray, P.; Marshall, S.; Javid, N.; Sefcik, J.; Boekhoven, J.; van Esch, J. H.; Santabarbara, S.; Hunt, N. T.; Ulijn, R. V. Biocatalytic induction of supramolecular order. Nat. Chem. 2010, 2, 1089−1094. (35) Korevaar, P. A.; George, S. J.; Markvoort, A. J.; Smulders, M. M. J.; Hilbers, P. A. J.; Schenning, A. P. H. J.; De Greef, T. F. A.; Meijer, E. W. Pathway complexity in supramolecular polymerization. Nature 2012, 481, 492−496. (36) Boekhoven, J.; Poolman, J. M.; Maity, C.; Li, F.; van der Mee, L.; Minkenberg, C. B.; Mendes, E.; van Esch, J. H.; Eelkema, R. Catalytic control over supramolecular gel formation. Nat. Chem. 2013, 5, 433− 437. (37) Bishop, K. J. M.; Wilmer, C. E.; Soh, S.; Grzybowski, B. A. Nanoscale Forces and Their Uses in Self-Assembly. Small 2009, 5, 1600−1630. (38) Dou, Y.; Dhatt-Gauthier, K.; Bishop, K. J. M. Thermodynamic costs of dynamic function in active soft matter. Curr. Opin. Solid State Mater. Sci. 2018, DOI: 10.1016/j.cossms.2018.11.002. (39) Hormoz, S.; Brenner, M. P. Design principles for self-assembly with short-range interactions. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 5193−5198. (40) Nguyen, M.; Vaikuntanathan, S. Design principles for nonequilibrium self-assembly. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 14231−14236. (41) van Esch, J. H.; Klajn, R.; Otto, S. Chemical systems out of equilibrium. Chem. Soc. Rev. 2017, 46, 5474−5475. (42) Gillespie, D. T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 1976, 22, 403−434. (43) Gillespie, D. T. A rigorous derivation of the chemical master equation. Phys. A 1992, 188, 404−425. (44) Husic, B. E.; Pande, V. S. Markov State Models: From an Art to a Science. J. Am. Chem. Soc. 2018, 140, 2386−2396. (45) Sorrenti, A.; Leira-Iglesias, J.; Markvoort, A. J.; de Greef, T. F.; Hermans, T. M. Non-equilibrium supramolecular polymerization. Chem. Soc. Rev. 2017, 46, 5476−5490. (46) Van den Broeck, C.; Esposito, M. Ensemble and trajectory thermodynamics: A brief introduction. Phys. A 2015, 418, 6−16. (47) Schnakenberg, J. Network theory of microscopic and macroscopic behavior of master equation systems. Rev. Mod. Phys. 1976, 48, 571. (48) Nicholson, S. B.; del Campo, A.; Green, J. R. Nonequilibrium uncertainty principle from information geometry. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2018, 98, 032106. (49) Jackson, E. A. Equilibrium statistical mechanics; Courier Corporation: 2000. (50) Algoet, P. H.; Cover, T. M. A Sandwich Proof of the ShannonMcMillan-Breiman Theorem. Ann. Probab. 1988, 16, 899−909. (51) Cover, T. M.; Thomas, J. A. Elements of Information theory; Wiley: 2006; Vol. 2. (52) Shannon, C. E. A Mathematical Theory of Communication. Bell Syst. Tech. J. 1948, 27, 623−656. (53) McMillan, B. The Basic Theorems Of Information Theory. Ann. Math. Stat. 1953, 24, 196−219. (54) Chung, K. L. A. Note on the Ergodic Theorem of Information Theory. Ann. Math. Stat. 1961, 32, 612−614. (55) Gaspard, P. Time-Reversed Dynamical Entropy and Irreversibility in Markovian Random Processes. J. Stat. Phys. 2004, 117, 599. (56) Zegers, P.; Fuentes, A.; Alarcón, C. Relative Entropy Derivative Bounds. Entropy 2013, 15, 2861−2873. (57) Schürmann, T.; Grassberger, P. Entropy estimation of symbol sequences. Chaos 1996, 6, 414−427. (58) Nicholson, S. B.; Greenberg, J. S.; Green, J. R. Entrance and escape dynamics for the typical set. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2018, 97, 012146. 4801

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802

Article

The Journal of Physical Chemistry B (59) Nicholson, S. B.; Alaghemandi, M.; Green, J. R. Learning the Mechanisms of Chemical Disequilibria. J. Chem. Phys. 2016, 145, 084112. (60) Nicholson, S. B.; Alaghemandi, M.; Green, J. R. Effects of temperature and mass conservation on the typical chemical sequences of hydrogen oxidation. J. Chem. Phys. 2018, 148, 044102.

4802

DOI: 10.1021/acs.jpcb.9b02811 J. Phys. Chem. B 2019, 123, 4792−4802