2.1 Inference for diffusion processes

Let us consider the stochastic differential equation (SDE)

with unknown drift A(x), noise intensity B(x) with Klimontovich interpretation [cf. the generator (2)] [23], and x∈Rd. The estimation of A(x) and B(x), or the corresponding diffusion matrix D(x):=B(x)B(x)T, is the general problem of inference for diffusion processes, which in the statistics literature is studied for instance in [24], [25].

For our illustrative purposes, working in terms of the infinitesimal generator of the process (1),

is particularly transparent. Indeed, applying the generator to the test functions fi(x)=xi and fij(x)=xixj, we see that the diffusion matrix may be computed by the formula

or, in matrix notation,

As a consequence, the drift is computed as^[1]

or

All formulas inspired by eq. (4) are referred to as Green–Kubo formulas and assume different forms depending on the application [5, Section 8.4]. We consider them as a way to infer the diffusion matrix from a sample of the time correlations of the second moments (cf. [24, Section 4.2.3]). Although we have considered only processes in Rd, the same conclusions can be drawn, at least formally, for infinite-dimensional systems; see [26] for a theory of measure-valued diffusion processes, and [27] for an application to a field theory.

2.2 Coarse-graining via the Green–Kubo formula

Let a physical system be described by the variables y=y1,y2,…,yn, with a large number n of degrees of freedom, and suppose that its dynamics is characterized by a “slow” and a “fast” time scale. The first goal of dynamic coarse-graining, in the sense used in this paper, is to identify a set of more macroscopic state variables X(y)∈Rd that resolve the slow dynamics, meaning that, in the limit of many degrees of freedom, they evolve only “slowly”. In addition, we suppose that the dynamics of the macroscopic variables, for large but finite n, can be well approximated by a stochastic process, which we also call X for simplicity: the most probable realizations of the stochastic process are regarded as paths of the slow dynamics, and the noise represents the idealized effect, at the macroscopic level, of the more microscopic, fast degrees of freedom. The identification of the coarse-grained variables (also called “collective variables”) is often a very challenging step [9].

In many cases the macroscopic variables X can be assumed to be the sum of short-time correlated interactions of many microscopic particles [14]. The short-time correlation implies that X can be assumed to be Markovian, in the presence of an infinitely large separation of time scales; in practice one has to deal with finite separations, which lead to memory effects. When the interactions are frequent and small, such that the process X can be considered to have continuous sample paths, then the process X is a diffusion and can be described by an SDE [28, Sec. 2.5]. In this case the noise is Gaussian.

We saw in (I) how the classical FDT suggests that the diffusion process Xt solves the SDE

with 2kBM(x)=D(x); the operator M(x) is called a friction matrix [29], [30], and the most probable path solves the gradient-flow equation

In contrast to (I), here we are assuming only dissipative components in the slow dynamics. This assumption is related to a condition of detailed balance [31, Section II.5] of the stochastic process Xt with respect to the distribution eS(x)/kB; see [32, Section 6.3.5] and [7].

The second goal of dynamic coarse-graining, for such systems, is to infer the dissipative structure, encoded by S(x) and M(x), by simulations. In this paper we always assume that the static simulations for the estimation of S(x) have already been performed, and the dynamic simulations for evaluating the friction matrix M(x) are based on the Green–Kubo formula (4).

Two remarks are in order. On the one hand, in the Green–Kubo formula, the limit τ→0 presupposes a continuous-time sample, while in concrete applications we always consider discrete-time observations. On the other hand, the limit τ→0 is not even desirable, since the diffusion process Xt is only an approximation of the true (in general non-Markovian) macroscopic dynamics, and this approximation is usually valid in a range of time scales that are large with respect to the microscopic, fast time scale t2. The limit τ→0 is thus replaced by τ being in the range

where t1 is the macroscopic, slow time scale. In other words, τ has to be macroscopically small, but microscopically large. Hence, the friction matrix may be estimated as

