Multidimensional characterization of stochastic dynamical systems based on multiple perturbations and measurements

Generalized nonlinear response theory is presented for stochastic dynamical systems. Experiments in which multiple measurements of dynamical quantities are used along with multiple perturbations of parameters of dynamical systems are described by generalized response functions (GRFs). These constitute a new type of multidimensional measures of stochastic dynamics either in the time or the frequency domains. Closed expressions for GRFs in stochastic dynamical systems are derived and compared with numerical non-equilibrium simulations. Several types of perturbations are considered: impulsive and periodic perturbations of temperature and impulsive perturbations of coordinates. The present approach can be used to study various types of stochastic processes ranging from single-molecule conformational dynamics to chemical kinetics of finite-size reactors such as biocells. C 2015 AIP Publishing LLC. [http: dx.doi.org 10.1063 1.4917527]


I. INTRODUCTION
Multidimensional spectroscopic techniques span a broad range of frequencies from NMR to the UV and use numerous types of pulse sequences and detection modes. 1 They all share one common thread: the molecule is excited by several pulses followed by a single detection event. nD techniques involve n perturbations followed by a single detection event and are represented by causal response functions. Thinking more broadly, one can envision other types of signals that involve n perturbations and m measurements. These measure both response and spontaneous fluctuations and are not causal. In the linear regime, there are two possible signals (excitation followed by measurement or two measurements) but they are connected by the fluctuation dissipation (FD) [2][3][4] relation and carry the same information. 5 This is not the case for the nonlinear response since there is no universal FD theorem. The basic ideas are not limited to optical pulses and hold for any type of perturbation (electrical, mechanical, chemical concentrations, etc.) In this paper, we discuss several examples of such new types of multidimensional signals.
The idea of using multiple perturbations and measurements has been discussed previously. Reference 6 considered generalized response functions (GRFs) corresponding to a series of perturbations of open quantum system followed by a number of measurements. That formalism allowed the calculation of van der Waals forces between molecules subject to time-dependent coupling. GRFs corresponding to impulsive perturbations of spin systems were considered in Ref. 7 and used the derived generalized multi-point correlation functions to study stochastic dynamics of disordered spin models. The idea of mixing multiple perturbations with measurements was extended in Ref. 8, to arbitrary type of time-dependent perturbation, and applied to examine classical dynamics of nonlinear Hamiltonian systems. That formalism a) Electronic address: mkryvohu@uci.edu was further extended to non-Hamiltonian systems, such as chemical reactions, and was used to study nonlinear kinetics and connectivity of species in chemical systems. 9 In the present paper, we introduce a generalized response theory of stochastic dynamical systems. The formalism is based on a perturbative expansion of the Fokker-Planck equation around steady state. The Fokker-Planck approach can be used for equilibrium 10 as well as non-equilibrium 11 systems; thus, the developed formalism is applicable to systems in thermal equilibrium as well as in non-equilibrium steady-state.
The paper is organized as follows. Expressions for nonlinear response functions of a general stochastic process are derived in Secs. II and III. Section IV considers the special case of impulsive perturbations. The concept of GRF is introduced in Sec. V. Section VI introduces generalized susceptibilities for a harmonic perturbation. GRFs of exactly solvable Ornstein-Uhlenbeck process (OUP) are calculated in Sec. VII and compared with analytical expressions. 2D GRFs corresponding to perturbations of temperature in bistable and nonlinear stochastic systems at thermal equilibrium are calculated in Sec. VIII, and the quantitative information contained in GRFs is analyzed. GRFs of a stochastic chemical network are calculated in Sec. IX and explored for signatures of chemical connectivity. We conclude with a discussion in Sec. X.

