**This page has been visited total ****since 11 June
2009**

*EROPHILE*

Eigenvector representation of observables and
probabilities in a high-dimensional Euclidean space.

or go to
JCP to see J. Chem. Phys. **130**, 044905 (2009);
DOI:10.1063/1.3063118

**Permalink:**
http://link.aip.org/link/?JCPSA6/130/044905/1

See also the "Dynamical Integretion over a
Markovian web "

See also the "Integretion over
a Marcovian web"

**EROPHILE****: A geometric representation for systems evolving through infrequent
transitions in a network of states**

Authors: Georgios C. Boulougouris^{1}
and Doros N. Theodorou^{2}

Affiliation:

*School
of Chemical Engineering, National Technical University of Athens, Zografou
Campus, GR-15780 Athens, Greece. *Fax: +30 210 772 3112

^{1}E-mail: gboul@chemeng.ntua.gr

^{2}E-mail: doros@central.ntua.gr.

**The dynamics of many
physical, chemical, and biological systems can be reduced to a succession of
infrequent transitions in a network of discrete states representing low energy
regions in configuration space. This enables
accessing long-time dynamics and predicting macroscopic properties. Here we develop a novel, perfectly general
statistical mechanical/geometric formulation that expresses both state
probabilities and all observables in the same Euclidean space, spanned by the
eigenvectors of the symmetrized time evolution
operator. Our formalism leads to simple expressions for nonequilibrium
and equilibrium ensemble averages, variances, and time correlation functions of
any observable and allows a rigorous decomposition of the dynamics into
relaxation modes. Applying it to sub-glass segmental relaxation in atactic
polystyrene up to times on the order of 10 ms, we probe the molecular mechanism
of the g and d processes and unequivocally identify the d process with rotation of a single phenyl
group around its stem. **

The dynamics of many physical,
chemical, and biological systems proceeds as a succession of transitions
between small regions in their configuration space, which we shall call
“states”. The states constitute
“basins” of low energy with respect to the generalized coordinates spanning
configuration space, or of low free energy with respect to a set of order
parameters providing a coarse-grained description of the system. Each state contains one or more local minima
of the (free) energy. Transitions
between states are infrequent events, in the sense that the mean waiting time
for a transition out of a state is long in comparison to the time required for
the system to establish a restricted equilibrium distribution among
configurations in the state. The entire configuration space can be tessellated into
states. Representing each state by a
point in configuration space and connecting all pairs of states between which a
transition is possible, one obtains a graph, or network of states. Examples of
phenomena that can be modelled as occurring
through a succession of transitions in a network of states include diffusion of
point defects and impurities in metals and semiconductors;^{1} of gas
molecules in amorphous polymers;^{2} of bulky hydrocarbons in microporous solids, such as zeolites;^{3}
structural relaxation and plastic deformation in glasses;^{4} phase
transitions in molecular and atomic clusters;^{5} surface diffusion;^{6}
protein folding;^{7} chemical reactions;^{8} and the function
of gene regulatory networks for protein expression.^{9}

The possibility of
coarse-graining dynamics into a sequence of transitions in a network of states
is of strategic importance for understanding and predicting macroscopic
properties from atomic-level structure and interactions. The longest times that can be simulated with
“brute force” molecular dynamics are microseconds, i.e., too short by many
orders of magnitude in comparison with the experimental time scales of most
phenomena of interest. A more efficient
strategy is to construct the network of states *i*
and compute the rate constants *k _{i}*

In this article we
devise a simple geometric interpretation for the time-dependent state
probabilities and for any observable in a system evolving through a sequence of
infrequent transitions in a network of states. This interpretation employs a Euclidean
space of dimensionality equal to the number of states, spanned by the
eigenvectors of the rate constant matrix, appropriately symmetrized. The time-dependent expectation values and
autocorrelation functions of all observables are obtained from dot products in
that space. Furthermore, the dynamics of
any observable is decomposed into contributions from different relaxation
modes, allowing a direct connection with spectroscopic measurements. As an
example application we study the reorientational
dynamics of pendant groups during physical ageing of glassy atactic
polystyrene.

