^{1}and Markos A. Katsoulakis

^{1,a)}

### Abstract

We propose a new sensitivity analysis methodology for complex stochastic dynamics based on the relative entropy rate. The method becomes computationally feasible at the stationary regime of the process and involves the calculation of suitable observables in path space for the relative entropy rate and the corresponding Fisher information matrix. The stationary regime is crucial for stochastic dynamics and here allows us to address the sensitivity analysis of complex systems, including examples of processes with complex landscapes that exhibit metastability, non-reversible systems from a statistical mechanics perspective, and high-dimensional, spatially distributed models. All these systems exhibit, typically non-Gaussian stationary probability distributions, while in the case of high-dimensionality, histograms are impossible to construct directly. Our proposed methods bypass these challenges relying on the direct Monte Carlo simulation of rigorously derived observables for the relative entropy rate and Fisher information in path space rather than on the stationary probability distribution itself. We demonstrate the capabilities of the proposed methodology by focusing here on two classes of problems: (a) Langevin particle systems with either reversible (gradient) or non-reversible (non-gradient) forcing, highlighting the ability of the method to carry out sensitivity analysis in non-equilibrium systems; and, (b) spatially extended kinetic Monte Carlo models, showing that the method can handle high-dimensional problems.

This work was supported, in part, by the Office of Advanced Scientific Computing Research, U.S. Department of Energy under Contract No. DE-SC0002339 and the EU Project No. FP7-REGPOT-2009-1 “Archimedes Center for Modeling, Analysis and Computation.” We would like to thank Andrew Majda, Petr Plecháč, Luc Rey-Bellet, and Dion Vlachos for many interesting and valuable discussions as well as Giorgos Arampatzis for providing us with the ZGB simulation algorithm. Finally, we would like to thank the reviewers for their helpful comments which resulted in the improvement of the present manuscript.

I. INTRODUCTION

II. DISCRETE TIME MARKOV CHAINS

III. CONTINUOUS-TIME MARKOV CHAINS

IV. FURTHER GENERALIZATIONS

V. NUMERICAL EXAMPLES

A. Statistical estimators for RER and FIM

B. Schlögl model

C. Reversible and non-reversible Langevin processes

D. Spatially extended kinetic Monte Carlomodels

VI. CONCLUSIONS

### Key Topics

- Entropy
- 59.0
- Field ion microscopy
- 33.0
- Markov processes
- 26.0
- Stochastic processes
- 25.0
- Monte Carlo methods
- 20.0

## Figures

Upper plot: The number of *X* molecules as a function of time. The stochastic process sequentially visits the two most probable states defined as the maxima of the PDF. Lower panel: RER as a function of time when *k* _{1} *A* is perturbed by 0.05 computed using (28) (dashed line) and using (30) (grey line). In both cases, the accuracy of the numerical estimators increase as the number of samples increases.

Upper plot: The number of *X* molecules as a function of time. The stochastic process sequentially visits the two most probable states defined as the maxima of the PDF. Lower panel: RER as a function of time when *k* _{1} *A* is perturbed by 0.05 computed using (28) (dashed line) and using (30) (grey line). In both cases, the accuracy of the numerical estimators increase as the number of samples increases.

Upper plot: Exact (black), estimated by (28) (grey), and estimated by (30) (white) RER for various directions. *k* _{2} is the most sensitive parameter followed by *k* _{1} *A* while the least sensitive parameters are *k* _{4} and *k* _{3} *B*. Lower plot: RER when parameter θ_{ k } is perturbed by +ε_{0} (black), when perturbed by −ε_{0} (white) and when computed by FIM (grey).

Upper plot: Exact (black), estimated by (28) (grey), and estimated by (30) (white) RER for various directions. *k* _{2} is the most sensitive parameter followed by *k* _{1} *A* while the least sensitive parameters are *k* _{4} and *k* _{3} *B*. Lower plot: RER when parameter θ_{ k } is perturbed by +ε_{0} (black), when perturbed by −ε_{0} (white) and when computed by FIM (grey).