By combining coarse-graining with the fluctuation-dissipation theorem, as described above, many systems have successfully been studied. On the side of theory, we refer to [1, Section 4.2], [33, Section 2.7], [5, Section 8.4], [6], [10, Chapter 4]; for numerical results, see, e. g., [34], [35], [27]. The power of the method stems from two important features: (i) short simulations, with respect to the slow time scale t1, are sufficient to sample the matrix M(x); (ii) a single numerical experiment provides us with both the evolution equations and their dissipative structure, without the need of resorting to a different experiment for each irreversible process.

However, as we have remarked in (I), some systems have fluctuations that are distinctly non-Gaussian in nature, suggesting that the method just described may not give correct results. Our main example in Section 4 below is of this type. In that section we apply and implement the generalized FDT of this paper to this example, but we also investigate to which extent the application of Gaussian-based methods would give incorrect results (see Section 4.2).

In Section 4.2 we study three Gaussian-based methods: “Green–Kubo,” based on the scheme of this section, the “chemical Langevin equation,” and the “log-mean equation.” Each of them is defined by a different choice of the pair S(x),M(x). It is possible to sample this pair correctly for the Green–Kubo method, but with the result of a wrong macroscopic phenomenological equation; for the other two methods the macroscopic evolution equation is correct, but a simulation would fail in sampling either S(x) or M(x).

Apparently, applying Gaussian-based methods to systems with non-Gaussian fluctuations is like forcing a square peg into a round hole: the results will be suboptimal. In order to deduce the dissipative structure of such non-Gaussian systems, in this series of papers, we generalize the FDT and propose an extended theory of coarse-graining.

3.1 Inference for Markov jump processes

Markov jump processes are completely defined by the transition rates between all of their states, and this information is encoded in the transition rate matrix, fully equivalent to the infinitesimal generator. Very often, a Markov jump process is observed during an experiment or a numerical simulation, and one wishes to determine the transition rates from the observations. This is the task of the statistical inference for Markov jump processes, which was studied for instance in [36] and reviewed in [37].

In the case of continuous-time observations, it is easy to check [36] that the best estimator for an element of the transition rate matrix is

where Nij is the number of transitions from the state i to j and Ri is the time spent in i.

The observations, however, are never continuous, but are made at discrete times. The main issue of this situation is the famous embedding problem [38], [36]: to the same discrete-time Markov chain, the transition rates of which are estimated in practice, there may correspond zero, one, or many continuous-time Markov processes that have the same finite-time transition rates. We will not focus on this issue in this paper, although we believe it to be important for future developments.

The typical physical setup where jump processes arise naturally is the dynamics of many particles in an energy landscape characterized by metastable states. By a metastable state one indicates the region of attraction of a local energy minimum such that the time scale for the system to equilibrate is much shorter than the time scale to escape from it [39], [40]. Since the system spends most of its time in these special states, the escape events are “rare,” and the study of these phenomena has given rise to the field of rare-event estimation and simulation [18], [19], [20], [21], [22].

The most famous, preliminary toy model in this class of systems was studied by Eyring [16] and Kramers [17], who gave the statistical-mechanical derivation of an explicit formula for the transition rates between the two minima of a double-well potential. We will study this system again from the standpoint of our generalized FDT in Section 4. In particular, our goal is not to compute the transition rates between metastable states, but to resolve the dissipative structure of the macroscopic phenomenological equations, as we explore in the next subsection.

3.2 Coarse-graining via the generalized FDT

By analogy with the Gaussian picture of Section 2.2, let us consider again a setup described by microscopic variables y, and with a separation of time scales. A set of more macroscopic variables X(y) is introduced to separate the time scales effectively. We now assume that X is not necessarily a diffusion process, but may be represented as a general Markov process, and the most probable paths are the solutions of the deterministic equation

We then aim to estimate the structure of this phenomenological equation by analyzing the noise that results, at the macroscopic level, from the neglected degrees of freedom. Among these two steps, the identification of the macroscopic variables and the estimation of the force-flux constitutive law, which together define what we call a coarse-graining procedure, this paper focuses entirely on the second one.