Consider a
system capable of adopting *N* discrete states, among which it undergoes
infrequent transitions. We will denote
by *P _{i}*(

_{} , or _{} (1)

with given initial distribution among
states **P**(0). *k _{i}*

_{} (2)

Two clarifications are in order
here. First, states are not necessarily
single basins around individual minima on the system’s energy hypersurface. Basins
communicating through fast transitions, relative to the in-basin residence
time, are lumped into the same state, such that the time scale separation
criterion holds for exit from the state. This lumping into “metabasins”
has been discussed in the context of glass-forming liquids.^{13,14} Secondly,
for many systems the full set of states is not known a priori. This is not prohibitive for invoking the
formulation to be developed herein. In
our Dynamic Integration on a Markovian Web (DIMW) approach, we have shown^{15} how one can start off from a small set of
states (those contributing nonzero elements to **P**(0)) and progressively
augment it “on the fly”, using a first passage time criterion. In DIMW a sequence
of master equations is solved
on an ever increasing network of states. In the solutions one can discern plateau
regions, over which **P**(*t*) remains
practically independent of time; these are indicative of restricted
equilibration of the system over a subset of states, before it moves on across
the high energy barriers which surround the subset. The formulation we present is applicable for
times up to the end of such a plateau region, with *N* being the number of
states in the subset and **P**(*∞*) the restricted equilibrium values attained by
the probabilities of these states at the plateau.

In
many applications, the master equation (1) is solved numerically, generating
stochastic trajectories by Kinetic Monte Carlo simulation.^{2,3,6} This is convenient, but may lead to
unacceptably large statistical error in multistate systems with rugged energy
hypersurfaces.^{15} A preferable
alternative is to solve the master equation analytically, via spectral
decomposition.

Here
we introduce an Eigenvalue Representation of Observables
and Probabilities in a HIgh-dimensionaL Euclidean
space (EROPHILE) for systems evolving according to
master equation (1). Our starting point
is to transform^{16} the state probability vector **P**(*t*)
into a reduced vector _{} with elements _{}. This satisfies the reduced master equation ** _{} **with

_{} (3)

where the normalization condition _{} has been used in separating out the
equilibrium contribution (_{}). The eigenvectors _{} form an orthonormal
basis set: _{}, 0 £ *m*, *n *£ *N*-1 (4)

with d* _{mn}* being the Kronecker
delta. They also satisfy

We
now consider the *N-*dimensional linear space from which the vector of
reduced probabilities _{} takes values.
From equations (3) and (4) it is clear that the eigenvectors _{} of the symmetrized
rate constant matrix _{}constitute
an orthonormal basis set for this space. Furthermore, by equation (4), at any time *t* ³ 0, the
tip of the reduced probability vector _{} lies on a (*N*-1)-dimensional hyperplane which passes through point _{} and is spanned by _{}.
This hyperplane is normal to the equilibrium
eigenvector _{}=_{}. It is described by the equation _{} and intersects each _{}-axis
at a distance _{}from
the origin. Note that
this hyperplane, on which _{}is bound to evolve, is dictated solely by the equilibrium state
probabilities, therefore by the evolution operator **K**, and is independent
of the initial probability distribution among states. A pictorial
representation of the hyperplane for the simple case
of a system with three states is given in figure 1a.

From
the solution, equation (3), it is clear that one can prepare the initial
conditions in such a way that the probability distribution among states decays
to equilibrium along a single exponential.
If one sets

_{}, or _{} (1 £ *i* £ *N* , 1 £ *m* £ *N*-1) (5)

then, by virtue of the orthonormality of eigenvectors, equation (4), the
time-dependent distribution among states will evolve as _{}, or _{}.

In this case we will say that the
system has been initially perturbed, and decays to equilibrium, along mode *m*. Each mode corresponds to a cascade of interstate
transitions grouped together in such a way as to bring about a redistribution
of the occupancy probabilities of all states along a single eigenvector.

In
the more general case, any initial probability distribution can be expressed as
a linear combination of *N-*1 perturbation vectors from the equilibrium
distribution: _{}, with _{} (6)

The resulting solution will be: _{}. (7)

In _{}-space,
the time-dependent probability distribution of the system will move along a
curved line lying on the hyperplane spanned by eigenvectors _{},
until it collapses on point _{} (see figure 1a).

Another
interesting special set of initial conditions is that wherein the system is
initiated in a single state *i*: _{}. (8)

This is readily cast in the form of
equation (6) with _{}. (9)

The time-dependent evolution of
state probabilities follows from equation (7).

We now take a significant step
forward, which will enable us to describe any observable as a vector “living”
in the same Euclidean space of reduced state probabilities depicted in figure
1.

Let
*A* be any observable depending on the configuration of the system. We symbolize by *A _{i}* the value of

_{} (1 £ *i* £ N) (10)

The Euclidean norm of _{} has the physical meaning of the root mean
square value of observable *A* at equilibrium:

_{} (11)

The time-dependent ensemble average
of *A*, starting from an arbitrary initial distribution among states, is
readily obtained using equation (8) as

_{} (12)

where we have set _{}.

The equilibrium expectation value of
the observable is simply

_{}. (13)

The equilibrium variance of *A*
is obtained as

_{} (14)

Expressing _{} in the basis of eigenvectors
_{}, _{}, we
can further write _{} (15)

The above properties have an
immediate geometric interpretation in our Euclidean _{}-space
(see figure 1b). We draw _{} as a vector in that space, using any
consistent units to measure the state values of the observable, *A** _{i}*. By equations (12)-(15),

Of particular importance in
understanding the dynamics and making contact with spectroscopic measurements
are time autocorrelation functions of the form _{}. In general, these are averages over the nonequilibrium ensemble representing the system under given
**P**(0), and therefore depend on both the time
origin *t*_{0} and the time interval *t*. Specializing to the case *t*_{0}=0,
we write

_{} (16)

where _{}stands for the conditional
probability of finding the system in state *j* provided it was in state *i* at *t*=0, and _{}.
From equation (7), _{} can be written as _{} (17)

with _{} given by Eq. (9).

Substituting Eq. (17) into Eq. (16)
one obtains:

_{} (18)

If the system starts
off at equilibrium, _{},
equation (18) simplifies further to

_{}
(19)

By equation (19), the equilibrium
time autocorrelation function of excursions from the mean admits a simple
decomposition into the modes, the amplitude of each mode being the squared
projection of the reduced observable vector onto the mode eigenvector. The dominance of a mode (or a set of modes)
in a relaxation measurement can be assessed through the changes in observable β_{n}_{ }occurring along the
mode(s).

For
the special case where the observable is a component of a three-dimensional
vector _{}** **characterizing the configuration of the
system, equation (18) gives:

_{}, (20)

where _{} symbolizes component *i*_{dim}
Î **v**^{(l)}
in state *i*, _{}, and
_{}.

For the autocorrelation function of the entire vector _{} at equilibrium, equation (20) gives:

_{} (21)

where _{}.

Subtracting equation (21) evaluated
at _{} from equation (21), one obtains:

_{} (22)

Equation
(18) informs us how the autocorrelation function behaves for a system perturbed
from equilibrium. We express the initial probability distribution in terms of
mode contributions via equation (6). We
also express the observable value in each state *i*
as _{} , with β* _{m}*
defined in equation (12). It is then
easy to recast equation (18) as

_{} (23)

In the limit of an equilibrium
initial distribution (_{}),
one recovers equation (19) from equation (23), as expected from linear response
theory. In general, the autocorrelation
function, equation (23), contains all relaxation modes on which the projection _{} of the observable is nonzero, even in cases
where the system is initially perturbed away from equilibrium along a single
mode.

EROPHILE is applicable to any stochastic
system, from economics to medicine.
Although developed here in a statistical mechanical context, it is as
general as statistics itself.

We apply EROPHILE
to analyze the sub-glass dynamics of a united atom model of atactic polystyrene
(a-PS) that we have recently simulated at 250 K and 1 bar by the DIMW method.^{15} DIMW enabled us to obtain an analytical
representation of the dynamics of the glassy polymer in its multidimensional
configuration space, based on solving the master equation recursively in an
ever-expanding network of states. Each
state is a basin constructed around an inherent structure. Rate constants for transitions from one state
to the other were computed via multidimensional transition state theory on the
rugged potential energy landscape of the polymer glass, whereas saddle points
around a given minimum were located by the “dimer” method.^{17} The
evolution of the probability distribution has been tracked^{15} over
more than 10 orders of magnitude on the time scale, reaching 10^{-5}s.

**Acknowledgements**

This work has been supported by the
Senate Committee on Basic Research of NTU Athens in
the form of a 2006 “Karathéodory” programme. Simulations have been performed in part on the Mare Nostrum
supercomputer at the

**References**

1 Vineyard, GH. Frequency factors and isotope effects in solid state
rate processes. *J. Phys. Chem. Solids* **3**, 121
(1957).

2 Gusev, A.A., Müller-Plathe,
F.; van Gunsteren, W.F.; Suter,
U.W. Dynamics of small molecules in bulk polymers. *Advan.** Polym. Sci.* **116**, 207-247 (1994).

3 Theodorou, D.N., Snurr, R.Q., *Comprehensive**
Supramolecular Chemistry* (Elsevier Science, Oxford, 1996), p. 507-548.

4 Stillinger, F.H. A topographic view of supercooled liquids and glass formation. *S**cience* **267**, 1935-1939 (1995).

5 *Advan.** Chem. Phys.* **115**, 1-111 (2000).

6 Voter, A.F.& Doll, J.D. Dynamical corrections to
transition-state theory for multistate systems – Surface self-diffusion in the
rare-event regime. *J. Chem. Phys. ***82**, 80-92
(1985).

7 Wolynes, P.G. Energy landscapes
and solved protein-folding problems, *Phil. Trans. Roy. Soc.
A* **363**, 453-464 (2005).

8 Wei, J. & Prater, C.D.
The structure and analysis of complex reaction systems. In *Advances in Catalysis*,
Vol. 13 (New York, Academic, 1962).
p. 204-390.

9 Kaznessis, Y.N. Multi-scale models for gene network engineering, *Chem.
*.

10 *J. Chem. Phys*. **68**, 2959 (1978).

11 *Classical and Quantum Dynamics in Condensed Phase Simulations* (World
Scientific, Singapore, 1998). p. 51-66.

12 van Kampen, N.G. *Stochastic
Processes in Physics and Chemistry*, North-Holland,

13 Doliwa, B.
& Heuer, A. What does the potential energy
landscape tell us about the dynamics of supercooled
liquids and glasses? *Phys. **Rev. Lett.* **91**, 235501
(2003).

14 Tsalikis, D., Lempesis,
N., Boulougouris, G.C., Theodorou, D.B. On the role of inherent structures in
glass-forming materials II. Reconstruction of the mean square displacement by
rigorous lifting of the inherent structure dynamics. *J. Phys. Chem. B*, in press.

15 Boulougouris,
G.C. & Theodorou, D.N. Dynamical
Integration of a Markovian Web: A first passage time approach. *J. Chem. Phys.* **127**,
084903 (2007).

16 Reichl, L.E. *A Modern Course in Statistical
Physics*, 2nd ed.

17 Henkelman, G. & Jónsson, H.
A dimer method for finding saddle points on high
dimensional potential energy surfaces using only first derivatives. *J. Chem. Phys.* **111**, 7010-7022 (1999).

18 Yano, O. &
Wada, Y. Dynamic
mechanical and dielectric relaxations of polystyrene below glass
temperature. *J. Polym.
Sci., Part **Α**-2: Polym.
Phys*. **9**, 669-686 (1971)

19 Crissman, J.M.
& McCammon, R.D. Apparatus for Measuring the
Dynamic Mechanical Properties of Polymeric Materials between 4° and 300°K *J. Acoust. Soc. Am*. **34**, 1703- 1706 (1962)

20 Hunt, B.I., Powles,
J.G., Woodward, A.E. Proton
spin-lattice relaxation measurements on some high polymers of differing
structure and morphology Polyethylene, poly-4-methyl-pentene-1, poly-l-leucine, poly-phenyl-l-alanine,
polystyrene and poly-α-methyl-styrene. *Polymer ***5**, 323- 342 (1964)

21 Illers, K.H. & Jenckel,
E. Dynamic mechanical behavior of polystyrene, poly(p-chlorostyrene), and
poly(p-bromostyrene) at low temperatures. *J. Polym. Sci., Part B: Polym.
Lett*. **41**, 528-531 (1959).

**Figure 1: Schematic of the EROPHILE representation
fo****r a simple three-state system. a: State
probabilities. **Any possible distribution of the
system among its states lies within a hyperplane
normal to the equilibrium unit vector _{} in the Euclidean space of reduced
state probabilities _{}. The hyperplane
contains point _{} and intersects each of the _{} axes at _{}. The hyperplane is spanned by the *N-*1 eigenvectors _{}(*n*=1, …, *N*-1) of the symmetrized
rate constant matrix. Any perturbation
of the probability distribution from equilibrium is written as a linear
combination of the latter *N-*1 eigenvectors _{}. The curved line on the hyperplane traces a trajectory of the system starting at
initial reduced probability distribution _{}and
ending at equilibrium. The quantities *a** _{n}* are projections of

From 23/9/2009, Thanks to revolvermaps http://www.revolvermaps.com/