II. THE FOKKER-PLANCK EQUATION
A general time-dependent stochastic process in N-dimensions, x(t), can be described either by stochastic differential equations or, equivalently, its time-dependent probability density W (x,t) which satisfies a Fokker-Planck equation 10 ∂W (x,t) ∂t where T + is the time ordering operation, which orders time to increase from right to left in the products ofL(x, τ) of Taylor-expansion of the exponential. As can be seen from Eqs. (3) and (4), the differential operatorL(x,t) acts on x, which is the value of the stochastic variable at a later time t, rather than on x ′ at earlier time t ′ . This representation is therefore known as the forward Kramers-Moyal expansion. In a similar way, it is possible to derive a backward representation of Eq. (3), with the operatorL + acting on x ′ rather than x, It has been shown 10 thatL + is actually the adjoint of the Fokker-Planck operatorL and has the following form: The solution of Eq. (6) will then be formally similar to Eq. (4), the difference is that the time ordering T − is now reversed and that the operatorL + acts on x ′ instead of x. Dual representation (4) and (7) of the same quantity P(x,t|x ′ ,t ′ ) will be useful in the following derivations since it allows to switch between the old and new variables x ′ and x. It is also useful to define the propagatorsÛ(x,t; t ′ ) and The operatorÛ + is the adjoint ofÛ. 10 Yet, one should keep in mind that unlike ordinary quantum or classical dynamics, the operatorÛ is non-Hermitian and thusÛ + Û . Equations (4) and (7) can be recast in terms of these propagators as If we measure a dynamical variable A, we can calculate the outcome either by propagating P forward withÛ or by propagating A backward withÛ + . So, physically we are moving forward in either case, but calculate a different intermediate object.

III. PERTURBATION IN THE INTERACTION PICTURE
Consider a stochastic process described with a Fokker-Planck operatorL 0 (x,t). We now impose a perturbation δL(x,t) on our stochastic process (small, time-dependent perturbation of control parameter, such as, for instance, variations of external temperature, or perturbation of kinetic rate coefficients), setting We wish to represent time evolution of the probability density of perturbed stochastic process in terms of the unperturbed operatorL 0 (x,t) and the perturbation operator δL(x,t). To that end, we introduce interaction operatorÛ I defined aŝ whereÛ 0 describes unperturbed stochastic process, Substituting Eqs. (20) and (21) into Eq. (18), we obtain SinceÛ 0 satisfies Eq. (18) for the unperturbed stochastic process, the first terms on each side of Eq. (23) cancel out, leading tô which can be rewritten as using the inverse operatorÛ −1 defined in Eq. (19). Equation (25) can be solved recursively forÛ I in exactly the same way as Eq. (4). Introducing compact notation which stands for the representation of operator δL(x,t) in the interaction picture, the solution of Eq. (25) takes the following form: Using the properties of the operatorÛ, given by Eqs. (12) and (19), one can further simplify Eq. (28), Equation (29) can be written in the final compact form