Upper plot: The stationary distributions for the unperturbed process (solid line), the most sensitive parameter *k* _{2} (dashed line), and the least sensitive parameter *k* _{3} *B* (dotted line). Lower plot: The stationary distributions for the unperturbed process (solid line), when the most sensitive parameter *k* _{2} is perturbed (dashed line), and when the most sensitive direction ε_{max } is perturbed (dotted line).

Upper plot: The stationary distributions for the unperturbed process (solid line), the most sensitive parameter *k* _{2} (dashed line), and the least sensitive parameter *k* _{3} *B* (dotted line). Lower plot: The stationary distributions for the unperturbed process (solid line), when the most sensitive parameter *k* _{2} is perturbed (dashed line), and when the most sensitive direction ε_{max } is perturbed (dotted line).

Upper plot: Relative entropy rate as a function of time for perturbations of *D* _{ e } (solid line), *a* (dashed line), and of *r* _{ e } (grey line) at the reversible regime (α = 0). The variance of the numerical RER is large, necessitating more samples for accurate estimation. Lower plot: RER for various directions where the most sensitive parameter is *a*. We observe that the FIM-based RER which is always positive has substantially less variance than the directly-computed numerical RER.

Upper plot: Relative entropy rate as a function of time for perturbations of *D* _{ e } (solid line), *a* (dashed line), and of *r* _{ e } (grey line) at the reversible regime (α = 0). The variance of the numerical RER is large, necessitating more samples for accurate estimation. Lower plot: RER for various directions where the most sensitive parameter is *a*. We observe that the FIM-based RER which is always positive has substantially less variance than the directly-computed numerical RER.

Upper plots: Level sets (or neutral spaces) for the reversible case (α = 0). Lower plots: Level sets for the irreversible case (α = 0.1).

Upper plots: Level sets (or neutral spaces) for the reversible case (α = 0). Lower plots: Level sets for the irreversible case (α = 0.1).

Upper plot: Relative entropy rate as a function of time for perturbations of both *k* _{1} (solid line) and of *k* _{2} (dashed line). An equilibration time until the process reach its metastable regime is evident. Lower plot: RER for various directions. The most sensitive parameter is *k* _{1}.

Upper plot: Relative entropy rate as a function of time for perturbations of both *k* _{1} (solid line) and of *k* _{2} (dashed line). An equilibration time until the process reach its metastable regime is evident. Lower plot: RER for various directions. The most sensitive parameter is *k* _{1}.

Typical configurations obtained by ε_{0}-perturbations of the most and least sensitive parameters. The comparison with the reference configuration reveals the differences between the most and least sensitive perturbation parameters.

Typical configurations obtained by ε_{0}-perturbations of the most and least sensitive parameters. The comparison with the reference configuration reveals the differences between the most and least sensitive perturbation parameters.

Vector field with the most (solid arrows) and least (dashed arrows) sensitive directions computed from eigenvalue analysis of FIM. The length of the arrows is proportional to the corresponding eigenvalue.

Vector field with the most (solid arrows) and least (dashed arrows) sensitive directions computed from eigenvalue analysis of FIM. The length of the arrows is proportional to the corresponding eigenvalue.

## Tables

The rate of the *k*th event when the number of *X* molecules is *x* is denoted by *c* _{ k }(*x*). Ω is the volume of the system.

The rate of the *k*th event when the number of *X* molecules is *x* is denoted by *c* _{ k }(*x*). Ω is the volume of the system.

Parameter values for the Schlögl model.

Parameter values for the Schlögl model.

Parameter values for the discretized Langevin system.

Parameter values for the discretized Langevin system.

The rate of the *k*th event of the *j*th site given that the current configuration is σ is denoted by *c* _{ k }(*j*; σ) where n.n. stands for nearest neighbors.

The rate of the *k*th event of the *j*th site given that the current configuration is σ is denoted by *c* _{ k }(*j*; σ) where n.n. stands for nearest neighbors.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content