In contrast to paper (I), here we restrict the form of the phenomenological equations to the generalized gradient flow

which corresponds to a purely dissipative dynamics governed by an entropy function S and a dissipation potential Ψ∗.

As elaborated in (I), following [7], the dissipation potential Ψ∗ may be found by studying the stochastic process Xt. In particular, we need to compute the following cumulant generating function (cf. eq. 28 in (I) and [41, Chapter 1]):

In this expression, the expectation is taken over all possible realizations of the stochastic process Xt starting from x and with time duration τ, and τ is far from the fast time scale t2 and the slow time scale t1:

The generalized FDT implies that

This correspondence between the left-hand side (a property of the stochastic process) and the right-hand side (the structure of the most probable evolution) has its origin in the connection between large deviations for Markov processes and the generalized gradient flows that was proven in [7], where it is explained how the connection, in the purely dissipative case, is based on a detailed-balance property of the stochastic process with respect to the distribution eS(x)/kB.

In this paper we suppose that the distribution eS(x)/kB has already been sampled to find the function S(x), so that the static properties of the system are fully known. The main focus, instead, rests upon the computation of the dissipation potential Ψ∗ by the formula (14), which completely characterizes the dynamics.

From the practical standpoint, it is convenient to evaluate the expression (14), for fixed x, at the values

Indeed, for ξ=0, we get the second term on the left-hand side of eq. (14) because Ψ∗(x,0)=0; with the other values of ξ we explore the dissipation potential in ξ-space.

It is not the purpose of this paper to construct efficient simulations, nor to pursue any statistical rigor, which we reserve for future work. For instance, in the same spirit of the inference of the infinitesimal generator of a continuous-time Markov process, a theory of inference for the nonlinear generator, the limit of the cumulant generating function (12) as τ→0, should be developed.

We have seen that the computation of the cumulant generating function for the stochastic process Xt, the left-hand side of eq. (14), together with the information on the static distribution, provides us with the dissipative structure, expressed in terms of an entropy function S and a dissipation potential Ψ∗, for the macroscopic phenomenological equation. Note, in particular, that we do not need to assume the nature of the process Xt, except that it is Markovian and it satisfies detailed balance: the Gaussian case gives rise to a quadratic dissipation potential, and thus to a friction matrix. The procedure based on the generalized FDT, thus, gives a unified way of dealing with both Gaussian and non-Gaussian Markovian fluctuations, and with both the Green–Kubo and the Kramers picture.

4.1 A multiscale view: The levels of description

4.1.1 The microscopic level: Diffusion in a double-well potential

Figure 1

The double-well potential V that we use in the numerical experiment. The regions A=(−∞,0) and B=(0,∞) are representative of the two chemical states, and the height of the barrier is ΔV:=V(0)−V(L).

We consider a large number n of reactive particles in a mixture. The constituents of the mixture do not directly enter our description of the system, which is based only on the state of the reactive particles. We assume that the reactive particles are independent, and the state of each particle is described by its position on the real line, which can be interpreted as a reaction coordinate [42, p. 1158] or collective variable or coordinate [46], [47]. We gather all positions in the array y1,y2,…,yn∈Rn.

Each particle follows an overdamped Langevin dynamics in the energy landscape of Figure 1; in this dynamics, the noise represents the effective interaction of each particle with all constituents in the mixture, which is physically described as a heat bath at temperature T. The two wells of the energy landscape correspond to the two chemical states A and B, and the motion is described by the SDEs

(16)dYti=−1γV′(Yti)dt+kBTγdWti(i=1,…,n),

where γ is a friction coefficient with dimensions of [mass]/[time], V′ is the derivative of the energy, and the Wi are independent Wiener processes.