IV. RESPONSE TO IMPULSIVE PERTURBATIONS
We wish to express the system's response to external perturbations. First, we consider the linear response of a stochastic process x(t) described by the Fokker-Planck operatorL 0 (x,t) and which, at time t = 0, is prepared in a steady state characterized by distribution function ρ ss (x). One can consider an experiment, in which the above stochastic system is subjected to a small perturbation at time t 1 : δL(x,t) = ϵV (x,t)δ(t − t 1 ) (where ϵ ≪ 1 and operatorV (x,t) depends on the problem of interest), and the value of some observable f (x), which is a function of stochastic variable x, is then measured at a later time t 2 . If such perturbation, measurement process is repeated several times, one can calculate the value of f (x(t 2 )) averaged over several realizations of the perturbed stochastic trajectories.
where P(x 2 ,t 2 |x 1 , 0) is the perturbed transition probability. Substituting Eq. (10) into Eq. (31) results in By inserting Eq. (30) into Eq. (32), we then express the perturbed operatorÛ in Eq. (32) using the properties of the unperturbed stochastic process, where we have used steady-state propertyÛ 0 (x 2 ; t 1 , 0)ρ ss (x 2 ) = ρ ss (x 2 ). We next define the impulsive linear response function and making use of Eq. (33), we obtain where we have defined the average of operator over unperturbed trajectories Â 0 ≡ Â (x)ρ ss (x)dx and for brevity omitted x in the arguments of the operatorsÛ 0 andV .
In a similar way, one can calculate the non-linear response functions of a stochastic process. The non-linear response function can be defined as the second derivative of response with respect to the perturbation. In case of impulsive perturbations, the second order term in δL in Eq. (30) survives only if we perturb the stochastic process twice (since impulsive perturbation cannot interact with itself and thus time-ordered integrals τ n > τ n−1 > · · · > τ 1 vanish). We thus consider an experiment with two impulsive perturbations: By substituting the latter into Eq. (30), we obtain the perturbative expansion of the operatorÛ, Defining the second-order response function as and substituting Eq. (36) into Eq. (32), we obtain the compact final expression where τ 2 = t − t 2 and τ 1 = t 2 − t 1 . Expressions (35) and (38) are very similar to the response functions in classical or quantum dynamics and have a clear physical meaning. For instance, from Eq. (38), it follows that the second-order response function is an average of the product of operators which act in the following sequence: perturbation operatorV at time t 1 , evolution operatorÛ 0 from time t 1 to time t 2 , perturbation operatorV at time t 2 , evolution operator U 0 from time t 2 to time t, and a measurement of x at time t. The perturbation operatorV is known and often can be controlled externally, while the evolution operatorÛ 0 depends entirely on the unperturbed Fokker-Planck operatorL 0 , and therefore, various (measurable) correlation functions which contain different occurrences ofÛ 0 operators should contain valuable information on the dynamics of the unperturbed stochastic process governed by the Fokker-Planck operator L 0 . The linear and second-order response functions, given by Eqs. (35) and (38), respectively, thus serve as time-dependent experimentally accessible measures of the unperturbed stochastic process x(t).
In the same way, one can obtain the nth order nonlinear response functions which represents experiments with n impulsive perturbations δL( where τ j = t − t j for j = n and τ j = t j+1 − t j for j = n − 1, . . . , 1.

V. GENERALIZED RESPONSE FUNCTIONS
In Sec. IV, we have derived expressions for nonlinear response functions, which can serve as multi-dimensional (i.e., multi-time) measures of the stochastic process x(t). These measures correspond to the standard protocol involving n perturbations followed by a single measurement. We now generalize these experiments to include m perturbations and k measurements. Consider, for instance, an experiment, in which we measure the value of an observable f (x) at time t 1 , then perturb our stochastic system with δL = ϵV (x,t)δ(t − t 2 ) at time t 2 , and finally measure the value of observable g(x) at time t 3 .
Having the values f (x(t 1 )) and g(x(t 3 )) of the stochastic variable x, one can construct the correlation function ⟨g(x(t 3 )) f (x(t 1 ))⟩, in which the average is taken over various stochastic realizations of the same perturbation protocol (i.e., measurement at t 1 , perturbation at t 2 , and another measurement at t 3 ). This correlation function should depend on the type of perturbation and the time t 2 of the perturbation. Such an experiment should thus provide non-vanishing (twodimensional) quantities The "+" index denotes measurement and "−" denotes perturbation. We shall now calculate the new measure R +−+ . The correlation function ⟨g( Substituting Eq. (29) into Eq. (41), we get Inserting Eq. (42) into Eq. (40), we finally obtain Equation (43) contains the following sequence of operations: measurement of f (x) at time t 1 , time evolutionÛ 0 (t 2 ,t 1 ) of unperturbed stochastic system from t 1 to t 2 , perturbationV at t 2 , time evolution ofÛ 0 (t 3 ,t 2 ) of unperturbed stochastic system from t 2 to t 3 , and a measurement g(x) at t 3 . This is a completely new (experimentally accessible) correlation function, which depends on two time intervals t 3 − t 2 and t 2 − t 1 , and thus should carry additional information about the unperturbed stochastic process.
In a similar way, we can consider various experiments with k perturbations and m = n − k − 1 measurements which are described by a new class of multidimensional measures, in which {+ · · · +, − · · · −} denote various permutations of m "+" (measurements) and k "−" (perturbations). Note that the leftmost symbol being always "+" due to causality. We call these new quantities GRFs. 8,9 GRFs simply extend the class of correlation and response functions. For instance, in the notation of Eq. (44), an ordinary two-time correlation function corresponds to R ++ , while the second-order nonlinear response function corresponds to R +−− . The GRF expression has a simple form and is a correlation function of appropriate ordering of measurements f (x) and perturbationsV with time evolution operatorsÛ 0 (t j ,t j−1 ) in between. It is also useful to recast the expressions for GRFs using Eq. (10) in terms of unperturbed transition probabilities P 0 (x,t|x ′ ,t ′ ). We demonstrate this for two-dimensional GRFs: (a) 3-point correlation function (b) second-order response function (c) GRF corresponding to 1 perturbation followed by 2 measurements, (d) GRF corresponding to a measurement followed by 1 perturbation and then 1 measurement, In (45)-(48),V (x) is perturbation operator specific to the type of perturbation used to disturb the stochastic system. Examples will be given in Secs. VII-IX.

VI. RESPONSE TO PERIODIC PERTURBATIONS; FREQUENCY-DOMAIN SIGNALS
Periodic perturbations are an alternative to impulsive perturbations that can be imposed on a stochastic system 12 via an oscillatory variation of one of its parameters, such as temperature, see also Sec. VII B. Rather than varying timeordered interactions, we now scan frequencies with no control over time ordering.
We consider the perturbed Fokker-Planck operator We first consider the linear response, n = 1. By substituting Eq. (49) into Eq. (30) and then into Eq. (32), we obtain the mean of observable f (x(t)), In case of a time-independent Fokker-Planck operator L 0 (x), the operatorÛ 0 (x; t 2 ,t 1 ) =Û 0 (x; t 2 − t 1 ) solely depends on the time difference t 2 − t 1 , so that this leads to a simple form of the Fourier transform of Eq. (50), The coefficient in front of the delta-function in Eq. (51), which we denote the susceptibility χ (1) (ω), is simply given by the Fourier transform of the time-domain linear response function in Eq. (35), From Eq. (51), one sees that the susceptibility can be measured as the spectral amplitude of the derivative ∂   f (ω)  /∂ϵ 1 at frequency ω 1 , i.e., In case of the second-order response, n = 2, one considers two simultaneous perturbations with two different frequencies δL(x,t) = [ϵ 1 e −ıω 1 tV (x) + ϵ 2 e −ıω 2 tV (x)]/2π. Again, substituting it to Eq. (30), we obtain The second-order susceptibility can then be defined as where the summation p is over the two permutations of frequencies ω 1 and ω 2 . The higher-order nonlinear susceptibilities can be obtained similarly, In analogy to Sec. V, one can also extend the class of nonlinear susceptibilities χ (m,k) to a general class of experiments in which a stochastic system is perturbed on k frequencies ω ′ 1 , . . . , ω ′ k , and m-point correlation function ⟨ f 1 (ω 1 ) . . . f m (ω m )⟩ is measured as a response to these perturbations, It follows from Eq. (57) that the m + k frequencies ω ′ i ,ω j are constrained via a delta-function, such that only m + k − 1 of them are independent. We therefore omit one frequency ω m from the arguments of χ (m,k) and consider (k + m − 1)-th order generalized susceptibility as a function of k + m − 1 frequencies χ (m,k) (ω m , . . . ,ω 2 ,ω ′ k , . . . ,ω ′ 1 ). The generalized nonlinear susceptibility can be expressed as a linear combination of Fourier transforms of the nth order (n = k + m − 1) GRFs R +···,−··· of Sec. V. Below, we present the frequency-domain versions of Eqs. (45)-(48), i.e., generalized 2D susceptibilities corresponding to various combinations of three perturbations and measurements: (a) 3-point correlation function (b) second-order susceptibility for periodic perturbations on frequencies ω ′ 1 and ω ′ 2 , (c) generalized susceptibility corresponding to perturbation on one frequency ω ′ 1 and measurements on 2 frequencies ω 1 and ω 2 such that ω ′ 1 = ω 1 + ω 2 , Note that while in the time domain we had 4 two-time GRFs (Eqs. (45) and (46)), in the frequency domain we have only 3 two-frequency generalized susceptibilities (Eqs. (58)-(60)). This is because we are losing time ordering when switching to the frequency domain, i.e., in the frequency domain, there is no difference between the + + − and + − + techniques since all processes involving of 1 perturbation and 2 measurements are occurring simultaneously.

VII. EXAMPLE 1: THE ORNSTEIN-UHLENBECK PROCESS
The OUP represents stochastic diffusion in a harmonic potential. It is frequently used as a simplified model of diffusion in external force fields, since it can be solved analytically. Its dynamics is governed by the stochastic Langevin equation where ξ(t) is a delta-correlated noise ⟨ξ(t)ξ(s)⟩ = δ(t − s). The corresponding Fokker-Planck equation reads 13 with the Fokker-Planck operator The transition probability of this model reads 10 which in the long time t − t ′ → ∞ limit leads to the stationary distribution

A. Impulsive perturbation of the coordinate
We now subject the OUP to some external perturbation and examine its response. We consider perturbation in a form of sudden displacement of the coordinate x by small ϵ at time t 1 . This perturbation can be described by an additional δ-function term in stochastic differential equation (62), Indeed, if we integrate both parts of equation (67) in the vicinity of t 1 , we get This also implies that in discrete stochastic simulations of Eq. (68), one needs to take such intermediate mesoscopic 23 values of ϵ so that ϵ is small compared to the macroscopic variation of x and at the same time ϵ ≫ √ D∆t. The Fokker-Planck equation with the perturbation in Eq. (68) reads This implies that the perturbation operator is Correspondingly, the infinitesimal perturbation operatorV (x) (defined in the beginning of section V as δL = ϵV (x)δ(t − t 1 )) appearing in Eqs. (45)-(48) readŝ We therefore consider the following simple observables f 1 (x 1 ) = x 1 , f 2 (x 2 ) = x 2 , and f 3 (x 3 ) = x 2 3 . We note, however, that other types of perturbations, like the temperature perturbation discussed in Sec. VII B, can result in non-zero 2D measures for any kind of observable. The exact results for 2D measures R +++ (τ 2 , τ 1 ) and R ++− (τ 2 , τ 1 ) and R +−+ (τ 2 , τ 1 ) and R +−− (τ 2 , τ 1 ) calculated with Eqs. (65), (66), and (71) are shown in Fig. 1. They compare well with the numerical results of Fig. 2, which were calculated using the non-equilibrium approach, Eq. (44). Both approaches are equivalent. We note that some signals look similar. This is to be expected for Gaussian processes, where these measures are related to each other via the generalized fluctuation-dissipation relations. 9

B. Periodic temperature modulation of the OUP
Harmonic perturbations can be realized by periodically modulating the strength of stochastic noise ξ(t) in Eq. (62), for instance, by changing the temperature. For a general stochastic process, of diffusion in potential U(x), periodic temperature modulation will result in periodic modulation of diffusion constant D = D 0 + ϵ cos(ωt), which depends linearly on temperature. The perturbation operator then reads The operator of infinitesimal perturbation appearing in Eqs.
(58)-(60) is then given bŷ By substituting Eq. (74) in Eqs. (70)-(60), one can calculate the exact 2D susceptibilities. Alternatively, these can be calculated using the non-equilibrium approach, Eq. (57), by taking derivatives of correlation functions with respect to the amplitudes of oscillatory perturbations. This mimics the experimental protocol. Results of both calculations are shown in Fig. 3. Since temperature perturbations correspond to the second-order (even) derivative, then χ (3,0) and χ (1,2) must vanish for symmetric potentials since they effectively involve averaging over the product of odd number of random variables. This can be one of the merits of the generalized 2D measures as compared to the regular 2D correlation and response functions for studying stochastic processes by temperature perturbations.

VIII. EXAMPLE 2: IMPULSIVE TEMPERATURE PERTURBATIONS OF DIFFUSION PROCESS IN ANHARMONIC AND DOUBLE-WELL POTENTIALS
Diffusion in anharmonic, such as double-well, potentials is an important process, which is often used to model biological processes. 14 For instance, conformational dynamics of proteins is often described as a diffusion in double well potential along a reaction coordinate. 15 One can study such stochastic processes using temperature perturbations since it is hard to control perturbations of stochastic variables and often the relevant stochastic variable (reaction coordinate) is not known. By applying temperature perturbations at different times, one can measure various multidimensional GRFs as described in Secs. IV-VII and obtain, e.g., information on conformational dynamics of single proteins. It is interesting to examine the GRFs for different potentials and what information they can provide. Such experiment can be done, for instance, by measuring the electron-transfer (ET) rate between donor and acceptor on different ends of a protein at different times. 16 The ET rate constant depends exponentially on the distance between donor and acceptor k(t) = k 0 exp(− βx(t)); thus, the stochastic variable X(t) ≡ − log[k(t)] can be a measure of conformational dynamics of the protein.
Temperature perturbations can be either impulsive ∆T δ(t − τ) (such as pulsed heating 17,18 ) or a step T-jump perturbation ∆T θ(t − τ), which in some cases may be easier to implement experimentally. For the purpose of consistency with the discussion in Secs. IV and V, here we consider impulsive temperature perturbations. Since δ(t − τ) = −dθ(t − τ)/ dτ, GRFs of T-jump perturbations can be easily obtained from impulsive GRFs by integration over the corresponding perturbation times.
In Fig. 4, we compare 2D generalized response functions obtained from numerical simulations of an overdamped diffusion process x(t) in non-equilibrium experiments with impulsive perturbations of temperature. We consider symmetric double well U(x) = x 2 + 3e −x 2 and asymmetric double well U(x) = x 2 + 3e −x 2 − 0.5x potentials and also harmonic U(x) = x 2 /2 and Morse U(x) = 10(1 − exp(−0.3x)) 2 potentials for comparison. The observable in the GRFs was taken to be the stochastic variable itself Fig. 4 shows that 2D GRFs are qualitatively different for harmonic, Morse, and double-well potentials. In particular, the harmonic system only has a single non-vanishing GRF, R ++− (τ 2 , τ 1 ). In fact, R +++ (τ 2 , τ 1 ) and R +−− (τ 2 , τ 1 ) vanish for both the harmonic and symmetric double-well potentials which is the consequence of symmetry. Second, Fig. 4 shows that both Morse and double well potentials have a strong signal R +−+ (τ 2 , τ 1 ), which should be an indicator of anharmonicity. Third, the decay of R +−+ (τ 2 , τ 1 ) along τ 1 in the symmetric and asymmetric double-well potential is significantly slower than the decay of R ++− (τ 2 , τ 1 ). This is related to the fact that the lowest non-zero eigenvalue of the Fokker-Planck operator (proportional to the Kramers over-the-barrier escape rate) is significantly smaller than the second non-zero eigenvalue in double-well potentials (see also Sec. VIII A). Fourth, 2D R +−+ (τ 2 , τ 1 ) of a symmetric potential at long times decays with τ 1 + τ 2 , while in asymmetric double-well potentials, the symmetry along the τ 1 = τ 2 direction is broken. This should provide an experimental measure of the asymmetry of the potential. In Sec. VIII A, we discuss quantitative aspects of these observations.
Third, from Eqs. (79), one can see that R ++− (τ 2 , τ 1 ) decays along τ 1 as exp(−λ 2 τ 1 ), while R +−+ (τ 2 , τ 1 ) decays as exp(−λ 1 τ 1 ). In double-well potentials, especially with high barriers, the first non-zero eigenvalue λ 1 is much smaller than the second non-zero eigenvalue λ 2 , which explains the difference in decays of R ++− (τ 2 , τ 1 ) and R +−+ (τ 2 , τ 1 ) observed numerically in Sec. VIII. λ 1 is due to the barrier interaction (or possibly "tunneling") splitting of the ground state in symmetric double well potential (see Ref. 10 on similarity of Schrödinger and Fokker-Planck equations) and is proportional to the barrier height, while the second non-zero eigenvalue λ 2 is the second "energy" level in the local well and is proportional to the width of the well. Thus, in double-well potentials with high energy barriers, λ 1 ≪ λ 2 . Numerically, λ 1 = B exp(−E b /κT) (when multiplied by probability to stay in the well it is exactly the Kramers thermally activated escape rate 14,19 ) and one can determine the barrier height E b by measuring the decay rate λ 1 . We also note that the different decay of R ++− (τ 2 , τ 1 ) and R +−+ (τ 2 , τ 1 ) is due to the fact that different GRFs preserve different terms in expression (76). Selectivity property of GRFs was also observed in nondissipative systems. 8 In weakly asymmetric potentials, we expect the longtime behavior of R ++− (τ 2 , τ 1 ) and R +−+ (τ 2 , τ 1 ) to be dominated by the terms in Eqs. (79). Thus, for instance, by measuring where ξ i (t)ξ j (t ′ ) = σδ i j δ(t − t ′ ), and numerical values of parameters are given in Ref. 9. We numerically solve Eqs.
(81) and calculate R l+, j+,i− (τ 2 , τ 1 ) = ∂ ∂∆N i N l (t 3 )N j (t 2 ) via the non-equilibrium approach for different i, j, l. The typical behavior of X 1 (t), X 2 (t), X 3 (t), and X 4 (t) is shown in Fig. 5(a), while the calculated GRFs R l+, j+,i− (τ 2 , τ 1 ) are shown in Fig.  5(b). The 2D plots in Fig. 5(b) look the same as in Fig. 2 of Ref. 9. As discussed in that reference, these 2D plots can suggest the order of X i , X l , and X j species in the chemical network, i.e., out of two permutations X i − X l − X j and X i − X j − X l , the one with longer correlation tail along the τ 2 axis corresponds to the correct order of species X i , X l , and X j . In particular, from Fig. 5(b), it follows that species X 1 , X 2 , X 3 , and X 4 are ordered as GRFs provide constraints with characteristic signatures that may allow to distinguish between possible reaction paths. An open question for future studies is the possibility to invert data to get the reaction topology. The same problem arises in multidimensional NMR (inverting distances between protons obtained from 2D data to get molecular structure), multidimensional IR to get secondary structure of proteins, and multidimensional electronic spectroscopy to get energy transfer paths in chromophore aggregates. Research on relation of GRFs to reaction topology is currently underway. Our preliminary results suggest that GRFs can be used to distinguish kinetic orders of chemical reactions, to determine kinetic parameters of elementary steps in complex chemical reactions, to identify feedback loops in reaction networks, etc.

X. DISCUSSION
Multidimensional techniques have now become standard tools in spectroscopy. However, these are not limited to optical pulses and electronic or vibrational transitions. One can use similar multidimensional ideas with different types of perturbations, to study arbitrary dynamical systems. In this paper, we considered stochastic dynamical systems and showed that their dynamics can be studied via perturbations of their coordinates or parameters, such as temperature. Moreover, perturbations and measurements can be mixed leading to new types of multidimensional measures of stochastic dynamics, which we call generalized response functions. Generalized response functions can be measured experimentally either in time or the frequency domains. Information contained in GRFs can provide additional insights on stochastic microscopic processes. In particular, the examples given in this paper showed that GRFs can reveal the symmetry or bistability of the mean field potential of stochastic process and can provide information on its anharmonicity and internal barrier height.
We also propose to use multiple impulsive perturbations of temperature to study conformational dynamics of single proteins. Conformational dynamics of proteins plays an important role in regulating their biological functions, but is not well understood. Temperature perturbations provide an accurate and controllable way to probe conformational dynamics. Multiple perturbations and measurements provide a wide range of new tools-experimentally measured GRFs, which entirely depend on equilibrium fluctuations of protein, and thus serve as multidimensional measures of its stochastic conformational dynamics.
GRFs may be also used to study chemical dynamics of small chemical systems such as biocells. Such finite chemical systems contain a small number of molecules of reacting species and thus cannot be described via macroscopic chemical kinetics. A stochastic description is more appropriate.
Because of the small number of molecules, fluctuations will be significant and could hide any response to external perturbations. Time-dependent averaged quantities, such as correlation and response functions, are thus required to extract information about the underlying microscopic chemical dynamics. The proposed GRF approach of mixed perturbations and measurements provides a wide range of new multi-point correlation functions, or measures, characterizing microscopic chemical dynamics. Information contained in these new measures should be explored for specific systems. In particular, in Sec. IX, we showed how 2D GRFs R ++− (τ 2 , τ 1 ) can be used to study connectivity of stochastic reaction networks.