Although eq. (1) and eq. (16) are formally equivalent, their physical roles should be clearly distinguished: eq. (1) represents the evolution equation of a macroscopic state variable with noise enhancement in the statistical-mechanical setting of this paper; eq. (16), instead, is an effective microscopic dynamics.

4.2 Green–Kubo formula and the shortcomings of diffusion approximations

The generalized FDT provides us with the dissipation potential and the entropy function of the RRE from the properties of the associated stochastic process. Specifically, it gives us (i) an entropy function that can be sampled from the static distribution and (ii) a dissipation potential that can be sampled from the cumulant generating function; moreover, (iii) the generalized gradient flow constructed from the entropy and the dissipation potential produces the correct form of the RRE. Can these three features be reproduced by assuming the stochastic process Xt to be a diffusion process? In this section we find an answer.

Note that the assumption of a diffusion process is not just an academic exercise since, in practical situations, we do not know the nature of the macroscopic process Xt, but only have access to the dynamics of a more microscopic model. Assuming a diffusion process means restricting ourselves to the classical scheme of Green–Kubo relations.

Let us assume that the “Green–Kubo” diffusion process

is an accurate dynamics for the macroscopic variables. Since it is a diffusion, the corresponding dissipation potential is quadratic, and the friction matrix may be computed by the Green–Kubo formula

Note that, in practice, the actual sampling is made with simulations at a more microscopic level and, since in this case we know that the CME is the correct dynamics for the macroscopic variables, from eq. (3) we expect the result

where f(x)=x2, g(x)=x, and Q is the generator (20) of the CME. We remark that MGK is, in general, a function of x, and it would be so if there were two distinct rate constants k−≠k+ for the backward and forward reactions. With the entropy (19), we get

which has (i) the correct stationary distribution, that is, the correct entropy (by construction), (ii) a Green–Kubo expression for the friction matrix that can be computed by simulations, but (iii) the wrong drift and, therefore, the wrong most probable evolution (cf. Figure 2). (Note that, while “spurious” drift terms sometimes appear as the result of the different possible interpretations of the noise in an SDE, in this case the noise is additive, and no such terms arise.)

Figure 2

Comparison between the correct most probable evolution, given by the RRE, and the most probable evolution of the “Green–Kubo” diffusion process. (a) Simulation of the true Kramers problem and average over 200 realizations. (b) Analytical solution of the RRE. (c) Simulation of the “Green–Kubo” diffusion process and average over 200 realizations (10 realizations have fallen out of the domain [0,1]). (d) Numerical solution of the most probable evolution corresponding to the “Green–Kubo” diffusion process. The values of the parameters are: ΔV/(kBT)=4, Δt=10−2, n=10, and k=k¯ according to the Kramers formula (34).

Other possibilities of constructing a diffusion process for chemical reactions have been proposed: their goal is to approximate the CME for a large number of particles. One of them is called the chemical Langevin equation [56], which can also be interpreted as the diffusion approximation^[3] of the CME, and reads

with the entropy

and the friction matrix (30). This process has (iii) the correct most probable evolution, (ii) a Green–Kubo expression for the friction matrix, but (i) the wrong entropy, viz., the wrong stationary distribution (cf. Figure 3).

A third choice is the log-mean equation [60]

with the entropy (19) and the friction matrix (27). It has (iii) the correct most probable evolution, (i) the correct stationary distribution, but (ii) the friction matrix cannot be sampled by the Green–Kubo formula (cf. later Figure 5).

Figure 3

Comparison between the entropy functions associated to the stationary distribution of the CME and the CLE. The stationary distribution eSCLE/kB implies a nonzero probability for states outside the domain [0,1].

A common feature of all three diffusion processes is that the macroscopic variable X has a positive probability of exiting the domain [0,1], and this phenomenon may be seen both in the stationary distribution (cf. Figure 3) and in the dynamics (in the simulation behind Figure 2, 10 out of 200 realizations of the “Green–Kubo” diffusion process have exited the domain [0,1]).

In conclusion, there is no way for a diffusion process to satisfy requirements (i)–(iii) simultaneously (cf. Table 1). In order to satisfy all three requirements, one should move to more general Markov processes and to nonquadratic dissipation potentials.

4.3 Numerical experiments

The goal of the coarse-graining scheme proposed here is to infer the structure of a more macroscopic level of description from a more microscopic one, which in this case is represented by the overdamped Langevin dynamics of many independent particles in a double-well potential. In such a simple situation, the macroscopic variable is easily chosen as the concentration of particles in one well. The dynamics of the concentration, when the thermal energy kBT is small with respect to the energy barrier ΔV, is well approximated by a Markov process. However, we stress the fact that, in general, the system over which we have full control is the microscopic one; even if, in this simple case, we know the form of the macroscopic Markov process and the dissipation potential (up to the parameter k that completely characterizes them), in the numerical calculations we pretend to have no macroscopic information.

Before describing our numerical method, and to highlight the free parameters of the system, let us make eq. (16) dimensionless through the friction coefficient γ, the energy barrier ΔV, and L, the distance between the maximum and the minimum points in the energy landscape. Defining by a˜ the dimensionless quantity relative to a, we find

From the numerical calculations and the values of the dimensional parameters γ, ΔV, and L, one may find the physically meaningful results. To simplify the notation, from here on we drop the tildes.

The method advanced in this paper aims at computing the dissipation potential numerically. When, like in this case, we know its functional form, we could also determine the value of its unknown parameter k and compare it with the reference one given by Kramers’ formula [17], [39]

In accordance with the approach of this paper, we suppose that the entropy function S, or its derivative ∂S/∂x [12], has already been found, and we concentrate on the dynamics by computing the cumulant generating function (12). The evaluation of the cumulant generating function depends on the energy landscape, which we choose as

and on five parameters:

1. The ratio r:=ΔV/(kBT), which controls the separation of time scales. Kramers mentions in his seminal paper [17] that r=2.5 is sufficient: the process becomes approximately Markov. We actually work with r=4, which implies k¯≈4·10−5. In real applications, however, this is not a parameter, but a model datum.

2. The time-step size Δt→0 of the numerical scheme used to simulate the SDE (16), which should resolve the microscopic dynamics, guarantee stability of the scheme, and be smaller than the local equilibration time t2.^[4] We consider the numerical value Δt=0.01.

3. The time interval τ. This time constant should be “macroscopically small,” i. e., much smaller than the typical jump time t1=k−1, but also larger than the equilibration time t2, in such a way that we retain only the macroscopic features of the process Xt. We take τ=k¯−1/50≈5·102, where k¯ is given by formula (34). In more general contexts, we may not know the values of the characteristic times in advance and should perform an appropriate estimation of them.

4. The number of particles n→∞. Since the particles are independent, the cumulant generating function factorizes into cumulant generating functions for the single particles. Its computation would then require the simulation of just one particle. However, since n controls the discretization of the space Xn, and to keep the generality of our framework, which should work for interacting particle systems as well, we choose n=10.

5. The sample size N→∞ over which we calculate the empirical expectation. To obtain good statistics, sufficiently many jumps should occur in the total observation time. Since the average jump time is t1, we need Nτ≫t1. With the choice N=105, the total average number of jumps is Nτ/t1≈2·103.

The numerical setup is the following.

1. For every x∈Xn, we select a series of values αj=(ξj−∂S/∂x)/(2kB), with the ξj logarithmically spaced^[5] in the interval [−ξmax,ξmax] with ξmax=8. The logarithmic spacing has the aim of resolving the region around ξ=0 in a sufficiently accurate way.

2. For every x, we run N simulations, of length τ, of the n independent SDEs (16). We use an Euler–Mayurama scheme with step size Δt. The starting point for nx particles is 1 and that for the others is −1, that is, the initial microscopic state is

   (1,1,…,1︸nx,−1,−1,…,−1︸n−nx).

   We need not care about equilibration in the wells because τ≫t2.

3. We index the simulations by ℓ and say that the random variable Xτ has xτℓ as its realization. After the ℓ-th simulation started from x, we compute the quantities

   (35)hjℓ(x):=eαjxτℓ−x.

   There is one such quantity for each x, j, and ℓ.

4. To estimate the expectation in eq. (12), we calculate

   (36)2kBτln(1N∑ℓ=1Nhjℓ(x)).

Figure 4

Numerical estimate of the dissipation potential Ψˆ∗/(nkB) compared with the reference result. The latter is represented by the smooth surface, and the former by the red dots. The values of the parameters are: r=4, Δt=10−2, τ=k¯−1/200, n=10, N=105, ξmax=8.

Figure 5

(a) The data points Ψˆ∗(x,ξ)/(nkB) (x=0.5,0.1) are fitted with (b) the function Acosh(ξ/2)−1 [cf. eq. (22)] and compared with (c) the reference dissipation potential with k=k¯. The confidence intervals are at level 99%. In addition, we compare (d) the dissipation potential (38) derived from the friction matrix (30), which can be computed by the Green–Kubo formula, with (e) the dissipation potential (37) derived from the friction matrix (27) with k=k¯. The two agree only near the equilibrium point x=0.5.

The results of the algorithm are displayed in Figure 4, where the solid surface is the reference dissipation potential (22) with k=k¯, and the red dots are its estimated values, which show good agreement. If we assume the functional form of the dissipation potential given by eq. (22), we may determine the reaction constant by fitting the simulation points with the function (22). We do so with the values of Ψˆ∗(x,ξ) at the points x=0.5,0.1. In Figure 5, we compare this fitting procedure to the reference dissipation potential with k=k¯. We observe a small discrepancy, which we do not fully understand, but we do not know a simple alternative to compute the dissipation potential numerically. We only mention that the Kramers formula overestimates the reaction constant k for any finite r; see [40] for a more sophisticated method to compute k.

In addition, we show how the diffusion setting does not work in this example. In the previous subsection we have seen how two diffusion approaches, “Green–Kubo” and the “chemical Langevin equation,” may be constructed from a friction matrix that can be computed by the Green–Kubo formula, but do not reproduce either the correct most probable evolution or the static distribution of the process Xt. We now demonstrate numerically how also the third diffusion process, the “log-mean equation,” is inconsistent, namely, the corresponding friction matrix (27) cannot be sampled by the Green–Kubo formula (30). The friction matrix (27) would give rise to the dissipation potential

whereas the Green–Kubo formula would suggest

The two agree in the region around the equilibrium point x=0.5, but differ significantly for more extreme values of x, as one observes in the second frame of Figure 5. From expressions (22), (37), and (38), we also note that, in the same region around x=0.5, all dissipation potentials agree for values of ξ close to 0: in particular, at x=0.5 and close to ξ=0, the dissipation potentials (37) and (38) are the second-order Taylor approximations of the dissipation potential (22).

We emphasize that our aim has been to directly resolve the dissipative structure of the macroscopic equations: this is not the aim of standard approaches [19], [21], [20]. We also expect that our viewpoint will constitute an advantageous tool for less simple systems, such as problems of homogenization [63] or systems where the structure of the macroscopic dynamics is not known in advance [64].

If we consider more general reaction networks, characteristic problems of standard approaches will remain: for instance, high local minima of the potential landscape are rarely explored and boundaries between macroscopic states are not easy to set [21]. The method proposed here presents the additional complication of being “stiff” because of the strong nonlinearity in eq. (35).

To tackle the issues of numerical efficiency, importance sampling has been established as the basis of many known numerical methods in statistical mechanics [65]: the probability distribution of a random variable is changed, often by “exponential tilting,” in such a way that rare events become less rare and can more easily be observed. Physically, this corresponds to biasing the system by an external force. Following the ideas in [20], we would like to improve our numerical algorithm by importance-sampling techniques.